A New Method for Determining Optimal Regularization Parameter in Near-Field Acoustic Holography

Tikhonov regularizationmethod is effective in stabilizing reconstruction process of the near-field acoustic holography (NAH) based on the equivalent sourcemethod (ESM), and the selection of the optimal regularization parameter is a key problem that determines the regularization effect. In this work, a newmethod for determining the optimal regularization parameter is proposed.The transfer matrix relating the source strengths of the equivalent sources to the measured pressures on the hologram surface is augmented by adding a fictitious point source with zero strength. The minimization of the norm of this fictitious point source strength is as the criterion for choosing the optimal regularization parameter since the reconstructed value should tend to zero. The original inverse problem in calculating the source strengths is converted into a univariate optimization problem which is solved by a onedimensional search technique. Two numerical simulations with a point driven simply supported plate and a pulsating sphere are investigated to validate the performance of the proposedmethod by comparison with the L-curve method.The results demonstrate that the proposed method can determine the regularization parameter correctly and effectively for the reconstruction in NAH.


Introduction
Near-field acoustic holography (NAH) is an effective technique for noise sources identification and acoustic field visualization.By measurements on a hologram surface near the sound source, the acoustic quantities such as pressure, particle velocity, and intensity anywhere including the source surface can be reconstructed, and the three-dimensional acoustic filed can be predicted.In the past years, many alternative methods of transform algorithm for realizing the NAH have been developed, for example, the spatial Fourier transform [1,2], the boundary element method [3][4][5], the Helmholtz equation least squares method [6,7], the statistically optimized method [8,9], and the equivalent source method (ESM) [10][11][12][13][14].Among the methods, the ESM has attracted more attention with the advantages of the simple principle, the high calculation efficiency, and good adaptability for the shapes of the sound sources.
However, the reconstruction of acoustic quantities of sound sources from near-field measurement data is an inverse problem, which is different from the traditional acoustic radiation calculation.Consequently, the reconstruction stability is a key problem in NAH technology due to the fact that the realization process is ill-posed and the reconstructed results are highly sensitive to signal-to-noise ratio (SNR) [15].If the conventional method for solving inverse matrix from the source to the hologram surface is used directly, the measurement error on the hologram surface as the input may be amplified sharply, and even the reconstructed results are completely unavailable.At present, the regularization approach has been widely used in the reconstruction process to restrict the influence of measurement errors on the reconstructed results and ensure the validity of the reconstruction performance.
Currently, the direct regularization methods and the iterative regularization methods are mostly used for dealing with ill-posed problem of NAH [16,17].The direct regularization methods include truncated singular value decomposition (TSVD) method [18] and Tikhonov regularization method [19][20][21].The iterative regularization methods include Landweber iterative regularization method [22] and conjugate gradient method [23].The essence of the regularization method is to find an optimal approximate solution to the actual solution by a certain control strategy.So far there is no single regularization that is absolutely superior to the rest.

Shock and Vibration
For each of different regularization methods, the selection criteria of the optimal regularization parameters are the critical factors that determine the regularization effect.
The methods for selection of optimal regularization parameter can be classified into prior [24] and posterior strategy.Compared with prior strategy, the posterior strategy does not require the prior knowledge of the noise variance in advance.The generalized cross validation (GCV) method [25] and the L-curve method [26,27] as the posterior strategy are the two most widely used methods for selecting the optimal regularization parameters, both of which have good adaptability in the engineering application.Although the convergence of the GCV criterion has been proven theoretically, the GCV function curve sometimes is too flat, which makes it difficult to determine the minimum value of the GCV function.The L-curve method describes the comparison between the solution norm and residual norm with the regularization parameters in a log-log scale coordinate and adopts the maximum curvature of L-curve as the location of the optimal parameter.The L-curve method is reliable for determination of the optimal regularization parameter in most cases, but it may fail to search the maximum curvature position when there is more than one inflexion on the curve or the curve is not L-shaped approximately [28].
This paper focuses on a new method for determining the optimal parameters of Tikhonov regularization in NAH technique based on the ESM.In the proposed method, the transfer matrix from the source to the hologram surface is updated by augmenting a fictitious point source with zero strength near the locations of the equivalent sources.The minimization of the norm of this fictitious point source strength can be as the criterion for selection of the optimal regularization parameter since the reconstructed value should tend to zero.Thus the original inverse problem is converted into a univariate optimization problem that can be solved by one-dimensional search technique.The numerical simulations are investigated to demonstrate the validity of the proposed method, and the results show that the novel proposed method is able to select the regularization parameter correctly and effectively.

