A Method for Solving Ill-Conditioned Nonlinear Least Squares Problems and Its Application in Image Distortion Correction Using Self-Calibration

: In this study, the ill-conditioning of the iterative method for nonlinear models is discussed. Due to the e ﬀ ectiveness of ridge estimation for ill-conditioned problems and the lack of a combination of the H-K formula with the iterative method, the improvement of the LM algorithm is studied in this paper. Considering the LM algorithm for ill-conditioned nonlinear least squares, an improved LM algorithm based on the H-K formula is proposed for image distortion correction using self-calibration. Three ﬁ nite di ﬀ erence methods are used to approximate the Jacobian matrix, and the H-K formula is used to calculate the damping factor in each iteration. The Brown model, quadratic polynomial model and Fourier model are applied to the self-calibration, and the improved LM algorithm is used to solve the model parameters. In the simulation experiment of space resection of a single image, we evaluate the performance of the LM algorithm based on the gain ratio (LM h ) and the improved LM algorithm based on the H-K formula (LM HK ), and the accuracy of di ﬀ erent models and algorithms is compared. A ridge trace analysis is carried out on the damping factor to illustrate the e ﬀ ects of the improved algorithm in handling ill-conditioning. In the second experiment, the improved algorithm is applied to measure the diameter of a coin using a single camera. The experimental results show that the improved LM algorithm can reach the same or higher accuracy as the LM h algorithm, and it can weaken the ill-conditioning to a certain extent and enhance the stability of the solution. Meanwhile, the applicability of the improved LM algorithm in self-calibration is veri ﬁ ed.


Introduction
Self-calibration is a type of analytical calibration method that describes the system error of a camera by the distortion models in the adjustment model.The self-calibration method is convenient and flexible, and it is widely used in camera calibration.According to the coordinates of the reference points, the additional parameters in the distortion model are solved to compensate for the influence of system error on the results.In order to improve the accuracy of the self-calibration, the camera parameters and distortion parameters need to be solved accurately [1].In the field of self-calibration, a lot of research about distortion models has been done.Based on many experiments and analyses, the Brown model was introduced into the analytical calibration method for distortion correction in 1971 [2].Currently, the Brown model and its improved model are still commonly used in photogrammetry [3].To address the problem of inconsistency between the optical center and the geometric center of the imaging system, self-calibration based on a simplified Brown model is proposed [4].David et al. show that the complex distortion model can improve the accuracy of self-calibration by an experiment of three-dimensional reconstruction [5].Based on the mathematical approximation theory, a distortion model based on the Fourier Axioms 2024, 13, 209 2 of 16 series was proposed for self-calibration purposes [6] and verified in an experiment of simulated distortion data fitting [7].The Brown model and Fourier model are applicable in image processing, pattern classification and scene analysis in the field of computer vision.With the development of the self-calibration method, it is worth considering whether the distortion model has a strong correlation or overparameterization, which leads to ill-conditioning.Unreasonable distortion models and parameters are one of the important causes of ill-conditioning.In the self-calibration method, the quality of feature points is another important cause of ill-conditioning.It depends on many factors, including the number, accuracy and distribution of feature points [8].Insufficient feature extraction will fail to fully reflect the role of each parameter, or unilaterally highlight the roles of some parameters.Furthermore, it is easy to cause ill-conditioning in the solution of parameters.Most of the above research avoids ill-conditioning from the aspect of model selection and feature extraction, but few research studies are carried out from the aspect of parameter iteration.In this work, the method for solving the self-calibration model is studied from the aspect of a numerical iterative method.
The solutions of ill-conditioned equations are very sensitive to small parameter value perturbations, showing poor numerical stability.For ill-conditioned equations, it is difficult to obtain accurate and reliable parameter estimates, which severely affects the accuracy and quality of data processing.For the solution of ill-conditioned problems, a variety of biased estimation methods and improved methods have been proposed to improve the quality of parameter estimates, such as regularization, truncated singular value decomposition and ridge estimation [9][10][11][12][13][14][15][16].The key to solving ill-conditioned problems by the regularization method is the selection of the stabilization functional and regularization parameter.The stabilization functional is constructed based on prior information of the model parameters, which can effectively improve the structure of the model and make solving it feasible [17].When prior information cannot be obtained, the stabilization functional is often expressed as a 2-norm constraint on model parameters [18] and then the ridge estimation of the parameters is derived.Ridge estimation is a special form of regularization that regards the identity matrix as a regularization matrix and then improves the reliability and stability of parameter estimation by reasonably selecting the regularization parameter.In ridge estimation, the regularization parameter is also called the ridge parameter.The common methods to determine the ridge parameters are the ridge trace method and Hoerl-Kennard (H-K) formula.After determining the appropriate ridge parameters, ridge estimation can effectively reduce the mean square error by properly modifying the ill-conditioned matrix [18].The adjustment criterion of the ill-conditioned uncertainty model based on ridge estimation can effectively suppress the influence of ill-conditioning [19].The common strategy for solving nonlinear parameters is to linearize the nonlinear model and solve it by an iterative method, and the Levenberg-Marquardt (LM) algorithm [20,21] is commonly used.An LM algorithm based on the gain ratio [22] can weaken the ill-conditioning by determining a damping factor, and it has been widely used in nonlinear optimization [23].However, there are some problems in the existing research: on the one hand, the damping factor determined by this algorithm is usually large, which greatly changes the structure of the Jacobian matrix; on the other hand, although ridge estimation can effectively improve the stability of the solution, there are few research studies on the combination of ridge estimation and the LM algorithm.In response to the shortcomings of the above research work, the combination of the method for determining the ridge parameters and the LM algorithm is studied in this paper.
In this work, with the LM algorithm steps as the basic framework, an improved LM algorithm is proposed.Firstly, the ill-conditioning of the iterative method for nonlinear models is discussed.Then, according to the parameter estimation criterion of the LM algorithm, the H-K formula is combined with the LM algorithm to calculate the damping factor in each iteration.Finite difference methods are used to approximate the Jacobian matrix, and finally, we present the algorithm steps of the improved LM algorithm based on the H-K formula and finite difference.The Brown model, quadratic polynomial model and Fourier model are discussed.In the numerical experiments, the three distortion models are applied to self-calibration, and the improved LM algorithm is used to solve the model parameters.In the simulation experiment of space resection of a single image, the accuracy and performance of different models and algorithms on the solution of parameters are evaluated, and a ridge trace analysis is carried out on the damping factor to illustrate the effects of the improved algorithm in handling ill-conditioning.The applicability of the improved algorithm in practical problems is verified by the measurement of a coin diameter using a single camera.The experimental results show that the improved LM algorithm can reach the same or higher accuracy as the LM algorithm based on the gain ratio, and it can weaken the ill-conditioning to a certain extent and enhance the stability of the solution under the condition of changing the matrix structure as little as possible.The improved LM algorithm proposed in this paper provides a new method for solving self-calibration model parameters, and it also provides a new idea for solving ill-conditioned nonlinear least squares.

