A coordinate transformation method based on Rodrigues' formula

Abstract

This invention belongs to the field of surveying and mapping engineering technology, specifically relating to a coordinate transformation method based on Rodrigues' formula. Addressing the difficulty in determining the a priori unit weight mean square error in robust estimation, this invention designs a new weighting approach and establishes a robust coordinate transformation model. First, the linearization model of coordinate transformation based on Rodrigues' formula is derived in detail, and the method for calculating the initial parameter values ​​is given. Second, the parameters are solved using a parameter adjustment method with additional constraints. Finally, the initial a priori unit weight mean square error is calculated using median estimation, and a corresponding weight function is designed and continuously updated through an iterative process. This invention's method is not limited by the size of the rotation angle, the model is simple and computationally convenient, and it solves the problem of the difficulty in reasonably determining the a priori unit weight mean square error. Compared with existing methods, the proposed method has higher coordinate transformation accuracy and stronger robustness.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of surveying and mapping engineering, and particularly relates to a coordinate conversion method based on a Rodrigues formula. BACKGROUND

[0002] The engineering control network refers to a special surveying control network laid for engineering construction, is a control network laid for meeting various engineering construction, construction lofting and safety monitoring, and is established by using special surveying instruments such as a GPS receiver, a theodolite, a total station, a level, a laser tracker and the like, is laid in different forms according to the size of the project and the topography of the location where the project is constructed, the type of the project construction and the like, and guides the installation and adjustment of the project construction and equipment. The control network is composed of limited points in geometric space, the control points are independent of each other (spatial independent points) and are connected with each other (geometric relationship between points, mainly the position relationship), the control points are composed of three coordinate vectors, and O-XYZ three directions are perpendicular to each other.

[0003] Due to the continuous change of time and space conditions, combined with the increasingly diversified information acquisition technology of the spatial position of the target point, the coordinates of the spatial points in different coordinate systems often need to be unified, which is an important factor affecting the accuracy of the engineering control network. The coordinate conversion is mainly used for the unification between two coordinate systems, and the number of common points is generally not less than three. The coordinates before and after the conversion are geometric quantities in three-dimensional space, the coordinate system can be regarded as an engineering surveying coordinate system (including a surveying station coordinate system, a surveying / global coordinate system, an engineering design coordinate system and the like), the coordinate system is generally set by a person, is convenient for construction, and the Z axis is parallel to the plumb line. The specific applications include the following cases:

[0004] ①The coordinate system conversion between the control networks established by different surveying instruments or systems, such as the coordinate system conversion between the GPS network and the control network established by the total station and the level. The origin, direction and scale of the coordinate system defined by different surveying systems are different, so if the unification is required, the coordinate conversion must be performed.

[0005] ②The coordinate system conversion of the point coordinates of the same engineering control network in two periods. After a period of time, the points in the control network inevitably deform and displace, so it is generally required to periodically re-survey the control network. The reference may also deform during re-surveying, so the point coordinates after re-surveying are generally obtained by using the barycenter reference or the quasi-stable reference. At this time, the two period coordinates of the same control network are obtained, and the coordinate systems are different due to the different references. Then, the coordinate system can be unified by using the common point conversion method.

[0006] 3. When using a laser tracker and other measuring instruments to measure the tunnel control network, the tunnel is long and narrow, so the network is also long and narrow, and single-station measurement can only obtain control points within part of the range of view. The overall control network generally adopts the method of free station setting and station transfer splicing. When splicing between two stations, coordinate conversion is inevitably involved, which is also a relatively classic scenario.

[0007] With the deepening of research, the model is more simple and efficient. The common four-parameter and seven-parameter models have a large model error when the linearization model is used for large angles because the rotation matrix is composed of trigonometric functions. In addition, the rotation matrix sometimes becomes singular, leading to ill-conditioned normal equations, which are generally suitable for small-angle conversion.

[0008] When the two periods of data to be converted do not contain abnormal observations, the least squares (LS) method is an optimal estimation method for parameters. However, in actual measurement, the quality of observation data is inconsistent due to the influence of factors such as the operator's level, instrument accuracy, and external environment, and even gross errors are generated. In addition, the common points may move and deform between the two periods, at which time LS is no longer applicable, and other measures need to be taken to resist the interference of abnormal data.

[0009] The current methods for processing gross errors are roughly divided into gross error detection and robust estimation. Due to the low efficiency of gross error detection, the method of robust estimation is generally used to process gross error data. The commonly used method in robust estimation is M-estimation, and the selected weight iteration method is widely used in the field of measurement as an M-estimation method. M-estimation can reduce the influence of large errors on parameter calculation to a certain extent and has become a common solution, but the accurate processing of gross errors is limited by specific conditions and is related to the structure of the control network. The pre-test unit weight error based on least squares is unreliable, and the pre-test accuracy is difficult to accurately grasp, resulting in low coordinate conversion accuracy. SUMMARY

