© J. Ma et al., Published by EDP Sciences, 2026

1 Introduction

Traffic noise has become the primary source of pollution affecting urban acoustic environment quality. Noise pollution has been identified by the World Health Organization as Europe's second-leading environmental health risk factor, and globally, there is growing consciousness of the health risks resulting from noise exposure [1]. The rapid pace of urbanization, coupled with the relentless expansion of transportation infrastructure, means that increasing amount of urban inhabitants are experiencing traffic noise levels beyond what is considered acceptable [2]. As key spaces where urban residents connect with nature and relieve stress, the acoustic environment quality of green spaces directly impacts users' physical recovery and psychological well-being [3]. Excessive exposure to road traffic noise not only undermines the restorative function of urban parks but also significantly reduces visitors' willingness to visit.

The field of interior architecture has seen significant progress in digital tools for tackling outdoor noise, yet similar advancements in landscape framework have been lacking. The use of topographical modifications to reduce outdoor noise holds promise as an alternative to traditional noise barrier solutions, offering noise reduction while preserving recreational value. Compared to man-made structures, topographical noise reduction has inherent advantages in terms of visual integration, ecological compatibility, and maintenance costs. However, research on topographical noise reduction faces two key challenges. First, the high cost of acoustic simulation calculations limits the systematic exploration of the parameter space. Although traditional acoustic simulation methods, such as the Boundary Element Method (BEM), can accurately solve sound propagation problems in complex terrain, they require enormous computational resources; even when simplified to two dimensions, their practical applicability remains severely limited [4].

Second, noise-reduction design through topography involves conflicting objectives and trade-offs. Striving for maximum noise reduction often requires extensive terrain modification, which not only significantly increases earthwork volumes but may also compromise the visual openness of the landscape [5]. Achieving a Pareto optimal balance among noise reduction, earthwork costs, and visual quality is a core issue that current research urgently needs to address.

In response to these challenges, scholars both domestically and internationally have conducted a series of exploratory studies. However, existing research still has significant shortcomings. There is a lack of an efficient mechanism linking high-fidelity acoustic simulation with terrain parameter optimization. The application of surrogate models in the field of landscape acoustics remains unexplored. Furthermore, a systematic optimization framework for terrain parameters under multi-objective constraints has yet to be established. This paper proposes, as a solution to these problems, an optimization method for acoustic simulation parameters of landscape terrain aimed at achieving noise reduction. The key benefits of this method may be outlined as follows.

• Develop a hybrid framework combining high-fidelity BEM simulations with Kriging surrogate models. By using Latin hypercube sampling to design a small number of high-precision BEM samples, we train a Kriging response surface model to replace time-consuming physical simulations, reducing the evaluation time for a single iteration in the optimization process from several hours to milliseconds.

• Establish a mathematical model for multi-objective optimization of topographic parameters. A constrained multi-objective optimization problem is formulated by introducing key parameters such as terrain slope, elevation, curvature, and vegetation impedance as decision variables, while pursuing the objectives of maximizing noise reduction, minimizing earthwork amount, and maximizing visual permeability.

• An enhanced multi-objective particle swarm optimization algorithm incorporating k-nearest neighbors, termed KMOPSO, was proposed and empirically shown to be superior. Addressing the limitation of traditional MOPSO, which is prone to premature convergence, this research put forward a guidance strategy that leverages KNN distance, normal adaptive weights, and an elite gene coverage scheme. Theoretical analysis demonstrated the algorithm's convergence and diversity. Experimental results indicate that KMOPSO significantly outperforms benchmark methods in terms of metrics such as Inverted Generational Distance (IGD) and Hypervolume (HV).

The following sections of this article are organized as described below. Section 2 provides an in-depth review of the current research status and highlights existing shortcomings. Section 3 presents the optimization of simulation parameters for the acoustic characteristics of garden space terrain under noise reduction guidance. Section 4 presents the specific design of the experiment. Section 5 evaluates the outcomes of the experiments. Section 6 offers an interpretation of the outcomes, a consideration of the study's restrictions, and guidance for future efforts. Section 7 summarizes the entire paper.

2 Literature review

2.1 Methods for optimizing acoustic parameters of landscape topography

Bar-Sinai et al. [6] provided a systematic review of research progress on reducing surface traffic noise using natural methods such as vegetation surfaces, gabions, earthen embankments, and ground-related effects. The benefit of these solutions is their seamless integration with landscape design, which not only reduces noise but also increases the aesthetic value of the site. Through a case study near Munich Airport, Chen et al. [7] presented a digital workflow for acoustic surface design and proved that terrain modification can feasibly reduce noise levels. Using approaches such as acoustic simulation and correction, spatial analysis, and statistics, Zhang et al. [8] explored how the morphology of urban green spaces influences regional environmental noise. Zhou et al. [9] analyzed field recordings from 67 urban parks and demonstrated that the density of understory vegetation has a positive impact on the loudness and diversity of bird songs in the park soundscape. Bai et al. [10] introduced fountain water features in urban parks to strike a balance between meeting users' acoustic preferences and reducing the audibility of traffic noise, thereby validating the effectiveness of water features as a means of enhancing the soundscape.

In the field of BEM acoustic simulation, Liu et al. [11] proposed a novel boundary element method capable of simulating outdoor sound propagation over uneven terrain. Similar wave-based approaches have been successfully applied to complex mountainous terrains, where both topography and meteorological conditions strongly influence sound propagation patterns [12]. Furthermore, experimental validations in valley-slope configurations have demonstrated the reliability of numerical models under inhomogeneous atmospheric conditions [13].

Accurate computation of the sound field under intricate boundary geometries and refraction conditions was realized by representing the total sound pressure as the sum of the incident pressure and the pressure scattered by obstacles, while also embedding meteorological and ground effects within the Green's functions. Kumar et al. [14] developed a two-dimensional BEM model for non-flat terrain. By incorporating an infinitely impedance plane into the Green's function, the boundary element model requires discretization only of the non-flat terrain sections. Liu et al. [15] used the BEM model to validate the effectiveness of outdoor sound propagation models under complex terrain configurations, such as those involving a combination of valleys and barriers.

2.2 Multi-objective optimization algorithms

Multi-objective improvement approaches are effective tools for solving problems involving multiple constraints and conflicting objectives. What lies at their core is the search for Pareto-optimal solutions in the presence of multiple conflicting goals. Relative to single-objective optimization algorithms, they are more responsive to the needs of engineering design in practice. At this time, multi-objective optimization algorithms applied to landscape design and acoustic optimization primarily feature the Multi-Objective Particle Swarm Optimization (MOPSO) [16], the Non-Dominance Sorting Genetic Algorithm (NSGA-II) [17], and Multi-Objective Simulated Annealing (MOSA) [18]. Yu et al. [19] combined a Kriging surrogate model with a multi-objective cooperative evolution algorithm (MOCA) to perform robust optimization of cavity geometry and sound absorption coefficients as stochastic design variables under vibroacoustic criteria. Liravi et al. [20] observed that multi-objective optimization provides a set of trade-off-driven Pareto-optimal solutions, which gives a wider variety of design choices for decision-making compared to a single optimum.