Ill-Conditioned Problems in the Iterative Method
For a nonlinear model, the Taylor formula is usually used to transform it into a linear form, and then the iterative method is used to solve it.We consider a nonlinear model expand it by the Taylor formula and take only the first-order term: where B is the Jacobian matrix consisting of the first-order partial derivatives of F(x) at x 0 .We rewrite it as an error equation, Under the condition of equal-precision independent observation, the iterative formula of the Gauss-Newton method can be obtained according to the least squares criterion: where k is the number of iterations, B = (B 1 , • • • , B n ) is a matrix of order m × n and its rank is r.In practical problems, usually m > n, and r B T B ≤ n < m can be obtained according to the properties of matrix rank.If B is rank-deficient, so is B T B. If it is of full rank but strong multicollinearity exists in its columns: where ω = (ω 1 , ω 2 , • • • , ω n ) T is the set of eigenvectors corresponding to the eigenvalues of B T B. The Jacobian matrix is ill-conditioned at this moment.Rank deficiency or illconditioning of the matrix occurs easily when the Gauss-Newton method is used.

LM Method Based on Hoerl-Kennard Formula and Finite Differences
According to the analysis in the previous section, rank deficiency or ill-conditioning occurs easily when the Gauss-Newton method is used, resulting in no unique or stable solution in the normal equation.Therefore, to avoid the ill-posedness of the normal equation, it is necessary to impose new constraints on the parameters.L2 regularization, namely, ridge regression, is commonly used.When L2 regularization is introduced, the parameter estimation criterion and its estimation formula are The iterative formula for the LM method is presented as Formula (4).This is an improved algorithm based on the Gauss-Newton method, and it is also an iterative method without a specific line search for calculating the step size, which is determined by the damping factor µ [22].The initial value of the damping factor can be set to µ = τ • max{B ii }.B ii denotes the diagonal elements of B T B, and τ is usually a small value.If the selected initial values of the parameters are believed to be good approximations to the estimated values, it can be set to τ = 10 −6 ; otherwise, it can be set to τ = 10 −3 or τ = 1.Then, the damping factor is increased or decreased according to the gain ratio (h) in the iteration.The algorithm steps are as follows: Step 1: Given the initial parameter value x k , the Jacobian matrix and initial value of the damping factor are calculated, and the convergence criteria are set: the threshold ε 1 of the gradient g, the threshold ε 2 of the error difference and the maximum number of iterations k max , along with k = 0 and υ = 2, are set.
The LM algorithm introduces the damping factor and identity matrix based on the Gauss-Newton method, which effectively resolves the rank deficiency of B T B and avoids the ill-conditioning caused by multicollinearity.When the Jacobian matrix is ill-conditioned, there is at least one eigenvalue close to 0 in B T B, while the degree of the eigenvalue close to 0 in B T B + µI will be improved, which will weaken the ill-conditioning.In the LM algorithm, the Jacobian matrix determines the descent direction of the function, and the damping factor affects both the direction and step size.Therefore, the key to solving the nonlinear model by the LM algorithm is to calculate the damping factor and the Jacobian matrix.Formula (4) shows that the calculation method of x in the iterative formula can be regarded as a ridge estimation derived from the regularization principle.Therefore, different from the above method based on the gain ratio, the damping factor can be determined by the method of selecting the ridge parameter.

