A fast Bayesian operational modal analysis method based on Newton iteration

Abstract

The present invention proposes a fast Bayesian operational modal analysis method based on Newton iteration, which can be used for on-site vibration testing and long-term health monitoring of engineering structures, including: S1, collecting vibration response data of the structure under environmental excitation and performing FFT transformation on it; S2, inputting frequency domain data; S3, giving the initial value of the iteration of the modal parameter; S4, regarding the modal response as a hidden variable, and starting the iterative operation of the EM algorithm; S5, using the convergence result of the EM algorithm in S4 as the initial value, starting the Newton iteration to update the parameters, and obtaining the modal parameter optimization result and its uncertainty. Compared with the prior art, the present invention can achieve higher recognition accuracy with a faster convergence speed, give full play to the advantages of the Newton iteration algorithm, and converge within a few seconds even for complex three dense mode problems. At the same time, the posterior covariance matrix can be output as a by-product, the calculation process is unified, and it is convenient for practical application.

Description

Technical Field

[0001] The present invention belongs to the technical field of data analysis for structural health monitoring and relates to a fast Bayesian operational modal analysis method based on Newton iteration. Background Art

[0002] Operational modal identification or analysis is to extract the dynamic characteristic parameters of the structure, including natural frequency, damping ratio and vibration mode, from the vibration response data of the structure under environmental excitation. These results can provide important basis for structural health monitoring, vibration control and finite element update. Since the excitation of the structure is not measured in this process, coupled with the existence of random noise and model errors, the identification results often show a certain degree of uncertainty. The introduction of the Bayesian principle can transform the operational modal analysis problem into a posterior inference problem. One of the benefits of this is that the algorithm can strictly quantify its identification uncertainty while optimizing the modal parameters. Although Bayesian operational modal analysis algorithms for different working conditions such as single-step testing, multi-step testing and asynchronous data have been developed in recent years, with the continuous improvement of acquisition technology and testing requirements, new challenges have been posed to the robustness, convergence accuracy and computational efficiency of these algorithms.

[0003] There are two main problems with the existing algorithms. First, as the number of parameters increases, the convergence accuracy is low or even fails to converge in some complex situations such as dense modes. Second, the identification of modal parameters and the quantification of uncertainty are usually two independent processes, ignoring their inherent connection, resulting in low computational efficiency and difficulty in applying to long-term health monitoring of structures. Summary of the invention

[0004] In view of the limitations of existing algorithms, the present invention proposes a fast Bayesian operational modal analysis method based on Newton iteration, which can greatly improve the calculation speed without losing the identification accuracy. At the same time, it can output the uncertainty of modal parameter identification as a by-product, unify the calculation process, and improve the calculation efficiency.

[0005] The technical solution of the fast Bayesian operational modal analysis method based on Newton iteration of the present invention can be given by the following steps:

[0006] S1. Collect the vibration response data of the structure under environmental excitation, and convert it into frequency domain data through fast Fourier transform;

[0007] S2. Obtain the sample power spectrum density matrix of the data in S1 at different frequencies, obtain the singular value spectrum after eigenvalue decomposition, select a frequency band near the resonance peak, and use the data within the frequency band as input;

[0008] S3. Obtain all peaks in the frequency band selected by S2 and set the iterative initial values ​​of the modal parameters; the modal parameters are the natural frequency, damping ratio, vibration mode, modal excitation power spectral density and prediction error power spectral density;

[0009] S4. Treat the modal response as a latent variable and perform iterative operations of the EM algorithm, repeating the iteration until convergence; in the EM algorithm, the E step is to calculate the expectation of the "complete data" negative log-likelihood function CDNLLF about the distribution of the latent variable, that is, to calculate the Q function, which is an upper bound of the negative log-likelihood function NLLF at the current step; the M step is to optimize the Q function and update the modal parameters;

[0010] S5. Using the convergence result of the EM algorithm in S4 as the initial value, start the Newton iteration to update the modal parameters; when finally converged, output the modal parameter optimization result and its uncertainty.

[0011] Compared with the prior art, the present invention has the following advantages:

[0012] 1. Fast convergence speed. It uses the quadratic convergence property of the Newton iteration algorithm, which is about an order of magnitude faster than the existing algorithms. Even for complex dense modal problems, it can converge within a few seconds or even a few tenths of a second.

[0013] 2. It has high convergence accuracy and can combine the advantages of the EM algorithm and Newton iteration. It has strong robustness and can give better optimization results than the current algorithm for highly dense modes.