The MOPSO algorithm is capable of efficiently finding global optimal solutions within the search space; it possesses strong global search capabilities and is frequently adopted to tackle multi-objective intelligent improvement issues. In Xu and Diao [21], a radial basis function neural network was combined with the MOPSO algorithm. A pre-trained radial basis function neural network served as the objective function within MOPSO to optimize landscape acoustic simulation parameters. Liu et al. [22] integrated a tabu table into the MOPSO algorithm to incorporate a tabu search mechanism, effectively preventing particles from repeatedly searching the same local regions and improving the algorithm's entire seek capability. Wang et al. [23] proposed a MOPSO algorithm that incorporates a competition mechanism. This algorithm selects elite particles based on dominance relationships and then updates the particles' velocities and positions through competition among these elite particles. Zhang et al. [24] grouped the population and introduced a dynamic mutation rate. They boosted the algorithm's global search performance via mutation-based updates of particle positions and velocities. Li et al. [25] proposed a crowding distance-based MOPSO approach. Using affinity propagation clustering along with a new crowding distance formula, they adjusted the crowding distance weights in the decision and objective spaces with greater precision, thus strengthening both convergence and diversity of the population. Feng et al. [26] applied MOPSO to the design of low-frequency noise suppression in landscape spatial acoustics. They achieved the dual objectives of reducing manufacturing costs while enhancing noise control effectiveness. The multi-objective optimization strategy they proposed is also widely applicable to other sound-absorbing materials.

To provide a more intuitive comparison of the applicability of different multi-objective optimization algorithms in optimizing acoustic parameters of landscape topography, Table 1 summarizes the strengths and limitations of these algorithms as described in representative literature.

2.3 Identified research gaps

Based on the above literature review, there are currently three major gaps in research on the optimization of acoustic parameters in landscape topography.

First, there is a lack of an efficient mechanism linking high-fidelity acoustic simulation with the optimization of terrain parameters. Although full-wave methods such as the BEM can accurately predict sound propagation in complex terrain, their enormous computational cost severely limits the systematic exploration of the parameter space. Most existing studies treat acoustic simulation as a tool for validating design proposals rather than as an optimization engine, resulting in terrain parameter adjustments that rely on empirical trial and error, making it difficult to achieve refined designs.

Second, the application of surrogate models in the field of landscape acoustics remains largely unexplored. Kriging surrogate models have already been successfully employed in fields such as materials acoustics and vibration acoustics to replace time-consuming simulations. However, there have been no reported studies on applying these models to the optimization of topographic acoustic parameters in landscape spaces. Given the highly nonlinear and multi-peaked complex mapping relationship between topographic parameters and acoustic responses, developing accurate surrogate models represents a key technical approach to overcoming computational bottlenecks.

Finally, a framework for optimizing terrain parameter systems under multi-objective constraints has yet to be established. Existing studies, while focusing on noise reduction, rarely take into systematic account engineering and landscape constraints such as earthwork volumes and visual permeability. The few studies that address multi-objective optimization have primarily focused on optimizing the shape or material parameters of noise barriers, rather than terrain modification.

There is an inherent trade-off among noise reduction, earthwork costs, and visual quality, and there is currently no unified optimization framework to balance these mutually constraining design objectives. How to efficiently search for a set of Pareto-optimal solutions within this multidimensional space of objectives is a core issue that urgently needs to be addressed in current research.

The proposed KMOPSO differs from existing MOPSO variants and surrogate-assisted methods in three aspects. Guidance mechanism: Unlike standard MOPSO that uses a single optimal particle, KMOPSO adopts a K-nearest neighbor distance strategy to select guiding particles from elite set, preserving Pareto front diversity. Adaptive perturbation: Instead of deterministic linearly decaying inertia weight, KMOPSO injects normally distributed random weights, enabling escape from local Pareto subsets. Surrogate integration. While prior work uses Kriging for single-objective or co-evolutionary frameworks, KMOPSO tightly couples a Kriging surrogate with the proposed multi-objective optimizer, exploiting predictive uncertainty to guide search.

3 Methodology

3.1 Overview of the proposed framework

The noise reduction-oriented terrain acoustic simulation parameter optimization method proposed in this paper adopts a three-layer progressive architecture composed of the high-fidelity physical simulation layer, the surrogate model layer, and the multi-objective optimization layer, forming a closed-loop technical process from data generation to model substitution and then global optimization. Figure 1 illustrates the data flow and control relationships within this architecture. The functional positioning of each layer is explained in detail below. The definitions of the main parameters involved in the proposed method are shown in Table 2.

High-fidelity acoustic simulation layer. This layer uses the BEM [27] as its core engine to establish a three-dimensional acoustic simulation model for garden spaces. Using standard terrain parameter sets (slope, elevation, curvature, and vegetation impedance coefficient), it solves the Helmholtz equation, thereby deriving the spatial variation of sound pressure levels from traffic noise on complex terrains for different frequency bands. The simulation results constitute the benchmark dataset for training subsequent surrogate models; their accuracy determines the reliability of the entire optimization framework.

Surrogate model layer. To overcome the high computational cost of BEM simulations and their difficulty in being directly embedded into optimization iterations, this layer employs Latin Hypercube Sampling (LHS) [28] to generate a limited number of highly representative sample points within the design space. Their response values are calculated using BEM. Based on these results, Kriging response surface surrogates are constructed to model the nonlinear correspondence between terrain parameters and acoustic responses. In milliseconds, the trained surrogate models produce predictions of noise reduction for any parameter combination, circumventing time-consuming physical simulations and furnishing the upper optimization layer with an efficient means of fitness evaluation.

Multi-objective optimization layer. Using Kriging surrogate models for fitness assessment, this layer applies the MOPSO to perform Pareto optimization with respect to three contradictory objectives—namely, maximizing noise reduction, minimizing earthwork volume, and maximizing visual transparency. Engineering constraints such as terrain slope and elevation variation range are also considered. Improvement results deliver a range of Pareto frontier solutions, presenting designers with trade-off options that are non-dominated and can be flexibly matched to distinct project budgets and landscape specifications.

Fig. 1 The proposed framework.

3.2 Acoustic simulation of landscape spaces using the boundary element method

BEM provides an efficient numerical means of solving acoustic wave issues. Whereas the finite element method requires discretization of the full computational domain, the core principle of BEM is to reduce a three-dimensional acoustic problem to a boundary integral equation, which calls for discretizing only the boundary of the region of interest. This dimensionality reduction significantly reduces the number of computational degrees of freedom. More importantly, for open semi-infinite domain problems such as garden spaces, BEM can naturally satisfy far-field radiation conditions without the need to set up absorbing boundary layers or perfectly matched layers as required by the finite element method, thus offering unique advantages. This section will elaborate on five aspects: governing equations, boundary integral forms, multi-domain coupling modeling for non-planar terrain, numerical discretization, and simulation parameter settings.

3.2.1 Acoustic control equations and boundary integral formulations