Methods for Selecting the Ridge Parameter
The damping factor is a number greater than 0, and different values lead to different ridge estimates.When µ → ∞ , x → 0 , and the solution is not relevant.Therefore, the value should be chosen to be as small as possible to solve the ill-posed problem using the proximal well-posed problem.The main methods of selecting the ridge parameter are ridge trace and the Hoerl-Kennard formula.
(1) Ridge trace: Each component xi (i = 1, 2, • • •) of the correction value of x is regarded as a function of the damping factor.When the damping factor changes between [0, ∞), several curves can be drawn in the rectangular coordinate system o − µ x, which are called ridge trace curves, as shown in Figure 1.With the increase in the ridge parameter, the model parameter estimates gradually stabilize and reach a stable state around µ * .The ridge trace method selects the smallest µ * that makes the change stable.To determine the appropriate value according to the ridge trace curves, the following three criteria must be generally satisfied: first, the ridge estimates of model parameters in the equation are roughly stable; second, the ridge estimation can make the parameter values reasonable; and third, the sum of the squared residuals does not increase substantially.The main disadvantage of this method is the lack of a strict theoretical basis.Determining the value of the damping factor using ridge trace curves is subjective, but this subjectivity can also organically combine qualitative analysis with quantitative analysis.However, the numbers of iterations may be large.If the ridge parameter in each iteration is selected by the ridge trace curves, the workload may be large.Therefore, this method is mainly used for ridge trace analysis rather than directly determining the ridge parameter in the iteration.
ridge trace and the Hoerl-Kennard formula.
(1) Ridge trace: Each component of the correction value of x is regarded as a function of the damping factor.When the damping factor changes between , several curves can be drawn in the rectangular coordinate system x ô μ − , which are called ridge trace curves, as shown in Figure 1.With the increase in the ridge parameter, the model parameter estimates gradually stabilize and reach a stable state around * μ .The ridge trace method selects the smallest * μ that makes the change stable.To determine the appropriate value according to the ridge trace curves, the following three criteria must be generally satisfied: first, the ridge estimates of model parameters in the equation are roughly stable; second, the ridge estimation can make the parameter values reasonable; and third, the sum of the squared residuals does not increase substantially.The main disadvantage of this method is the lack of a strict theoretical basis.Determining the value of the damping factor using ridge trace curves is subjective, but this subjectivity can also organically combine qualitative analysis with quantitative analysis.However, the numbers of iterations may be large.If the ridge parameter in each iteration is selected by the ridge trace curves, the workload may be large.Therefore, this method is mainly used for ridge trace analysis rather than directly determining the ridge parameter in the iteration.(2) Hoerl-Kennard (H-K) formula: The H-K formula is a common method for determining the ridge parameter according to the canonical form of the error equation.In order to analyze the problem, the canonical form of the error equation is introduced.
Considering the symmetry of matrix B B T , the Formula ( 5) can be obtained according to matrix theory where Ω is an orthogonal matrix; Λ is diagonal matrix and the diagonal element is the eigenvalue of B B T .Formula (6) can be obtained by identity transformation of the original error equation ( ) (2) Hoerl-Kennard (H-K) formula: The H-K formula is a common method for determining the ridge parameter according to the canonical form of the error equation.In order to analyze the problem, the canonical form of the error equation is introduced.
Considering the symmetry of matrix B T B, the Formula ( 5) can be obtained according to matrix theory where Ω is an orthogonal matrix; Λ is diagonal matrix and the diagonal element is the eigenvalue of B T B. Formula ( 6) can be obtained by identity transformation of the original error equation Formula ( 7) is the canonical form of the error equation; Y is the canonical parameter, and its estimation is According to x = ΩY, the parameter estimation of the original error equation can be obtained.In canonical form, the mean square error (MSE) of the canonical parameter is Axioms 2024, 13, 209 6 of 16 Hoerl and Kennard proposed the H-K formula [14] based on the canonical form of the error equation.(ω 1 , • • • , ω n ) are assumed to be the eigenvectors corresponding to the eigenvalues (λ Then µ is calculated by the H-K formula The mean square error of parameter estimation is It can be seen that, when µ is large, although this estimation method can overcome the multicollinearity, it will also introduce more bias.The determination of µ according to the H-K formula is related to the eigenvalues.An ill-conditioned matrix contains very small eigenvalues, so a smaller value can be obtained according to the H-K formula, so that the parameter estimation can be stable as soon as possible with less bias.

Finite Difference Form of the Jacobian Matrix
The Jacobian matrix can be calculated according to the differentiation rule of functions of several variables, as expressed in Formula (1).However, in practical problems, due to the large number of observation equations and the associated large-scale matrices, this method may require a huge amount of calculation.In a common numerical method, the finite difference approximates the differential of the function by the corresponding function value after parameter discretization.This method of approximating the derivative by the difference quotient can effectively reduce the calculational burden.In this paper, the forward difference method (FDM), backwards difference method (BDM) and central difference method (CDM) are adopted to approximate the Jacobian matrix.According to the definition of forward difference, ∆F i = F i x j + δ − F i x j , and δ is the difference step.F i x j + δ is expanded by the Taylor formula, and only the first-order term is considered: Similarly, the backwards difference and central difference are Compared with the forward difference and backwards difference, the accuracy of approximating the Jacobian matrix by the central difference is higher, and the calculational burden is multiple times larger.
The steps of the improved LM algorithm based on the H-K formula and finite differences are as follows: Step 1: Given an initial parameter value x k , set convergence criteria: the threshold ε 1 of gradient g, the threshold ε 2 of the error difference and the maximum number of iterations k max , along with k = 0, are set.
Step 2: Given the difference step, ∂F i ∂x j ) is calculated by the finite difference method, the matrix B k is obtained, and Step 3: The equation , and the (k + 1)th parameter estimate x k+1 = x k + xk is obtained.
Step 4: is the optimal parameter estimate, and the iteration is terminated; otherwise, we set k = k + 1 and return to Step 2.