NAH Based on ESM
The basic idea of the NAH based on ESM is that the actual radiated sound field can be replaced by superposing the acoustic field produced by a number of equivalent sources distributed inside the source surface.The source strengths of the equivalent sources can be evaluated by matching the acoustic pressures measured on the hologram surface, and the acoustic quantities in the acoustic field can be obtained by the source strengths and the transfer matrices constructed by the equivalent sources.
As shown in Figure 1, the hologram surface S h is placed near the outside of the source surface S, and the equivalent source surface S q is located on the other side of the source surface S. Suppose that there are  equivalent sources on the surface S q , and the acoustic pressure at the field point r can be expressed as where (r) is the acoustic pressure at the field point r,  is the density of the medium,  is the speed of sound,  = / = 2/ is the wave number,  is the wavelength,  is the angular frequency, r q and (r q ) are the location vector and the strength of the th equivalent source, respectively, and (r, r q ) is the free field Green's function obtained from the expression Given that  measurement points (r h1 , r h2 , . . ., r h ) are selected on the hologram surface S h and by applying (1) to each point, the pressure column vector at the measurement positions can be represented in matrix form as where P h = [(r h1 ), (r h2 ), . . ., (r h )] T is the pressure column vector at  measurement points on the hologram surface, Q = [(r q1 ), (r q2 ), . . ., (r q )] T is the source strength column vector of the equivalent sources, G h (, ) = i(r h , r q ) is the transfer matrix relating the source strengths of the equivalent sources on the surface S q to the measured pressures on the surface S h , and r h is the location vector of the mth measurement point on the hologram surface S h .If the number of measurement points is not less than the number of equivalent sources, that is,  ≥ , the source strength column vector Q can be uniquely obtained by solving inverse of G h as where the superscript "+" denotes the generalized inverse.
Shock and Vibration 3 By using the source strength column vector Q, the pressure and particle velocity on any field point r can be calculated as where G fp and G fv are the transfer matrices relating the source strengths of the equivalent sources on the surface S q to the pressure and particle velocity at the field point r, respectively, and Similarly, the pressure and normal velocity on the surface of the acoustic source can be reconstructed in the matrix form as where G sp and G sv are the transfer matrices relating the source strengths of the equivalent sources on the surface S q to the pressure and normal velocity on the surface of the source, respectively, and in which r s is the location vector of the lth point on the surface of the source and n is the outward normal of the source surface S.

Tikhonov Regularization
In practice, there always exist errors in the measured pressures on the hologram surface.If (4) is used to calculate the source strengths of the equivalent sources Q directly, the reconstructed results are not reliable.So it is very critical to take measures in stabilizing the reconstruction procedure.Now the Tikhonov regularization is one of the most effective and popular methods to deal with the inverse problem.This method transfers the solution of (3) into the minimization problem as where  > 0 is the regularization parameter, which controls the weight between the residual norm and the solution norm.The first term of the right side is used to make the approximate solution satisfy (3) as far as possible but cannot be strictly true, because the measured errors in P h will cause the solution to deviate from the true value in the form of the oscillation with high frequency.The role of the second term is to suppress the oscillation of the solution and make the solution as smooth as possible.
The normal equation associated with ( 9) is given by where Ι is the identity matrix.The Tikhonov regularization solution can be expressed as The regularization parameter  determines the compromise between the residual norm ‖G h Q − P h ‖ 2 2 and the solution norm ‖Q‖ 2  2 , and it plays an indispensable role in designing stable inversion algorithms and in obtaining physically plausible solutions.More precisely, when  is too large, the residual norm becomes huge, since the solution is overly smooth, which often leads to undesirable effects.Conversely, it may be unstable and plagued with spurious and nonphysical details when  is too small.Therefore, the determination of the optimal regularization parameter constitutes one of the major inconveniences and difficulties in applying existing regularization method to NAH technique.

