Strutt Help

Harmonoise Excess Attenuation    1/1,1/3

Strutt|Environmental Noise|Harmonoise Excess Attenuation inserts the combined excess attenuation for a terrain profile into the active row of the worksheet calculated using the Harmonoise EU noise propagation model.

Unlike other environmental calculations in Strutt, Harmonoise excess attenuation is calculated as the cumulative excess attenuation for a full terrain profile (which may include diffraction edges such as noise barriers or natural terrain edges) - i.e. Harmonoise calculates the overall excess attenuation due to ground reflection, diffraction and atmospheric scattering/refraction. The user only has to add the distance loss and any atmospheric absorption loss.

The Harmonoise model is rather theoretically-complicated, but can be summarised as follows:

The Harmonoise model is described in detail in the various published papers describing the model (see References section at the end). The following sections are a summary of key features of the model.

Diffraction losses in Harmonoise are calculated using the Degout barrier formula, which is described in detail in the Barrier Attenuation calculation.

The path difference `delta` for diffraction is calculated automatically from the terrain as follows based on the magnitude of the diffraction angle `theta = theta_S+theta_R`

Ground attenuation is calculated using a complicated model, which accounts for several propagation factors.

The ground effect is calculated over each terrain segment, and a local coordinate system is used for each segment with the source/receiver heights `h_s` and `h_r` calculated by Strutt relative to each segment. These local source heights are used for image source/receiver calculations.

Each segment is characterised as either 'concave' (both local source height `h_s` and local receiver height `h_r` are positive), 'hull' (local source and receiver heights are both zero, e.g. segment `P_5-P_6` in the example above), or 'convex' (one or both of the local source/receiver heights `h_s`,`h_r` are negative).

Two overall models are used, the Concave model, which is used when diffraction is not significant, and the Transition model, which 'blends' ground reflection and illuminated-zone diffraction effects.

Concave Model

The concave model applies for terrain that contain no convex segments, and consists of two calculations - the ground reflection from flat ground, and ground reflection from deep 'valley'-shaped terrain, where the flat ground model becomes unstable. A transition factor `F_G` is used to 'blend' the two calculation regimes, and also helps to compensate for the omission of higher-order reflection paths from the model for valley-shaped terrain

The ground model uses the concept of Fresnel weightings, which are essentially a geometric weighting of how 'efficient' a ground segment is in contributing to the source-receiver reflection. A frequency-modified version of the weighting is used in the Harmonoise model. Refer below for calculation of Fresnel weightings.

The ground reflection from flat ground is the weighted sum of the reflection from each ground segment, weighted using the modified Fresnel weighting `w_k` of each segment:

`Delta L_(G,flat) = sum w_k Delta L_k`

The ground reflection from a single ground segment `k` is calculated using the following equation:

`Delta L_k = 10log_10(|1+C_k D_k Q_k|^2 + (1-C^2)|D_k Q_k|^2)`

The first term is for coherent reflection from the ground (i.e. pressure doubling), the second term is for incoherent reflection (i.e. energetic reflection).

`C_k` is a coherence factor for each ground segment `k` - i.e. whether the reflection is coherent (`C_k=1`), incoherent (`C_k=0`) or intermediate.

`D_k` is a geometric weighting factor for the reflection strength from segment `k`, which takes into account the additional path difference, any diffraction loss (e.g. from image sources/receivers) and phase differences associated with the reflection from segment `k`.

`Q_k` is the spherical wave reflection coefficient for each ground segment `k`

The ground reflection from valley-shaped terrain is similar, but is calculated as an overall term for the terrain segment as a whole, rather than as the sum of individual elements:

`Delta L_(G,"valley") = 10log_10(|1 + sum w_k C_k D_k Q_k|^2+sum w_k(1-C_k^2)|D_k Q_k|^2)`

The overall concave terrain ground reflection is calculated using a transition factor `F_G`:

`Delta L_(G,C) = F_G DeltaL_(G,flat) + (1-F_G)Delta L_(G,"valley")`

`F_G = 1-e^(-1//x_G^2)` is an empirically-tuned transition factor for flat/valley shaped terrain.

`x_G = N_w/sqrt(1+(f//f_c)^2)` is a non-dimensional frequency transition factor

`f_c` is a transition frequency (refer below for additional detail)

`N_w = sum w_k` is effectively a measure of the 'curvature' of the terrain:

Transition model The transition model uses a mixture of ground reflection, which is calculated using a modified version of the Concave ground model, and illuminated-zone diffraction, which is calculated using the diffraction model.

If `P_i...P_j` is the ground segment of points containing at least one convex surface, and `P_k` is the vertex which has the smallest negative path difference (i.e. closest to 0) - effectively the 'highest point' under the line-of-sight from `P_i` to `P_j`, the ground attenuation for the segment `P_i...P_j` under the transition model `Delta L_(G,T)` is calculated as follows:

`Delta L_(G,T) = chi Delta L_1 + (1-chi)Delta L_2`

`Delta L_1 = Delta L_D(P_i,P_k,P_j) + Delta L_(G,C)(P_i,P_k)+Delta L_(G,C)(P_k,P_j)` is the ground attenuation including diffraction:

`Delta L_D(P_i,P_k,P_j)` is the illuminated zone diffraction loss for the diffraction edge `P_k`
`Delta L_(G,C)(P_i,P_k)` is the ground loss (calculated using the Concave ground model) for the terrain from `P_i` to the diffraction edge `P_k`
`Delta L_(G,C)(P_k,P_j)` is the ground loss (calculated using the Concave ground model) for the terrain from the diffraction edge `P_k` to the end of the terrain interval at `P_j`

`Delta L_2 = Delta L_(G,C)(P_i,P_j)` is the ground loss (no diffraction) for the entire terrain between `P_i` and `P_j`

`chi` is a transition factor accounting for the balance between diffraction and reflection

`chi = chi_2 + (1-chi_1)(1-chi_2)`, where:

`chi_1={(1-e^(-1//tau_1^2) " for " tau_1 > 0),(1 " for " tau_1 <= 0) :}`

`tau_1 = (delta_("diff")-delta_(spek,avg))/(lambda//8)` compares the path difference for diffraction to the average terrain height
`delta_(spek,avg) = (sum w_k delta_(spek,k))/(sum w_k)` is the weighted average of the path difference for reflection `delta_(spek,k)` for each terrain segment

(Sign convention is that path differences are negative in the illuminated zone, hence `delta_("diff") > delta_(spek,avg)` because diffraction edge is the 'highest point' of the terrain.)

`chi_1` quantifies whether the diffraction is significant relative to the surface roughness of the terrain - by Rayleigh's criterion, if the size of the obstacle causing diffraction is less than `lambda//8` compared to the average terrain, then diffraction should be neglected since the object is better characterised as "roughness" of the ground. `chi_1 rarr 1` indicates the diffraction edge being insignificantly-small compared to the terrain.

`chi_2 = {(1-e^(-1//tau_2^2) " for " tau_2 > 0),(1 " for " tau_2 <= 0) :}`

`tau_2 = |delta_"diff"|/(lambda//64)`

`chi_2` is a 'transition' factor for continuity between the illuminated zone and the shadow zone. `chi_2 rarr 1` as the diffraction edge approaches the line-of-sight.

Coherence factor `C_k` accounts for the coherence of the ground-reflected component.

`C_k = C_A C_B`, where:

`C_A = e^(-0.5 sigma_phi^2)` accounts for phase fluctuations due to frequency-averaging and uncertainty of the source/receiver position

`sigma_phi` is the standard deviation of the phase fluctuations between direct and reflected sound
`(sigma_phi / phi)^2 = (sigma_f/f)^2[+(sigma_(h_s)/h_s)^2][+(sigma_(h_r)/h_r)^2]`

`sigma_f/f` accounts for integration of the calculation across a frequency band; `sigma_f/f = 1/3 (2^(0.5B)-2^(-0.5B))` where `B=1` for octave-bands, `B=1/3` for 1/3-octave bands.

`sigma_(h_s)` and `sigma_(h_r)` are the standard deviations of the source/receiver heights (input parameters to the model). The `sigma_(h_s)` term is only considered for the primary source (i.e. not included when calculating from a secondary source at a diffraction edge); similarly the `sigma_(h_r)` term is only included for the primary receiver

`phi` is the phase difference for the reflection.
Where the path difference of the reflection `delta_(spek,k)` is known, phase is calculated as `phi = k delta_(spek,k)`
Otherwise, the phase is approximated by `phi ~~ (2 pi f)/c_0 (2 h_s h_r)/d_(SR)`

`C_B = e^(-0.375 D_T gamma_T k^2 rho^(5//3) d_(SR))` accounts for phase differences due to fluctuations of the sound speed and atmospheric turbulence

`D_T` is a constant (0.364)
`gamma_T` is a turbulence structure parameter for wind speed and temperature fluctuations:
`gamma_T = (C_T/T)^2 +22/3 (C_W/c_0)^2`; however in practice the turbulence parameters `C_T` and `C_W` are not usually known, and it is more common to specify values of `gamma_T` as a tuneable input to the model.
Strutt uses a default value of `gamma_T = 5*10^-6` as recommended by Harmonoise for "moderate turbulence".
`rho = (h_s h_r)/(h_s+h_r)` is half of the mean separation of the direct and reflected ray paths. `rho=0` for the case of a hull segment.

For convex segments, `C_k` is calculated as above, except using the absolute value `|h_s|`,`|h_r|` of the source receiver heights.

