Measuring Point Planning and Fitting Optimization of the Flange and Spigot Structures of Aeroengine Rotors

: An optimized measuring point planning and ﬁtting method for rotor ﬂange and spigot structures was proposed to achieve precise measurement of position and pose of the aeroengine rotors during docking processes. Firstly, the impact of circumferential phase angle, distribution range angle, total number of measuring points, and number of distribution rings on measurement uncertainty was analyzed. The measuring point planning schemes for ﬂange and spigot were proposed. Secondly, the Gauss Newton iterative solution principle considering damping factors was clariﬁed. Subsequently, an optimized iterative reweighting method consisting of weight iterative estimation, singular value detection under the Chauvenet criterion, and clustering detection was proposed for ﬁtting the ﬂange annular end face. A mapping point total least squares method with practical geometric signiﬁcance was proposed for ﬁtting the spigot cylinder face. Finally, measuring and ﬁtting experiments were performed. The singular measuring point detection methods were veriﬁed. Under the optimized ﬁtting methods, the goodness of ﬁt and average orthogonal distance of ﬂange and spigot structures are 0.756 and 0.089 mm, respectively, which have higher ﬁtting accuracy than the other traditional methods.


Introduction
Flange and spigot are the typical positioning structures at the end of an aeroengine rotor, realizing axis orientation, axial positioning, and radial positioning between adjacent rotors [1][2][3]. A large amount of measurement work is required during the rotor docking assembly processes, to obtain the initial position and pose of each rotor and to provide data for docking trajectories [4][5][6]. Therefore, the measuring and fitting accuracy is crucial, directly determines the assembly accuracy of the rotors, and has a significant impact on the dynamic performance of the entire engine [7][8][9][10].
Automated systems have gradually been applied in the measurement of aeroengine rotors [11][12][13][14][15]. Among these, the multi-axis precision mechanical contact measurement system and laser measurement system have the advantages of high reliability and automation [16][17][18]. A large number of point cloud coordinates are generated after measurement. The question of how to arrange the number and position of measuring points is a key issue in balancing measurement accuracy and efficiency [19][20][21][22]. Based on the measuring point planning, optimizing the corresponding fitting methods for different structures is an effective way to further reduce the uncertainty of model parameters [23][24][25]. In summary, reasonable measuring point planning and optimized fitting methods of the flange and spigot structures are important prerequisites for obtaining the rotor position and pose information quickly and accurately.
In recent years, some scholars have conducted relevant research on the measuring point planning strategy of aeroengine rotors [26,27]. Zhang et al. have proposed a redundant measuring point detection method for the measurement of rotor end faces, and the optimized algorithm improved the accuracy by 21 µm compared with the traditional least squares (LS) fitting method [28]. Zhu et al. constructed a heuristic floating forward search algorithm to efficiently find a near-optimal solution for measuring point layout [29]. Ding et al. used the eigenvalue method to process the point cloud data of the rotor end face, which has the function of excluding gross error points in the measurement points [30]. However, the existing measuring point planning methods mainly deal with integral surfaces with simple geometric shapes, and their applicability to complex spatial feature structures such as flange and spigot is limited.
With regard to the research regarding the fitting methods for typical axes and circular contours: Xie et al. established a spatial axis LS fitting model with non-cylindrical data filtered out [31]. The Limacon model is widely used in circular profile measurement to characterize eccentricity parameters [32,33]. Chetwynd et al. improved the parameter estimation method of LS fitting radius in the Limacon model [34]. Sun et al. proposed the vector projection method in cylindrical contour fitting. Compared with the traditional Limacon model, the accuracy of rotor roundness has been improved by 2.2 µm, and the accuracy of coaxiality has been improved by 9.39 µm [35,36]. Zhang et al. constructed a particle swarm optimization algorithm for fitting aeroengine rotor models based on a two-dimension laser displacement sensor, and the measurement accuracy of cylindricity error was improved by 1.768 µm compared with the traditional method [37]. Wang et al. used the equivalent homogenization processing and 3D LS method for fitting aeroengine rotor annular plane [38]. For different geometric features, researchers also carried out comparative studies on various fitting methods based on the same set of point cloud data. Song et al. compared the fitting performances of six normal line fitting methods on different types of geometries, e.g., sphere, cylinder, and ellipsoid [39]. Nouira et al. compared the fitting accuracy of different models such as LS circle, minimum area circle, maximum inscribed circle and minimum circumscribed circle, and proposed a fitting method based on a micro displacement model [40]. The above research constructed the corresponding fitting algorithms for the points with specific distribution forms, achieving higher fitting accuracy compared with traditional methods. However, few targeted combinatorial optimization algorithms have been established to meet the requirements for fitting the spatial position and pose of the rotor flange and spigot structures. Moreover, it is no longer possible to measure the flange and spigot again after the rotors are docked, so higher fitting accuracy is required in a single fitting process.
In this paper, a systematic study on measuring point planning and fitting methods for typical rotor flange and spigot structures was carried out. First, the impact of different measuring point distribution parameters on the model uncertainty was analyzed and a measuring point planning strategy was proposed. Then, the optimized iterative reweighting method for fitting flange and the mapping point total least squares method for fitting spigot were proposed, respectively. Finally, the measurement and fitting experiments were performed. The singular value point detection methods were compared and the fitting accuracy of the optimized fitting methods was verified.
This paper is organized as follows. In Section 2, the measuring point planning method is proposed; in Section 3, the optimized fitting algorithm is proposed; in Section 4, the measurement experiments are performed; and in Section 5, the conclusions are presented.