Method for Determining Optimal Regularization Parameter
Many researchers have developed a number of selection criteria and algorithms from different points of view, and each of them has their merit, shortcoming, and applicable condition [29].Most of the methods endeavor to excavate the prior knowledge of the original problem or set some assumptions to derive some criteria for determining the regularization parameters.Mathematically, the regularization solution is expected to be located in the neighborhood of the true solution, but in fact the current criteria are usually focused on the information of outputs, and the neighborhood is out of the question due to the absence of information of the true solution.For the reconstruction problem in NAH based on ESM, the following gives a new criterion, which can make use of the true solution information and has a clear physical meaning.
After determining the positions of all the equivalent sources, the corresponding source strengths will be obtained by matching the acoustic pressures measured on the hologram surface.In this case, there is no sound source or source strength is zero except for the positions of the equivalent sources.Here, the proposed method is to add a fictitious point source to the location where there is no equivalent source near the surface S q ; then the transfer matrix can be reconstructed to obtain a new augmented transfer matrix.Thus, the source strength column vector of the equivalent sources Q can be expressed as where Q 1 and Q 2 are the source strength column vector of the equivalent sources existing originally and the fictitious point source added later, respectively.
Substituting ( 12) and ( 13) into (10), it can be transformed as From ( 14), the unknown source strength of the fictitious point source can be obtained.Furthermore, the true value of Q 2 is known (zero), which provides available information for the solution of the problem; that is, a good regularization solution for the updated augmented system corresponding to the part of the fictitious point source Q ()  2 should tend to zero.Therefore, the minimization of the norm of this fictitious source strength Q ()  2 can be the criterion for selection of the optimal regularization parameter: Qualitatively, the change trend of the norm of strength of the fictitious source ‖Q 2 ‖ 2 with the increase of regularization parameter  can be deduced as follows: (a) When  = 0, regularization does not work, and ‖Q 2 ‖ 2 is much larger than its true value (true value of ‖Q 2 ‖ 2 is zero) because it is contaminated by the measurement errors.
(b) With the increase of , the regularization starts to work and ‖Q 2 ‖ 2 tends to its true value gradually.In this process, ‖Q 2 ‖ 2 decreases with the increase of ; that is, ‖Q 2 ‖ 2 goes from the large oscillation to its true value, which is zero.
(c) With the decrease of ‖Q 2 ‖ 2 , when ‖Q 2 ‖ 2 reaches the point where the position is closest to the true value, it starts to deviate from its true value as  gradually increases.So in this process, ‖Q 2 ‖ 2 increases with the increase of .
(d) Note that the norm of strengths of the equivalent sources ‖Q‖ 2 decreases monotonically with ; that is, ‖Q‖ 2 → 0 as  → ∞, and thus ‖Q 2 ‖ 2 → 0. When ‖Q 2 ‖ 2 increases to a certain level with , it begins to tend to zero with the increase of .When  is too large, the solving process produces an overregularization with too large regularization error.
Therefore, it can be seen that the optimal regularization parameter should be at the point corresponding to the first local minimum where ‖Q 2 ‖ 2 goes from decreasing to increasing.
The minimization problem related to (15) can be solved by one-dimensional search technique, which is exactly the same as the L-curve method in the Matlab package developed by Hansen [30].In order to capture the optimal solution more easily in the search interval, the logarithmic scale coordinate is also used to ensure that the curve of ‖Q 2 ‖ 2 with  can be flat and broad in the vicinity of the optimal solution.
Furthermore, considering the initial situation without the fictitious point source, the corresponding solution equation can be expressed as According to (14), if Q 2 ̸ = 0, Q 1 will contain the additional interference generated by Q 2 .In order to obtain a better solution, after determining the optimal regularization parameter  opt , ( 16) is used to solve for Q 1 as the final solution, which can be expressed as Thus, by using the reconstructed source strength column vector Q , the acoustic quantities anywhere in the sound field can be predicted.
In addition, the essence of the proposed method to select the regularization parameter can be illustrated by the physical process of the equivalent source method.If the original  equivalent sources can match the sound pressure on the hologram surface properly, then a good approximation for the sound field can be obtained by using the  equivalent sources.To date, many literatures have been published and several satisfactory results have been obtained for setting the appropriate .In general, if the number of measurement points is greater than the number of equivalent sources ( ≤ ), the strengths of the  equivalent sources can be uniquely determined by the singular value decomposition (SVD) technique.Thus the sound field can be replaced directly by using the  equivalent sources, and the source strength column vector of the original equivalent sources can be obtained by using (10) as where Q 1 is the source strength column vector of the original equivalent sources before the fictitious source is added.When  + 1 equivalent sources are used to approximate the sound field and the positions of  equivalent sources are unchanged, Q 1 can be obtained by using (14) as By comparing ( 18) with (19), it is possible to obtain According to the properties of 2-norm, there is an inequality as where cond(G 0] T can be regarded as a good special solution to approximate the sound field when using  + 1 equivalent sources.That is, constructing such a good special solution of the  + 1 dimensional source strengths can be a suitable criterion to determine the optimal regularization parameter.When this criterion is satisfied, a good regularization parameter can be determined because the sound field can be approximated using the  + 1 equivalent sources (the source strengths column vector is [Q T 1 0] T ) and using N equivalent sources.