Assuming the sound field is a time-harmonic field, with the temporal factor taken as e^−jωt, where j is the imaginary element, ω = 2πf is the angular frequency, f is the frequency, and t is time, then the sound pressure $p(r,t)=\phi (r)e^{-j\omega t}$. The spatial complex amplitude φ(r) satisfies the Helmholtz equation, shown below.

$${\nabla }^2\phi (r)+k^2\phi (r)=0$$(1)

where ∇² is the Laplace operator, k = ω/c is the wavenumber, c is the speed of sound in air, r ∈ Ω, Ω is the acoustic field domain, referring to the semi-infinite air domain above garden spaces here.

For acoustic scattering problems in infinite or semi-infinite domains, the direct boundary element method is adopted, introducing a boundary integral equation. Let S be the boundary of the computational domain. The following integral equation holds for the sound pressure at any point P located on the boundary. By integrating the sound pressures and their normal derivatives over the boundary, it obtains the sound pressure at any point inside the domain.

$$c(P)\phi (P)={\int }_S(\phi (Q)\frac{\partial G(P,Q)}{\partial n_Q}-G(P,Q)\frac{\partial \phi (Q)}{\partial n_Q})dS(Q)+{\phi }_{\text{inc}}(P)$$(2)

where φ(P) is the total sound pressure at boundary point P; φ(Q) is the total sound pressure at another source point Q on the boundary; ∂φ(Q)/∂n_Q is the sound pressure derivative in the outward normal direction at point Q, G(P,Q) is the free-space Green's operation, representing the sound pressure produced through a unit-strength point source located at Q and measured at P. $R=|r_P-r_Q|$ is the distance between P and Q; $\partial G(P,Q)/\partial n_Q$ is the directional derivative of the Green's function in the outward normal direction at point Q. φ_inc(P) represents the incident sound pressure, directly generated by traffic noise sources. c(P) is a coefficient related to boundary geometry. If P lies on a smooth boundary and the outward normal is continuous, then c(P) = 1/2. If P lies within the domain, then c(P) = 1. If point P is located at a corner or along an edge, it needs to be calculated according to the solid angle at that point.

3.2.2 Multiphysics coupling modeling

The terrain of garden spaces is often not simply flat ground, but contains complex features such as rolling hills, depressions, water surfaces, and hard paving. To accurately simulate the reflection, diffraction, and scattering of sound waves on these terrains, directly discretizing the entire terrain boundary requires a huge computational effort. For this reason, this paper adopts the multi-domain coupling strategy proposed by Defrance and Gabillet [29]. Its key concept is to integrate the effect of an infinite impedance plane into the fundamental solution of the Green's function via the image method, allowing numerical discretization to concentrate solely on non-flat terrain features, thereby substantially lowering computational costs. For general terrain, the total sound field can be decomposed into three parts, where φ_inc is the direct sound field from traffic noise sources, φ_ref is the reflected sound field produced by a flat reference ground surface, and φ_scat is the additional scattered sound field caused by terrain undulations.

$$\phi ={\phi }_{inc}+{\phi }_{ref}+{\phi }_{scat}$$(3)

Substituting equation (3) into the boundary condition yields that, on the non-flat terrain surface S_irr, the total sound pressure obeys the impedance boundary condition.

$$\partial \phi /\partial n=jk\beta \phi$$(4)

where $\beta ={\rho }_{0}c/Z$ is the normalized admittance, and Z is the ground acoustic impedance. For commonly vegetated ground surfaces, the impedance Z can be calculated using the Delany-Bazley empirical model. Substituting equation (3) into equation (4), we obtain the boundary condition that the scattered field must satisfy, as shown in equation (5). The scattered field's normal derivative on the boundary is specified. Equation (5) together with the radiation condition constitutes a well-posed problem for φ_scat.

$$\left(\partial {\phi }_{scat}/\partial n\right)=-(\partial ({\phi }_{inc}+{\phi }_{ref}\left)\right)/\partial n+jk\beta \phi$$(5)

3.2.3 Discretization and numerical solutions

To solve the integral equation (2), it is necessary to discretize the continuous boundary integral equation into an algebraic equation system. This article adopts the constant element boundary element method, dividing the terrain surface into many small elements, and assuming that sound pressure and normal velocity are constants on each element. The specific steps are as follows.

• Boundary discretization. Divide the terrain surface S into N triangular or quadrilateral elements. For garden spaces, the typical element size should be less than 1/6 of the shortest wavelength to ensure calculation accuracy.

• Integral equation discretization. Apply equation (2) to the center point P_i(i = 1,…,N) of each element and assume thatand ∂φ/∂n are constants on element e_j. The integral is then approximated as the sum of integrals over individual elements.

• Matrix assembly and solution. Combine all boundary element discrete equations to form a linear algebraic equation system Hp = Gq + p_inc, where H and G are coefficient matrices containing the integration of the Green's function and its normal derivative, p and q represent vectors of sound pressure and normal velocity at boundary nodes, respectively. Combining with impedance boundary conditions, express q as a function of p, substitute it into the original equation system and rearrange to obtain Ap = b. The complex linear system is solved using the iterative method to obtain the scattered sound pressure distribution on undulating terrain surfaces.

• After obtaining the unknown boundary quantities, use the Helmholtz integral formula to back-calculate the total sound pressure p(R) at any receiving point R in space. Afterwards, determine the A-weighted sound pressure level L_Aeq. Extract the sound pressure level attenuation ΔL of specific evaluation points or regional grids as acoustic performance indicators. To ensure computational stability, a semi-analytical approach combining analytical and numerical integration is adopted for singular integrals, and discrete accuracy is verified via mesh convergence tests. When the change in receiving point sound pressure level is less than 0.5 dB after adjacent mesh refinement, it is judged that the mesh density meets engineering calculation requirements.

The following parameters are fixed across all BEM simulations: frequency range from 63 Hz to 4000 Hz in 1/3-octave bands (center frequencies: 63, 80, 100, 125, 160, 200, 250, 315, 400, 500, 630, 800, 1000, 1250, 1600, 2000, 2500, 3150, 4000 Hz). The A-weighted sound pressure level is computed using standard weighting factors. The terrain surface is discretized into constant triangular elements with maximum edge length ${\lambda }_{\text{min}}/6=0.0142\text{m}$ (where λ_min = c / f_max = 340/4000 = 0.85 m). The iterative solver tolerance for the linear system $Ap=b$ is set to 10⁻⁶. The impedance boundary uses the Delany-Bazley model for vegetated ground with flow resistivity σ  = 20000 kPa∙s/m² and porosity . The incident field is modeled as a line source at 1.5 m height, representing a highway traffic noise source with sound power level 85 dB per meter.

3.3 Kriging response surface models based on latin hypercube sampling

Although BEM simulations are highly accurate, they are time-consuming for single calculations. For a model containing thousands of elements, one frequency sweep may require several hours or even dozens of minutes. Multi-objective optimization algorithms usually need thousands to tens of thousands of fitness evaluations; if the BEM is used directly, calculation times can extend to months, which is unacceptable in engineering. Therefore, it is essential to construct a response surface model. This uses a small number of high-precision simulation samples as training data to build a fast approximation mapping from terrain parameters to noise reduction responses. After training, the single-prediction time for the surrogate model requires only milliseconds, making large-scale optimization possible. This section employs Kriging models as the surrogate model. Kriging is an interpolation method based on Gaussian processes. Compared to polynomial response surfaces, Kriging can accurately fit highly nonlinear functions and is especially suitable for the complex, non-convex relationship between terrain parameters and acoustic responses discussed in this paper.