Distortion Models
In photogrammetry, since the beam of the imaging system within the field of view does not strictly meet the ideal center projection, the actual image points produce a position error, which is called distortion.Considering the distortion of image points, the actual image coordinates can be regarded as the sum of the ideal image coordinates and distortion, which is expressed by a collinearity equation with additional parameters: where (u, v) are the actual image coordinates, (X, Y, Z) are the corresponding ground control point coordinates, (u 0 , v 0 , f ) are the elements of interior orientation, (X S , Y S , Z S ) are the translation elements of exterior orientation, R ij (i, j = 1, 2, 3) are the direction cosine composed of the angle elements of exterior orientation (φ 1 , φ 2 , φ 3 ) and (∆u, ∆v) are the distortions, which are generally functions of the image coordinates.Formula ( 16) is also called a self-calibration model.Distortion models mainly include physical models and mathematical models.The Brown model, polynomial model and Fourier model are mainly described in this paper.

Brown Model
The Brown model is a commonly used physical model, which was originally designed for large-area film cameras.Various forms of distortion that occur in camera imaging are considered.In practical applications, the parameters of the Brown model need to be selected according to the camera imaging characteristics.With the application of a digital camera in aerial photography, the model has been simplified into radial distortion and decentering distortion.The Brown model is expressed as where s = u 2 + v 2 , and K i (i = 1, 2, 3) and P i (i = 1, 2) are radial distortion parameters and decentering distortion parameters.The distortion of the image edges can be described only by the radial distortion and decentering distortion when using a digital camera.In practical applications, radial distortion usually remains K 1 , K 2 .

Polynomial Model and Fourier Model
In addition to the physical model, the distortion model can also be established from the mathematical point of view (in this section, a is used to represent the distortion coefficient of different mathematical models).The orthogonal polynomial in the mathematical model can effectively reduce the correlation between parameters and improve the stability of the solution.The polynomial model is a commonly used distortion correction model, which is expressed mathematically as where t is the order of the polynomial and a ij are the parameters to be determined.However, a higher-order polynomial is also prone to overparameterization, which will not only increase the calculational work required to obtain the numerical solution but also lead to the instability of the solution.Therefore, the order should not be too high in practical applications.The quadratic polynomial (QP) model is often selected as a distortion model.Considering the correlation of the coefficients in u and v, the quadratic orthogonal polynomial model is According to the Weierstrass second approximation theory, a binary Fourier series orthogonal polynomial model represented by the image coordinates can be obtained [6]: where a i are the parameters to be determined and COS i,j = cos(iu + jv), SIN i,j = sin(iu + jv), u = (u−width/2) height π (width and height are the image width and height).

Numerical Experiments and Analysis
The first experiment involves a simulation experiment of the space resection of a single image based on the collinearity equation with additional parameters, and the second involves measuring the diameters of coins using a single camera.The improved LM algorithm is used in the numerical experiments, and the convergence conditions are set as follows: The Jacobian matrix in each iteration is approximated by the finite difference methods introduced in Section 2. The H-K formula is used to calculate the damping factor, which is compared with the method according to the gain ratio.The experimental environment is MATLAB R2021a running on a 1.80 GHz PC with Windows 7.

