Shock and Vibration
Volume 2020 (2020), Article ID 5491957, 10 pages
https://doi.org/10.1155/2020/5491957

Research Article

Eccentricity Parameters Identification for a Motorized Spindle System Based on Improved Maximum Likelihood Method

Wengui Mao,¹ Chaoliang Hu,¹ Jianhua Li,¹ Zhonghua Huang,¹ and Guiping Liu²

¹Hunan Provincial Key Laboratory of Wind Generator and Its Control, College of Mechanical Engineering, Hunan Institute of Engineering, Xiangtan 411101, China
²State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, College of Mechanical and Vehicle Engineering, Hunan University, Changsha 410082, China

Correspondence should be addressed to Wengui Mao; mwglikai@163.com

Received 22 November 2019; Revised 18 June 2020; Accepted 30 June 2020; Published 29 July 2020

Academic Editor: Jie Zhang

Copyright © 2020 Wengui Mao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

As a kind of rotor system, the electric spindle system is the core component of the precision grinding machine. The vibration caused by the mass imbalance is the main factor that causes the vibration of the grinding machine. Identifying the eccentricity parameters in an electric spindle system is a key issue in eliminating mass imbalances. It is difficult for engineers to understand the approximate range of eccentricity by experience; that is, it is difficult to obtain a priori information about eccentricity. At the same time, due to the geometric characteristics of the electrospindle system, the material factors and the randomness of the measurement response, these uncertain factors, even in a small case, are likely to cause large deviations in the eccentricity recognition results. The search algorithm used in the maximum likelihood method to identify the eccentricity parameters of the electrospindle system is computationally intensive, and the sensitivity in the iterative process brings some numerical problems. This paper introduces an Advance-Retreat Method (ARM) of the search interval to the maximum likelihood method, the unknown parameter increment obtained by the maximum likelihood method is used as the step size in the iteration, and the Advance-Retreat Method of the search interval is used to adjust the next design point so that the objective function value is gradually decreasing. The recognition results under the three kinds of measurement errors show that the improved maximum likelihood method improves the recognition effect of the maximum likelihood method and can reduce the influence of uncertainty factors on the recognition results, and the robustness is satisfactory.

1. Introduction

The high-speed motorized spindle is supported by bearings, which integrates the functions of the machine tool spindle and the high-speed motor in structure, and realizes the “zero transmission” between the variable frequency motor and the machine tool spindle. It is a product that combines the motor with the spindle. The rotor of the motor is the rotating part of the spindle. In theory, the motorized spindle can be regarded as a high-speed motor. It focuses on many key technologies, such as the design, manufacture, assembly and control of high-speed motor, and high-performance machine tool spindle, and has the characteristics of high-speed and zero-transmission chain. The motorized spindle is the core component of optical grinding machine tools, and its performance has a significant effect on the machining accuracy because of the combination of spindle and tool. In recent years, with the increase of grinding machine speed, machine tool failures frequently occur, and the research on safe operation and fault diagnosis of high-speed grinding machines has become an important research content in the development of machine tool technology [1, 2]. Unbalance is one of the most common faults in high-speed grinding machine, which causes equipment vibration, coupling deflection, bearing wear and tear, and other hazards [3, 4]. The fault will cause the electric spindle system to bear the eccentricity parameter, which will cause changes in the dynamic performance of the spindle system, such as eigenvalues, thus aggravating the operation fault. Motorized spindles widely used for high-speed precision machine tools are very sensitive to the mass unbalance of rotors; thus, their balancing problem is always a research hotspot. Although many significant studies were done regarding the theory and application of various rotor balancing technologies for motorized spindles. Dynamic balance technology including influence coefficient method, cross-correlation method, and mode method [5, 6] is commonly used to balance machine tool spindle faults at present. The basic process involves installing the acceleration and speed sensors, collecting the spindle vibration signal, using a suitable data processing method to the extract amplitude and phase of the signal, and outputting the correction masses and mounting positions. However, the particularity of motorized spindles is not considered carefully in the existing balancing approaches. When the rotor unbalance of a motorized spindle occurs in operation, it is subject to both the mass unbalance-induced inertia force and air gap unbalance-induced electromagnetic force, which is an important feature that distinguishes the motorized spindle from a mechanical spindle [7]. In addition, dynamic balance technology of machine tool depends on whether the unbalance parameters of the machine tool are obtained accurately or not. The eccentricity of grinding wheel and motor is the main factor of unbalanced fault in high-speed motorized spindle system. The eccentricity parameters should be accurately obtained before using dynamic balance technology. Using parameter identification technology to identify the eccentricity parameter of motorized spindle system is a necessary work in diagnosing unbalance fault of rotor system. However, for the electric spindle system, in engineering practice, the eccentricity parameters are mainly caused by the disk parts of the motor and the grinding wheel. It is affected by many factors, and it is difficult for the engineer to understand the approximate range of the eccentricity parameters based on experience; that is, it is difficult to obtain prior information of the eccentricity parameter. At the same time, due to the randomness of rotor geometry, material factors, and measurement response, these uncertainties, even in a relatively small case, may also lead to large deviations in the identification results of eccentricity parameter. How to reduce the influence of uncertainties such as randomness and identify quickly eccentricity parameter is a core of researches on engineering inverse problem under uncertainty.