Based on existing research analysis, it is known that the main slope θ, relative elevation h, terrain curvature k, vegetation normalized impedance real part $Z_{veg}$, and vegetation impedance frequency-dependent exponent $f_{veg}$ have significant effects on acoustic performance. Vegetation impedance employs a simplified single-parameter model $Z=Z_{\text{veg}}{(k/k_0)}^{-f_{\text{veg}}}$, $k_0=2\pi f_0/c$. This model can capture the acoustic characteristics of typical vegetation. The response variable is defined as noise reduction performance indicators, as shown below.

$$y=\Delta L_{avg}=\frac{1}{M}\sum _{m=1}^M\Delta L\left(r_m\right).$$(6)

The average reduction in A-weighted sound pressure level (dB) at all receiving points within the garden space. The higher this value, the better the overall noise reduction effect. This approach simplifies both the input and output dimensions of the surrogate model while avoiding introducing additional uncertainties.

To train the Kriging model, several sample points need to be generated within a 5-dimensional design space, and the BEM response value y corresponding to each sample point should be calculated. The selection of sample points is critical. If a simple full factorial design is used, the number of samples increases exponentially with dimensions. If random sampling is adopted, it may cause clustering of samples in certain areas while missing other regions. As a space-filling approach, LHS can achieve even coverage of the full parameter space with fewer samples. Each interval for each variable is exactly covered by one sample point; thus, on any projection along a single dimension, the sample points are uniformly distributed. Meanwhile, because the permutation is randomly generated, samples also have good space-filling properties in high-dimensional spaces.

After generating the sample points, BEM simulations $x^j$ are run for each sample point to obtain response values $y^{\left(j\right)}=\Delta L_{avg}\left(x^{\left(j\right)}\right)$, forming a training set D = {(x(j), $y^{\left(j\right)}){\}}_{j=1}^N$. The Kriging model treats an unknown function $f\left(x\right)$ as a Gaussian random process. The essential form of the model is presented below.

$$f\left(x\right)=\mu +\epsilon \left(x\right)$$(7)

where μ represents the global trend term, indicating the average value of the function across the entire design space. $\epsilon \left(x\right)$ is a zero-mean stationary Gaussian random process, and its covariance function characterizes the correlation between responses at different points, as shown below, where ${\sigma }^2$ is the process variance, R is the correlation function, and θ is the correlation length hyperparameter.

$$Cov\left[\epsilon \right(x^{\left(i\right)}),\epsilon (x^{\left(j\right)}\left)\right]={\sigma }^{2}R(x^{\left(i\right)},x^{\left(j\right)};\theta ).$$(8)

Correlation function R determines the smoothness of the response values with respect to design variable changes. This paper adopts a Gaussian-type correlation function as follows.

$$R(x^{\left(i\right)},x^{\left(j\right)};\theta )=exp(-\sum _{k=1}^d{\theta }_k{(x_k^{\left(i\right)}-x_k^{\left(j\right)})}^2)$$(9)

Parameter ${\theta }_k>0$ controls the sensitivity of the k th variable to the response. The larger ${\theta }_k$ is, the more rapidly the correlation decreases due to small changes in this variable, indicating that the response is highly sensitive to it. Conversely, if ${\theta }_k$ is very small, the response varies smoothly along this dimension.

After completing the BEM simulations of training samples and extracting response values, hyperparameters of the Kriging model must be determined based on maximum likelihood estimation. The concentrated log-likelihood function is jointly optimized using genetic algorithms and quasi-Newton methods to obtain optimal θ, enabling fast prediction for any new parameter combination while quantifying predictive uncertainty. Subsequently, an independent Latin hypercube validation set is used to evaluate the model accuracy. When the accuracy proves insufficient, or when the predicted variance of solutions on the Pareto front in the following multi-objective optimization goes above a threshold, high-uncertainty points are appended to the training set, prompting retraining of the Kriging model. This continuously improves the local accuracy of the surrogate model near optimal regions within limited BEM computation budgets, providing reliable and efficient fitness evaluation functions for the MOPSO algorithm.

3.4 Optimization of acoustic simulation parameters for landscape topography

3.4.1 Mathematical model

To seek Pareto optimality among noise reduction effect, earthwork volume and visual transparency, the optimization problem related to garden-space terrain acoustic simulation parameters is presented in this work as a constrained multi-objective optimization problem. Listed below are the decision variables, objective functions, and constraints. According to the above analysis, the key parameters affecting the acoustic performance of garden space terrain include terrain geometric parameters and ground vegetation acoustic parameters. Denote the decision vector as $x=[\theta ,h,k,$ $Z_{veg},f_{veg}{]}^T$, where $0^{{̂}^{}}\le \theta \le {30}^{{}^{̂}}$,$0m\le h\le 5m$,$-0.5\le k\le 0.5$,$0.2\le Z_{veg}\le 2$, $0.5\le f_{veg}\le 1$.

This optimization problem includes three conflicting objective functions, corresponding to noise reduction effect, earthwork volume and visual transparency. Noise reduction effect objective$f_1\left(x\right)$. This objective uses the noise reduction performance indicatoras a measure of noise reduction. The higher this value is, the more significant the noise reduction. In multi-objective optimization, all objectives are usually unified into minimization form, and $f_1\left(x\right)$ is defined as follows, $\Delta L\left(r_m\right)$ stands for the reduction in A-weighted sound pressure level observed at receiving point $r_m$, M is the total number of receiving points.

$$f_1\left(x\right)=-\Delta L_{avg}\left(x\right)=-\frac{1}{M}\sum _{m=1}^M\Delta L\left(r_m\right)$$(10)

Earthwork volume objective$f_2\left(x\right)$. The earthwork volume caused by terrain modification directly determines engineering costs and ecological disturbance. This paper adopts the absolute value of the volume change before and after terrain modification as a measure. For a given terrain surface $z=\zeta (x,y)$ and benchmark plane $z=0$, the earthwork volume V can be approximately calculated by integrating the absolute values of terrain elevation, as shown below.

$$V\left(x\right)={\int }_A\left|\zeta \right(x,y;x)-{\zeta }_0(x,y\left)\right|Adxdy$$(11)

where A is the horizontal projection area of the garden space, and ${\zeta }_0(x,y)$ is the original terrain elevation. Since the original terrain can be considered as known during optimization, $\zeta (x,y;x)$ is uniquely determined by slope θ, elevation h, and curvature k through a parametric terrain generation function. To simplify calculations and maintain compatibility with the surrogate model, an empirical formula based on decision variables is adopted.

$$f_2\left(x\right)=V\left(x\right)\approx {\alpha }_{1}h+{\alpha }_2\left|k\right|+{\alpha }_3\theta$$(12)