Measuring Point Planning Method
The schematic that illustrates the measurement of the flange and spigot structures of the aeroengine rotors is shown in Figure 1. With regard to the measuring point planning of the flange annular end face. The end face description equation can be represented by z M = a 1 + a 2 x M + a 3 y M in the measurement coordinate system O M -xyz. Without loss of generality, allow N measuring points M i (x Mi , y Mi , z Mi ) distributed on the plane z M = 0 to be fitted. To simplify the calculation, the case of N = 3 is analyzed first. The measuring points M 1 (0, .
Machines 2023, 11, x FOR PEER REVIEW 4 of 21  It can be calculated that when x M3 has the minimum value, both δ B and δ C have the minimum values simultaneously. The uncertainty of attitude column matrix a representing the fitted normal vector of the rotor end face is shown in Equations (4)- (6).
The increase in x M2 and y M3 leads to a decrease in uncertainty δ ax and δ ay , i.e., the relative distance of the measuring points M 1 , M 2 , and M 3 should be as far as possible. Then δ ax = δ ay and x M2 = ± √ 3R C and the angle among O M -M i is 2π/3. When the calculation method is generalized to N ≥ 3, it can be concluded that the uncertainty of the end face parameters is minimized when the measuring points are uniformly distributed along the circumference.  For situations where the measurement trajectory is limited, measuring a large arc range can suffice instead of measuring the entire circle. The analysis of the impacts of the number of measuring points and the uniform distribution range on the uncertainty of the end face model parameters are necessary. The measuring points are arranged symmetrically as shown in Figure 2b. Parameter a1 determines the third-row element az in the workpiece pose column matrix a, which is an important parameter for determining the spatial angle of the end face axis (the angle between the end face normal vector and the coordinate axis OM-z). The uncertainty δa1 of parameter a1 using the LS method is calculated as Equation (7). For situations where the measurement trajectory is limited, measuring a large arc range can suffice instead of measuring the entire circle. The analysis of the impacts of the number of measuring points and the uniform distribution range on the uncertainty of the end face model parameters are necessary. The measuring points are arranged symmetrically as shown in Figure 2b. Parameter a 1 determines the third-row element a z in the workpiece pose column matrix a, which is an important parameter for determining the spatial angle of the end face axis (the angle between the end face normal vector and the coordinate axis O M -z). The uncertainty δ a1 of parameter a 1 using the LS method is calculated as Equation (7).
where N is the number of measuring points, R C is the radius of the circular path, θ i is the circumferential angle of the measuring point M i , and σ M represents measurement systematic error. Taking the flange end face of an aeroengine low pressure turbine (LPT) mid shaft as an example, the measuring path radius R C = 149 mm. The relationship between the calculated uncertainty δ a1 and the number of measuring points N/the angle of uniform distribution range ϕ is shown in Figure 3a. δ a1 reduces significantly when ϕ ≥ 180 • . Therefore, the measurement range should cover at least half the arc area. The reduction rate of the minimum value of the model parameter uncertainty min(δ a1 /σ M ) gradually decreases when N increases step by step in the arithmetic progression of 4, 12, 20, 28, 36. The reduction rate of min(δ a1 /σ M ) is less than 10% when N increases from 28 to 36. Taking into account the impact on measurement accuracy and efficiency, N = 28 is a reasonable value. Due to the radial width on the circular end face, the influence of the number of measuring rings distributed along multiple concentric circles/arcs is analyzed while N = 28. The uncertainty is shown in Figure 3b when R C = 149, 137 and 125 mm, respectively. When 30 • ≤ ϕ ≤ 80 • , δ a1 under a multi-ring distribution reduces significantly compared with the single ring scheme and δ a1 under a double-ring distribution is lower than under a three-ring distribution. Due to the limitation of the flange edge width (24 mm) in this example, when the total number of measuring points remains unchanged, the excessive rings result in a small radial spacing and a large circumferential spacing of measuring points. Therefore, the arrangement of measuring points along the double concentric rings is the optimal method and has min(δ a1 /σ M ) = 0.188. In summary, the N = 28 measuring points evenly distributed in double rings form a reasonable measuring point planning scheme, as shown in Figure 2c. Similarly, the optimized measuring point planning strategies for end faces at different positions (shown in Figure 4) in an aeroengine LPT rotor assembly are analyzed, as listed in Table 1.
Taking the flange end face of an aeroengine low pressure turbine (LPT) mid shaft as an example, the measuring path radius RC = 149 mm. The relationship between the calculated uncertainty δa1 and the number of measuring points N/the angle of uniform distribution range φ is shown in Figure 3a. δa1 reduces significantly when φ ≥ 180°. Therefore, the measurement range should cover at least half the arc area. The reduction rate of the minimum value of the model parameter uncertainty min(δa1/σM) gradually decreases when N increases step by step in the arithmetic progression of 4, 12, 20, 28, 36. The reduction rate of min(δa1/σM) is less than 10% when N increases from 28 to 36. Taking into account the impact on measurement accuracy and efficiency, N = 28 is a reasonable value. Due to the radial width on the circular end face, the influence of the number of measuring rings distributed along multiple concentric circles/arcs is analyzed while N = 28. The uncertainty is shown in Figure 3b when RC = 149, 137 and 125 mm, respectively. When 30° ≤ φ ≤ 80°, δa1 under a multi-ring distribution reduces significantly compared with the single ring scheme and δa1 under a double-ring distribution is lower than under a three-ring distribution. Due to the limitation of the flange edge width (24 mm) in this example, when the total number of measuring points remains unchanged, the excessive rings result in a small radial spacing and a large circumferential spacing of measuring points. Therefore, the arrangement of measuring points along the double concentric rings is the optimal method and has min(δa1/σM) = 0.188. In summary, the N = 28 measuring points evenly distributed in double rings form a reasonable measuring point planning scheme, as shown in Figure 2c. Similarly, the optimized measuring point planning strategies for end faces at different positions (shown in Figure 4) in an aeroengine LPT rotor assembly are analyzed, as listed in Table 1.   (c) (d)     With regard to the measuring point planning of the spigot cylindrical face. Without loss of generality, let N points M proj_i (x proj_Mi , y proj_Mi , z proj_Mi ) projected by the measuring points be distributed on the assumed plane z M = 0. To simplify the calculation, the case of N = 3 is analyzed first. The mapping points M proj_1 (0, R proj_M , 0), M proj_2 (R proj_M cosθ 2 , R proj_M sinθ 2 , 0), and M proj_3 (R proj_M cosθ 3 , R proj_M sinθ 3 , 0) are shown in Figure 2d. The coordinate uncertainty of x proj_MO and y proj_MO in the model parameters of the standard circle is shown in Equations (8) and (9). Thus, the uncertainty δ O of the center coordinates of the mapped circle is calculated as: When ∂δ O /∂θ 2 = 0, the uncertainty δ O has a minimum value if θ 2 = θ 3 /2. The relationship between δ O and circumferential phase angle θ 3 is shown in Figure 3c. When θ 2 = 2π/3, θ 3 = 4π/3, δ O has the minimum value 1.155σ M . When the calculation method is generalized to N ≥ 3, it can be concluded that the uncertainty of the center coordinate has a minimum value when the measuring points are uniformly distributed along the circumference.
The number of points and the uniformly distributed range angle also have an impact on the uncertainty δ O . The mapping points are arranged symmetrically as shown in Figure 2e. δ O is calculated as Equation (11) when using the LS method.
The relationship between the uncertainty δ O and the number of measuring points N/the angle of uniform distribution range ϕ is shown in Figure 3d. δ O reduces significantly when ϕ ≥ 180 • . Therefore, the measurement range should cover at least half the arc area. The reduction rate of the minimum value of uncertainty min(δ O /σ M ) gradually decreases when N increases step by step in the arithmetic progression of 4, 12, 20, 28, 36, 44. The reduction rate of min(δ O /σ M ) is less than 10% when N increases from 36 to 44. Taking into account the impact on measurement accuracy and efficiency, N = 36 is a reasonable value. Similarly, the optimized measuring point planning strategies for cylindrical faces at different positions (shown in Figure 4) in an aeroengine LPT rotor assembly are analyzed, as listed in Table 2.
The solution process of the Gauss Newton method for searching the minimum value χ 2 is shown in Figure 5. The gradient of ∆a = 0 when F(a + ∆a) is minimized. The model parameter adjustment matrix ∆a = −H −1 g. The gradient vector g and Hessian matrix H are decomposed into the Jacobian matrix J LS , as shown in Equation (13).

LS principle and Gauss Newton Iterative Solution Method Considering Damping Factors
Measuring point fitting can be considered a minimization process. ( ) The solution process of the Gauss Newton method for searching the minimum value χ 2 is shown in Figure 5. The gradient of Δa = 0 when F(a + Δa) is minimized. The model parameter adjustment matrix Δa = −H −1 g. The gradient vector g and Hessian matrix H are decomposed into the Jacobian matrix JLS, as shown in Equation (13).
Model parameter iteration adjustment matrix Δa is shown in Equation (14). The iteration termination target is shown in Equation (16). Model parameter iteration adjustment matrix ∆a is shown in Equation (14). The iteration termination target is shown in Equation (16).
where W is the weight diagonal matrix; r is the residual column vector; and ε is the iteration termination value. The traditional Gauss Newton algorithm may fall into saddle points or iterate in wrong directions. The damping factor µ is introduced in the Levenberg-Marquardt method. Thus, the iterative adjustment matrix ∆a is revised to: The damping factor µ is determined by the adaptive coefficient λ µ and the maximum diagonal element n max of the quadratic normal matrix N = J T LS WJ LS , as shown in Equation (18). The empirical value of λ µ can be set to 0.001 in the first iteration. The direc- tion is correct if χ 2 (a + ∆a) < χ 2 (a) after iteration. Then update a + ∆a→a and decrease λ µ in the next iteration. Otherwise, update a and increase λ µ .

Optimized Fitting Method for Flange Annular End Face
The weight w i is theoretically independent of each observation value z Mi and is inversely proportional to the square of the standard uncertainty of the probability distribution function of z Mi . However, the uncertainty cannot be directly calculated in practical applications. Therefore, iterative operations can be used to obtain a more accurate estimated weight. The iterative mathematical model makes the singular points deviate more significantly. Thus, singular points can be removed accurately before weighted fitting is performed again. The optimized iterative reweighted LS process for the flange annular end face is shown in Figure 6. method. Thus, the iterative adjustment matrix Δa is revised to: The damping factor µ is determined by the adaptive coefficient λµ and the maximum diagonal element nmax of the quadratic normal matrix , as shown in Equation (18). The empirical value of λµ can be set to 0.001 in the first iteration. The direction is correct if χ 2 (a + Δa) < χ 2 (a) after iteration. Then update a + Δa→a and decrease λµ in the next iteration. Otherwise, update a and increase λµ.

Optimized Fitting Method for Flange Annular End Face
The weight wi is theoretically independent of each observation value zMi and is inversely proportional to the square of the standard uncertainty of the probability distribution function of zMi. However, the uncertainty cannot be directly calculated in practical applications. Therefore, iterative operations can be used to obtain a more accurate estimated weight. The iterative mathematical model makes the singular points deviate more significantly. Thus, singular points can be removed accurately before weighted fitting is performed again. The optimized iterative reweighted LS process for the flange annular end face is shown in Figure 6.

Iterative Estimation of Weights Based on Residuals
The residual Δi is approximately mapped to the weight wi. When the observed value zMi under (xMi, yMi) is very close to the model function, it causes 1/Δi to become too high. It is therefore necessary to set a reasonable upper limit value for wi: Figure 6. The process of optimized iterative reweighted LS fitting method with singular values removed.

Iterative Estimation of Weights Based on Residuals
The residual ∆ i is approximately mapped to the weight w i . When the observed value z Mi under (x Mi , y Mi ) is very close to the model function, it causes 1/∆ i to become too high. It is therefore necessary to set a reasonable upper limit value for w i : The residual ∆ i obeys a normal distribution due to the influence of independent random factors. Its lower bound value λ L is determined by Equation (20). To simplify operations, λ L can be set as the median of absolute residual |∆ i |, causing 50% of the observed values to have the same weight, 1/λ 2 L .
where ∧ σ wz M is the estimated value of the standard deviation of the weighted observation value z Mi and w i /N is the normalized mean weight value.
For the first round of fitting, i.e., the initial model parameter estimation method with equal weight, the reader can refer to Appendix A. Then, w i can be related to ∆ i and ∧ σ wz M in the second round. The number of measuring points N for each flange is less than 100 to ensure measurement efficiency. Then, the high value of ∧ σ wz M causes λ L to become too high. It is difficult to distinguish observations with low |∆ i | from the measuring points. Therefore, it is necessary to re-estimate the weight values after obtaining the mathematical model in the last round of fitting. One should then repeat the iterations between point cloud fitting and weight estimation until the changes in the model parameters are small enough.

Standard Residual Singular Value Detection Based on the Chauvenet Criterion
If there is a singular value in the measuring points, its weight needs to be zeroed: Singular points are found by calculating thresholds. Under the LS framework, the observed mean z Mi in the standard deviation formula corresponds to the estimated value. There are three elements in the parameter vector a of the flange end face. The standard deviation ∧ σ z M is estimated by Equation (23).
The Chauvenet criterion defines the probability of observed values occurring within the mean bandwidth range. The probability P in of ∆ i occurring within the interval where erf is the Gauss error function. Set the nominal number of occurrences of singular values to v 0 = NP out . The relationship between the coefficient к 0 and the number of measuring points N is shown in Equation (26).
where erfinv is the Gauss inverse error function.

Clustering Detection
The principle of clustering detection requires that qualified observations be classified into a one-dimensional category. Their absolute deviations Λ i = |∆ i | are relatively concentrated. Simultaneously, the deviations of singular values are far away from this category. Thus, the specified interval value is searched for in the deviation sorting and is used as the criterion for singular value determination. Cluster detection can avoid misjudgment caused by the sparsity of observation distribution.
The absolute deviations Λ i are arranged in ascending order and renumbered ordinal number n: · · · ≤ Λ j ≤ · · · ≤ Λ k ≤ · · · ≤ Λ l ≤ · · · The difference of absolute deviations between adjacent orders is: There is a significant gap between the absolute deviations of singular values and the absolute deviations before their orders. The global criterion is that the difference of absolute deviations is higher than the typical threshold, i.e., d[n] ≥ к 1 d glob [n]. The local criterion is that the difference in absolute deviations is higher than the threshold generated by partial fore correlation values, i.e., d[n] ≥ к 2 d loc [n]. The threshold λ 0 is set as Equation (30).
The global threshold d glob [n] and local threshold d loc [n] are defined as the weighted average of specific differences, shown as Equations (31) and (32). The weight decreases as the absolute deviation of the current measuring point increases and the slope is related to the number of measurement points N. The d loc [n] mainly depends on the influence of the differences of the adjacent absolute deviations. Local changes in the differences are overestimated when the absolute deviations are equal. d loc [n] should be adjusted using Equation (33

Optimized Fitting Method for Spigot Cylindrical Face
The coordinates of the projection points M proj_i on the fitted end face of the spigot measuring point M i are calculated as Equation (34).
To obtain the spatial circle center O proj_M (b 1 ,b 2 ,b 3 ), the fitting equation after mapping is shown as Equation (35).
Compared with the traditional LS method, an optimized total least squares (TLS) method is proposed to minimize the square of the distance between the mapping of observation points and the tangent points. The tangent point is on the tangent line of the model and is orthogonal to the minimum distance vector of the measuring point and the model. Thus, TLS has more geometric significance for fitting the mapped spatial circle of the spigot, as shown in Figure 7. Due to the symmetry of the circular curve, weight estimation and singular value detection are unnecessary, and observations with relatively higher errors can also be used as conditional inputs. Referring to the definition of χ 2 (b), the fitted model is transformed into f x proj_Mi , y proj_Mi , z proj_Mi | b ≡ 0 by shifting terms and the end face also provides a limiting condition function. The objective function of TLS is shown in Equation (36).
Compared with the traditional LS method, an optimized total least squares (TLS) method is proposed to minimize the square of the distance between the mapping of observation points and the tangent points. The tangent point is on the tangent line of the model and is orthogonal to the minimum distance vector of the measuring point and the model. Thus, TLS has more geometric significance for fitting the mapped spatial circle of the spigot, as shown in Figure 7. Due to the symmetry of the circular curve, weight estimation and singular value detection are unnecessary, and observations with relatively higher errors can also be used as conditional inputs. Referring to the definition of χ 2

Measurement Experiment
Measurement and fitting experiments for an aeroengine low-pressure turbine shaft are performed. The three-axis measurement system includes: measuring head RENISHAW OMP 40-2 (a high-precision three-dimensional contact probe that can be used on coordinate measuring machines, with unidirectional repeatability of 1 µm and can be equipped with various types of measuring needles) and optical interface OMI-2T (an anti-interference infrared optical receiver for probe signals, with a detection range of 6 m). The initial position and pose of the workpiece in each experiment is randomized using an adjustment fixture. The measuring points are arranged according to the proposed planning schemes. The traditional fitting methods and the proposed optimized fitting methods are used to process the same set of measurement data, respectively. The singular point detection methods are verified and the results, such as goodness of fit and model parameter uncertainty, are compared.

Measurement and Fitting Experiment of Flange end Face
The flange end face measurement experiment is shown in Figure 8. In each measurement experiment of the rotor with random position and pose, five fitting schemes are used to fit the same coordinate data of the measurement points. The processes of different fitting schemes are shown in Figure 9.  Taking one experiment as an example, the coordinates of the measuring points are listed in Table 3. The initial estimated model parameter vector ainit is obtained after equal weight fitting and the compensation step size ∆a calculated based on the weighting times in different fitting schemes are listed in Table 4. All Δai ≤ 0.01 in the fourth iteration of weight estimation. Thus, the iteration is terminated.   Taking one experiment as an example, the coordinates of the measuring points are listed in Table 3. The initial estimated model parameter vector ainit is obtained after equal weight fitting and the compensation step size ∆a calculated based on the weighting times in different fitting schemes are listed in Table 4. All Δai ≤ 0.01 in the fourth iteration of weight estimation. Thus, the iteration is terminated.  Taking one experiment as an example, the coordinates of the measuring points are listed in Table 3. The initial estimated model parameter vector a init is obtained after equal weight fitting and the compensation step size ∆a calculated based on the weighting times in different fitting schemes are listed in Table 4. All ∆a i ≤ 0.01 in the fourth iteration of weight estimation. Thus, the iteration is terminated.

Iterations Initial Values of Model Parameters and Iterative Compensation Values
Equal Due to iterative reweighting, the weight of potential singular points that deviate significantly decreases. The fitting end face is closer to the measuring points with smaller deviation. Thus, the singular points can be easier to detect. The numbers of singular points of five measurement experiments with rotor random positions and poses are listed in Table 5. Whether the measuring points have better weights has a smaller impact on the standard residual method, but a greater impact on the clustering detection method. Thus, it is necessary to perform iterative estimation of weights and model parameters before detecting singular points. The two methods have different judgment criteria, but both have achieved good results in detecting singular points and can achieve better fitting accuracy. Thus, both are reasonable methods for detecting singular points and can be used in parallel when the numerical control system has sufficient computing power. When the singular points are removed, the weight of the model parameter vector a is re-evaluated after four iterations of reweighting. ∆a i ≤ 0.01 is achieved after one compensation calculation. Thus, the iteration is terminated. The fitting accuracy results of the experimental example are listed in Table 6. Because the standard residual method and clustering detection method have removed the same singular points in this experimental example, the corresponding fitting accuracy results are equal. It can be seen that the weighted fitting obtains a goodness of fit, g fit , that is less than 1.000 and closer to 1.000 than the equal weight fitting. Thus, the fitting model has a better consistency with the measuring point data. Simultaneously, the uncertainties of model parameters δ a1 , δ a2 and δ a3 decrease to varying degrees and iterative reweighting further improves fitting accuracy results compared with the ordinary weighting. The fitting after singular value detection further reduces the uncertainty of the model parameters in a small range and all uncertainties δ a are less than 0.027. The fitting accuracy results of other experiments show similar patterns to this experiment example. In summary, "iterative reweighting + singular point detection (using residual method or cluster detection method) + iterative reweighting" is the optimal fitting scheme for the annular end face of the rotor flange.

Measurement and Fitting Experiment of the Spigot Cylindrical Face
The experiment in which we measured the spigot cylindrical face is shown in Figure 10. Each spigot measurement experiment is performed after the corresponding flange measurement experiment. The same coordinate data of each group of measuring points is fitted through three fitting schemes, including "OLS fitting", "TLS fitting", and "Iterative TLS fitting". experimental example are listed in Table 6. Because the standard residual method and clustering detection method have removed the same singular points in this experimental example, the corresponding fitting accuracy results are equal. It can be seen that the weighted fitting obtains a goodness of fit, gfit, that is less than 1.000 and closer to 1.000 than the equal weight fitting. Thus, the fitting model has a better consistency with the measuring point data. Simultaneously, the uncertainties of model parameters δa1, δa2 and δa3 decrease to varying degrees and iterative reweighting further improves fitting accuracy results compared with the ordinary weighting. The fitting after singular value detection further reduces the uncertainty of the model parameters in a small range and all uncertainties δa are less than 0.027. The fitting accuracy results of other experiments show similar patterns to this experiment example. In summary, "iterative reweighting + singular point detection (using residual method or cluster detection method) + iterative reweighting" is the optimal fitting scheme for the annular end face of the rotor flange.

Measurement and Fitting Experiment of the Spigot Cylindrical Face
The experiment in which we measured the spigot cylindrical face is shown in Figure 10. Each spigot measurement experiment is performed after the corresponding flange measurement experiment. The same coordinate data of each group of measuring points is fitted through three fitting schemes, including "OLS fitting", "TLS fitting", and "Iterative TLS fitting".
The coordinates of the measuring points are listed in Table 7. The measuring points are projected towards the fitted flange end face zM = 758.4004 − 0.0243xM + 0.0300yM. The iterative compensation ∆b using LS and TLS fitting methods are calculated by the initial estimate model parameter vector binit, which are listed in Table 8. In the iterative TLS fitting, the fourth iteration is terminated for all Δbi ≤ 0.01.  The coordinates of the measuring points are listed in Table 7. The measuring points are projected towards the fitted flange end face z M = 758.4004 − 0.0243x M + 0.0300y M . The iterative compensation ∆b using LS and TLS fitting methods are calculated by the initial estimate model parameter vector b init , which are listed in Table 8. In the iterative TLS fitting, the fourth iteration is terminated for all ∆b i ≤ 0.01. TLS is a nonlinear fitting, and the Jacobian matrix J LS contains the model parameter vector b itself. The orders of magnitude of the diagonal element calculation result of covariance matrix C are too large and distort model the parameter uncertainty δ b . With regard to the goodness of fit g fit , b is gathered to the same side of the model equation in TLS method and the observation value defaults to 0. This causes the calculation results of g fit to become too small to have significance in their evaluation. Therefore, the average orthogonal distance from the mapping point to the fitting model d TLS = f TLS (b)/N is used as the main accuracy evaluation indicator. The fitting accuracy results of different fitting methods of this experiment example are listed in Table 9. TLS fitting achieves smaller average orthogonal distance d TLS compared with the traditional LS fitting and d TLS can be further reduced after iterative calculation of ∆b. The fitting accuracy results of other experiments show similar patterns to this experiment example. In summary, the iterative TLS fitting method is the optimal scheme for the rotor spigot. The final fitting equation is calculated as Equation (42) Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The estimated initial model parameters of the flange end face can be calculated as follows: The estimated initial model parameters of the spigot cylindrical face can be calculated as follows: x proj_Mi − b 1, init 2 + y proj_Mi − b 2, init 2 + a 1 + a 2 x proj_Mi + a 3 y proj_Mi − b 3, init