[0014] 3. The calculation process is unified, and each iteration step can give the result of the inverse operation of the Hessian matrix. The last iteration can output it as the posterior covariance matrix without the need to additionally solve the identification uncertainty. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 The figure is a schematic diagram of the process of the method provided by the present invention.

[0016] Figure 2 Schematic diagram of a model of an embodiment of the present invention.

[0017] Figure 3 The power spectral density and singular value spectrogram of the simulated data for the embodiment.

[0018] Figure 4 This is the modal recognition result of the method proposed in the present invention in the embodiment. DETAILED DESCRIPTION

[0019] The present invention is further described and illustrated below in conjunction with specific embodiments. The embodiments are merely exemplary of the present disclosure and do not define the scope of limitation. The technical features of each embodiment of the present invention may be combined accordingly without conflicting with each other.

[0020] like Figure 1 As shown, the fast Bayesian operation mode identification method based on Newton iteration of the present invention comprises the following steps:

[0021] S1. Collect the vibration response data of the structure under environmental excitation, denoted as Here n represents the total number of degrees of freedom measured, and N represents the number of time domain sample points collected by each data channel. The vibration response data is time domain data, and it is regarded as a stationary random process, so its fast Fourier transform (FFT) can be defined as:

[0022]

[0023] Where i is the imaginary unit, Δt is the sampling interval, and k is the index of the frequency domain point.

[0024] According to the modal superposition method, the theoretical response of the structure can be considered as the sum of the contributions of each mode. Adding modeling errors and data noise, the FFT data can be modeled as:

[0025]

[0026] In the formula is the vibration mode matrix, m is the modal order in the frequency band, is the modal excitation, is the prediction error, describing the discrepancy between the theoretical response of the structure and the measured data. is a diagonal matrix whose diagonal elements contain the frequency response functions of various orders, given by:

[0027]

[0028] In the formula, f i and i are the i-th order natural frequency and damping ratio respectively. When q is 0, 1 and 2, it corresponds to acceleration, velocity and displacement data respectively.

[0029] S2, obtain the sample power spectrum density matrix of the data in S1 at different frequencies, obtain the singular value spectrum after eigendecomposition, select a frequency band near the resonance peak, and use the FFT data within the frequency band as input. Assuming that the modal excitation and prediction error are independent stationary random processes, it can be obtained that they obey a circularly symmetric complex Gaussian distribution with zero mean, and their covariance matrices can be expressed as and As a result, the FFT data It also obeys a complex Gaussian distribution with zero mean, and its covariance matrix is ​​given by:

[0030]

[0031] Where 'E[·]' represents the expected operation, is the modal response η k =h k p k The power spectral density (PSD), I n is the identity matrix, symbol '[·] T '、'[·] H 'and'[·] * ' represent transpose, conjugate and conjugate transpose operations respectively.

[0032] Without loss of generality, the unknown parameters to be analyzed in the present invention can be expressed as a real number vector θ = [f; ζ; vec(Φ); S ii ; a; b; S e ], where 'vec(·)' means arranging the matrix vertically into column vectors, represents the diagonal elements of S, and are the independent real and imaginary parts of the off-diagonal elements of S. Thus, the parameter quantity within a frequency band is

[0033] Assuming that the acquired data is long enough, that is, the number of FFT points in the selected frequency band is much greater than 1, it can be proved that They are independent of each other at different frequencies, so the negative log-likelihood function (NLLF) can be expressed as L(θ) = ∑ k L k (θ), where:

[0034]

[0035] Wherein, 'ln[·]' represents the natural logarithm operation, and '|·|' represents the determinant operation.

[0036] Assume that the prior distribution of θ is uniform. According to the Laplace approximation, its posterior probability density function can be approximated by a Gaussian distribution, and the mean and covariance matrix are given by the following formula:

[0037]

[0038] In the formula represents the parameter value that minimizes NLLF, represents the Hessian matrix of NLLF at the optimal point.

[0039] At this time, the Bayesian model does not have a global solution, because the vibration mode is essentially a eigenvector and can be normalized arbitrarily. In this invention, we stipulate that the modulus of each vibration mode is 1 and introduce the constraint function g i (θ), which can be expressed as:

[0040]

[0041] This makes the model globally solvable in a physically correct sense.