where ${\alpha }_1$, ${\alpha }_2$, ${\alpha }_3$ are positive coefficients fitted by finite volume element computation. Their specific values depend on site area and unit conversion. This approximate form reflects the positive contribution of elevation change, terrain curvature, and slope to earthwork volume.

Visual transparency objective$f_3\left(x\right)$. Several observation points are evenly distributed within the site, and the proportion of visible site area for each observation point is calculated. The average value over the entire site is taken as a measure of visual transparency T(x). Higher transparency indicates better landscape quality. $f_3\left(x\right)$ is defined as follows.

$$f_3\left(x\right)=1-T\left(x\right)=1-\frac{1}{N_{obs}}\sum _{i=1}^{N_{obs}}\frac{A_{visible}^{\left(i\right)}}{A_{total}}$$(13)

where N_obs is the number of observation points, $A_{visible}^{\left(i\right)}$ is the visible area from the i th observation point, and $A_{total}$ is the total site area. $f_3\left(x\right)\in [0,1]$, lower values indicate better visual transparency. In summary, the multi-objective optimization problem in this study can be expressed as follows.

$${min}_{x\in \chi }F\left(x\right)={\left[f_1\right(x),f_2(x),f_3(x\left)\right]}^T.$$(14)

3.4.2 Optimizing the solution of mathematical models

To address the drawback of traditional multi-objective particle swarm algorithms being easily guided by a single optimal particle, this paper proposes a MOPSO algorithm based on KNN distance. The KMOPSO algorithm is used to solve the optimal parameters for landscape space terrain acoustic simulation. First, by calculating the distance relationship between candidate particles and K elite particles, a competitive optimal particle selection strategy is designed to replace the traditional global optimal selection method, significantly enhancing the spatial distribution diversity of the particle swarm. Second, an adaptive weight control strategy based on normal distribution is constructed to dynamically adjust the scale of the population learning from guiding particles. Finally, superior particle gene random coverage operations are introduced to enhance the exploration ability of the population. The process of optimizing and solving the simulation parameters for the acoustic conditions of the garden space is shown in Figure 2.

The KNN algorithm evaluates similarity between two samples by calculating their Euclidean distance; shorter distances indicate closer proximity between the two samples. For N-dimensional vectors X and Y, their Euclidean distance d is expressed as follows.

$$d(X,Y)=\sqrt{\sum _{i=1}^N{\left(X_i-Y_i\right)}^2}$$(15)

where X_i andare the Y i th elements of an N-dimensional space vector.

Inspired by the KNN idea, this paper proposes a K-nearest distance learning strategy (K-ND) to select leading particles. In K-ND, an elite particle set is used to provide candidate competitive particles to guide the population's search. To select elite particles that can maintain a balance between convergence and diversity, this paper uses non-dominated sorting and crowding distance-based sorting from NSGAII to choose the elite particles. For every competition among elite particles, first K elite particles are randomly selected from the elite particle set, and the distance between each elite particle and a specific particle in the current population is calculated according to equation (15). By comparing the distances of each pair of particles, the closest particle to the specified one is chosen as the best particle. Once the global optimal particle has been identified, the position and velocity of the current particle p can be updated by learning from the optimal particle.

Let $v_i$ and $p_i$ denote the speed and location of the i th particle in the population. $1\le i\le N$, for the i th particle, the updated velocity $v_i^{\text{'}}$ and updated position $p_i^{\text{'}}$ are calculated according to the following equations.

$$v_i^{\prime }=w{v}_i+r_1(p_w-p_i)$$(16)

$$p_i^{\prime }=p_i+v_i^{\prime }\times H(\epsilon -lr)$$(17)

where $r_1\in [0,1]$ is a randomly generated vector; l_r is the learning rate; w is the inertia weight, is a group of vectors randomly generated from a normal distribution in which the mean value is controlled by the weight . The values for w_r are calculated according to the following equations, where w_max denotes the maximum weight and w_min the minimum weight.

$$w_r=w_{max}-\frac{(w_{max}-w_{min})iter}{maxiter}.$$(18)

During the iterations of a multi-objective optimization algorithm, an external archive set is utilized to record the non-dominated solutions that have been found. Due to the limited size of the archive set, its convergence and diversity must be guaranteed simultaneously. The following strategy is adopted in this paper. After each generation iteration, the current population is merged with all existing solutions in the archive set for non-dominated sorting. Non-dominated level 1 solutions are retained as new archive candidates. If the size of the archive exceeds the preset upper limit N_arch, each solution's crowding distance in the objective space is obtained. Defined as the cumulative Euclidean distance between adjacent solutions on the Pareto front, crowding distance indicates the local solution density. The calculation formula is as follows.

$$C{D}_j=\sum _{m=1}^M\frac{f_m^{next}-f_m^{prev}}{f_m^{max}-f_m^{min}}$$(19)

where M is the amount of targets, $f_m^{next}$ and $f_m^{prev}$ are the predecessor and successor objective values for the j th solution at the m th objective after sorting in ascending order. Until the archive size reaches its maximum allowed capacity, solutions characterized by the minimal crowding distance are eliminated. This strategy prioritizes solutions located in sparse regions of the front, thus maintaining population diversity. By integrating the above mechanisms, the complete pseudocode for KMOPSO algorithm to solve optimal acoustic simulation parameters for garden space terrain is shown in Table 3.

Fig. 2 The process of optimizing and solving the simulation parameters.

3.5 Theoretical analysis

Standard MOPSO is prone to premature convergence caused by attraction toward a single global optimal particle. The KNN competition guidance, normal adaptive weighting and elite gene coverage operations introduced in this paper improve the population's ergodicity and Pareto front coverage from the perspective of stochastic processes. KNN-guided local manifold tracking. Let the guiding particle $p_{g,i}^{\left(t\right)}$ for particle i at generation t be selected from external archive A using k th Nearest Neighbor Distance (K-ND) strategy. If non-dominated solutions in A are considered as discrete sampling points on true Pareto front, then KNN strategy essentially allows particles to move along local tangent directions of the Pareto manifold. The expected velocity update can be expressed as follows.

$$E\left[v_i^{(t+1)}\right]=wE\left[v_i^{\left(t\right)}\right]+c_{1}E[p_{w,i}^{\left(t\right)}-p_i^{\left(t\right)}]+c_{2}E[p_{g,i}^{\left(t\right)}-p_i^{\left(t\right)}].$$(20)

Since $p_{g,i}^{\left(t\right)}$ is always a non-dominated solution in neighborhood, its expected gradient direction is proportional to the local Pareto density. Theoretical analysis shows that under assumptions where objective functions satisfy Lipschitz continuity and decision space is bounded, this local guidance strategy enables swarm dispersion radius to decay at rate $O(1/\sqrt{t})$ while maintaining monotonic non-decreasing property of hypervolume indicator, thereby achieving rigorous trade-off between convergence speed and distribution uniformity.