[0010] The present application provides a coordinate conversion method based on the Rodrigues formula to solve the problem of low coordinate conversion accuracy caused by the unreliable pre-test unit weight error of the prior art method in the process of engineering control network data coordinate conversion.

[0011] To solve the above technical problems, the technical solutions included in the present application and the beneficial effects corresponding to the technical solutions are as follows:

[0012] The coordinate conversion method based on the Rodrigues formula of the present application comprises the following steps:

[0013] 1) express the coordinate transformation model based on the Rodrigues formula, the coordinate transformation model contains 8 unknown parameters, respectively, 3 translation parameters Δx, Δy, Δz, 3 elements of unit vector a, b, c, 1 rotation parameter θ, and 1 scale parameter λ;

[0014] 2) determine the initial value of the 8 unknown parameters, the column matrix composed of the 8 unknown parameters is parameter Then, according to the initial value of the 8 unknown parameters, the coefficient matrix B and the constant term l are determined, and the error equation is listed;

[0015] 3) adopt the parameter adjustment method with the restriction condition to perform initial calculation, obtain the parameter And the cofactor matrix Q xx of the parameter ;

[0016] 4) according to the parameter obtained in step 3) and the listed error equation, determine the observation correction value V; estimate the pre-test unit weight error σ0, and according to the designed weight function, re-calculate the parameter And the cofactor matrix Q xx of the parameter ; judge whether the iteration termination condition is met, if not, update the coefficient matrix B, the constant term l and the error equation using the latest parameter , and re-execute step 4) until the iteration termination condition is met, and the final parameter is obtained, wherein the designed weight function is:

[0017]

[0018] In the formula, p i is the prior weight, and the numerical value is 1; ε1 and ε2 are coefficients, and ε1 < ε2; ω i is the weight factor, and v i is the observation correction value of the i point;

[0019] 5) determine the coordinate transformation model according to the final parameter ; and perform coordinate transformation according to the determined coordinate transformation model.

[0020] The beneficial effects of the above technical solutions are that the application proposes a new weight setting idea, and establishes a robust model for coordinate transformation based on the new weight function, the method is not limited by the size of the rotation angle, the model is simple and convenient to calculate, solves the problem of difficult determination of the pre-test unit weight error, and designs the corresponding weight factor to resist the interference of abnormal data, thereby improving the accuracy, reliability and robustness of the coordinate transformation.

[0021] Further, in order to obtain reliable scale parameters to improve the accuracy of the first adjustment result, the initial value of the scale parameter λ in step 2) is determined by the following method:

[0022] The coordinates of any two common points in the two sets of coordinate systems are calculated by the following distance calculation:

[0023]

[0024] In the formula, (x i ,y i ,z i ) is the coordinate of the common point i in the first set of coordinate systems; (x j ,y j ,z j ) is the coordinate of the common point j in the first set of coordinate systems; (X i ,Y i ,Z i ) is the coordinate of the common point i in the second set of coordinate systems; (X j ,Y j ,Z j ) is the coordinate of the common point j in the second set of coordinate systems.

[0025] The average value of the distance calculation results of any two common points is obtained, and the average value is the initial value of the scale parameter λ.

[0026] Further, in order to obtain reliable rotation parameters to improve the accuracy of the first adjustment result, the initial values of the rotation parameter θ and the three elements a, b, and c of the unit vector are calculated by the following formula:

[0027]

[0028] In the formula, θ0 is the initial value of the rotation parameter θ; tr(·) represents the trace operator; R is a rotation matrix, and a i is a parameter in the rotation matrix, i=1,2,…,9; a0, b0, and c0 are the initial values of the three elements a, b, and c of the unit vector, respectively.

[0029] Further, in order to obtain reliable translation parameters to improve the accuracy of the first adjustment result, the initial values of the three translation parameters Δx, Δy, and Δz are calculated by the following formula:

[0030]

[0031] In the formula, ΔX0 is the initial value of the translation matrix composed of the three translation parameters Δx, Δy, and Δz; λ0 is the initial value of the scale parameter λ; R0 is the initial value of the rotation matrix composed of the rotation parameter and the three translation parameters Δx, Δy, and Δz; X1 is the coordinate of the first phase; X2 = Hξ, ξ = [ΔxΔyΔz a1 a2 a3 a4 a5 a6 a7 a8 a9] T , i X2 represents the coordinates of point i in the second phase.

[0032] Furthermore, to reduce the interference of outlier data, in step 4), the median estimation method is used to estimate the prior unit weight standard error σ0:

[0033]