Space Resection of a Single Image Based on the Collinearity Equation with Additional Parameters
In this experiment, the image points in a single aerial image are simulated for space resection.It is assumed that the local coordinate system is a North-East-Down (NED) coordinate system, the flight altitude is 50 m, the design focal length of the camera is 9 mm, the pixel size is 2.4 µm/px and the image width and height are 5472 pixels and 3648 pixels.At the moment of exposure, the focal length is 8.9 mm, the coordinates of the principal points are u 0 = 2737.2px, v 0 = 1827.4px and the elements of exterior orientation are . This experiment is carried out according to the following steps: Step 1: Simulate the data.The values of X and Y of the ground points are centered on the plane coordinates of the camera station, which are uniformly distributed in the range of [X S − 10, X S + 10] and [Y S − 10, Y S + 10].And the values of Z are in the range of [−1 , 1]  since the origin of the NED coordinate system is set on the ground.A total of 120 image coordinates, and the corresponding ground point coordinates, are simulated, and Gaussian noises are added to the image coordinates as observations.
Step 2: Select the distortion models.The Brown model, quadratic polynomial (QP) model and Fourier model are regarded as distortion models, which are added to the collinearity equation to form the self-calibration models.
Step 3: Initialization.The initial values of the angle elements of the exterior orientation are set to φ 1 = φ 2 = φ 3 = 0 • , the initial value of Z S is Z 0 S = −50 m and the initial values of X S and Y S are calculated according to the following formula where m is the number of ground points.The initial values of the elements of the interior orientation are u 0 0 = width/2, v 0 0 = height/2, f 0 = 9mm, and the initial values of the additional parameters are set to 0. The normalized image coordinates are substituted into the self-calibration models composed of the three distortion models.
Step 4: Solve the parameters.The LM algorithm based on the forward difference, backwards difference and central difference methods with the gain ratio and the H-K formula (LM FDM+h , LM BDM+h , LM CDM+h , LM FDM+HK , LM BDM+HK and LM CDM+HK ) are used to solve the parameters.The elements of the interior orientation, and additional parameters, are determined while solving for the elements of the exterior orientation.
The experiment evaluates the performance of the algorithms from the following aspects: the accuracy is compared using the sum of squared residuals (SSR), the maximum residuals of the image points, the reprojection errors (REs) and the true errors of the parameters; the efficiency is compared using the number of iterations and the running time; the influence on the ill-conditioning is compared using the condition number; the stability of the solution is analyzed using the ridge trace curve.Table 1 presents the SSR and the maximum residuals of each algorithm at the optimal solution.It can be seen from Table 1 that, for the same algorithm, the fitting accuracy of different models for image coordinate observations is high, and the SSR of the Fourier model is generally the smallest, reaching an accuracy of 10 −6 .Because the three difference methods are consistent in the approximation of the Jacobian matrix, the difference between the results obtained by the three difference methods is small.In contrast, the performance of the LM algorithm based on CDM is better.For the maximum residuals obtained by the LM algorithm based on the gain ratio (LM h ) and the improved LM algorithm based on the H-K formula (LM HK ), there is no significant law, which means that the LM HK algorithm has a poor fitting effect on image points with a large error.However, the SSR corresponding to the LM HK algorithm is generally smaller, reaching an accuracy of 10 −6 , indicating that the LM HK algorithm has a higher fitting accuracy for the observations.According to the introduction of the H-K formula in Section 2, the H-K formula can determine a small damping factor and less bias is introduced, so that the solution of the parameters can reach a higher accuracy.This is proved by the comparison results of the LM h algorithm and LM HK algorithm.Through the above analysis, another discovery obtained from the table is explained: for the same model, the LM CDM+HK algorithm can reach the same or higher accuracy as other algorithms.The LM algorithm based on CDM shows a better performance in Table 1.Therefore, in order to show the accuracy of the improved LM algorithm more clearly and intuitively, Figure 2 and Table 2 present the distribution and the root mean square error (RMSE) of reprojection errors (REs) obtained by the LM CDM+h and LM CDM+HK corresponding to each model.In Figure 2, the dotted line indicates the position where the RE is 1 pixel.It can be seen from Figure 2 that the maximum of the REs for all methods is less than 2 pixels.For the same model, the number of image points with an RE greater than 1 pixel in the results of the LM CDM+HK algorithm is generally less than that of the LM CDM+h algorithm, and a consistent conclusion can be obtained from Table 2: the RMSE corresponding to the LM CDM+HK algorithm is the same as or smaller than that of the LM CDM+h algorithm, which indicates that the LM CDM+HK algorithm can reach the same or higher accuracy as the LM CDM+h algorithm.For the same algorithm, the RMSE of the Brown model is the closest to 1 pixel, which is the worst accuracy of all the methods.The Brown model considers radial distortion and decentering distortion without considering other forms of distortion.The QP model and Fourier model, established from the perspective of function approximation, can accurately fit the unknown distortion in the image.The accuracy of the Fourier model solved by the LM CDM+HK algorithm is less than 0.8 pixels, which is the highest accuracy of all methods.Table 3 presents the true errors of the parameters.It can be seen from Table 3 that the parameter estimates of LM CDM+h and LM CDM+HK have an equivalent accuracy.From the true error of the exterior orientation elements, it can be found that, compared with the LM CDM+h , the true error obtained by the LM CDM+HK is generally smaller and the result is closer to the true value.The LM algorithm based on CDM shows a better performance in Table 1.Therefore, in order to show the accuracy of the improved LM algorithm more clearly and intuitively, Figure 2 and Table 2 present the distribution and the root mean square error (RMSE) of reprojection errors (REs) obtained by the LMCDM+h and LMCDM+HK corresponding to each model.In Figure 2, the dotted line indicates the position where the RE is 1 pixel.It can be seen from Figure 2 that the maximum of the REs for all methods is less than 2 pixels.For the same model, the number of image points with an RE greater than 1 pixel in the results of the LMCDM+HK algorithm is generally less than that of the LMCDM+h algorithm, and a consistent conclusion can be obtained from Table 2: the RMSE corresponding to the LMCDM+HK algorithm is the same as or smaller than that of the LMCDM+h algorithm, which indicates that the LMCDM+HK algorithm can reach the same or higher accuracy as the LMCDM+h algorithm.For the same algorithm, the RMSE of the Brown model is the closest to 1 pixel, which is the worst accuracy of all the methods.The Brown model considers radial distortion and decentering distortion without considering other forms of distortion.The QP model and Fourier model, established from the perspective of function approximation, can accurately fit the unknown distortion in the image.The accuracy of the Fourier model solved by the LMCDM+HK algorithm is less than 0.8 pixels, which is the highest accuracy of all methods.Table 3 presents the true errors of the parameters.It can be seen from Table 3 that the parameter estimates of LMCDM+h and LMCDM+HK have an equivalent accuracy.From the true error of the exterior orientation elements, it can be found that, compared with the LMCDM+h, the true error obtained by the LMCDM+HK is generally smaller and the result is closer to the true value.Figure 3 shows the iterative changes in the SSR for different models solved by LM CDM+h and LM CDM+HK .In order to clearly show the difference, the first iteration has been removed.Table 4 presents the number of iterations and the running time of each algorithm and model.According to Table 4, the number of iterations and the running time of the LM HK algorithm are generally less than those of the LM h algorithm, indicating that a high fitting accuracy can be obtained by the improved algorithm, with a higher iterative efficiency.According to the running time, compared with the LM CDM+h algorithm, the efficiency of the LM CDM+HK algorithm is improved by 64%, 55% and 33% corresponding to the three models.The algorithms using CDM to approximate the Jacobian matrix require more time, which is consistent with the principle of central difference.Compared with FDM and BDM, CDM needs to calculate one more approximate partial derivative of each variable in the iteration within the difference step range, and the calculational burden is the largest.As shown in Figure 3, the LM CDM+h algorithm makes the SSR of the Brown model reach a stable state after five iterations, which indicates that the descending speed is slower.However, the LM CDM+HK algorithm makes the SSR of all models stable after three iterations.Figure 3 shows the iterative changes in the SSR for different models solved by LMCDM+h and LMCDM+HK.In order to clearly show the difference, the first iteration has been removed.Table 4 presents the number of iterations and the running time of each algorithm and model.According to Table 4, the number of iterations and the running time of the LMHK algorithm are generally less than those of the LMh algorithm, indicating that a high fitting accuracy can be obtained by the improved algorithm, with a higher iterative efficiency.According to the running time, compared with the LMCDM+h algorithm, the efficiency of the LMCDM+HK algorithm is improved by 64%, 55% and 33% corresponding to the three models.The algorithms using CDM to approximate the Jacobian matrix require more time, which is consistent with the principle of central difference.Compared with FDM and BDM, CDM needs to calculate one more approximate partial derivative of each variable in the iteration within the difference step range, and the calculational burden is the largest.As shown in Figure 3, the LMCDM+h algorithm makes the SSR of the Brown model reach a stable state after five iterations, which indicates that the descending speed is slower.However, the LMCDM+HK algorithm makes the SSR of all models stable after three iterations.