Ergodicity enhancement with normal perturbation and elite coverage. Inertial weight perturbation $w~N(w_r,{\sigma }_w^2)$ injects controllable randomness into velocity updates. When the population falls into a local Pareto subset, variance components provided by σ_w allow particles to escape attraction domains with probability greater than 0. Combined with elite gene coverage operations, the algorithm introduces exploration operators similar to crossover and mutation in decision space. Markov chain theory tells us that when the state transition matrices are irreducible and positive recurrent, the algorithm converges asymptotically to the global Pareto optimal solution set with probability one.

4 Experimental design

4.1 Dataset

The dataset collected in reference [30] is used as field data source in this paper. Nationwide field surveys within the RESTORE project framework enabled the Swiss Federal Institute for Forest, Snow and Landscape Research (WSL) to build this dataset., recording road traffic noise exposure levels ($L_{den}$ ), green space landscape characteristics parameters and users' subjective evaluations of sound environments across multiple sample plots in Swiss green areas. Based on spatial locations of sample plots, elevation, slope and curvature at each site were extracted using 30 m resolution digital elevation models. Vegetation types were mapped to normalized impedance real parts $Z_{veg}$ and frequency-dependent exponents $f_{veg}$, constructing complete five-dimensional input vectors. Take $L_{den}$ as the reference sound pressure level and infer the measured equivalent noise reduction $\Delta L_{field}$ by comparing the difference of sound levels at the same location under conditions with and without terrain undulations. Combined with $\Delta L_{avg}$ simulated by BEM, they form the target responses for multi-fidelity modeling. The original dataset contains 120 valid sample points. After feature extraction and outlier removal, 100 samples are obtained. According to a stratified sampling strategy, the samples are divided into training set (70%), validation set (10%), and test set (20%). Among them, the training set is used for constructing Kriging surrogate models, the validation set for hyperparameter tuning, and the test set for final generalization performance evaluation.

The original dataset from [27] contains 120 sample points. The following preprocessing steps are applied. Outlier removal: Samples with standardized residual > 3 standard deviations from the regression line between BEM-simulated and measured noise reduction are removed. Feature scaling. All five input variables $(\theta ,h,k,Z_{\text{veg}\prime },f_{\text{veg}})$ are normalized to the range [0,1] using min-max scaling: $x_{\text{norm}}=(x-x_{\text{min}})/(x_{\text{max}}-x_{\text{min}})$. Stratified sampling: The remaining 115 samples are split into training (70%, 80 samples), validation (10%, 12 samples), and test (20%, 23 samples). Stratification is based on the noise reduction response quantiles to ensure distributional similarity across subsets. Handling missing values: No missing values exist in the final dataset. The preprocessed dataset and the train/validation/test split indices are provided as supplementary material.

4.2 Benchmark methods

For a comprehensive evaluation of the KMOPSO algorithm's performance when optimizing parameters for acoustic simulations within garden spatial terrain, five typical multi-objective optimization algorithms are selected as baselines for comparison.

NMGA [17]: Utilizes fast non-dominated sorting and crowding distance selection mechanism to approach the Pareto front while maintaining population diversity.

MOSA [18]: Controls solution set evolution through Metropolis acceptance criteria and Pareto dominance relations, achieving a transition from global exploration to local refinement during temperature descent.

MOCA [19]: Embeds Kriging surrogate models into a multi-objective coevolutionary framework for robust optimization design of cavity geometries and sound absorption coefficients under vibration acoustics criteria.

CMOPSO [24]: Combines Affinity Propagation clustering with a new crowding distance formula to dynamically adjust the crowding weights in decision and objective spaces, enhancing population convergence and diversity.

MOPSO [25]: Applies MOPSO to low-frequency noise suppression design in garden space acoustics. Achieves dual objectives of reducing manufacturing costs while improving noise control effectiveness.

4.3 Evaluation metrics

To quantitatively evaluate the quality of Pareto fronts obtained by each algorithm, this paper selects IGD, HV, spacing metric (SP), and Maximum Pareto Front Error (MPE) as evaluation metrics. IGD measures the average distance between the algorithm's front and the true Pareto front, comprehensively reflecting convergence and diversity; smaller values are better. HV measures the volume of the region enclosed by the algorithm's front and a reference point, with larger values indicating broader coverage in objective space and stronger convergence. SP evaluates the uniformity of solution distribution on the front by calculating the standard deviation of Euclidean distances between adjacent solutions. MPE measures the maximum distance between the algorithm's front and the true front. All metrics are calculated based on normalized objective values, and report average values and standard deviations from 30 independent runs.

4.4 Implementation details

The experiments are conducted on an Intel Xeon Gold 6248 workstation (20 cores, 256 GB RAM) with Ubuntu 20.04 system, using Python 3.8 combined with Scikit-learn, PlatEMO and MATLAB/OpenBEM for surrogate model training and optimization calculations. The Kriging model adopts a constant mean and Gaussian correlation function, and estimates hyperparameters through genetic algorithm BFGS joint estimation, using a training set of 80 samples and a validation set with coefficient of determination R²> 0.95 before being used in subsequent optimizations. All optimization algorithms are uniformly configured with population size 50, maximum iterations 200, and external archive upper limit of 100, where KMOPSO's dedicated parameters take K = 3, inertia weight is 0.9, learning factor is 2, and gene coverage probability is 0.1. The five benchmark algorithms are all configured according to the original literature recommendations. All the random algorithms were run independently for 30 times, with different random seeds (1 to 30) used each time. The mean and standard deviation of each indicator were calculated, and pairwise comparisons were conducted using the Wilcoxon rank sum test with a significance level of α  = 0.05.

5 Experimental results analysis

5.1 Multi-objective optimization performance comparison

To comprehensively evaluate the effectiveness of KMOPSO algorithm in parameter optimization for acoustic simulation of garden space topography, KMOPSO is compared with benchmark methods under identical experimental conditions through 30 runs each, calculating four metrics IGD, HV, SP and MPE. The results are shown in Table 4.

To assess the statistical reliability of the performance indicators of each algorithm, Table 4 reports the mean and standard deviation of 30 independent runs. It can be seen that the standard deviation of KMOPSO in each indicator is smaller than that of the comparison methods (for example, the standard deviation of IGD is 0.024, while that of NMGA is 0.042), indicating its better robustness. The Wilcoxon rank sum test results show that the differences between KMOPSO and all benchmark methods are statistically significant at the p < 0.05 level. For the ablation experiments in Table 6, the same 30 independent runs were adopted. The differences between each variant and the complete KMOPSO were all significant (p < 0.05), verifying the effectiveness of each improved module.

Compared to CMOPSO and MOPSO, KMOPSO achieves 39.9% lower IGD and 3.0% higher HV due to its KNN guidance and normal adaptive weights. Relative to surrogate-assisted MOCA, the proposed framework obtains better SP, demonstrating the advantage of combining Kriging with KNN-guided particle swarm search for landscape acoustic optimization.

To compare the convergence speed of different algorithms on test functions, with evaluation times as the x-axis and IGD as the y-axis, we display the convergence speeds of different algorithms on multi-objective optimization problems. Figure 3 shows the trend in the IGD metric for different algorithms.