[0034] In the formula, Let Q be the cofactor of the correction for the observation at point i, and Q vv =Q LL -BQ xx B T Q LL Q is the cofactor matrix of the observations; xx The parameter is the cofactor matrix.

[0035] Furthermore, in step 3), the limiting condition is:

[0036]

[0037] In the formula, C = (0 0 0 0 2a0 2b0 2c0 0), where a0, b0, and c0 are the initial values ​​of the three elements a, b, and c of the unit vector, respectively.

[0038] Furthermore, the parameters are obtained using the following formula.

[0039]

[0040] In the formula, N aa =B T PB, where P is the weight matrix; W = B T Pl;

[0041] Furthermore, the parameters cofactor array Q xx :

[0042]

[0043] Further, the order of determining the initial value of the scale parameter, the initial value of the translation parameter and the initial value of the rotation parameter is: first determining the scale parameter, then determining the rotation parameter and the three elements of the unit vector, and finally determining the translation parameter; and the initial value of the parameter that has been solved will be used as the known value for the subsequent parameter initial value calculation.

[0044] Further, in step 4), the iteration termination condition is that the difference between the parameters before and after is less than the corresponding threshold value, and the translation parameter threshold value is 1×10 -3 mm, and the scale parameter and rotation parameter threshold values are set to 1×10 -8 . BRIEF DESCRIPTION OF DRAWINGS

[0045] Figure 1-1 is a schematic diagram of a single set of control point layout of the present application;

[0046] Figure 1-2 is a schematic diagram of the overall distribution of control points of the present application;

[0047] Figure 2 is a schematic diagram of the position distribution of control points of the present application;

[0048] Figure 3 is a schematic diagram of the point position residual of the present application. DETAILED DESCRIPTION

[0049] The present application derives a coordinate conversion linearization model based on the Rodrigues formula, and gives a calculation method of the parameter initial value. When weighting, a new strategy is adopted to determine the pre-test unit weight error to overcome the difficulty in determining the pre-test accuracy, and a corresponding weight factor is designed to resist the interference of abnormal data, thereby improving the accuracy and reliability of coordinate conversion.

[0050] The following will first introduce the various aspects of the coordinate conversion model based on the Rodrigues formula used in the present application.

[0051] 1. Mathematical model.

[0052] The mathematical model of coordinate conversion can be expressed as:

[0053] X2 = ΔX + λRX1 (1)

[0054] In the formula, X1 and X2 respectively represent the three-dimensional coordinates of a batch of points in two sets of coordinate systems; ΔX is the translation parameter; λ is the scale parameter; and R is the rotation matrix.

[0055] For any common point in the coordinate system, expand formula (1) to obtain:

[0056]

[0057] In the formula, a iare the parameters in the rotation matrix R, i = 1, 2, …, 9.

[0058] A point A(x, y, z) in three-dimensional space rotates around a unit vector by an angle θ to get A'(X, Y, Z), the mathematical relationship can be represented by the Rodrigues formula, that is:

[0059]

[0060] where A' = (X, Y, Z) T ; A = (x, y, z) T ; and

[0061]

[0062] Expanding equation (3) and expressing it in matrix form, we have:

[0063]

[0064] Equation (5) represents that the two coordinate systems in which point A is located have the same origin and scale, and the coefficient matrix is the rotation matrix R, which is composed of the unit vector and the rotation angle θ, containing 4 parameters, i.e. (a, b, c, θ), and a, b, c satisfy the following constraint conditions:

[0065] a 2 +b 2 +c 2 = 1 (6)

[0066] Without loss of generality, equation (5) can be expressed as a general coordinate transformation model, that is:

[0067] A' = ΔX + λRA (7)

[0068] There are 8 unknown parameters in equation (7), including 3 translation parameters Δx, Δy, Δz, 3 elements a, b, c of the unit vector, 1 rotation parameter θ, and 1 scale parameter λ. Considering equation (6) as a constraint condition for the 3 elements, there are 7 independent parameters, and at least 3 points are needed to solve the parameters.

[0069] 2. Error equation.

[0070] As can be seen from equation (5), the 4 unknown parameters contained in the rotation matrix R are in a nonlinear relationship in the equation, so equation (7) is linearized, the Taylor series expansion is carried out at its approximate value and the first order term is retained, then:

[0071]

[0072] In the formula, (λ0, Δx0, Δy0, Δz0, a0, b0, c0, θ0) is the approximate value of unknown parameters; (dλ, dΔx, dΔy, dΔz, da, db, dc, dθ) represents the parameter correction number. The error equation can be listed from formula (8), that is:

[0073]

[0074] In the formula, V is the coordinate correction number under the target coordinate system; the column matrix composed of 8 unknown parameters is parameter B is the coefficient matrix; l is the constant term, which is expressed as:

[0075]

[0076] Based on the least square principle, 3 or more common points can be used to solve the unknown parameters. Formula (6) is used as the constraint condition of the parameters, and is linearized and expressed as:

[0077]

[0078] In the formula, C=(0 0 0 0 2a0 2b0 2c0 0);

[0079] Combined with formula (9), the parameter adjustment method with the attached constraint condition can be used to solve the unknown parameters. The expression of parameter can be solved by the following formula, that is:

[0080]

[0081] In the formula, N aa =B T PB, P is the weight matrix; W=B T Pl; Therefore, the coordinate conversion correction number V can be obtained according to formula (9). Further, the post-test unit weight variance estimate can be expressed as:

[0082]

[0083] In the formula, f represents the degree of freedom; n represents the number of common points; t is the number of parameters; s is the number of constraint condition equations; and o is the number of observation values with zero weight. According to the propagation law of variance and covariance, the cofactor matrix of parameter is:

[0084]

[0085] 3. Solution of parameter initial value.

[0086] In the process of parameter solving, the initial value directly affects the result of the first adjustment, and the weight iteration method always weights according to the residual error of the previous adjustment result. Improper initial value may even lead to non-convergence of the result, so the reliability and accuracy of the initial value of the parameter is crucial. To obtain a more reliable initial value of the scale parameter, translation parameter and rotation parameter, the following strategies are adopted: first, calculate the scale parameter; then, take measures to eliminate the translation parameter, and temporarily ignore the correlation between parameters to solve the rotation parameter; finally, obtain the translation parameter. In the process of parameter solving, the solved parameters will be used as known values for subsequent parameter calculation.

[0087] The scale parameter can be calculated using the distance between corresponding coordinate points in two sets of coordinate systems, i.e.:

[0088]

[0089] In the formula, (x i ,y i ,z i ) are the coordinates of common point i in the first set of coordinate systems; (x j ,y j ,z j ) are the coordinates of common point j in the first set of coordinate systems; (X i ,Y i ,Z i ) are the coordinates of common point i in the second set of coordinate systems; (X j ,Y j ,Z j ) are the coordinates of common point j in the second set of coordinate systems.

[0090] For two sets of data with n common points, according to formula (15), we can obtain scale parameters, and take the average value as the initial value λ0 of the scale parameter λ, which is:

[0091]

[0092] For any common point, formula (2) can also be expressed as:

[0093]

[0094] In the formula, the parameter vector ξ = [Δx Δy Δz a1 a2 a3 a4 a5 a6 a7 a8 a9] T .

[0095] The matrix form of formula (17) is:

[0096] X2 = Hξ (18)

[0097] In the formula,

[0098] From equation (17), it is known that there are 13 parameters in the equation. λ is calculated from equation (16). Without considering the correlation of the 9 elements in the rotation matrix, at least 4 points are needed to solve the remaining parameters.

[0099] The first three columns of the H matrix are the same. After eliminating the first three columns of the equation, H' is obtained. The specific steps are as follows: ① Form 3n equations according to equation (17) using n points; ② Choose 1 point, and subtract the corresponding 3 equations of the point from the corresponding equations of the remaining points; ③ After arranging, 3(n-1) equations are obtained. At this time, the 9 parameters contained in the equation correspond to the 9 elements of the rotation matrix R one by one. The first 3 elements of the vector ξ are eliminated at this time. The remaining parameters form a vector ξ', which can be calculated by the following formula:

[0100] ξ' = (H' T H' -1 (H' T X2 (19)

[0101] The ξ' obtained from equation (19) contains the 9 elements of R. The initial value of the rotation parameter is extracted from R by the following formula:

[0102]

[0103] In the formula, tr(·) represents the trace operator.

[0104] The initial value of the rotation matrix R0 is composed of the initial value (θ0, a0, b0, c0). Combined with equation (1), the initial value of the translation parameter ΔX0 can be calculated by the following formula:

[0105]

[0106] In the formula, X1 is the first period coordinate.

[0107] At this point, the parameter initial value calculation is completed. The above is the complete process of parameter initial value calculation. According to the order of scale, rotation and translation parameters, the parameter initial value is calculated based on the full use of all observation data. On the one hand, the parameters that have been solved are used as known information, which can increase the accuracy and stability of the subsequent parameter calculation results; on the other hand, taking the average of all data will greatly reduce the interference of abnormal data, thereby increasing the reliability of the parameter results.

[0108] 4. Pre-test accuracy estimation.

[0109] The M estimation method solves the parameters by iteration. In abnormal cases, the empirical value or theoretical value of the pre-test unit weight error σ0 cannot reflect the accuracy of the actual observation value, making it difficult to determine the pre-test accuracy. To reduce the interference of abnormal data, σ0 is obtained by the median estimation method, and med(·) represents the median operator, i.e.:

[0110]

[0111] In the formula, Q vv = Q LL -BQ xx B T ; is the observation value correction factor of the i point, Q LL is the observation value correlation matrix; Q xx is the parameter correlation matrix; v i is the observation value correction value of the i point.

[0112] According to formula (22), the median estimation method can resist the influence of abnormal data to a certain extent, but the residual error and its correlation matrix are equal-weight processing results, and the influence of abnormal values cannot be avoided, which leads to inaccurate σ0. To make the estimation more robust and greatly suppress the influence of abnormal values on the residual error and its correlation matrix, σ0 is updated according to formula (22) in the iteration process until convergence, so as to obtain the pre-test accuracy without the interference of abnormal data. When calculating, the weight factor is constructed according to σ0 and the residual error, that is:

[0113]

[0114] The weight determination rule is: when |v|≤ε1·σ0, it is regarded as normal observation data; when ε1·σ0<|v|≤ε2·σ0, it is regarded as low-quality observation data; and when |v|>ε2·σ0, it is regarded as a gross error. Wherein, the value range of ε1 and ε2 is 1.5-2.0 and 3.0-5.0 respectively. The corresponding equivalent weight function is represented as:

[0115]

[0116] In the formula, p i is the prior weight, and the value is 1.

[0117] The above weight determination considers: 1) through iterative calculation, the weight distribution tends to be reasonable, the influence of abnormal data on the result gradually decreases, and σ0 becomes smaller and converges. 2) As can be seen from formula (24), the critical value of abnormal data determination is in a multiple relationship with σ0, and the critical value also decreases as σ0 decreases. On the one hand, the examination and screening of abnormal data are more strict, which can effectively exclude harmful information. On the other hand, since σ0 is constantly decreasing in iteration, to avoid the situation that the critical value is too low to cause the effective information to be unable to be fully utilized, the critical value of abnormal data determination can be adjusted appropriately, so as to reasonably distribute the weight and increase the reliability of the result.

[0118] Based on the above introduction, a coordinate conversion method based on the Rodrigues formula will be described in detail in combination with the drawings and embodiments.

[0119] Method embodiment:

[0120] An embodiment of the coordinate conversion method based on the Rodrigues formula of the present application is aimed at the coordinate conversion of the point coordinates before and after the two periods of the tunnel control network of a particle accelerator, and the whole method process is as follows:

[0121] Step one, express the coordinate conversion model based on the Rodrigues formula.

[0122] Step two, calculate the initial value of the parameter. The unknown parameters are three translation parameters (Δx, Δy, Δz), three elements (a, b, c) of the unit vector, one rotation parameter θ, and one scale parameter λ. The initial value of the parameter is obtained according to formulas (15)-(20), and is denoted as (λ0, a0, b0, c0, θ0, Δx0, Δy0, Δz0).

[0123] Step three, list the error equation. Substitute the initial value of the parameter determined in step one into formulas (8) and (10) to obtain the coefficient matrix B and the constant term I, and list the error equation according to formula (9).

[0124] Step four, first solution of the parameter. In combination with the constraint condition equation (11), the initial weight matrix is set as I 3n , the parameter adjustment method with the constraint condition is used to solve the conversion parameter, the parameter is calculated by formula (12), and the cofactor matrix Q of the parameter xx is calculated by formula (14).

[0125] Step five, select the weight iteration method to solve the parameter. Substitute the parameter obtained in step three into formula (9) to obtain the observation correction V, calculate σ0 by formula (22), weight according to formula (24), and then use the parameter adjustment method again, i.e., use formulas (12) and (14) to solve the parameter and the cofactor matrix Q of the parameter xx . The pre-test unit weight error σ0 is recalculated according to formula (22) in each iteration process, and the latest parameter result is used to update formulas (9)-(11). The iteration termination condition in this embodiment is that the parameter difference before and after step four is less than the corresponding threshold value, and the iteration is stopped; wherein the threshold value of the translation parameter is set to 1×10 -3 mm, and the threshold values of the scale parameter and the rotation parameter are set to 1×10 -8 .

[0126] The effectiveness of the method of the present application is verified by the simulation test and the example analysis as follows.

[0127] 1. Simulation test.

[0128] Taking the linear accelerator tunnel control network as an example, a network consisting of P1 to P2 was designed. 24 A three-dimensional control network consisting of 24 points was constructed. Following the actual measurement scenario, the instrument center was located near the geometric center of the control network. A laser tracker was used to simulate measurements of the 24 points, providing the first set of data. The point coordinates are shown in Table 1. The layout of a single set of control points is as follows: Figure 1-1 As shown, the overall distribution of control points is as follows: Figure 1-2 As shown.