[0042] Therefore, Bayesian operational modal analysis includes two aspects: one is to minimize NLLF to find the optimal parameter value; the other is to obtain the second-order derivative of NLLF to synthesize the Hessian matrix and calculate its inverse. In theory, both can be solved by directly processing formula (5), but this formula is very complicated, and all modal parameters are coupled to E k In, E k It exists in the form of an inverse and a determinant, so this is not a good choice.

[0043] Since the excitation is not measured in the operational modal analysis, the physical quantities related to it conform to the concept of hidden variables in statistics. In practice, the introduction of appropriate hidden variables can often simplify the calculation. Based on this, the convenient and efficient expectation maximization (EM) algorithm can be used to solve the problem, which can make the calculation more efficient. Before this, the initial value of the algorithm iteration must be given, as shown in the following steps.

[0044] S3, pick up all the peaks in the frequency band selected by S2 (one corresponds to sparse mode, and multiple corresponds to dense mode), and use the corresponding frequency points as the iterative initial value of the natural frequency; for civil engineering structures, the damping ratio is usually taken as 1% according to experience; the vibration mode takes the eigenvector corresponding to the maximum eigenvalue of the sample PSD matrix at the initial value of the natural frequency point as the initial value. The initial value of the modal excitation PSD is given by the following formula:

[0045]

[0046] Where N f The number of FFT data points for the selected frequency band.

[0047] The initial value of the prediction error PSD is given by:

[0048]

[0049] In the formula, The columns contain The eigenvectors corresponding to the first m largest eigenvalues ​​of .

[0050] It should be noted that the PSD matrix and covariance matrix mentioned in the present invention are essentially the same in meaning, but the former is a physical concept and the latter is a mathematical concept.

[0051] S4, the modal response As a hidden variable, the EM algorithm iterative operation begins. The E step is to calculate the expectation of the "complete data" negative log-likelihood function (CDNLLF) about the distribution of the hidden variable, hereinafter referred to as the Q function, which is an upper bound of the NLLF at the current step; the M step is to optimize the Q function and update the parameters.

[0052] Specifically, based on the long data assumption, CDNLLF can be expressed as in:

[0053]

[0054] Where c is a constant term that is independent of θ.

[0055] According to the model assumptions, we can get the given θ and Under the condition of k It obeys a Gaussian distribution, and its mean and covariance matrix are given by:

[0056]

[0057] In the formula

[0058] The expected Q function of formula (10) with respect to the distribution described in formula (11) can be expressed as Q(θ)=∑ k Q k (θ), where:

[0059]

[0060] Where 'tr[·]' represents the trace operation of the matrix, 'Re(·)' represents taking the real part, represents the latent variable η k The second-order origin moment of .

[0061] Obtain the first-order derivative of formula (12) with respect to the modal parameters and organize them into the following table:

[0062] Table 3: First-order derivative of the Q function

[0063]

[0064] Note: f 2 =f⊙f,f 3 =f⊙f⊙fandζ 2 =ζ⊙ζ.

[0065] In the table, we define Where 'x' represents a parameter in θ'. In order to temporarily avoid considering S as a complex symmetric matrix constraint, a new parameter vector θ' = [f; ζ; vec(Φ); S e; vec(S)].

[0066] Let the above derivatives be 0 to obtain the analytical update form of some parameters, which are given by the following formulas:

[0067]

[0068] Because the natural frequency and damping ratio are coupled in the frequency response function in a nonlinear manner, their analytical update form cannot be obtained. The present invention uses the Newton iteration algorithm to update these two parameters, and the iteration equation is given by the following formula:

[0069]

[0070] In the formula, v = [f; ζ], and They represent the Hessian matrix and gradient vector of the Q function with respect to v. The second-order derivative of the Q function will be given in the next step.

[0071] In summary, the EM algorithm can be expressed as follows: E step, the first-order and second-order origin moments of the latent variables are calculated based on the modal parameters of the current step based on formula (11); M step, the modal parameters of the next step are updated according to formulas (13) to (16). These two steps are repeated until convergence. In this stage, the convergence criterion can be defined as S e The relative change between the two steps is less than a set threshold, such as 10 -4 .

[0072] The EM algorithm is simple, intuitive, and easy to apply, especially because it can give an analytical update form of a vibration mode with a huge amount of parameters (proportional to n), so its single-step iteration speed is extremely fast. But on the other hand, because it is essentially a gradient-based optimization algorithm, it is easily affected by the local properties of the objective function such as saddle points, resulting in the effective number of movement steps of the parameters being only reflected in the first few steps. After that, even if it iterates a hundred steps, the objective function value hardly changes. Taking this into account, the present invention introduces the Newton iteration algorithm, fully combines the advantages of the two, and achieves fast and accurate solutions, which are specifically given by the following steps.