The KMOPSO algorithm demonstrates good convergence performance: It converges within 2500 evaluations on multi-objective optimization functions and reaches the minimum IGD value. The MOPSO algorithm does not converge within the first 10,000 evaluations. Although the MOSA algorithm shows a convergence trend, its convergence speed and accuracy are inferior to that of KMOPSO. NMGA, MOCA, and CMOPSO algorithms have slower convergence speeds than KMOPSO but their final convergence values also reflect good performance for these three algorithms, which is consistent with the analysis of experimental results, thus confirming that the KMOPSO algorithm can achieve a better Pareto front with fewer evaluation times and is suitable for optimizing parameters in acoustic simulation of garden space topography.

Fig. 3 The trend in the IGD metric for different algorithms: (a) NMGA; (b) MOSA; (c) MOCA; (d) CMOPSO; (e) MOPSO; (f) KMOPSO.

5.2 Noise reduction results analysis

For evaluating how effectively the optimized scheme reduces noise in an actual garden setting, we choose from the KMOPSO-derived Pareto front the best solution for noise reduction, the best solution for earthwork, and a compromise solution. The BEM simulation receives the topography parameters of these three solutions to compute the spatial distribution of A-weighted sound pressure levels at 100 receiver locations within the site. The results are compared against a flat reference ground surface and summarized in Table 5. The optimal noise reduction solution achieves the highest average noise reduction of 5.82 dB but at the cost of high earthwork volume reaching up to 385.6 m³ and low visual transparency of only 0.31. The optimal earthwork solution requires merely 42.3 m³ of earthwork with a transparency level of 0.89, yet its noise reduction effect is only 1.53 dB. The balanced solution achieves a good balance among the three aspects, reducing noise by 3.71 dB, earthwork volume of 186.4 m³, and transparency level of 0.62. These results confirm that KMOPSO can provide a Pareto-optimal solution for parameter optimization in acoustic simulation of garden space topography.

Figure 4 visually reveals the spatial coupling mechanism between the terrain topological structure and the noise reduction effect through a 2D top-down perspective: By comparing the three optimization schemes, it can be seen that (a) the optimal noise reduction scheme constructs a single large uplift terrain (100 m × 80 m), forming a significant sound shadow area behind the sound source, resulting in a high-value noise reduction band of 4-6 dB (the dark red area), but at the cost of 385.6 m of earthwork; (b) the optimal earthwork scheme only performs minor terrain modification (80 m × 40 m), with the noise reduction amount mainly <2 dB (the blue area), although the earthwork volume is only 42.3 m, but the acoustic benefit is limited; (c) the balanced solution scheme innovatively adopts a double-peak terrain layout (100 m × 60 m), through "8" shape contour lines to construct a gradient noise reduction band, achieving a continuous noise reduction zone of 2-4 dB (the orange-yellow area) with 186.4 m³ of earthwork volume, verifying the effectiveness of the multi-objective optimization strategy proposed in this paper. That is, under the constraint of limited engineering cost, by optimizing the terrain shape rather than simply increasing the elevation, the Pareto optimal balance between noise reduction performance and visual transparency can be achieved, providing a quantifiable design basis for actual landscape engineering.

Figure 5 illustrates the distribution of sound pressure levels around the sound barrier before and after optimization using KMOPSO. One can see that shape optimization has greatly reduced the sound pressure level in the shadow region of the Y-shaped sound barrier. thus, the proposed KMOPSO method provides an effective means to improve the noise reduction performance of sound barriers.

Figure 6 presents the evaluation of acoustic scene signal-to-noise ratio (SNR) under various traffic sounds. For a fixed SNR, the noise reduction performance of the acoustic scene deteriorates as the sound pressure level of the traffic noise rises. Given the same traffic sound pressure level, the acoustic scene's noise reduction performance grows with increasing SNR. For a traffic sound level of 50 dB, the noise reduction performance of the acoustic scene under a 15 dB signal-to-noise ratio is markedly better than under 5 dB, 0 dB, or −5 dB. at traffic sound levels of 55 dB and 60 dB, the noise reduction effect of the acoustic scene with a 15 dB signal-to-noise ratio is significantly higher than that of 5 dB, 0 dB, and −5 dB, and the noise reduction effect of the acoustic scene with a 5 dB signal-to-noise ratio is significantly higher than that of −5 dB and 0 dB.

Fig. 4 2D top-down perspective visualization.

Fig. 5 The distribution of sound pressure levels around the noise barriers before and after the treatment: (a) Before optimization; (b) After optimization.

Fig. 6 Comparison of soundscapes noise reduction effects: (a) Soundscapes evaluation; (b) Sound scene noise reduction effect.

5.3 Ablation experiment

To verify the effectiveness of each module in the KMOPSO algorithm, four comparative experiments were designed: full KMOPSO, a variant without KNN guidance (w/o KNN), a variant without normal adaptive weights (w/o NAW), and a variant without elite gene coverage (w/o EGC). Each variant was independently run 30 times under identical experimental conditions, with IGD and HV as evaluation metrics. The results are shown in Table 6.

The three improved modules have positive contributions to the performance of the KMOPSO algorithm, with their impact decreasing gradually. Specifically, after removing KNN guidance, IGD increases from 0.385 to 0.493, and HV decreases from 0.801 to 0.764, showing the most significant performance degradation, indicating that the KNN competitive guidance strategy is a core mechanism for maintaining population diversity and front uniformity. After removing normal adaptive weights, IGD increases to 0.467, and HV decreases to 0.775, demonstrating that this strategy plays an important role in balancing global exploration and local exploitation. After removing elite gene coverage, IGD is 0.428 and HV is 0.788. Although the impact is relatively small, it still has statistical significance, proving that it can effectively avoid premature convergence as an auxiliary exploration mechanism. The complete KMOPSO significantly outperforms all variants, verifying the effectiveness of the synergistic action of the three modules.

5.4 Case study

To verify the engineering practicality of the KMOPSO method in actual garden space terrain noise reduction design, an urban park in Zurich, Switzerland is selected as the case site. The park covers an area of approximately 1.2 hectares and borders the A1 highway with a measured daytime equivalent traffic noise $L_{Aeq}$ of about 68 dB. The current condition of the site is flat grassland with sparse vegetation and prominent noise problems. The optimization objectives are set as: noise reduction effect $\Delta L_{avg}\ge 4\text{dB}$, earthwork volume $V_{earth}\le 200m^3$, and visual permeability $T_{visual}\le 0.5$. The KMOPSO algorithm is used for parameter optimization, yielding a Pareto front containing 76 non-dominated solutions. A total of 23 solutions that meet the constraint conditions are selected from them. The method employs fuzzy membership to calculate the aggregate satisfaction of each solution, and the solution with the highest satisfaction value is taken as the recommended design. This paper sets three alternative schemes for comparison. The first is flat terrain, maintaining the status quo without any modifications. The second is a traditional sound barrier. The third is random terrain. BEM simulations are used to calculate the noise reduction effect and earthwork volume of each scheme respectively, with the outcome implied in Table 7.