On engineering inverse problem under uncertainty, when the sample information of uncertain parameters is sufficient, the probability density can be used to describe the uncertainty of parameters. When the prior information of the parameters identified is unknown, the maximum likelihood method is often used to identify structural parameters or load parameters. In the maximum likelihood method, the parameter value when the maximum value is taken by the likelihood function is used as the identification result, and the corresponding confidence interval is calculated [8–10]. Liu et al. [11] proposed the maximum likelihood method based on the sensitivity matrix method (SMM), which is suitable for dealing with the uncertainty inverse problems with insufficiency and imprecision in the input and/or output parameters. This proposed method can transform some known parameters into explicit solvable form and identify unknown parameters by iteration method. However, due to the unknown prior information of the parameters identified, that is, the interval of the parameters identified is unknown, and the sensitivity will bring some numerical problems in the iterative process, which affects the accuracy and reliability of the parameter identification result, and sometimes causes the identification process to diverge and the effective identification result is not obtained [12]. In this paper, Advance-Retreat Method (ARM) of the search interval is introduced to the SMM, the strict requirements of the traditional iterative numerical method for the search space are avoided, and the SMM is improved to make the iteration suitable for complex engineering optimization problems. It is possible to improve the identification effect of the SMM.

2. Improved SMM

2.1. The Basic Idea of SMM to Identify Eccentricity Parameters

The motion equation (1) of the motorized spindle system with unbalance is expressed as follows [13]:

where , , , and are the mass, damping, stiffness, and gyroscope matrices which can be obtained by the method used in [13]; , , and are the unbalanced displacement, velocity, and acceleration vectors of the motorized spindle system, respectively; and is an unbalanced force.

In the paper, the unbalanced force of the motorized spindle system is caused by the eccentricity parameters of the motor and grinding wheel. The eccentricity parameters are unknown; the unbalanced response measured by experiment is known. The identification problem is to determine the unknown eccentricity parameter from the observed unbalanced response. However, the input part parameters, such as motor length and motor radius, and output parameters such as unbalanced response are all partly known; that refers to the measurement errors. The identification problem belongs to the uncertainty inverse problem with both insufficiency and randomness in the input and output of structural systems. The unbalanced response of equation (1) is expressed as a forward solver as follows:

where is vector of input parameters, is vector of output parameters, and is the forward solver representing the translation process from input to output. The subscript u denotes unknown and the subscript k denotes known. The input vector consists of two parts; one part is the known parameter vector , which obeys the probability density distribution, such as the length and radius of the motor. The other part is the unknown parameter vector , such as the eccentricity parameters of the motorized spindle system. The output vector also contains two parts: the unknown and the known . is a known vector derived from experimental observations and subject to a certain probability density distribution, such as the unbalanced displacement.

The maximum likelihood method is a reverse method for finding the mean and confidence interval of an unknown parameter by the probability density distribution of the known parameters in the input and/or output of the system. Introducing the sensitive matrix [14], such as equation (3), the incomplete information of input and output known parameters in equation (2) can be transformed into explicit solvable form of equation (4):

where is the change in output parameters, is the change in input parameters, and m and n are the number of output and input parameters, respectively.

where , .

According to the maximum likelihood principle, equation (4) is transformed into a similar least squares form, as in

where and are variances, , ,, and are deviations of the ith parameter in the input and output known vector , .

The iteration method is used to solve the increment of unknown parameters. According to probability theory, the least square of equation (5) is the maximum likelihood solution:

Considering the nonlinearity of the structural system, the confidence interval is usually calculated by the Markov Chain Monte Carlo method (MCMC). 95% confidence interval of unknown parameters [15] is

where is the covariance of unknown parameters, is an rows and m columns matrix, representing the distance between the calculated value of parameter and the average value . is the number of calculations.

2.2. ARM for Determining Search Interval

The basic idea of sensitive matrix method and maximum likelihood method (SMM) is to transform the inverse problem of stochastic uncertainty into an explicit iterative optimization problem similar to the least squares form. The unknown parameters are identified by minimizing the deviation between the experimental and calculated values. The mean and confidence intervals of unknown parameters are gained by using input and output known parameters with normal distribution characteristics. That is, the unknown eccentricity parameters are identified by solving the optimization problem that minimizes the error between the unbalanced displacement calculated and measured at the corresponding measuring point by adjusting repeatedly the eccentricity parameters of the motorized spindle system, as in (8)

where and are the input unknowns parameters corresponding to current step i and next step i + 1, respectively. and are the unbalanced displacement numerically and accurately calculated based on equation (2), is the output known parameters obtained via measurement. and are allowable errors of iteration stopping.

Optimization methods such as Gradient Descendent [16] require step size and search direction as iterative guidance in the search interval. The SMM can obtain the increment of unknown parameters in each iteration process, but the prior information of unknown parameters is not considered, so the search interval in the iteration process is unknown. It will lead to blindness of search and affect the accuracy and reliability of parameter identification results. In this paper, the maximum likelihood function is improved by introducing the ARM for determining search interval [17, 18], the trend of “high -low-high” of the value of the objective function formed by the design points is explored, and the value of the objective function is towed toward the direction of gradual descent. Equation (8) shows that the SMM should continuously reduce the error between the calculated results and the measured results in the iteration process. This results in the comparison of the objective function values of the next step with those of the previous step, as well as the comparison of the calculated objective function values with the test results. In this way, the value of the objective function will be pulled down gradually during the iteration process, and the error between the calculated value of the objective function and the test results will be reduced or increased. Considering the unknown search interval and the problems arising from the iteration process, the ARM is introduced to improve equation (8). The unknown parameter increment obtained by the SMM is used as the step size in the iteration process and the ARM of the search interval to adjust the next design point, so that the objective function value is towed in the direction of gradual descent. When the error between the objective function value and the test result becomes larger, the design point of the previous step is selected as the initial design point and the search is restarted. The selection of design points is shown infd9

3. Solution Procedure

To sum up the above arguments, the flow chart of the present method is shown in Figure 1. The main steps are as follows:

 Step 1. Set the error constant and the probability density distribution of known parameters, such as normal distribution and , to test and obtain the unbalanced displacement at the sensitive points of sensors that are easy to arrange Step 2. Assume the unknown parameter value , combining them together with the known parameters into the forward solver to calculate the corresponding output  Step 3. Calculate the sensitivity matrix centering on the present parameters based on equation (3), and the increment of unknown parameters in equation (6) is obtained by SMM. Update , , , and  Step 4. If all these errors are within the predefined toleration, this set of parameters is considered to be the total solution of the problem, and the solution procedure ends. Otherwise, the unknown parameter is updated by the ARM equation (9), and return to Step 2. Namely, Step Calculate , . Compare and , if , move right. By increasing the step size , get . If , retreat to the left, ; then is obtained Step 5. Use equation (7) to calculate the confidence interval of unknown parameters.

Figure 1: Identification flow chart.

4. Eccentricity Parameter Identification

4.1. Description of the Numerical Example

The motorized spindle system that referred to an inverse problem of identification for the eccentricity parameters is given in the section to validate the proposed method. The model parameters of the motorized spindle system are shown in Figure 2. The diameters of the grinding wheel, the motor, and shaft are 1.4 m, 2.4 m, and 0.6 m, respectively, their lengths are 0.04 m, 0.5 m, and 1.3 m, respectively, and the front and rear bearings are installed at a distance of 0.5 m and 0.2 m from both ends of the shaft. The performance parameters of the motorized spindle system are listed in Table 1. The average air gap length is 8 mm without eccentric, the air permeability coefficient is 1.25610⁻⁶, the air gap fundamental wave magnetomotive force coefficient is 5.2, and the rotating frequency is 100 Hz.

Figure 2: Structural parameter model of spindle system.

Table 1: The performance parameters of the motorized spindle system.