Models
Brown QP Fourier  To further analyze the performance of the improved LM algorithm based on the H-K formula, a ridge trace analysis of the different methods is carried out.Figure 4 shows the ridge trace curves of x changing with the damping factor at the optimal solution.Table 5 shows the condition number (C/C 0 ) of the normal matrix with or without the damping factor.It can be seen from Table 5 that both the LM CDM+h algorithm and the LM CDM+HK algorithm can reduce the condition number of the matrix and weaken the ill-conditioning, and the effect of the LM CDM+h algorithm is more significant.Especially for the Brown model and the QP model, the condition number is reduced to 26.823 and 32.576, which can be considered to mean that the matrix is well-posed.However, it can be seen from Figure 4 that the order of magnitude of the damping factor for the Brown model, QP model and Fourier model should be 10 −9 , 10 −14 and 10 −11 .According to the selection principle of the damping factor and Table 5, the actual value determined by the LM CDM+h algorithm is generally too large, while the values of damping factor determined by the LM CDM+HK algorithm are relatively consistent with the trend of the ridge trace curves, which is closer to the result of the ridge trace analysis.Therefore, it is found that the LM CDM+h algorithm weakens the ill-conditioning of the normal matrix more significantly, since this algorithm changes the structure of the normal matrix by selecting a larger damping factor.Upon the premise of changing the structure of the normal matrix as little as possible, the LM CDM+HK algorithm can weaken the ill-conditioning to a certain extent and make the solution stable by selecting a smaller damping factor.To further analyze the performance of the improved LM algorithm based on the H-K formula, a ridge trace analysis of the different methods is carried out.Figure 4 shows the ridge trace curves of x ˆ changing with the damping factor at the optimal solution.
Table 5 shows the condition number (C/C0) of the normal matrix with or without the damping factor.It can be seen from Table 5 that both the LMCDM+h algorithm and the LMCDM+HK algorithm can reduce the condition number of the matrix and weaken the illconditioning, and the effect of the LMCDM+h algorithm is more significant.Especially for the Brown model and the QP model, the condition number is reduced to 26.823 and 32.576, which can be considered to mean that the matrix is well-posed.However, it can be seen from Figure 4 that the order of magnitude of the damping factor for the Brown model, QP model and Fourier model should be  5, the actual value determined by the LMCDM+h algorithm is generally too large, while the values of damping factor determined by the LMCDM+HK algorithm are relatively consistent with the trend of the ridge trace curves, which is closer to the result of the ridge trace analysis.Therefore, it is found that the LMCDM+h algorithm weakens the ill-conditioning of the normal matrix more significantly, since this algorithm changes the structure of the normal matrix by selecting a larger damping factor.Upon the premise of changing the structure of the normal matrix as little as possible, the LMCDM+HK algorithm can weaken the ill-conditioning to a certain extent and make the solution stable by selecting a smaller damping factor.