The flat terrain scheme has no earthwork cost and high permeability, but its noise reduction is only 0.8 dB, which cannot meet basic requirements. The traditional sound barrier scheme achieves the highest noise reduction, but with an earthwork volume of up to 420 m and visual permeability as low as 0.22, it has poor integration into garden landscapes and a comprehensive satisfaction of merely 0.61. The random terrain scheme yields unstable noise reduction effects that fail to achieve the expected goals. In contrast, the KMOPSO-optimized scheme achieves optimal balance among noise reduction, earthwork volume, and permeability, with a comprehensive satisfaction reaching 0.91, significantly outperforming other schemes.

6 Discussion

6.1 Interpretation of results

The KMOPSO algorithm performs excellently in three aspects: algorithm performance, acoustic physics, and engineering applicability. Algorithmically, its IGD, HV, and other indicators comprehensively outperform traditional methods, converging within 2500 evaluations, which is 40% faster than conventional MOPSO. The KNN guidance and normal adaptive weights are key improvements. Acoustically, noise reduction increases logarithmically with elevation and curvature, showing significant early gains in earthwork but diminishing marginal returns later on. By reducing surface reflection coefficients, a Pareto improvement between acoustics and visual aesthetics is achieved. Engineering-wise, the optimized scheme achieves 4.53 dB of noise reduction using only 195 m³ of earthwork, with comprehensive satisfaction far exceeding that of traditional sound barriers. The terrain naturally integrates into the site, validating its practicality and green advantages.

6.2 Limitations

Although the proposed framework has performed well in theoretical verification and case testing, its physical modeling and optimization process still rely on the following core assumptions and engineering simplifications, which need to be clearly defined when applied.

• Core physical assumptions. Steady-state harmonic sound field assumption. The BEM solution is based on time-harmonic excitation, discretizing the broadband traffic noise into multiple single-frequency components for superposition, ignoring the transient pulse noise of vehicles and the high-frequency nonlinear radiation characteristics of tire-road friction. Uniform medium and static meteorological assumption: The model assumes that the air density and sound speed are constant, and does not couple the wind shear, temperature gradient, and atmospheric turbulence effects on the refraction of long-distance sound propagation. Terrain-vegetation impedance decoupling assumption. The surface absorption is modeled using the Delany-Bazley empirical model, simplifying the vegetation layer into a uniform porous medium, without considering three-dimensional canopy scattering, root-soil coupling absorption, and the dynamic changes in plant moisture content.

• (2) Simulation - Measurement Deviation Mechanism. In open flat terrain and low vegetation coverage scenarios, the model prediction error is usually controlled within ± 0.8 dB; however, in complex micro-topography or high-vegetation density forest areas, the error may expand to 1.2–1.5 dB. This is mainly due to the insufficient resolution of the boundary element grid discretization scale for high-frequency diffraction, as well as the aliasing effect of non-sound source interference (such as bird calls, human voices) in the measured background noise. Additionally, the generalization ability of the Kriging surrogate model highly depends on the design of the LHS sampling space, and local overfitting may occur in the extreme combination regions of decision variables.

• (3) Engineering application boundary: This method is applicable to urban road traffic noise scenarios dominated by medium and low-frequency (125–1000 Hz) energy. In strong winds (>5 m/s), inversion layers, or complex weather conditions, it is recommended to introduce parabolic equation (PE) or ray-wave hybrid algorithms for multi-physics field correction. The target function for earthwork volume adopts linear fitting of geometric parameters. In actual engineering, it is necessary to conduct dynamic verification in combination with geological investigation reports, underground pipeline distribution, and compaction techniques. Visual transparency calculation is based on the intersection of discrete viewpoint rays. Factors such as dynamic eye scanning habits and night lighting occlusion are not included. In the future, multi-sensory coupling optimization can be achieved by combining eye movement tracking and field entropy model.

6.3 Future work

Addressing the limitations stated above, future investigations might be expanded in the following ways.

• Coupled meteorological-topographic multiphysics models. Introduce parabolic equation (PE) or ray-wave hybrid algorithms, superimpose wind temperature profiles and atmospheric attenuation coefficients to build an all-weather sound propagation simulation engine, enhancing predictive robustness under complex meteorological conditions.

• Dynamic vegetation acoustic impedance and 3D canopy acoustics. Combine Biot's porous medium theory with measured forest structure parameters to establish a seasonally adaptive impedance model. Introduce an equivalent volume scattering term into the BEM boundary conditions to achieve "topography-vegetation" collaborative soundscape optimization.

• Full lifecycle multi-objective expansion. Incorporate earthwork carbon footprint, construction energy consumption, vegetation maintenance costs and biodiversity indices into objective functions to establish an optimization framework integrating sound environment, ecology and economy that supports low-carbon garden design decisions.

7 Conclusion

Urban traffic noise significantly impacts the acoustic comfort of landscaped spaces. While terrain modification offers both visual integration and ecological benefits for natural noise reduction, existing research faces two major challenges: high-precision acoustic simulation is costly, making it difficult to support large-scale parameter exploration; and there is an inherent trade-off between noise reduction effectiveness, earthwork volume, and visual openness, with a lack of a unified multi-objective optimization framework.

In response to the aforementioned challenges, this study introduces a method oriented toward noise reduction for the optimization of acoustic simulation parameters of landscape topography. By constructing a Kriging surrogate model using Latin hypercubic sampling, the time required for a single acoustic evaluation is reduced from several hours to milliseconds. A K-MOPSO algorithm based on KNN distance guidance, normal adaptive weights, and elite gene coverage was designed, and its convergence and diversity were theoretically proven. Experimental results indicate that K-MOPSO significantly outperforms the baseline method in metrics such as IGD and HV. Case study results show that the optimized scheme achieves a noise reduction of 4.53 dB with 195 m of earthwork, a visual permeability of 0.58, and a comprehensive satisfaction score of 0.91, far surpassing traditional sound barriers and random terrain configurations. The method presented in this paper provides a scientific and quantifiable decision-making tool for green noise reduction design in landscape spaces. Future research will focus on multi-physics coupling and multi-objective optimization across the entire life cycle.

Funding

This research received no external funding.

Conflicts of interest

All authors declare that they have no conflicts of interest.

Data availability statement

This article has no associated data generated and/or analyzed.

Author contribution statement

Supervision, Writing—original draft, Methodology, Validation, Bihua Zou; Formal Analysis, Jiefeng Liu; Data Curation, Rivas-Rivera Vaneshca.

References

All Tables

All Figures

	Fig. 1 The proposed framework.
In the text

	Fig. 2 The process of optimizing and solving the simulation parameters.
In the text

	Fig. 3 The trend in the IGD metric for different algorithms: (a) NMGA; (b) MOSA; (c) MOCA; (d) CMOPSO; (e) MOPSO; (f) KMOPSO.
In the text

	Fig. 4 2D top-down perspective visualization.
In the text

	Fig. 5 The distribution of sound pressure levels around the noise barriers before and after the treatment: (a) Before optimization; (b) After optimization.
In the text

	Fig. 6 Comparison of soundscapes noise reduction effects: (a) Soundscapes evaluation; (b) Sound scene noise reduction effect.
In the text