[0129] Table 1 Control point coordinates / mm

[0130]

[0131] Figure 1-2 This demonstrates the spatial layout of a single set of control points within the tunnel and the distribution of six sets of control points in the coordinate system of the tracking instrument station. The interval between adjacent sets of points is 7m, and the longest side does not exceed 18m, ensuring the instrument can perform its measurement accuracy effectively. The scale remains unchanged between the second and first phases of data, with a translation of (5.0, 8.0, 0.3)m and a rotation vector... The rotation angle is 50°. According to the Leica AT901-B laser tracker, the angular measurement accuracy within a 20m measurement range is ±(15μm + 6×10). -6 Random errors were added to the original observations of the two sets of data using the nominal accuracy of the ranging accuracy ±(0.5μm / m·D) and the distance measurement accuracy ±(0.5μm / m·D) to obtain the simulation test data.

[0132] Using the LS method to process the experimental data, the maximum absolute value of the coordinate component residuals reached 0.185 mm, and the observation accuracy σ of the coordinate components was considered to be ±0.185 mm. Based on the simulation data, errors were added to the coordinate components to conduct experiments, and the following four schemes were designed to generate simulation data.

[0133] Option 1: The data is left unprocessed and contains only random error.

[0134] Option ②: Considering that 2σ < 0.4 mm < 3σ, at P1, P 10 P 15 P 20 Adding an error of 0.400 mm to the X-coordinate component of each point is considered as low-quality observation data.

[0135] Option ③: P1, P 10 P 15 P 20 An error of 1.0 mm is added to the X-coordinate component of each point, which is considered as gross error data.

[0136] Option 4: P1, P 10 Add an error of 0.400mm to the X-coordinate component of each point, P 15 P20 The error of 1.0 mm is added to the X coordinate component of the point, respectively, which is considered as the case of coexistence of gross error data and low quality observation data.

[0137] The method of the present application is compared with LS, Huber scheme, Tukey scheme and IGG3 scheme (the IGG3 scheme refers to the M estimation method with the equivalent weight function being the IGG3 weight function), and the weight functions corresponding to each scheme are shown in Table 2.

[0138] Table 2 Equivalent weight function

[0139]

[0140] From Table 2, the critical value e = 4.685σ0 in the Tukey function and the critical value e = 4.685σ0 in the IGG3 function For the standardized residual, the value ranges of k0 and k1 are 1.0-1.5 and 2.5-8.0, respectively, and the critical value of the equivalent weight function is adjusted according to the weight setting result during calculation, so as to make the result reasonable. The theoretical value is taken as σ0 during weight setting, and formula (13) is commonly used in practice. Since the method of the present application is relatively strict in abnormal data screening, in order to ensure the reliability of the result, ε1 and ε2 are taken as 2.0 and 5.0, respectively, during calculation. In order to analyze the accuracy of the coordinate conversion parameters corresponding to several methods, the difference between the parameter calculation result and the true value is shown in Table 3.

[0141] Table 3 Difference between parameter result and true value of different methods

[0142]

[0143]

[0144] From Table 3, the result is analyzed as follows.

[0145] Analysis of scheme ①: under the condition of containing only random error, the parameter results of the five methods are close to the true value, dλ does not exceed -0.2×10 -6 , the difference between the translation parameters in each direction and the true value is not more than 0.100 mm, and dθ is about 2″.

[0146] Analysis of schemes ②-④: after adding gross error, the deviation of the parameter results of the five methods from the true value increases. Especially in scheme ③, the deviation of dλ corresponding to the Huber scheme from the true value is the largest, reaching 1.2×10 -6 ; for the translation and rotation parameters, the deviations of all results from the true value are about 0.100 mm and 10×10 -6 , respectively, and the deviations of dz and da of the IGG3 scheme are 0.122 mm and -10.1×10 -6 , respectively, which are the largest among all results. However, the largest dλ of the method of the present application in the four schemes is only -0.2×10-6 and dΔx is not more than 0.010 mm, and the parameters are stable.

[0147] To verify the robustness of the several methods, the corresponding results are shown in Table 4.

[0148] Table 4 different methods Comparative / mm

[0149]

[0150] From Table 4, it can be seen that, since LS does not have robustness, when there are low-quality observation data or gross errors in the data, the LS result is obviously larger than the other four methods. The results of the four schemes in the table are analyzed as follows.

[0151] Analysis of scheme ①: when the data only contains random errors, the Huber scheme is equivalent to the LS result, while the Tukey scheme, the IGG3 scheme and the method of the application are smaller. Through analysis, the weight of some observation values is reduced when the weight is determined in the latter three schemes.