Parts
Parameters
Specifications

Rotor
Elastic modulus (GPa)
211

Shear modulus (GPa)
81.2

Density (kg/m³)
7810

Grinding wheel
Mass (kg)
11000

Motor
Mass (kg)
15000

The front bearing
Stiffness coefficient (MN·m⁻¹⁾
85

The rear bearing
Stiffness coefficient (MN·m⁻¹⁾
65

According to the solution procedure described in Figure 2, the unknown input parameters are eccentricity parameters ; the known input parameters are motor length and diameter . The output parameters are the radial displacement .

If the radial displacement of this measured point is very sensitive to eccentric parameters, but it is not convenient to arrange sensors at this point, the radial displacement of this measured point is attributed to unknown parameters. Considering the influence of random errors in measurement results, the known parameters in input and output are random variables subject to a normal distribution. Taking 3% measurement deviation as an example, the known parameters are expressed as and .

4.2. Sensitivity Analysis

The optimization method implementation criterion requires that the unbalanced response of the measuring point is sensitive to the eccentricity parameter. It is necessary to evaluate the influence of the eccentricity parameter on the response data. This is accomplished using a sensitivity analysis combined with forward solver based on Finite Element Simulation Analysis to spindle system. The results of the sensitivity analysis in the x-direction and the y-direction of each point on the spindle are shown in Figure 3. It can be seen that the eccentricity is sensitive to the unbalanced response. It is feasible to identify eccentric parameters based on optimization method through minimizing the error squared of the unbalance response between the experiment results and the computational ones. The experiment unbalance response is measured by displacement sensors. Due to the electromagnetic effects of eddy current sensor, it should be installed in a suitable position. In this paper, there are two unknown parameters and two known parameters in the input parameters. In order to avoid the morbidity of the sensitive matrix in formula (3), the number of known parameters and unknown parameters in the output parameters should be the same. Figure 3 shows that the unbalanced response in x-direction at both ends of the spindle is large, and the sensor can be easily arranged, so it can be used as a measurable point. The y-direction unbalanced response of the first and second bearings is of great concern. It is not easy to arrange sensors and treat them as unmeasurable points.

(a)

(b)

Figure 3: The sensitivity analysis.

4.3. Recognition Process

Table 2 gives the corresponding solution for each iteration. According to the identification process, the parameter vector is written as follows: , , , , , and . Assume the unknown eccentricity parameter is . The mean values of motor length and diameter are . The unbalanced response analysis of the motorized spindle system is used as the forward solver. With the input information, the results of radial displacement at both ends of axis and the first and second bearings are calculated to be , respectively, as shown in the first row of Table 2. Then, an initial increment is obtained from this set of initial parameters based on SMM; based on equation (8), the unknown parameters are computed to be  μm (see the 2nd row of Table 2). And subsequently, substituting the originals of eccentricity parameters with the newly computed values of 0.0639 and 0.1676, and combining them with the known into the forward solver FEM once again, a new set of radial displacement is calculated to be  μm, respectively, as shown in the second row of Table 2. Then, the error defined in equation (8) is examined to decide if the solution procedure ends. When the stopping criterion in equation (8) is not satisfied, the unknown parameters are updated by the ARM from equation (9), the updated eccentricity parameters are 0.1916 μm and 0.5027 μm, and the solving process continues to the next iteration. After 981 times of such iterations, the maximal error of the sought parameters with respect to their measured values decreases to 0.00075352, which satisfies the stopping criterion and the calculation is completed. Table 2 gives the corresponding solutions at each of the iterations. Figure 4 shows that eccentricity parameter identification iterative process converges; the displacement error is decreasing at every iteration, which makes the eccentricity parameter to be identified approach the most advantageous. The eccentric parameters are substituted into improved maximum likelihood method from zero. With the increase of iterative steps, the simulated unbalanced response at the grinding wheel and shaft end approaches the experimental displacement measured by the displacement sensor step by step, and the displacement error decreases, about 200 steps. The change of eccentric parameters is small. The given convergence criterion is fulfilled after 981 iterations of training. The iterative evolution process of parameter identification in Figure 4 shows that the improvement measures made in this paper to the maximum likelihood method are effective.

Table 2: Eccentricity parameter identification procedure using improved maximum likelihood method.

No. of iterations
Method
Input parameters
Output parameters (radial displacement)

Unknown parameter
Known parameters
Unknown parameter
Known parameters

e₁ (μm)
e₂ (μm)
L (m)
D (m)
d₃ (μm)
d₄ (μm)
d₁ (μm)
d₂ (μm)

Initial 1
FEM
0
0
0.5
2.4
0
0
0
0

SMM
0.0639
0.1676

0.0366
0.0459
0.1055
0.0735

2
ARM
0.1916
0.5027
0.5
2.4
0.1098
0.1378
0.3164
0.2205

SMM
0.2047
0.7212

0.1174
0.1473
0.3382
0.2356

3
ARM
0.2309
1.1582
0.5
2.4
0.1325
0.1663
0.3817
0.2660

SMM
0.2874
1.2521

0.1649
0.2069
0.4750
0.3310

…
…
…
…
…
…
…
…
…
…

979
ARM
19.9342
100.1121
0.5
2.4
11.4363
14.3565
32.9509
22.9610

SMM
20.0208
100.5260

11.4860
14.4188
33.0941
23.0607

980
ARM
19.9342
100.1121
0.5
2.4
11.4363
14.3565
32.9509
22.9610

SMM
19.9707
100.4100

11.4574
14.3829
33.0115
23.0032

981
ARM
19.9342
100.1121
0.5
2.4
11.4363
14.3565
32.9509
22.9610

SMM
19.9997
100.5923

11.4740
14.4037
33.0594
23.0365

Target

20
100
0.5
2.4
11.4740
14.4037
33.0594
23.0366

(a)

(b)

(c)

(d)

Figure 4: Eccentricity parameter identification iterative process.

4.4. Identification and Discussion for Different Measurement Deviations

The influence of the three measured deviations of 1%, 3%, and 5% of the known parameters on the identification results of eccentricity parameters is considered at the motor and grinding wheel. The identification results are given in Table 3. Comparing to the different measurement deviations, the improved maximum likelihood method can reduce the influence of uncertainty factors on the recognition results, and the robustness is satisfactory. In each iteration step, equation (7) is used to calculate 95% confidence interval. Comparisons of the confidence interval of the eccentricity parameters calculated in each iteration step are shown in Figure 5. Figure 5 shows the 95% confidence interval iterative process for parameters under three errors. The eccentric parameters change a little at about 500 steps and float near the mean. It can be seen that the confidence interval increases with the increase of the measurement error.

Table 3: Parameter identification means and 95% confidence interval.

Identification parameter
Different uncertainties (%)
Mean
Confidence interval
Standard deviation

Motor eccentricity (μm)
1
19.9997
[19.9106, 20.0888]
1.5e − 5

3
19.9995
[19.7059, 20.2932]
2.5e − 5

5
20.0000
[19.4178, 20.5822]
0

Grinding wheel eccentricity (μm)
1
100.5923
[100.1375, 101.0471]
5.923e − 3

3
100.9543
[98.4023, 103.5063]
9.543e − 3

5
100.0314
[99.4962, 100.5665]
3.14e − 4

(a)

(b)

Figure 5: 95% confidence interval iterative process for parameters under three errors. (a) Motor eccentricity. (b) Grinding wheel eccentricity.

5. Conclusions

The unbalanced fault of the motorized spindle system of the grinder is caused by the coupling of several mass eccentricity parameters, so it is difficult for engineers to understand the approximate range of mass eccentricity parameters based on experience. Due to the existence of uncertain parameters such as the structure parameters and test conditions of the grinder motorized spindle system and the unknown prior information of the mass eccentricity parameters, the search calculation used in the maximum likelihood method to identify the mass eccentricity parameters of rotor unbalanced faults is faced with a large amount of calculation, and the sensitivity in the iterative process will bring some numerical problems. The Advance-Retreat Method of the search interval is introduced into the sensitive matrix method and the maximum likelihood method to improve the maximum likelihood method. The iteration in the improved maximum likelihood method is suitable for dealing with complex engineering optimization problems and avoids the stringent requirements of the traditional iterative numerical method on the search space. In this paper, the sensitive matrix method and maximum likelihood method are developed based on Advance-Retreat Method to identify the eccentricity parameters in a motorized spindle system with input size and output test response with random measurement error. The feasibility and validity of developed algorithm has been tested with numerical example. Three measurement deviations are considered in the numerical application. The comparison of parameter identification means and 95% confidence interval indicates that it is robust to the insufficient observations.

Data Availability

The authors confirm that the data supporting the conclusion of the article are shown in the relevant figures and tables in the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was a project supported by the National Natural Science Foundation of China (Grant nos. 51775180 and 51875193) and the Hunan Provincial Natural Science Foundation of China (Grant no. 2019JJ60036).

References

W. L. Ming, S. Y. Min, Y. Lei, Y. Y. Lin, and L. Jing, “Fault diagnosis for centre wear fault of roll grinder based on a resonance demodulation scheme,” Journal of Physics: Conference Series, vol. 842, no. 1, pp. 12–57, 2017.

X. Guo, B. Z. Li, and J. G. Yang, “Dynamic performance test and analysis of spindle system of high-speed grinding machine tools,” Applied Mechanics and Materials, vol. 226–228, pp. 720–724, 2012.

M. H. Tao, M. S. Fan, J. R. Galantu, and Y. Ying, “Research on the heavy block grinding wheel dynamic balancing system based on fuzzy control,” in Proceedings of the International Conference on Advanced Computer Science and Electronics Information, Beijing, China, August 2013.

F. H. Wei, J. M. Qing, and L. Heng, “Program design of an online dynamic balancing system for grinding-wheel and spindle,” in Proceedings of the 2011 IEEE International Conference on Computer Science and Automation Engineering, Shanghai, China, June 2011.

L. F. Zhang, J. Zha, C. Zou, X. Y. Chen, and Y. L. Chen, “A new method for field dynamic balancing of rigid motorized spindles based on real-time position data of CNC machine tools,” The International Journal of Advanced Manufacturing Technology, vol. 102, no. 5–8, pp. 1181–1191, 2019.

Zhang and Y. J. Cai, “A new double-face online dynamic balance device and its control system for high speed machine tool spindle,” Journal of Vibration and Control, vol. 22, no. 4, pp. 1037–1048, 2016.

J. F. Liu, F. Li, J. H. Yong, L. Tao, and P. Zhang, “Investigation of spindle-tool assembly dynamics for optical grinding motorized spindles,” Optik, vol. 216, pp. 1–10, 2020.

X. L. Gang, “Study of constraint and impact of a nuisance parameter in maximum likelihood method,” Journal of Physics G: Nuclear and Particle Physics, vol. 46, no. 8, pp. 1–22, 2019.

F. H. Huan, W. Y. Rong, and J. X. Hua, “A maximum likelihood method for estimating probabilistic strain amplitude-fatigue life curves,” Acta Mechanica Solida Sinica, vol. 31, no. 1, pp. 80–93, 2018.

R. V. Ermakov, A. A. Seranova, A. A. L’vov, and D. M. Kalikhman, “Optimal estimation of the motion parameters of a precision rotating stand by maximum likelihood method,” Measurement Techniques, vol. 62, no. 2, pp. 139–146, 2019.

G. R. Liu, Y. G. Xu, and Z. P. Wu, “Total solution for structure mechanics problems,” Computer Methods in Applied Mechanics and Engineering, vol. 191, pp. 989–1012, 2011.

Z. Zhou, L. L. Zhou, Y. Ren, D. Z. Ma, and S. C. Fan, “A new type of solving method for reliability model parameter optimization,” Applied Mechanics and Materials, vol. 291–294, pp. 1895–1900, 2013.

F. Michael, E. T. P. John, D. G. Seamus, and W. L. Arthur, Dynamics of Rotating Machines, Cambridge University Press, Cambridge, UK, 2015.

W. G. William, “Review of parameter identification procedures in groundwater hydrology,” Water Resources Research, vol. 22, pp. 95–103, 1986.

R. C. Aster, B. Borchers, and C. H. Thurber, Parameter Estimation and Inverse Problems, Academic Press, Cambridge, MA, USA, 2005.

Q. H. Zhou, Y. R. Zhang, F. X. Xu, Y. Geng, and X. D. Sun, “An improved trust region method for unconstrained optimization,” Science China Mathematics, vol. 56, no. 2, pp. 425–434, 2013.

D. Auroux and J. Blum, “A nudging-based data assimilation method: the back and forth nudging (BFN) algorithm,” Nonlinear Processes in Geophysics, vol. 13, no. 1, pp. 305–319, 2008.

S. Q. Lu, K. X. Ying, and K. Jian, “The application of a method of advance and retreat in multidimensional nonlinear constrained optimization,” Journal of Taiyuan Heavy Machinery Institute, vol. 22, no. 3, pp. 200–203, 2001.