Geometric weighting factor `D_k` accounts for the strength of the reflected sound relative to the direct sound. At its simplest (a ground reflection with no diffraction), it simply accounts for the additional propagation distance associated with the reflection.

`D_k` is calculated for the following geometric scenarios:

The geometric terms `D_k` given above strictly only apply to concave or hull ground segments - i.e. for a convex segment (because either the source or receiver height below the segment is negative) an image reflection is not geometrically possible.
However, ground reflection from convex surfaces is still included in the Harmonoise model, by multiplying the `D_k` terms for convex segments by `(p_d (S,X,R prime))/(p_d (S,X,R))`, where `X` is the point where the line of sight from `S`-`R` intersects the ground plane of the convex element (i.e. the 'imaginary' specular reflection point for image source `S prime`):

In effect, for convex segments, Harmonoise calculates the ground reflection from a diffracted ray that 'bends' around point `X` to reflect off the convex segment, and the `D_k` term is reduced in strength to account for the additional diffraction.

Spherical reflection coefficient `Q_k` is calculated using a modified version of the classical theory for reflection of sound above a ground plane (the Weyl-Van der Pol equation) based on the work of Chien and Soroka:

`Q_k = R_p + (1-R_p)(F_Q)^(n_G)`, where:

`R_p` is the plane-wave reflection coefficient `(Zcos(theta)-Z_0)/(Zcos(theta)+Z_0)`

`Z` is the impedance of the ground plane, calculated from the flow resistivity `sigma` using the Delaney and Bazley theory for oblique incidence as described in the Porous Absorber calculation)

`Z_0` is the characteristic acoustic impedance (`rho_0 c_0`) of air

The second term `(1-R_p)(F_Q)^(n_G)` in the spherical wave coefficient is known as the ground wave term, and corrects the plane wave coefficient `R_p` for the interaction of a curved wavefront and the ground plane

`F_Q` is known as the boundary loss factor and is given (analytically) by:

`F_Q (w) = 1+j sqrt(pi) w e^(-w^2)"erfc"(-jw)`, where `"erfc"(x)` is the complimentary error function

`w = 0.5 (1+j) sqrt(k d_(S prime R))(cos(theta)+1/Z)` is a complex numerical distance term: `w = w_r+jw_i`

`d_(S prime R)` is the distance of the reflected path (image source-receiver)

Evaluating `F_Q` analytically is too mathematically-intensive, and therefore Strutt uses an approximate form to calculate `F_Q`:

Harmonoise modifies the boundary loss factor `F_Q` to be more-realistic for situations where the sound is close to grazing incidence via term `n_G`:
`n_G = 1-0.7e^(-h_m/(lambda//32))`

Note that `Q_k` is calculated as if for an infinite-sized reflecting surface. Finite-sized surface effects are taken into account via the use of Fresnel weightings `w_k` for each segment.

For convex segments, `Q_k` is calculated as above, except using the absolute value `|h_s|`,`|h_r|` of the source receiver heights.

The Fresnel weightings, `w_k` are used to account for the "efficiency" with which a surface can reflect from the source to the receiver.

The Fresnel weightings are based on the Fresnel ellipsoid, which has the image source `S prime` and the receiver `R` as its two focal points:

The ellipsoid is also defined by the Fresnel parameter `n_F` which is proportional to the wavelength of sound, and effectively defines the locus of points that can reflect efficiently to the receiver point. A frequency-dependent value of `n_F` is used in Harmonoise (based on validation studies), which helps to account for phenomena such as surface roughness/scattering that decrease the strength of the ground-reflected component at high frequency.

The intersection between the Fresnel ellipsoid and the ground plane of the segment describes an elliptical region, the Fresnel ellipse. A reflection in Harmonoise is only considered significant if it occurs within the Fresnel ellipse.
Note that the centre of the Fresnel ellipse is not coincident with the specular reflection point.

The centre of the ellipse occurs at a distance `d_F` from the projection of the source to the ground plane, and the semi-major axis of the ellipse is `a`:
`d_F = d_(SR)/2(1+(h_s^2-h_r^2)/(D^2-d_(SR)^2))`

`a = 1/2 sqrt((D^4+(D_S^2-D_R^2)^2-2D^2(D_S^2+D_R^2))/(D^2-d_(SR)^2))`

`D = lambda/n_F + sqrt((h_S+h_R)^2+d_(SR)^2)`



These paraneters are used to calculate the unmodified Fresnel weightings `w_(F,k)` for each ground segment, using `n_F=lambda//8` for the unmodified weighting. The unmodified Fresnel weighting `w_(F,k)` is the proportion of the surface that lies within the Fresnel ellipse (i.e. the proportion of the surface that can efficiently-reflect to the receiver):

`w_(F,k) = F_w (xi_2)-F_w (xi_1)`, where:

`F_w (x)= {(0 " for " x <= -1),(1-1/pi (cos^-1(x)-x sqrt(1-x^2)) " for " -1 < x < 1),(1 " for " x >= 1) :}`
`xi` is a non-dimensional distance within the Fresnel ellipse: `xi_m = (d_m-d_F)/a`, i.e. `|xi_m| > 1` indicates a distance outside of the Fresnel ellipse.
`xi_1` and `xi_2` are the `xi`-coordinates of the beginning and end of the terrain segment, respectively.
For continuity, `F_w = 0` at the source and `F_w = 1` at the receiver.

The unmodified Fresnel weightings do not provide accurate results at high frequency, and therefore Harmonoise uses the modified Fresnel weighting, `w_k` in calculations. (The unmodified weightings are, however, used as intermediate steps in some calculations.)
The modified Fresnel weighting is centred around the specular reflection point, and uses a frequency-dependent value of `n_F`:

`f_c` is the transition frequency (see below).

The modified weightings use a modified version of the non-dimensional distance:

`w_k = F_w (xi prime_2) - F_w (xi prime_1)`, where:

`xi prime_m = (xi_m-xi_C)/(1-xi_m xi_c)`

`xi_C = (d_C-d_F)/a`

`d_C = alpha(f)d_F+[1-alpha(f)]d_(SP)`

`d_(SP) = d_(SR) h_S/(h_S+h_R)` is the distance to the specular reflection point


Transition frequency `f_c` is used to divide the calculation into high-frequency and low-frequency regimes, with frequency-dependent parameters such as `F_G`, `n_F` and `alpha (f)` being calculated based on the ratio of the band centre frequency `f` and `f_c`.
`f_c = sqrt(f_(min) f_(max))`

`f_(min)` and `f_(max)` are based on the phase difference between the direct and reflected sound and are defined as follows:

`phi_(max) (f_(min))=pi/2`

`phi_(max) (f_(max))=pi`

where, for each frequency:

`phi_(max)` is the maximum value of the phase difference `phi_k` of all the terrain segments `k` over the calculation interval.

`phi_k = "Im"(Q_k)+k*delta_(spek,k)`

i.e. the phase change associated with the ground reflection (the imaginary part of the reflection coefficient `Q_k`) and the phase change associated with the additional path difference `delta_(spek,k)` between direct and reflected sound).
In evaluating `phi_k`, segments with zero unmodified Fresnel weighting `w_(F,k)` and/or convex segments are excluded.

In general, `phi_(max)` equal to `pi` or `pi/2` will not occur precisely at the centre band frequency `f` used for calculation. This means that `f_(max)` and `f_(min)` are determined as follows:

Atmospheric scattering is calculated as follows:

`Delta L_(scat) = 25 + 10log_10 gamma_T + 3 log_10 f/1000 + 10 log_10 d_(SR)/100`

The effect of atmospheric scattering is to limit the practical maximum excess attenuation via the ground profile (i.e. screening or ground absorption) by scattered sound "short-circuiting" the diffraction/reflection paths.

Refraction of sound waves is accounted for by a coordinate transformation of the terrain profile - in effect 'straightening out' the curved ray paths to an equivalent geometry with straight ray paths. The Harmonoise excess attenuation is then calculated using the "neutral-atmosphere" method described above using the new coordinate system.

The coordinate transformation is made in the complex plane. If the original set of coordinates are (in Cartesian) form `y = x + jz`, the mapping transforms to an equivalent coordinate system `dot y= dot x+ j dot z`

Harmonoise assumes a linear vertical sound speed profile `c(z) = c_0(1+z/R_c)`, where `R_c` is the approximate radius of curvature of the curved ray paths.

The coordinate transformation results in a transformed system with approximately-constant sound speed. The assumptions behind the coordinate transformation and the approximation of constant sound speed are valid only for `|R_c|>5d_(SR)`. The new coordinates `dot x` and `dot z` are given by the following equations:

`dot x = (C_0^2 ddot x)/(ddot x^2+(C_0+ddot z)^2`

`dot z = (C_0(ddot x^2+ddot z^2+ddot z C_0))/(ddot x^2+(C_0+ddot z)^2)`

`ddot x = x - (x_0+x_N)/2`

`ddot z = z-(z_0+z_s+z_n+z_r)/2`

`C_0=2 ((z_s+z_r)/2+R_C)`

Recall that the source and receiver coordinates in the global coordinate system are `(x_0,z_s)` and `(x_n,z_n+z_r)` and that Strutt assumes `z_0=0` (all heights are relative to terrain at the source point).

The overall Harmonoise excess attenuation is calculated as:

`A_(excess) = 10log_10(10^((sum Delta L_D + sum Delta L_G)//10)+10^(Delta L_"scat"//10))`
(i.e. the sum of all diffraction and ground reflection terms, logarithmically added with the scattered sound level.


Comments or suggestions to