[0152] Analysis of schemes ②-④: after adding errors to the data, the Huber scheme is closest to the LS result, and the Tukey scheme is second, indicating that the resistance ability of the two schemes to abnormal data is poor. In contrast, the IGG3 scheme and the method of the application are close in results, and are smaller than the Huber and Tukey schemes, indicating that the IGG3 scheme and the method of the application have strong robustness. In contrast, the result of the method of the application is the smallest among the four robustness schemes, indicating that the method of the application is optimal and has stronger robustness.

[0153] To verify the reliability of the pre-test precision obtained by the method of the application, the results of σ0 before and after iteration are shown in Table 5.

[0154] Table 5 results of σ0 before and after iteration / mm

[0155]

[0156] From Table 5, it can be seen that, whether or not to add errors to the data, the value of σ0 becomes smaller after iteration. Observing schemes ②-④, σ0 before iteration becomes larger as the order of magnitude of the added error increases, indicating that the median estimate is only resistant to the influence of abnormal data to a certain extent; the σ0 result after iteration is relatively stable, and is about 0.048 mm, which is obviously less affected by abnormal data, verifying the effectiveness of the method of the application in determining the pre-test precision.

[0157] 2. Example analysis.

[0158] ​​​Two period data of a ring tunnel control network are selected, and a total of 44 homonymic points are contained, and the control point distribution is as shown in the figure Figure 2 . As can be seen from Figure 2 , there is obvious position offset and rotation between the two period data, the LS is used to process the data, and the common point position residual is obtained, and the result is as shown in the figure Figure 3 . As can be seen from Figure 3 , the point position residual of BD512 point is the largest, and the value reaches 2.873mm, so it can be regarded as a gross error point, and the point position residual of the remaining points is within 1.000mm, and further processing is needed. After removing the BD512 point, the data is processed again by using the LS, and the obtained is reduced from 0.402mm to 0.309mm, and the maximum value of the point position residual is 0.919mm, which is within 1.000mm. After removing the BD512 point, the coordinate conversion precision is obviously improved, which indicates that the point is a gross error point.

[0159] The data is processed by using four kinds of robust methods, the LS processing result is combined, and corresponding parameters are adjusted, so that the result is reasonable, and the obtained is as shown in Table 6.

[0160] Table 6 different methods

[0161]

[0162] As can be seen from Table 6, the Huber scheme, the Tukey scheme and the IGG3 scheme are all greater than the method, which preliminarily indicates that the method has the highest post-test precision. In the calculation process, the values of ε1 and ε2 are 2.0 and 4.0 respectively, and the test result shows that the weight values of the three coordinate components of the BD512 point are all 0, and the weight value distribution result ensures the contribution of most normal observation values to the parameter result, and reduces the influence of low-quality observation data. In addition, the values of k0 and k1 of the IGG3 scheme are 1.5 and 3.0 respectively, and the treatment of the gross error point is consistent with the result of the method, and the treatment of other points is similar to the method. In comparison, the Huber scheme and the Tukey scheme only reduce the weight of the corresponding coordinate component of the BD512 point when treating the BD512 point, and do not completely remove it, which indicates that the screening of abnormal data is not strict.

[0163] In summary, the test result reflects that the weight value distribution result of the method conforms to the actual situation, and can accurately remove the gross error and reasonably reduce the weight of low-quality observation data. Compared with other methods, the post-test unit weight error obtained by the method is the smallest, and the robustness is stronger.

[0164] In summary, the application deduces a linearization model of coordinate transformation based on the Rodrigues formula and gives a calculation method of the initial value of parameters. In view of the problem that the pre-test unit weight error is difficult to determine reasonably, a new weight design idea is designed and a robust model of coordinate transformation is established. Taking the particle accelerator tunnel control network as an example, simulation test and example analysis are carried out, and comparison is made with the existing method. The results show that the method of the application is not limited by the rotation angle, the model is simple and efficient, the pre-test accuracy does not need to be known during calculation, and the problem that the pre-test unit weight error is difficult to determine reasonably is solved. Compared with the existing method, the weight design method of the application is simple, the pre-test accuracy result is reliable, the weight distribution result is more reasonable, and the robustness is stronger.

[0165] In this embodiment, the point coordinates of the same engineering control network (particle accelerator tunnel control network) measured in two periods are used for coordinate transformation. As other embodiments, the coordinate transformation between control networks established by different measuring instruments or systems, or the coordinate transformation when two stations are spliced can also be used.

Claims