[0073] S5. Using the convergence result of the EM algorithm in S4 as the initial value, start Newton iteration to update the parameters. The iteration equation is given by the following formula:

[0074]

[0075] Where λ is the step size factor, which is used to control the amplitude of the parameter movement in each iterative step, and the column of U contains the constraint function A set of bases of the null space of the gradient is introduced to consider the constraints described in formula (7), so that the matrix in the brackets is a full-rank matrix, which is convenient for direct inversion. The sum of E1, E2 and E3 gives the Hessian matrix of NLLF, represents the gradient of the Q function.

[0076] As mentioned above, the NLLF given by formula (5) is very complex in form. It would be very cumbersome to directly use it to obtain the first-order and second-order derivatives of the modal parameters to obtain the gradient and Hessian matrix. In contrast, the modal parameters in formula (10) are partially decoupled and can be divided into two groups: S e and Φ are coupled in the third term, and f, ζ, and S are coupled in H k In the equation, the two sets of parameters are no longer coupled. This property is conducive to the determination of derivatives, especially for the second-order derivatives, the cross terms generated by the two sets of parameters are zero, which makes the Hessian matrix more sparse and convenient for storage and calculation. To this end, the present invention introduces the Fisher equation and the Louis equation to establish the relationship between the first-order derivative and the second-order derivative of NLLF and CDNLLF, which are given by the following equations respectively:

[0077]

[0078] Based on this, the problem is transformed into finding the first and second derivatives of CDNLLF. For the former, since the initial value of the Newton iteration is the convergence result of the EM algorithm, and in the M step, the first-order derivative of the Q function has been given, see Table 1, so this part can directly use the previous results; for the latter, the Louis equation actually divides the Hessian matrix of NLLF into three parts, E1 is the product of the gradient of the Q function. Similarly, this part does not need to be calculated additionally; E2 is the Hessian matrix of the Q function. Compared with NLLF, this part is more sparse, and only one more derivative needs to be calculated based on the first-order derivative. The results are shown in the following table:

[0079] Table 4: Second-order derivative of the Q function

[0080]

[0081] Note: In the formula, we define where 'x' and 'y' both represent a parameter in θ'.

[0082] E3 is the expectation of the CDNLLF gradient product. For the gradient part, since the Q function is calculated based on the CDNLLF to calculate the expectation of the hidden variable, and the derivative and expectation are both linear operators and can be exchanged, the result of the Q function gradient can be used to replace the origin moment in the expression with the corresponding hidden variable. Specifically, Replace Replace with η k As for the expectation, given the current step parameter θ (t) and data Under the condition of kThe distribution of is a Gaussian distribution, and its moments have been mathematically concluded and can be expressed by the mean and covariance described in formula (11). The first to fourth order moments of the latent variables used in the present invention can be calculated by the following formulas:

[0083]

[0084] Where a is a constant vector, A, B and C are constant matrices,

[0085] Through the calculation of the above three parts, we can get the gradient vector of NLLF with respect to θ′ and the Hessian matrix However, this is not the ultimate goal, because the parameter vector we define is θ, which differs from θ′ only in the expression of S. By introducing some constant matrices, S can be decomposed into:

[0086]

[0087] In the formula, the constant matrix D h and D s Respectively with S ii , a and b are multiplied to form the diagonal elements of S, the real and imaginary parts of the off-diagonal elements, and they are converted into vector form. Thus, the relationship between θ and θ′ is given by the following formula:

[0088] θ′=Dθ (50)In the formula, 'blkdiag(·)' represents a block diagonal matrix.

[0089] Therefore, the gradient vector and Hessian matrix of NLLF with respect to θ can be expressed as:

[0090]

[0091] In the iterative process, in order to ensure that the parameters updated in each step can meet the constraints (such as the range of the natural frequency is within the selected frequency band, the damping ratio is within the range of 0 to 1, etc.) and can reduce the value of NLLF, the step size factor λ is introduced. Its value can be initially set to 1, and the iterative equation (17) is executed. If the obtained parameters do not meet the requirements, they are halved until they meet the requirements. The convergence criterion for this stage can be defined as the maximum value of the relative change of the parameters in the two-step iteration is less than the set threshold, such as 10 -3 It is worth mentioning that since the inverse of the NLLF Hessian matrix is ​​calculated at each iteration, when it finally converges, the result can be output as the posterior covariance matrix to characterize the uncertainty of parameter analysis without additional calculation.