Numerical Simulations
To examine the performance of the proposed method to determine the optimal regularization parameter in NAH, two numerical simulations were investigated, and the reconstruction results by using the proposed method were compared with those by using the L-curve method.

Simulation 1:
A Vibrating Plate.This simulation case is a point driven simply supported aluminum plate with the thickness of 3 mm mounted in an infinite baffle.The dimension of the plate is 0.3 m × 0.3 m, and a harmonic force with an amplitude of 1 N is exerted vertically to the plate at the center.The radiated sound field is calculated by using Rayleigh's integral formulation [31].The center of the plate is set as the coordinate origin, and the -axis and -axis are parallel to the edges of the plate, respectively.The hologram plane and the equivalent source plane are placed at z = 0.05 m and  = −0.03m, respectively, with the same dimension as the plate and a measurement gird of 25 × 25.The random noise is added to the measurement data at a level of 30 dB.A fictitious point source with zero strength is located at 0 m, 0 m, and −0.03 m.The pressure and normal surface velocity on a plane 0.03 m away from the plate (z = 0.03 m) were reconstructed based on the "measured" pressures on the hologram surface, and the L-curve method and the proposed method are used to determine the optimal parameters of Tikhonov regularization.To quantify the reconstruction accuracy, the relative error is defined in the percentage form as where A re is the reconstructed pressure or normal surface velocity and A th is the corresponding theoretical value. determined by the L-curve method regularization parameter are determined as 1.1933 × 10 −2 and 3.7901 × 10 −3 at the corner of the L-curve, respectively, and the relative errors between the reconstructed pressure and the theoretical value are 6.68% and 6.37% by using the parameters for regularization, respectively.Figure 3 shows the change curve of the object function ‖Q 2 ‖ 2 with the regularization parameter  by using the proposed method at 500 Hz and 800 Hz.The trends of the curves are consistent with the speculation mentioned above.The regularization parameters can be chosen as 1.0010 × 10 −2 and 3.5481 × 10 −3 at the minimum of the curve, and the relative errors are 6.63% and 6.44%, respectively.For comparison, the locations of the parameter determined by using the Lcurve method are also marked in Figure 3.It can be found that the values of the regularization parameter determined by using the L-curve method and the proposed method are very close, and both methods can give an appropriate regularization parameter for the pressure reconstruction in NAH.In addition, as shown in Figure 3, if a linear scale is used for searching the regularization parameter, the optimal solution may be crossed in the very narrow V-shaped region containing it and stray into the right side of the monotonically decreasing interval, which will result in a search failure.In the proposed method, the logarithmic scale is used to extend the search interval of regularization parameter, so it is easier to search for the optimal solution.
Figure 4 presents the comparisons of the reconstructed pressure with the corresponding theoretical value by using the two methods at 500 Hz and 800 Hz.It can be seen that the reconstructed results by using these two methods show good agreements with the theoretical values.The L-curve of the pressure reconstruction for the vibrating plate at 1000 Hz is shown in Figure 5.It can be seen that there is more than one inflexion on the L-curve, which causes the L-curve method to have a wrong detection for the regularization parameter.
The regularization parameter is determined as 6.9374, and the reconstruction error is 53.46%.
As a comparison, Figure 6 shows the curve of the proposed method for choosing the regularization parameter.The regularization parameter by using the proposed method is 2.1877 × 10 −3 , and the reconstructed relative error is 6.68%, which is much smaller than that by using the L-curve method.By the comparison of the two errors, it can be seen that the proposed method is valid for the pressure reconstruction when the L-curve method fails to determine the regularization parameter.The comparison of the reconstructed pressure with the corresponding theoretical value at 1000 Hz is shown in Figure 7.It can be seen that the reconstructed pressure by using the L-curve method differs largely from the theoretical value due to the overregularization, while the reconstructed pressure by using the proposed method can give a better reconstruction.
Figure 8 shows the L-curve of the normal surface velocity reconstruction for choosing regularization parameter at 600 Hz.It can be seen that the position of the optimal parameter is not located at the corner of the L-curve, and there is a wrong detection problem in the L-curve method.The regularization parameter is determined as 233.27, which seems too large.Using the parameter determined by the Lcurve method for regularization, the reconstruction relative error is 43.68%.The curve of the proposed method used to determine the optimal regularization parameter is shown in Figure 9.The regularization parameter can be determined as 7.7618 × 10 −2 by detecting the minimum of the curve.The relative error of the normal surface velocity in the proposed method is 5.97%, which is much smaller than that in the Lcurve method.
Figure 10 shows the comparison of the amplitude distributions of the reconstructed normal surface velocities with the corresponding theoretical values by using the two methods at 600 Hz.It can be seen that there is overregularization in ) reconstructed by using the L-curve method at 500 Hz, (c) reconstructed by using the proposed method at 500 Hz, (d) theoretical value at 800 Hz, (e) reconstructed by using the L-curve method at 800 Hz, and (f) reconstructed by using the proposed method at 800 Hz. the L-curve method due to the fact that the optimal parameter determined by this method is too large, which results in losing too much information of the real holographic pressures.Meanwhile the reconstructed normal surface velocity by using the proposed method shows good agreement with the theoretical value, which indicates that the proposed method can choose a good regularization parameter when the Lcurve method has wrong detection.The influence of the position of the additional fictitious source on the reconstruction accuracy is also investigated to  determined by the L-curve method  point source coincides with the position of the measurement point on the hologram surface), which means that the proposed method is not sensitive to the position of the additional fictitious point source.Actually, as long as the position of the additional fictitious source is near the equivalent source surface, the proposed method can get a good regularization parameter to ensure stability of the reconstruction.