1. A method for coordinate transformation based on the Rodrigues' formula, characterized in that, The method comprises the following steps: 1) expressing a coordinate transformation model based on a Rodrigues formula, wherein the coordinate transformation model comprises eight unknown parameters, i.e., three translation parameters Δx, Δy, Δz, three elements a, b, c of a unit vector, one rotation parameter θ, and one scale parameter λ; 2) Determine the initial value of 8 unknown parameters, which form a column matrix of parameters Further, according to the initial value of the 8 unknown parameters, determine the coefficient matrix B and the constant term l, and list the error equation; 3) using the parameter adjustment method with limited conditions to make the initial calculation, and obtaining the parameters and the cofactor matrix Q xx of parameters ; 4) the parameters obtained according to step 3) and the listed error equation to determine the observation correction V; estimate the pre-test unit weight error σ0, and according to the designed weight function, the parameter adjustment method is used to re-calculate to obtain the parameters and the parameter co-factor matrix Q xx ; determine whether the iteration termination condition is met, if not, update the coefficient matrix B, the constant term l, and the error equation with the latest parameters , and re-execute step 4) until the iteration termination condition is met, and the final parameters are obtained , wherein the set weight function is: where p i is a priori weight, and the values are all 1; ε1 and ε2 are coefficients, and ε1 < ε2, the value ranges of ε1 and ε2 are 1.5-2.0 and 3.0-5.0 respectively; ω i is a weight factor, and v i is an observation value correction value of the i point. 5) in accordance with the final parameters determining a coordinate conversion model; and performing coordinate conversion of the coordinates of the points before and after the control network of the particle accelerator tunnel in accordance with the determined coordinate conversion model.

2. The coordinate conversion method based on the Rodrigues' formula according to claim 1, characterized in that, In step 2), the initial value of the scale parameter λ is determined by the following method: The coordinates of any two common points in the two sets of coordinate systems are calculated by the following distance calculation: wherein (x i ,y i ,z i ) is the coordinate of common point i in the first period coordinate system; (x j ,y j ,z j ) is the coordinate of common point j in the first period coordinate system; (X i ,Y i ,Z i ) is the coordinate of common point i in the second period coordinate system; (X j ,Y j ,Z j ) is the coordinate of common point j in the second period coordinate system; The average value of the distance calculation results of any two common points is obtained, and the average value is the initial value of the scale parameter λ.

3. The coordinate conversion method based on the Rodrigues' formula according to claim 1, wherein, The initial values of the rotation parameter θ and the three elements a, b, c of the unit vector are calculated by the following formula: where θ0is an initial value of the rotation parameter θ; tr(·) denotes a trace operator; R is a rotation matrix, and a i are parameters in the rotation matrix, i = 1, 2, …, 9; a0, b0, c0are initial values of the three elements a, b, c of the unit vector, respectively.

4. The coordinate conversion method based on the Rodrigues' formula according to claim 1, wherein, The initial values of the three translation parameters Δx, Δy, Δz are calculated by the following formula: wherein ΔX0is the initial value of the translation matrix composed of 3 translation parameters Δx, Δy, Δz; λ0is the initial value of the scale parameter λ; R0is the initial value of the rotation matrix composed of the rotation parameters and 3 translation parameters Δx, Δy, Δz; X1is the first period coordinate; X2=Hξ, ξ=[ΔxΔyΔz a1a2a3a4a5a6a7a8a9] T , i X2is the second period corresponding i point coordinate.

5. The coordinate conversion method based on the Rodrigues' formula according to claim 1, wherein, In step 4), the initial value of the scale parameter λ is determined by the following method: where is the i observation correction factor for the covariates, and Q vv = Q LL -BQ xx B T , Q LL is the observation covariates matrix; Q xx is the parameter covariates matrix.

6. The coordinate conversion method based on the Rodrigues' formula according to claim 1, wherein, In step 3), the restriction condition is: In the formula, C=(0 0 0 0 2a0 2b0 2c0 0), a0, b0, and c0 are initial values of the 3 elements a, b, and c of the unit vector, respectively; 7. The coordinate conversion method based on the Rodrigues' formula according to claim 6, wherein, The parameters are obtained using the following formula where N aa = B T PB, P is a weighting matrix; W = B T Pl; 8. The coordinate conversion method based on the Rodrigues' formula according to claim 7, wherein, The parameters The cofactor matrix Q xx :

9. The method of claim 1, wherein, The determination order of the initial value of the scale parameter, the initial value of the translation parameter, and the initial value of the rotation parameter is: first determining the scale parameter, then determining the rotation parameter and the three elements of the unit vector, and finally determining the translation parameter; and the initial value of the determined parameter will be used as the known value for the subsequent parameter initial value calculation.

10. The coordinate conversion method based on the Rodrigues' formula according to any one of claims 1 to 9, characterized in that, In step 4), the iteration termination condition is that the difference between the parameters of the previous and the current iteration is smaller than the corresponding threshold, and the threshold for the translation parameter is 1 x 10 - 3 mm, and the thresholds for the scale and rotation parameters are set to 1 x 10 -8 .

Images