Measurement of a Coin Diameter Using a Single Camera
Experiment 3.1 shows that the LM algorithm based on CDM presents a better performance in solving the self-calibration model.Therefore, LMCDM+h and LMCDM+HK are used to measure the diameters of coins using a single camera in this experiment.Nine images of a calibration pattern are taken from different angles to calibrate the camera by Zhang's calibration method.Taking the calibration result as the initial value and using the detected point of the last image, LMCDM+h and LMCDM+HK are combined with the Brown model, QP model and Fourier model to estimate the model parameters and then measure the diameters of the coins.The detected points and the coins are shown in Figure 5.The SSR and the maximum residuals in pixels of each algorithm and model are shown in Table 6.It can be seen from Table 6 that the accuracy of all the models and algorithms is high, since good initial values are obtained.The SSR and maximum residuals of the LMCDM+h algorithm and LMCDM+HK algorithm are consistent, which indicates that the fitting accuracy is equivalent.The orders of magnitude of SSR and the maximum for the Brown model are

Measurement of a Coin Diameter Using a Single Camera
Experiment 3.1 shows that the LM algorithm based on CDM presents a better performance in solving the self-calibration model.Therefore, LMCDM+h and LMCDM+HK are used to measure the diameters of coins using a single camera in this experiment.Nine images of a calibration pattern are taken from different angles to calibrate the camera by Zhang's calibration method.Taking the calibration result as the initial value and using the detected point of the last image, LMCDM+h and LMCDM+HK are combined with the Brown model, QP model and Fourier model to estimate the model parameters and then measure the diameters of the coins.The detected points and the coins are shown in Figure 5.The SSR and the maximum residuals in pixels of each algorithm and model are shown in Table 6.It can be seen from Table 6 that the accuracy of all the models and algorithms is high, since good initial values are obtained.The SSR and maximum residuals of the LMCDM+h algorithm and LMCDM+HK algorithm are consistent, which indicates that the fitting accuracy is equivalent.The orders of magnitude of SSR and the maximum for the Brown model are 3 10 − and 2 10 − .For the same algorithm, the accuracy of the Brown model is lower, while the accuracy of the QP model and Fourier model are higher, which indicates that the image has other forms of distortion except radial distortion and decentering distortion, and the mathematical model can compensate for these distortions more effectively.Table 7 presents the number of iterations and running time.It can be seen from Table 7 that the LMCDM+HK algorithm has a higher iteration efficiency than LMCDM+h.This is consistent with the conclusion in Section 3.1.According to the running time, the efficiency of the LMCDM+h algorithm is improved by 12%, 28% and 30% corresponding to the three models.For the same algorithm, the QP model and Fourier model need fewer iterations, while the Brown  It can be seen from Table 6 that the accuracy of all the models and algorithms is high, since good initial values are obtained.The SSR and maximum residuals of the LM CDM+h algorithm and LM CDM+HK algorithm are consistent, which indicates that the fitting accuracy is equivalent.The orders of magnitude of SSR and the maximum for the Brown model are 10 −3 and 10 −2 .For the same algorithm, the accuracy of the Brown model is lower, while the accuracy of the QP model and Fourier model are higher, which indicates that the image has other forms of distortion except radial distortion and decentering distortion, and the mathematical model can compensate for these distortions more effectively.Table 7 presents the number of iterations and running time.It can be seen from Table 7 that the LM CDM+HK algorithm has a higher iteration efficiency than LM CDM+h .This is consistent with the conclusion in Section 3.1.According to the running time, the efficiency of the LM CDM+h algorithm is improved by 12%, 28% and 30% corresponding to the three models.For the same algorithm, the QP model and Fourier model need fewer iterations, while the Brown model requires more iterations.However, the running time of the Brown model is the least, since the model has the least number of distortion parameters.