Simulation 2:
A Pulsating Sphere.A pulsating sphere with a radius of 0.1 m is investigated as depicted in Figure 12.
There are 62 nodes distributed on the source surface with the equal angle interval /6.The analytical pressure of a pulsating sphere with a radius of  at the point  in the sound field can be expressed as where V 0 is the radial velocity,  is the sphere radius, and  is the vibration frequency of sound source.
In the simulation, the radial velocity V 0 is 0.01 m/s, the air density  is 1.2 kg/m 3 , and the speed of sound  is 344 m/s.The  determined by the L-curve method equivalent source surface is a spherical surface with a radius of 0.08 m and has the same distributed form as the source surface with 62 equivalent point sources.The hologram surface with the size of 1 m × 1 m is located 0.2 m away from the center of the sphere, and there are 16 × 16 measurement points distributed on the surface with a uniform interspacing.All the acoustic pressures on the hologram surface are given by (23).A fictitious point source with zero strength is located at 0.0739 m, 0.0173 m, and 0.040 m and the field point located at 0 m, 0 m, and 1 m is chosen as the reconstructed point.The random noise is added to the measurement data at a level of 30 dB.
Figure 13 shows the L-curves of the pressure reconstruction for the pulsating sphere at 300 Hz, 500 Hz, 700 Hz, and 900 Hz.It can be seen that all the curves have a quite clear L-shape.The regularization parameters are determined as 1.0128 × 10 −4 , 5.003 × 10 −4 , 8.1372 × 10 −4 , and 1.1166 × 10 −3 at the corner of the L-curve, respectively, and the relative errors between the reconstructed pressure and the theoretical   By using the proposed method, the curves for choosing regularization parameter at 300 Hz, 500 Hz, 700 Hz, and 900 Hz are shown in Figure 14.The optimal regularization parameters are chosen as 2.7541 × 10 −3 , 6.1662 × 10 −3 , 1.3832 × 10 −2 , and 2.6923 × 10 −2 corresponding to the bottom of the curve, respectively, and the reconstruction relative errors are 9.89%, 3.22%, 11.86%, and 13.38%, respectively, all of which are much smaller than those by using the L-curve method.It can be seen that the proposed method can give a good regularization parameter in NAH technique based on the ESM.
Figure 15 shows the relative errors between the reconstructed pressure and the theoretical value by using the two methods for choosing the regularization parameter over the frequency range from 50 Hz to 1000 Hz.It can be seen that the reconstruction errors by using the proposed method are always smaller than those determined by using the L-curve method, and the relative errors of pressure are below 20% over most of the whole frequency range, which demonstrates the validity of the proposed method for the selection of regularization parameter in NAH.
The influence of the SNR on the reconstructed pressure at different frequencies is also examined.Figure 16 shows the reconstruction relative error versus SNR in the range of 10 dB to 50 dB by using the proposed method at 200 Hz, 300 Hz, 450 Hz, 600 Hz, 800 Hz, and 1000 Hz, respectively.It can be seen that the relative errors calculated by using proposed method are insensitive to the change of SNR, which indicates that the proposed method is of good stability for choosing regularization parameter.

