Ground-borne vibration due to railway traffic is an increasingly important environmental concern. Following earlier developments of metamaterials in the fields of electromagnetics, photonics, and phononics, this research aims to design an array of resonators installed between a railway track and buildings for vibration reduction. We develop an optimization method to enhance the vibration reduction performance of the metamaterials in a wide frequency range between 1 and 80 Hz and an extended volume of soil. The optimization relies on a 3D coupled finite element – boundary element method. We consider the wavefield generated by a point load at the soil's surface to evaluate the vibration reduction. The resonators are modeled as single-degree-of-freedom systems on top of square concrete foundations and the soil is represented as a horizontally layered elastic halfspace. In the optimization problem, the dynamic properties of each row of the resonators are design variables. Three alternative objective functions are considered: (1) the average displacement of a grid of receiver points is minimized; (2) the energy dissipated by the metasurface is maximized; and (3) the power flow through a vertical surface in the soil behind the metasurface is minimized. In layered soils, vibration propagates as a superposition of multiple dispersive surface waves resulting in a non-uniform response behind the metasurface. Therefore, we introduce power flow analysis as a performance metric. A local optimization method with a gradient-based algorithm is used to solve the optimization problem imposing an upper limit for the total mass of the resonators. Results obtained for the three formulations of the optimization problem are compared. In all three cases, the optimizer converges to a solution with the upper limit of mass, as expected. In cases (1) and (2), no specific pattern is observed for the resonance frequencies of the different rows of the resonators. In case (3), however, the optimizer converges to a classic metawedge where vibration reduction relies on rainbow trapping. When the performance of different solutions is evaluated according to the metrics not considered in the particular optimization problem, it is found that case (1) shows high vibration reduction and case (2) shows high power dissipation at some frequencies. case (3), however, shows a uniform reduction in power flow through a vertical surface in the frequency range of interest.