[0092] The specific effects of the method of the present invention are demonstrated below in conjunction with the embodiments and drawings.

[0093] Example:

[0094] The structure of this embodiment is an imaginary shear frame structure. Figure 2 As shown, the two side lengths are 0.4m and 0.45m respectively, which leads to very similar stiffness and mass distribution of the structure in the x and y directions, and it can be expected that obvious dense modes will be generated. Assuming that the measurement degrees of freedom are 8, distributed in the four corners, there are three modes contributing to the structural response, namely the first order of translation in the x and y directions and the first order of torsion. The natural frequencies are set to 4.20Hz, 4.25Hz and 4.40Hz; the damping ratios are 1.00%, 1.50% and 2.00%; the diagonal element of the modal excitation PSD is 1 (μg) 2 / Hz, non-diagonal elements S 21 =0.5e iπ / 4 , conjugated at symmetric positions, and the rest are 0; the prediction error PSD is 10 (μg) 2 / Hz. Assume that the measurement is performed at a sampling frequency of 100 Hz for 20 minutes.

[0095] The PSD and singular value spectrum (SV) of the simulated data are shown in Figure 3 As shown, the SV bracket is the selected frequency band, and the circle represents the iterative initial value of the third-order modal natural frequency. Figure 4 The recognition results of the method proposed in the present invention are well in line with the theoretical preset values, and the corresponding uncertainties are in brackets. The entire recognition process takes about 5 seconds, which is about an order of magnitude faster than the existing algorithm. And its NLLF value is also smaller than that of the existing algorithm, reflecting its higher recognition accuracy.

[0096] The above-mentioned embodiments only express several implementation methods of the present invention, and the description is relatively specific and detailed, but it cannot be understood as limiting the scope of the present invention. For ordinary technicians in this field, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. A fast Bayesian operational modal analysis method based on Newton iteration, characterized in that The steps include: S1. Collect the vibration response data of the structure under environmental excitation, and convert it into frequency domain data through fast Fourier transform; The step S1 is: The vibration response data is recorded as n represents the total number of degrees of freedom measured, and N represents the number of time domain sample points collected by each data channel; The Fast Fourier Transform FFT is defined as: In the formula, i is the imaginary unit, Δt is the sampling interval, and k is the index of the frequency domain point; The FFT data is modeled as: In the formula is the vibration mode matrix, m is the modal order in the frequency band, is the modal excitation, is the prediction error, is a diagonal matrix whose diagonal elements contain the frequency response functions of various orders, given by: In the formula, f i and i are the i-th order natural frequency and damping ratio respectively; when q is 0, 1 and 2, they correspond to acceleration, velocity and displacement data respectively; S2. Obtain the sample power spectrum density matrix of the data in S1 at different frequencies, obtain the singular value spectrum after eigenvalue decomposition, select a frequency band near the resonance peak, and use the data within the frequency band as input; Assuming that the modal excitation and prediction error are independent stationary random processes, they obey a circularly symmetric complex Gaussian distribution with zero mean, and their covariance matrices are expressed as and Therefore, the FFT data It also obeys a complex Gaussian distribution with zero mean, and its covariance matrix is ​​given by: Where E[·] represents the expected operation, is the modal response η k =h k p k The power spectral density PSD, I n is the identity matrix, symbol [·] T ,[·] H and[·] * denote transpose, conjugate, and conjugate transpose operations respectively; The modal parameters can be expressed as a real vector θ = [f; ζ; vec(Φ); S ii ; a; b; S e ], where vec(·) means arranging the matrix vertically into column vectors, represents the diagonal elements of S, and are the independent real and imaginary parts of the off-diagonal elements of S; Independent of each other at different frequencies, the second negative log-likelihood function NLLF is expressed as L(θ)=∑ k L k (θ), where: In the formula, ln[·] represents the natural logarithm operation, and |·| represents the determinant operation; Assume that the prior distribution of θ is uniform. According to the Laplace approximation, its posterior probability density function can be approximated by a Gaussian distribution, and the mean and covariance matrix are given by the following formula: In the formula represents the parameter value that minimizes NLLF, represents the Hessian matrix of NLLF at the optimal point; S3. Obtain all peaks in the frequency band selected by S2 and set the iterative initial values ​​of the modal parameters; the modal parameters are the natural frequency, damping ratio, vibration mode, modal excitation power spectral density and prediction error power spectral density; S4. Treat the modal response as a latent variable, perform iterative operations of the EM algorithm, and repeat the iteration until convergence; in the EM algorithm, the E step is to calculate the expectation of the first negative log-likelihood function CDNLLF about the distribution of the latent variable, that is, to calculate the Q function, which is an upper bound of the second negative log-likelihood function NLLF at the current step; the M step is to optimize the Q function and update the modal parameters; The S4 is specifically: Based on the long data assumption, CDNLLF is expressed as in: Where c represents a constant term that is independent of θ; According to the model assumptions, we can get the given θ and Under the condition of k It obeys a Gaussian distribution, and its mean and covariance matrix are given by: In the formula The expected value of formula (9) with respect to the distribution described in formula (10) is obtained as Q function, which is expressed as Q(θ)=∑ k Q k (θ), where: Where tr[·] represents the trace operation of the matrix, Re(·) represents taking the real part, represents the latent variable η k The second-order origin moment of The E step is to calculate the first-order and second-order origin moments of the latent variables based on the modal parameters of the current step based on formula (10); The M step is to update the modal parameters of the next step according to formulas (12) to (15); Among them, the natural frequency and damping ratio are coupled in the frequency response function in a nonlinear way, and their analytical update forms cannot be obtained. Therefore, the Newton iteration algorithm is used to update these two parameters; S5. Using the convergence result of the EM algorithm in S4 as the initial value, start the Newton iteration to update the modal parameters; when finally converged, output the modal parameter optimization result and its uncertainty.