Conclusion
A new optimal parameter selecting method for the Tikhonov regularization was proposed in NAH based on the ESM.In the proposed method, the criterion for determining the optimal regularization parameters was given by adding a fictitious point source with zero strength to make a part of the solution of the augmented problem known, which has a clear physical meaning.A one-dimensional search technique was used to minimize the norm of the fictitious point source strength for determination of the optimal parameter in the Tikhonov regularization method.The validity of the proposed method has been demonstrated by numerical simulations with kinds of sound sources including a point driven simply supported plate and a pulsating sphere.It was shown that the proposed method can get a better regularization parameter for the reconstruction with good accuracy and stability in NAH compared with the L-curve method.Future research will emphasize the existence and uniqueness of the minimum point on the curve of the norm of the fictitious point source strength with the regularization parameter in mathematics.The present work also can provide a new idea for solving other inverse problems, that is, try to add new known information  determined by the L-curve method The proposed method The L-curve method

Figure 1 :
Figure 1: Schematic diagram of NAH based on ESM.

Figure 2 :
Figure 2: L-Curve of the pressure reconstruction for the vibrating plate at (a) 500 Hz and (b) 800 Hz.

Figure 2
shows the L-curves of the pressure reconstruction for the vibrating plate using the Matlab package developed by Hansen at 500 Hz and 800 Hz.The values of

Figure 3 :
Figure 3: Object function versus the regularization parameter of the pressure reconstruction for the vibrating plate at (a) 500 Hz and (b) 800 Hz.

Figure 4 :
Figure 4: Comparison of reconstructed pressures.(a) Theoretical value at 500 Hz, (b)reconstructed by using the L-curve method at 500 Hz, (c) reconstructed by using the proposed method at 500 Hz, (d) theoretical value at 800 Hz, (e) reconstructed by using the L-curve method at 800 Hz, and (f) reconstructed by using the proposed method at 800 Hz.

2 Figure 5 :
Figure 5: L-Curve of the pressure reconstruction for the vibrating plate at 1000 Hz.

Figure 6 :Figure 7 :
Figure 6: Object function versus the regularization parameter of the pressure reconstruction for the vibrating plate at 1000 Hz.

2 Figure 8 :
Figure 8: L-Curve of the normal surface velocity reconstruction for the vibrating plate at 600 Hz.

Figure 9 :
Figure 9: Object function versus the regularization parameter of the normal surface velocity reconstruction for the vibrating plate at 600 Hz.

Figure 10 :Figure 11 :Figure 12 :
Figure 10: Comparison of reconstructed normal surface velocities at 600 Hz.(a) Theoretical value, (b) reconstructed by using the L-curve method, (c) reconstructed by using the proposed method, and (d) normal surface velocity along the middle row.

Figure 14 :
Figure 14: Object function versus the regularization parameter of the pressure reconstruction for the pulsating sphere at (a) 300 Hz, (b) 500 Hz, (c) 700 Hz, and (d) 900 Hz.

Figure 15 :
Figure 15: Reconstruction errors of pressure by using the two methods versus frequency.

Figure 16 :
Figure 16: Reconstruction errors of pressures by using the proposed method versus SNR.