Discussion
In this work, we study the combination of the H-K formula and the LM algorithm, and an analysis of the LM algorithm and its improvement are carried out in the numerical experiment.Consistent with the previous study [20][21][22][23], it is found that the LM algorithm based on the gain ratio (LM h ) is effective in handling nonlinear least squares problems, and it can weaken the ill-conditioning of the Jacobian matrix by determining the damping factor.However, from the perspective of ridge estimation, the damping factor determined by this algorithm is generally too large, and it greatly changes the structure of the matrix.In the existing literature about the methods for solving ill-conditioning [9][10][11][12][13][14][15][16][17][18][19], there is no research on the combination of the H-K formula and the LM algorithm.Therefore, an improved LM algorithm based on the H-K formula (LM HK ) is proposed and used to solve the self-calibration model composed of the Brown model, QP model and Fourier model.
From the numerical experiment and analysis, it is found that, compared with the LM h algorithm, the LM HK algorithm can reach the same or higher accuracy, and it makes the SSR stable after fewer calculations (three iterations) in the simulation experiment, indicating that a high fitting accuracy can be reached by the improved algorithm with a higher iterative efficiency.In the ridge trace analysis, it is found that the damping factor determined by the improved LM algorithms is smaller, and it is consistent with the trend of the ridge trace curves, which means that this algorithm can weaken the ill-conditioning to a certain extent and make the solution stable by selecting a smaller damping factor and changing the structure of the Jacobian matrix as little as possible.
From the aspect of the distortion model, we consider the fitting effects of the Brown model, QP model and Fourier model on image distortion.In this paper, the study of the Brown model is different from the previous study [2,4].In the simulation experiment, the specific distortion form is uncertain.Under this condition, the fitting effect of the Brown model on distortion is studied.It is found that if there are unknown distortions in the image, the mathematical model, such as the QP model and Fourier model, can effectively compensate for these distortions.For the same algorithm, the fitting accuracy of the Fourier models for image point observations is higher, which is consistent with the previous study [7].
Although the study reveals the above discoveries, there are also limitations in this paper.In this paper, the LM algorithm is improved from the perspective of ridge estimation.However, in addition to the ridge trace method and H-K formula, the method for determining the ridge parameters also includes the minimum mean square error and cross validation, which are not considered.Therefore, an important future direction for improving the LM algorithm is a combination of the calculation of the damping factor and other methods for determining the ridge parameter.

Conclusions
Ill-conditioning generally exists in nonlinear least squares problems.An LM algorithm based on the gain ratio is the commonly used method, and it is discussed in this paper.This algorithm linearizes the nonlinear model and weakens the ill-conditioning of the Jacobian matrix by determining a damping factor.However, the damping factor determined by this algorithm is usually large, which greatly changes the structure of the matrix.Since the H-K formula can weaken the ill-conditioning by determining a smaller damping factor, an improved LM algorithm based on the H-K formula is proposed for solving nonlinear least squares problems.From the perspective of the ridge estimation of nonlinear parameters, the damping factor is calculated by the H-K formula.For the ill-conditioned problem in image distortion, the Brown model, quadratic polynomial model and Fourier model are discussed for the space resection of a single image.The improved LM algorithm based on the H-K formula is used to solve the self-calibration model, and it is applied to measure the diameter of coins using a single camera.The applicability of the improved LM algorithm in self-calibration is verified.Through numerical experiments, it is found that the improved LM algorithm can stabilize the parameter values under the condition of changing the matrix structure as little as possible, and it can reach the same or a higher accuracy than the LM algorithm based on the gain ratio.The improved LM algorithm based on the H-K formula provides a new method for solving the self-calibration model, and it also provides a new idea for solving ill-conditioned nonlinear least squares.

Figure 2 .
Figure 2. The distribution of reprojection errors.

Figure 2 .
Figure 2. The distribution of reprojection errors.

Figure 3 .
Figure 3. Iterative changes of SSR in different methods.

Figure 3 .
Figure 3. Iterative changes of SSR in different methods.

10 −
. According to the selection principle of the damping factor and Table

Figure 4 .
Figure 4. Ridge trace curves of x ˆ at the optimal solution.

Figure 4 .
Figure 4. Ridge trace curves of x at the optimal solution.

Figure 5 .
Figure 5.The detected points and the coins.

3. 2 .
Measurement of a Coin Diameter Using a Single Camera Experiment 3.1 shows that the LM algorithm based on CDM presents a better performance in solving the self-calibration model.Therefore, LM CDM+h and LM CDM+HK are used to measure the diameters of coins using a single camera in this experiment.Nine images of a calibration pattern are taken from different angles to calibrate the camera by Zhang's calibration method.Taking the calibration result as the initial value and using the detected point of the last image, LM CDM+h and LM CDM+HK are combined with the Brown model, QP model and Fourier model to estimate the model parameters and then measure the diameters of the coins.The detected points and the coins are shown in Figure 5.The SSR and the maximum residuals in pixels of each algorithm and model are shown in

Figure 5 .
Figure 5.The detected points and the coins.

Figure 5 .
Figure 5.The detected points and the coins.

Table 1 .
SSR and maximum residuals of each algorithm.

Table 4 .
The number of iterations k and the running time (s).

Table 4 .
The number of iterations k and the running time (s).

Table 5 .
The condition number and damping factor at the optimal solution.

Table 5 .
The condition number and damping factor at the optimal solution.

Table 5 .
The condition number and damping factor at the optimal solution.

Table 7 .
The number of iterations k and the running time (s).To measure the coins, the top-left and top-right corners of the bounding box are converted into world coordinates.Then the Euclidean distance between them is calculated in millimeters.The actual diameter is 19.05 mm.The diameters of the coins calculated by different methods are shown in Table8.It can be seen from Table8that the numerical results calculated by the different methods are consistent.The measurements of the first coin and the second coin are accurate to within 0.004 mm and 0.164 mm.All the methods have a sufficiently high accuracy.

Table 8 .
The diameter of coins calculated by different methods (mm).