2. The fast Bayesian operational modal analysis method based on Newton iteration according to claim 1 is characterized in that: The S3 is specifically: The frequency points corresponding to all peaks in the selected frequency band are used as the iterative initial values ​​of the natural frequency; the initial value of the damping ratio is set according to the structural properties, and the vibration mode takes the eigenvector corresponding to the maximum eigenvalue of the sample PSD matrix at the initial value of the natural frequency point as the initial value; The initial value of the modal excitation power spectral density is given by the following formula: Where N f The number of FFT data points for the selected frequency band; The initial value of the prediction error power spectral density is given by: In the formula, The columns contain The eigenvectors corresponding to the first m largest eigenvalues ​​of .

3. The fast Bayesian operational modal analysis method based on Newton iteration according to claim 1 is characterized in that: The iterative equation for updating the natural frequency and damping ratio using the Newton iteration algorithm is: In the formula, v = [f; ζ], and They represent the Hessian matrix and gradient vector of the Q function with respect to v respectively.

4. The fast Bayesian operational modal analysis method based on Newton iteration according to claim 1 is characterized in that: In S5, the iterative equation of Newton iteration is given by the following formula: Where λ is the step size factor, which is used to control the amplitude of the parameter movement in each iterative step, and the column of U contains the constraint function A basis for the null space of the gradient, the sum of E1, E2 and E3 gives the Hessian matrix of the NLLF, E1 is the product of the gradient of the Q function, E2 is the Hessian matrix of the Q function, E3 is the expectation of the product of the CDNLLF gradients, represents the gradient of the Q function.

5. The fast Bayesian operational modal analysis method based on Newton iteration according to claim 4 is characterized in that: Introducing the constant matrix D h and D s , S is decomposed into: In the formula, the constant matrix D h and D s Respectively with S ii , a and b are multiplied to form the diagonal elements of S, the real and imaginary parts of the off-diagonal elements, and they are converted into vector form; then, the relationship between θ and θ′ is given by the following formula: θ′=Dθ (18) In the formula, θ′=[f; ζ; vec(Φ); S e ; vec(S)] is the newly defined parameter vector, blkdiag(·) represents a block diagonal matrix; Therefore, the gradient vector and Hessian matrix of NLLF with respect to θ are expressed as: In the Newton iteration process of S5, in order to ensure that the updated parameters in each step can meet the constraints and reduce the value of NLLF, the step size factor λ is introduced. Its value can be initially set to 1. The iterative equation (17) is executed. If the obtained parameters do not meet the requirements, they are halved until they meet the requirements. Since the inverse of the NLLF Hessian matrix is ​​calculated at each iteration, when it finally converges, the result is output as the posterior covariance matrix to characterize the uncertainty of parameter analysis without additional calculation.

Images