Matrix sparse array optimization for 3D ultrasound imaging often uses simulated annealing (SA) algorithms to reduce the number of elements and to find their best position, intensity and size. The aim of optimization is to obtain the best possible radiated beam pattern in terms of main lobe width and side lobe level. This has a direct impact on the image resolution and contrast. However, to limit the convergence time, the calculation of the cost function that drives the optimization was usually based on analytically computed radiated beam patterns. In this work, realistic acoustic simulations are based on the FIELD II software, and a novel multi-depth cost function is introduced. Such a cost function allows taking into account the radiated pressure field at several depths, and not only at the focal depth. The multi-depth cost function was inspired by the recently introduced Blackman-tapered 256-element spiral array [1]. The positions of 256 elements of two 7 MHz sparse arrays were obtained by optimizing the beam patterns at: a) the focal depth only (opti1d); b) three different depths (15, 25 and 35 mm) (opti3d). Both the arrays had a 6 mm aperture radius, λ-size square elements and are focused at 25 mm depth. The results show that using the multi-depth cost function significantly improves the sensitivity (+1.2dB) and the SLL (−9.7 dB and −10.7 dB at 15mm and 35mm depths, respectively). The optimization results were compared to the beam patterns produced by a 256-element spiral array and by a 2D dense array (2292 elements), respectively.