Method for implementing high-precision orientation and evaluating orientation precision of large-scale dynamic photogrammetry system

Abstract

The present invention provides a method for implementing high-precision orientation and evaluating orientation precision of a large-scale dynamic photogrammetry system, including steps: a) selecting a scale bar, arranging code points at two ends of the scale bar, and performing length measurement on the scale bar; b) evenly dividing a measurement space into multiple regions, sequentially placing the scale bar in each region, and photographing the scale bar by using left and right cameras; d) limiting self-calibration bundle adjustment by using multiple length constraints, adjustment parameters including principal point, principal distance, radial distortion, eccentric distortion, in-plane distortion, exterior orientation parameter and spatial point coordinate; and e) performing traceable evaluation of orientation precision of the photogrammetry system. The present invention can effectively reduce the relative error in length measurement.

Claims

1. A method for implementing high-precision orientation and evaluating orientation precision of a large-scale dynamic photogrammetry system, comprising steps: a) selecting a scale bar, arranging code points at two ends of the scale bar, and performing length measurement on the scale bar; b) evenly dividing a measurement space into multiple regions, sequentially placing the scale bar in each region, and photographing the scale bar by using left and right cameras, wherein continuous photographing is carried out in the following manner: b1) adjusting left and right cameras so that fields of view of the left and right cameras fully cover the scale bar, synchronously moving the left and right cameras along a perpendicular direction of a line connecting the code points at the two ends of the scale bar, and at the same time rotating the scale bar by using the perpendicular direction of the line connecting the code points at the two ends of the scale bar as an axis of rotation; and b2) adjusting the left and right cameras, synchronously moving the left and right cameras along the perpendicular direction of the line connecting the code points at the two ends of the scale bar in a direction opposite to that in the step b1), and at the same time rotating the scale bar in a direction opposite to that in the step b1); c) performing preliminary relative positioning and orientation of the photogrammetry system according to images shot by the left and right cameras; d) limiting self-calibration bundle adjustment by using multiple length constraints, adjustment parameters comprising principal point, principal distance, radial distortion, eccentric distortion, in-plane distortion, exterior orientation parameter and spatial point coordinate; e) performing traceable evaluation of orientation precision of the photogrammetry system.

2. The method according to claim 1, wherein a distance between the left and right cameras is the same as a region length according to which the measurement space is evenly divided.

3. The method according to claim 1, wherein the photographing process of the step b) is repeated multiple times.

4. The method according to claim 1, wherein the step c comprises the following steps: c1) determining a solution space of a fundamental matrix by using a five-point method: automatically selecting five imaging points in a position of the scale bar, wherein positions of the imaging points are a center and four vertexes of the shot image, and the automatic selection is performed using the following algorithm:
abs(xl.sub.i)+abs(yl.sub.i)+abs(xr.sub.i)+abs(yr.sub.i).fwdarw.min
(xl.sub.i)+(yl.sub.i)+(xr.sub.i)+(yr.sub.i).fwdarw.max
(xl.sub.i)+(yl.sub.i)+(xr.sub.i)+(yr.sub.i).fwdarw.max
(xl.sub.i)+(yl.sub.i)+(xr.sub.i)+(yr.sub.i).fwdarw.max
(xl.sub.i)+(yl.sub.i)+(xr.sub.i)+(yr.sub.i).fwdarw.max wherein [xli yli] represents image plane coordinates of the i.sup.th point captured by the left camera, and [xri yri] represents image point coordinates of the i.sup.th point captured by the right camera; and c2) obtaining a solution of an essential matrix by using an elimination method.

5. The method according to claim 4, wherein the obtaining of the solution of the essential matrix in the step c2) comprises the following key steps: c21) building a tenth-degree polynomial of an unknown number w; c22) solving a real root of the tenth-degree polynomial in the step c21; and c23) judging the solution of the essential matrix.

6. The method according to claim 5, wherein in a method of building the tenth-degree polynomial of w in the step c21, a Toeplitz matrix is used: c211) describing all polynomials as an array:
a.sub.mx.sup.m+A.sub.m-1x.sup.m-1+ . . . +a.sub.2x.sup.2+a.sub.1x+a.sub.0 A=[a.sub.ma.sub.m-1 . . . a.sub.2a.sub.1a.sub.0].sup.T
b.sub.nx.sup.n+b.sub.n-1x.sup.n-1+ . . . +b.sub.2x.sup.2+b.sub.1x+b.sub.0 B=[b.sub.nb.sub.n-1 . . . b.sub.2b.sub.1b.sub.0].sup.T; and c212) performing programming by using the Toeplitz matrix to compute polynomial multiplication, with a formula being: $C = T * B = [\begin{matrix} a_{m} & 0 & 0 & .Math. & 0 \\ a_{m - 1} & a_{m} & 0 & .Math. & 0 \\ a_{m - 2} & a_{m - 1} & a_{m} & .Math. & 0 \\ .Math. & .Math. & .Math. & .Math. & .Math. \\ 0 & a_{0} & a_{1} & .Math. & a_{m} \\ .Math. & .Math. & .Math. & .Math. & .Math. \\ 0 & 0 & 0 & 0 & a_{0} \end{matrix}] [\begin{matrix} b_{n} \\ b_{n - 1} \\ .Math. \\ b_{1} \\ b_{0} \end{matrix}]$ wherein T is a Toeplitz matrix corresponding to a polynomial A, a quantity of columns in T is equal to a quantity of elements in B, and a quantity of rows in T is equal to (m+1)+(n+1)1.

7. The method according to claim 6, wherein the step c212 further comprises: performing programming by using the Toeplitz matrix to compute polynomial division, specifically comprising the following steps: c213) describing n as n=T(d)q+r by using a Toeplitz matrix of d, wherein it is set that a numerator polynomial is n, a denominator polynomial is d, a quotient is q, and a remainder is r; and c214) calculating an optimal solution of a quotient of the polynomial division: q=(T.sup.TT).sup.1T.sup.Tn.

8. The method according to claim 5, wherein a method for solving the real root of the tenth-degree polynomial in the step c22 is: c221) building a Sturm sequence of the tenth-degree polynomial, with a formula being:
p.sub.0(x)=p(x)
p.sub.1(x)=p(x)
p.sub.2(x)=rem(p.sub.0,p.sub.1)
p.sub.3(x)=rem(p.sub.1,p.sub.2)
.
.
.
0=rem(p.sub.m-1,p.sub.m) wherein rem represents calculating the remainder of the polynomial, and P(x) is a known tenth-degree polynomial; c222) searching all single intervals of the polynomial by recursive dichotomy in combination with the Sturm's Theorem; and c223) after the single intervals are obtained, calculating a numerical solution of a real root of |C(w)|=0 in each single interval by dichotomy.

9. The method according to claim 1, wherein the step d comprises the following steps: d1) obtaining a camera imaging model; d2) building an error equation of the photogrammetry system; d3) performing adaptive proportional adjustment of error equation terms; d4) performing a partitioning operation on the normal equation; d5) obtaining self-calibration bundle adjustment with length constraints; and d6) evaluating precision of the photogrammetry system.

10. The method according to claim 1, wherein the step e of performing traceable evaluation of orientation precision of the photogrammetry system is performed by using direct linear transformation, and comprises the following steps: e1) obtaining an error equation of a single imaging point; e2) when the photogrammetry system comprises two cameras, obtaining an error equation set corresponding to a spatial point [X Y Z]; e3) obtaining a corresponding normal equation and coordinate correction amount, and optimizing spatial point coordinates through multiple iterations; e4) obtaining a measured value of a corresponding length of the scale bar through the step e3; e5) obtaining an average value of measurement results of lengths of the scale bar in all orientations; e6) performing scaling by using a ratio of an average length to a gage length as an orientation result, and performing an uncertainty analysis on the length measurement result, to obtain an orientation precision evaluation; e7) obtaining a quality parameter of a measurement instrument; e8) calculating a relative error of each length; and e9) calculating an uncertainty of the relative error in length measurement.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) To describe the technical solutions in the embodiments of the present invention or the prior art more clearly, the drawings that need to be used in the embodiments or the prior art are briefly introduced. It would be obvious that the drawings in the following description are merely some embodiments of the present invention, and those of ordinary skill in the art may further obtain other drawings according to these drawings without creative efforts.

(2) FIG. 1 is an experimental assembly of a method for high-precision positioning and orientation of a large-scale multi-camera photogrammetry system according to the present invention;

(3) FIG. 2 is a schematic experimental diagram of the system;

(4) FIG. 3 shows images of a scale bar in two cameras in the system; and

(5) FIG. 4 is a diagram showing the relative positions between scale bars at different positions and attitudes in space and two cameras.

DETAILED DESCRIPTION OF THE INVENTION

(6) The objectives and functions of the present invention and methods for achieving such objectives and functions will be illustrated through exemplary embodiments. However, the present invention is not limited to the exemplary embodiments disclosed below, but may be implemented in different forms. The essence of this specification is merely for the purpose of helping those skilled in the art to comprehensively understand the specific details of the present invention.

(7) The present invention provides a large-scale dynamic photogrammetry system. The system includes a camera 101, a power source 102, a network cable 103, a signal line 104 and a host computer 105. The power source 102 supplies power to the camera 1, the camera 2 and the host computer 105. A network cable 103a is connected between the camera 1 and the host computer 105, and a network cable 103b is connected between the camera 2 and the host computer 105, for providing communication between the cameras and the host computer. A signal line 104a is connected to the camera 1 and the host computer 105, and a signal line 104b is connected to the camera 2 and the host computer 105, for transmitting signals between the host computer 105 and the camera 101. The host computer 105 is further equipped with a USB signal adapter. The host computer 105 sends a synchronization signal via the USB signal adapter 106, and controls the two cameras 101 via the signal line 104 to synchronously capture images.

(8) FIG. 2 is a schematic experimental diagram of the system. The measurement space 201 mainly includes two cameras 202 and a scale bar 203. The measurement space 201 provides an experimental verification environment for the method for implementing high-precision orientation and evaluating orientation precision of the system.

(9) The present invention provides a method for implementing high-precision orientation and evaluating orientation precision of a large-scale dynamic photogrammetry system, including steps:

(10) a) selecting a scale bar, arranging code points at two ends of the scale bar, and performing length measurement on the scale bar.

(11) b) evenly dividing a measurement space into multiple regions, sequentially placing the scale bar in each region, and photographing the scale bar by using left and right cameras, where continuous photographing is carried out in the following manner:

(12) b1) adjusting left and right cameras so that fields of view of the left and right cameras fully cover the scale bar, synchronously moving the left and right cameras along a perpendicular direction of a line connecting the code points at the two ends of the scale bar, and at the same time rotating the scale bar by using the perpendicular direction of the line connecting the code points at the two ends of the scale bar as an axis of rotation; and

(13) b2) adjusting the left and right cameras, synchronously moving the left and right cameras along the perpendicular direction of the line connecting the code points at the two ends of the scale bar in a direction opposite to that in the step b1), and at the same time rotating the scale bar in a direction opposite to that in the step b1);

(14) c) performing preliminary relative positioning and orientation of the photogrammetry system according to images shot by the left and right cameras, where in the photographing process, a distance between the left and right cameras is the same as a region length according to which the measurement space is evenly divided, and preferably, in some embodiments, the photographing process of the step b) is repeated multiple times.

(15) d) limiting self-calibration bundle adjustment by using multiple length constraints, adjustment parameters including principal point, principal distance, radial distortion, eccentric distortion, in-plane distortion, exterior orientation parameter and spatial point coordinate.

(16) e) performing traceable evaluation of orientation precision of the photogrammetry system.

(17) Based on the method of the step b, the scale bar of the present invention is designed and calibrated through the following operations:

(18) The main body of the scale bar is a carbon fiber bar, and two backlight reflection code points are fixed at two ends of the bar as endpoints for length measurement. During calibration, the scale bar at a particular position is imaged in two cameras. As shown in FIG. 3, white point areas are backlight reflection points 301 adhered at two ends of the scale bar.

(19) As shown in FIG. 4, the three-dimensional space includes coordinates 402 of the two cameras and multiple scale bars 401 at different positions, and backlight reflection points 403 are adhered at two ends of the scale bar and used for calibrating the length of the scale bar. The number of images shot should be not less than 60.

(20) Step c:

(21) According to an embodiment of the present invention, the step c includes the following steps:

(22) c1) determining a solution space of a fundamental matrix by using a five-point method: automatically selecting five imaging points in a position of the scale bar, where positions of the imaging points are a center and four vertexes of the shot image, and the automatic selection is performed using the following algorithm:
abs(xl.sub.i)+abs(yl.sub.i)+abs(xr.sub.i)+abs(yr.sub.i).fwdarw.min
(xl.sub.i)+(yl.sub.i)+(xr.sub.i)+(yr.sub.i).fwdarw.max
(xl.sub.i)+(yl.sub.i)+(xr.sub.i)+(yr.sub.i).fwdarw.max
(xl.sub.i)+(yl.sub.i)+(xr.sub.i)+(yr.sub.i).fwdarw.max
(xl.sub.i)+(yl.sub.i)+(xr.sub.i)+(yr.sub.i).fwdarw.max(3.1)

(23) where [xli yli] represents image plane coordinates of the i.sup.th point captured by the left camera, and [xri yri] represents image point coordinates of the i.sup.th point captured by the right camera;

(24) The point satisfying the above formulas is a point closest to the expected position.

(25) The corresponding image point satisfies:
X.sup.TEX=0(3.2)

(26) where E is an essential matrix, and X and X are respectively coordinates obtained by normalizing respective intrinsic parameters of corresponding points on left and right images.

(27) Five points are selected and substituted into the above formula, to obtain the following equation set:

(28) $\begin{matrix} [\begin{matrix} X_{1}^{} X_{1} & X_{1}^{} Y_{1} & X_{1}^{} Z_{1} & Y_{1}^{} X_{1} & Y_{1}^{} Y_{1} & Y_{1}^{} Z_{1} & Z_{1}^{} X_{1} & Z_{1}^{} Y_{1} & Z_{1}^{} Z_{1} \\ X_{2}^{} X_{2} & X_{2}^{} Y_{2} & X_{2}^{} Z_{2} & Y_{2}^{} X_{2} & Y_{2}^{} Y_{2} & Y_{2}^{} Z_{2} & Z_{2}^{} X_{2} & Z_{2}^{} Y_{2} & Z_{2}^{} Z_{2} \\ .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. \\ X_{5}^{} X_{5} & X_{5}^{} Y_{5} & X_{5}^{} Z_{5} & Y_{5}^{} X_{5} & Y_{5}^{} Y_{5} & Y_{5}^{} Z_{5} & Z_{5}^{} X_{5} & Z_{5}^{} Y_{5} & Z_{5}^{} Z_{5} \end{matrix}] [\begin{matrix} e_{11} \\ e_{12} \\ e_{13} \\ e_{21} \\ e_{22} \\ e_{23} \\ e_{31} \\ e_{32} \\ e_{33} \end{matrix}] = [\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}] & (3.3) \end{matrix}$

(29) A fundamental system of solutions of the above formula is obtained by singular value decomposition (SVD), thus extending a solution space of the essential matrix:
E=wE.sub.w+xE.sub.x+yE.sub.y+E.sub.z(3.4)

(30) c2) obtaining a solution of an essential matrix by using an elimination method;

(31) The essential matrix has the following two natures:
|E|=0
2EE.sup.TEtr(EE.sup.T)E=0(3.5)

(32) Ten equations of unknown numbers w, x, y may be obtained by using the above formulas. The ten equations are rewritten by using x, y, x.sup.2, xy, y.sup.2, x.sup.3, x.sup.2y, xy.sup.2 and y.sup.3 as unknown variables and was a known variable, to eliminate x and y. The following equation set is obtained:

(33) $\begin{matrix} \underset{10 10}{C (w)} [\begin{matrix} 1 \\ x \\ y \\ x^{2} \\ x y \\ y^{2} \\ x^{3} \\ x^{2} y \\ x y^{2} \\ y^{3} \end{matrix}] = 0 & (3.6) \end{matrix}$

(34) where C(w) is a polynomial matrix of w.

(35) The homogeneous linear equation set described by the formula (3.6) has a non-zero solution, and therefore the determinant of C(w) is zero. Theoretically, a polynomial expressed by the determinant of C(w) is described, and a root of the polynomial is solved, thus obtaining numerical solutions of w and C(w); further, values of x and y can be obtained by calculating a non-zero solution of the formula (3.6). This inevitably involves the description of the polynomial of a 10-order symbolic matrix determinant, and how to solve a real root of a higher-degree polynomial.

(36) According to an embodiment of the present invention, the obtaining of the solution of the essential matrix in the step c2) includes the following key steps:

(37) c21) building a tenth-degree polynomial of an unknown number w;

(38) A solution of the determinant is calculated by programming using the following method which is suitable for both symbolic polynomial and numerical operations:

(39) m=1;

(40) for k=1: n1 for I=k+1: n for j=k+1: n A(i, j)=(A(k, k)A(i, j)A(I, k)A(k, j))/m;

(41) end end m=A(k, k);

(42) end return A(n, n);

(43) The above algorithm involves polynomial multiplication and division.

(44) According to an embodiment of the present invention, in a method of building the tenth-degree polynomial of w in the step c21, programming and computation are performed by using a Toeplitz matrix, with a method being:

(45) c211) describing all polynomials as an array:
a.sub.mx.sup.m+A.sub.m-1x.sup.m-1+ . . . +a.sub.2x.sup.2+a.sub.1x+a.sub.0 A=[a.sub.ma.sub.m-1 . . . a.sub.2a.sub.1a.sub.0].sup.T
b.sub.nx.sup.n+b.sub.n-1x.sup.n-1+ . . . +b.sub.2x.sup.2+b.sub.1x+b.sub.0 B=[b.sub.nb.sub.n-1 . . . b.sub.2b.sub.1b.sub.0].sup.T (3.7)

(46) c212) performing programming by using the Toeplitz matrix to compute polynomial multiplication, with a formula being:

(47) $\begin{matrix} C = T * B = [\begin{matrix} a_{m} & 0 & 0 & .Math. & 0 \\ a_{m - 1} & a_{m} & 0 & .Math. & 0 \\ a_{m - 2} & a_{m - 1} & a_{m} & .Math. & 0 \\ .Math. & .Math. & .Math. & .Math. & .Math. \\ 0 & a_{0} & a_{1} & .Math. & a_{m} \\ .Math. & .Math. & .Math. & .Math. & .Math. \\ 0 & 0 & 0 & 0 & a_{0} \end{matrix}] [\begin{matrix} b_{n} \\ b_{n - 1} \\ .Math. \\ b_{1} \\ b_{0} \end{matrix}] & (3.8) \end{matrix}$

(48) where T is a Toeplitz matrix corresponding to a polynomial A, a quantity of columns in T is equal to a quantity of elements in B, and a quantity of rows in T is equal to (m+1)+(n+1)1.

(49) According to an embodiment of the present invention, the step c212 may further include: performing programming by using the Toeplitz matrix to compute polynomial division, specifically including the following steps:

(50) c213) describing n as n=T(d)q+r by using a Toeplitz matrix of d,

(51) where it is set that a numerator polynomial is n, a denominator polynomial is d, a quotient is q, and a remainder is r; and

(52) c214) calculating an optimal solution of a quotient of the polynomial division: q=(T.sup.TT).sup.1T.sup.Tn.

(53) The method of solving the polynomial division by using the Toeplitz matrix does not have the error accumulation phenomenon of a long division method and can avoid the ill-conditioned problem of division of a large number by a small number. The obtained result is an optimal solution in the least squares sense.

(54) c22) solving a real root of the tenth-degree polynomial in the step c21;

(55) Solving the real root of the higher-degree polynomial can rely only on a numerical method such as Newton's method and dichotomy, but the key problem lies in determining a single interval.

(56) According to an embodiment of the present invention, single intervals are searched by using the Sturm's Theorem in the present invention. For a known polynomial P(x), its Sturm sequence may be determined:
p.sub.0(x)=p(x)
p.sub.1(x)=p(x)
p.sub.2(x)=rem(p.sub.0,p.sub.1)
p.sub.3(x)=rem(p.sub.1,p.sub.2)
.
.
.
0=rem(p.sub.m-1,p.sub.m)(3.9)

(57) where rem represents calculating the remainder of the polynomial.

(58) The Sturm's Theorem describes the number of real roots of the polynomial P(x) within an interval (a, b]. Assuming that the sign change frequency of the function value of the Sturm sequence at the endpoint a is c(a) and the sign change frequency of the function value of the Sturm sequence at the endpoint b is c(b), the number of real roots of the polynomial P(x) within this half-open interval is c(a)c(b). In the present technology, a recursive dichotomy technique is designed, and all single intervals are searched in a wide numerical range (a, b], with pseudo-code of the corresponding function being shown as follows:

(59) determining the Sturm sequence sturmSequence of the polynomial P(x);

(60) function [rootIntervals]=rootIntervalsDetect(a, b, sturmSequence)

(61) determining the number of real roots, numRoot, within (a, b] by using the Sturm's Theorem;

(62) rootRegion=[ ];

(63) if numRoot<1

(64) elseif numRoot==1 rootRegion=[rootRegion; a b];

(65) elseif numRoot>1 a1=a; b1=(a+b)/2; a2=(a+b)/2; b2=b; o1=rootIntervalsDetect(a1, b1, sturmSequence); rootRegion=[rootRegion; o1]; o2=rootIntervalsDetect(a2, b2, sturmSequence); rootRegion=[rootRegion; o2];

(66) end

(67) After single intervals are obtained, the real root of |C(w)|=0 is solved in each single interval by dichotomy. The solved w is substituted into the expression of C(w) in the formula (3.6), unknown terms about x and y are solved by SVD, and values of x and y are further obtained. The solved w, x and y are substituted into the formula (3.4), to obtain the essential matrix E, the number of solutions of which is the same as the number of real roots of w.

(68) c23) judging the solution of the essential matrix.

(69) Each w has a camera position relationship and corresponding reconstructed spatial point coordinates. In special cases, for example, imaging of a planar scene, a feasible solution cannot be judged according to the conventional criterion of smallest image plane error, because there are two solution having similar image plane errors. To fundamentally solve the root judgment problem, the present technology uses spatial information as a constraint and a length ratio of two scale bars at orthogonal positions as a criterion of judgment. An essential matrix whose reconstructed space length ratio is closest to 1 is the final feasible solution.

(70) Further, the translation amount and the spatial coordinates are scaled according to the ratio of the space length to the reconstructed length. The obtained exterior orientation parameter and spatial coordinates may be used as initial parameter values for subsequent self-calibration bundle adjustment.

(71) Step d:

(72) According to an embodiment of the present invention, the step d includes the following steps:

(73) d1) obtaining a camera imaging model;

(74) The camera orientation is described by using [X.sub.0 Y.sub.0 Z.sub.0 Az El Ro].sup.T, coordinates of a spatial point are [X Y Z] .sup.T, and coordinates [X Y Z] .sup.T in a camera coordinate system are described as:

(75) $\begin{matrix} [\begin{matrix} \overline{X} \\ \overline{Y} \\ \overline{Z} \end{matrix}] = R ([\begin{matrix} X \\ Y \\ Z \end{matrix}] - [\begin{matrix} X_{0} \\ Y_{0} \\ Z_{0} \end{matrix}]) = [\begin{matrix} a_{1} & b_{1} & c_{1} \\ a_{2} & b_{2} & c_{2} \\ a_{3} & b_{3} & c_{3} \end{matrix}] ([\begin{matrix} X \\ Y \\ Z \end{matrix}] - [\begin{matrix} X_{0} \\ Y_{0} \\ Z_{0} \end{matrix}]) & (4.1) \\ where \\ R = [\begin{matrix} \begin{matrix} \cos (Ro) \cos (Az) + \\ \sin (Ro) \sin (El) \sin (Az) \end{matrix} & \begin{matrix} \cos (Ro) \sin (Az) - \\ \sin (Ro) \sin (El) \cos (Az) \end{matrix} & \sin (Ro) \cos (El) \\ \begin{matrix} - \sin (Ro) \cos (Az) + \\ \cos (Ro) \sin (El) \sin (Az) \end{matrix} & \begin{matrix} - \sin (Ro) \sin (Az) - \\ \cos (Ro) \sin (El) \cos (Az) \end{matrix} & \cos (Ro) \cos (El) \\ \cos (El) \sin (Az) & - \cos (El) \cos (Az) & - \sin (El) \end{matrix}] & (4.2) \end{matrix}$

(76) Linear projection coordinates of this spatial point are:

(77) $\begin{matrix} \begin{matrix} x_{l} = - f \frac{\overline{X}}{\overline{Z}} = - f \frac{a_{1} (X - X_{0}) + b_{1} (Y - Y_{0}) + c_{1} (Z - Z_{0})}{a_{3} (X - X_{0}) + b_{3} (Y - Y_{0}) + c_{3} (Z - Z_{0})} \\ y_{l} = - f \frac{\overline{Y}}{\overline{Z}} = - f \frac{a_{2} (X - X_{0}) + b_{2} (Y - Y_{0}) + c_{2} (Z - Z_{0})}{a_{3} (X - X_{0}) + b_{3} (Y - Y_{0}) + c_{3} (Z - Z_{0})} \end{matrix} & (4.3) \end{matrix}$

(78) An object distance of it relative to the imaging system is:
s=Z(4.4)

(79) Researches show that a radial distortion parameter of a spatial point with an object distance of s on the imaging plane of the imaging system is:

(80) $\begin{matrix} k_{1 {ss}^{}} =_{s^{}} \frac{C_{s^{}}^{2}}{C_{s_{1}}^{2}} k_{1 {ss}_{1}} + (1 -_{s^{}}) \frac{C_{s^{}}^{2}}{C_{s_{2}}^{2}} k_{1 {ss}_{2}} & (4.5) \\ k_{2 {ss}^{}} =_{s^{}} \frac{C_{s^{}}^{4}}{C_{s_{1}}^{4}} k_{2 {ss}_{1}} + (1 -_{s^{}}) \frac{C_{s^{}}^{4}}{C_{s_{2}}^{4}} k_{2 {ss}_{2}} \\ k_{3 {ss}^{}} =_{s^{}} \frac{C_{s^{}}^{6}}{C_{s_{1}}^{6}} k_{3 {ss}_{1}} + (1 -_{s^{}}) \frac{C_{s^{}}^{6}}{C_{s_{2}}^{6}} k_{3 {ss}_{2}} \\ _{s^{}} = \frac{S_{2} - S^{}}{S_{2} - S_{1}} .Math. \frac{S_{1} - f}{S^{} - f} \end{matrix}$

(81) where s is the focus length of the imaging system, k.sub.1ss.sub.1, k.sub.2ss.sub.1, k.sub.3ss.sub.1 and k.sub.1ss.sub.2, k.sub.2ss.sub.2, k.sub.3ss.sub.2 are calibration results of the radial distortion parameter over two distances, and C.sub.s.sub.1 and C.sub.s.sub.2 are respectively image distances of Gaussian imaging formulas corresponding to the two distances.

(82) Assume that principal point offsets of the camera is x.sub.p and y.sub.p, and eccentric distortion and in-plane distortion parameters corresponding to this spatial point are P.sub.1, P.sub.2, B.sub.1, B.sub.2. The image point distortion amount at (x.sub.l, y.sub.l) is:

(83) $\begin{matrix} x = \overline{x} (K_{1 {ss}^{}} r^{2} + K_{2 {ss}^{}} r^{4} + K_{3 {ss}^{}} r^{6}) + P_{1} (2 {\overline{x}}^{2} + r^{2}) + 2 P_{2} \overline{xy} + B_{1} \overline{x} + B_{2} \overline{y} y = \overline{y} (K_{1 {ss}^{}} r^{2} + K_{2 {ss}^{}} r^{4} + K_{3 {ss}^{}} r^{6}) + P_{2} (2 {\overline{y}}^{2} + r^{2}) + 2 P_{1} \overline{xy} \begin{matrix} \overline{x} = x - x_{p} & \overline{y} = y - y_{p} & r = \sqrt{{\overline{x}}^{2} + {\overline{y}}^{2}} \end{matrix} & (4.6) \end{matrix}$

(84) Therefore, final image point coordinates of the spatial point on the camera image plane at this position is:
x=x.sub.l+x.sub.px
y=y.sub.l+y.sub.py custom character {tilde over (L)}=({tilde over (X)})(4.7)

(85) The exterior orientation parameter and the spatial coordinates are correlated to the radial distortion parameter through s.

(86) d2) building an error equation of the photogrammetry system;

(87) for a dual-camera photogrammetry system, using (xij, yij) to represent imaging coordinates of the j.sup.th point in space for the i.sup.th image, the error equation is:

(88) $\begin{matrix} {\hat{l}}_{ij} = l_{ij} + v_{ij} = K_{ij} {\hat{x}}_{ij} .Math. [\begin{matrix} x_{ij} -_{x} (X_{ij}^{0}) \\ y_{ij} -_{y} (X_{ij}^{0}) \end{matrix}] + [\begin{matrix} v_{xij} \\ y_{yij} \end{matrix}] = A_{ij} [\begin{matrix} X_{0 i} \\ Y_{0 i} \\ Z_{0 i} \\ A z_{i} \\ E l_{i} \\ R o_{i} \\ x_{0 i} \\ y_{0 i} \\ f_{i} \\ K_{1 i} \\ K_{2 i} \\ K_{3 i} \\ P_{1 i} \\ P_{2 i} \\ B_{1 i} \\ B_{2 i} \end{matrix}] + B_{ij} [\begin{matrix} X_{j} \\ Y_{j} \\ Z_{j} \end{matrix}] = A_{ij}_{i} + B_{ij} {\dot{}}_{j} i = 1, 2; j = 1, 2, .Math. n & (4.8) \end{matrix}$

(89) where X.sub.ij.sup.0 is all parameters related to the imaging point (x.sub.ij, y.sub.ij), (v.sub.xij, v.sub.yij) is a residual error, and A.sub.ij and B.sub.ij are Jacobian matrices of an imaging model (4.7) for the exterior orientation parameter of the i.sup.th image, coordinates of the j.sup.th spatial point, and the camera imaging parameter.

(90) 0 $\begin{matrix} A_{ij} = [\begin{matrix} \frac{x_{ij}}{X_{0 i}} & \frac{x_{ij}}{Y_{0 i}} & \frac{x_{ij}}{Z_{0 i}} & \frac{x_{ij}}{A z_{i}} & \frac{x_{ij}}{E l_{i}} & \frac{x_{ij}}{R o_{i}} & \frac{x_{ij}}{x_{p}} & \frac{x_{ij}}{y_{p}} & \frac{x_{ij}}{f} & \frac{x_{ij}}{K_{1 {ss}_{2}}} & \frac{x_{ij}}{K_{2 {ss}_{2}}} & \frac{x_{ij}}{K_{3 {ss}_{2}}} & \frac{x_{ij}}{P_{1}} & \frac{x_{ij}}{P_{2}} & \frac{x_{ij}}{B_{1}} & \frac{x_{ij}}{B_{2}} \\ \frac{y_{ij}}{X_{0 i}} & \frac{y_{ij}}{Y_{0 i}} & \frac{y_{ij}}{Z_{0 i}} & \frac{y_{ij}}{A z_{i}} & \frac{y_{ij}}{E l_{i}} & \frac{y_{ij}}{R o_{i}} & \frac{y_{ij}}{x_{p}} & \frac{y_{ij}}{y_{p}} & \frac{y_{ij}}{f} & \frac{y_{ij}}{K_{1 {ss}_{2}}} & \frac{y_{ij}}{K_{2 {ss}_{2}}} & \frac{y_{ij}}{K_{3 {ss}_{2}}} & \frac{y_{ij}}{P_{1}} & \frac{y_{ij}}{P_{2}} & \frac{y_{ij}}{B_{1}} & \frac{y_{ij}}{B_{2}} \end{matrix}] & (4.9) \\ B_{ij} = [\begin{matrix} \frac{x_{ij}}{X_{j}} & \frac{x_{ij}}{Y_{j}} & \frac{x_{ij}}{Z_{j}} \\ \frac{y_{ij}}{X_{j}} & \frac{y_{ij}}{Y_{j}} & \frac{y_{ij}}{Z_{j}} \end{matrix}] \end{matrix}$

(91) Each term in A.sub.ij is solved through the following process:

(92) $\begin{matrix} \frac{C_{s^{}}}{s^{}} = - \frac{f^{2}}{{(s^{} - f)}^{2}} \frac{_{s^{}}}{s^{}} = - \frac{(S_{1} - f) (S_{2} - f)}{(S_{2} - S_{1}) {(S^{} - f)}^{2}} \frac{k_{1 {ss}^{}}}{s^{}} = \frac{_{s^{}}}{s^{}} \frac{C_{s^{}}^{2}}{C_{s_{1}}^{2}} k_{1 {ss}_{1}} + 2 \frac{C_{s^{}}}{s^{}} \frac{_{s^{}} C_{s^{}}}{C_{s_{1}}^{2}} k_{1 {ss}_{1}} + (- \frac{_{s^{}}}{s^{}}) \frac{C_{s^{}}^{2}}{C_{s_{2}}^{2}} k_{1 {ss}_{2}} + 2 \frac{C_{s^{}}}{s^{}} \frac{(1 -_{s^{}}) C_{s^{}}}{C_{s_{2}}^{2}} k_{1 {ss}_{2}} \frac{k_{2 {ss}^{}}}{s^{}} = \frac{_{s^{}}}{s^{}} \frac{C_{s^{}}^{4}}{C_{s_{1}}^{4}} k_{2 {ss}_{1}} + 4 \frac{C_{s^{}}}{s^{}} \frac{_{s^{}} C_{s^{}}^{3}}{C_{s_{1}}^{4}} k_{2 {ss}_{1}} + (- \frac{_{s^{}}}{s^{}}) \frac{C_{s^{}}^{4}}{C_{s_{2}}^{4}} k_{2 s s_{2}} + 4 \frac{C_{s^{}}}{s^{}} \frac{(1 -_{s^{}}) C_{s^{}}^{3}}{C_{s_{2}}^{4}} k_{2 {ss}_{2}} \frac{k_{3 {ss}^{}}}{s^{}} = \frac{_{s^{}}}{s^{}} \frac{C_{s^{}}^{6}}{C_{s_{1}}^{6}} k_{3 {ss}_{1}} + 6 \frac{C_{s^{}}}{s^{}} \frac{_{s^{}} C_{s^{}}^{5}}{C_{s_{1}}^{6}} k_{3 {ss}_{1}} + (- \frac{_{s^{}}}{s^{}}) \frac{C_{s^{}}^{6}}{C_{s_{2}}^{6}} k_{3 s s_{2}} + 6 \frac{C_{s^{}}}{s^{}} \frac{(1 -_{s^{}}) C_{s^{}}^{5}}{C_{s_{2}}^{6}} k_{3 {ss}_{2}} & (4.10) \end{matrix}$

(93) It can be seen from the relationship described by the formula (4.4) that:

(94) $\begin{matrix} \frac{s^{}}{X_{0}} = a_{3} \frac{x}{s^{}} = \overline{x} (\frac{k_{1 {ss}^{}}}{s^{}} r^{2} + \frac{k_{2 {ss}^{}}}{s^{}} r^{4} + \frac{k_{3 {ss}^{}}}{s^{}} r^{6}) \frac{y}{s^{}} = \overline{y} (\frac{k_{1 {ss}^{}}}{s^{}} r^{2} + \frac{k_{2 {ss}^{}}}{s^{}} r^{4} + \frac{k_{3 {ss}^{}}}{s^{}} r^{6}) \frac{x}{X_{0}} = \frac{x}{s^{}} \frac{s^{}}{X_{0}} \frac{y}{X_{0}} = \frac{y}{s^{}} \frac{s^{}}{X_{0}} & (4.11) \end{matrix}$

(95) Through the description of the formula (4.3), partial derivatives of x.sub.l and y.sub.l relative to X.sub.0 can be solved:

(96) $\begin{matrix} \begin{matrix} \frac{x_{l}}{X_{0}} = \frac{f (a_{1} \overline{Z} - a_{3} \overline{X})}{{\overline{Z}}^{2}} & \frac{y_{l}}{X_{0}} = \frac{f (a_{2} \overline{Z} - a_{3} \overline{Y})}{{\overline{Z}}^{2}} \end{matrix} & (4.12) \end{matrix}$

(97) Therefore, with reference to the formula (4.7), a partial derivative of an observed value relative to the term X.sub.0 in the exterior orientation parameters is:

(98) $\begin{matrix} \begin{matrix} \frac{x}{X_{0}} = \frac{x_{l}}{X_{0}} - \frac{x}{X_{0}} & \frac{y}{X_{0}} = \frac{y_{l}}{X_{0}} - \frac{y}{X_{0}} \end{matrix} & (4.13) \end{matrix}$

(99) Likewise, partial derivatives of the observed value relative to other exterior orientation parameters can be obtained:

(100) $\begin{matrix} \begin{matrix} \frac{s^{}}{Y_{0}} = b_{3} & \frac{x}{X_{0}} = \frac{x}{s^{}} \frac{s^{}}{Y_{0}} & \frac{y}{Y_{0}} = \frac{y}{s^{}} \frac{s^{}}{Y_{0}} \end{matrix} \begin{matrix} \frac{x_{l}}{Y_{0}} = \frac{f (b_{1} \overline{Z} - b_{3} \overline{X})}{{\overline{Z}}^{2}} & \frac{y_{l}}{Y_{0}} = \frac{f (b_{2} \overline{Z} - b_{3} \overline{Y})}{{\overline{Z}}^{2}} \end{matrix} \begin{matrix} \frac{x}{Y_{0}} = \frac{x_{l}}{Y_{0}} - \frac{x}{Y_{0}} & \frac{y}{X_{0}} = \frac{y_{l}}{X_{0}} - \frac{y}{X_{0}} \end{matrix} & (4.14) \end{matrix}$

(101) $\begin{matrix} \begin{matrix} \frac{s^{}}{Z_{0}} = c_{3} & \frac{x}{Z_{0}} = \frac{x}{s^{}} \frac{s^{}}{Z_{0}} & \frac{y}{Z_{0}} = \frac{y}{s^{}} \frac{s^{}}{Z_{0}} \end{matrix} \begin{matrix} \frac{x_{l}}{Z_{0}} = \frac{f (c_{1} \overline{Z} - c_{3} \overline{X})}{{\overline{Z}}^{2}} & \frac{y_{l}}{Z_{0}} = \frac{f (c_{2} \overline{Z} - c_{3} \overline{Y})}{{\overline{Z}}^{2}} \end{matrix} \begin{matrix} \frac{x}{Z_{0}} = \frac{x_{l}}{Z_{0}} - \frac{x}{Z_{0}} & \frac{y}{Z_{0}} = \frac{y_{l}}{Z_{0}} - \frac{y}{Z_{0}} \end{matrix} & (4.15) \end{matrix}$

(102) The description of partial derivatives of linear terms x.sub.l and y.sub.l relative to different angles is complex, and reference can be made to relevant literatures. Herein, only partial derivatives of x and y relative to different angles are analyzed.

(103) $\begin{matrix} \frac{s^{}}{A z} = - (\cos (El) \cos (Az) (X - X_{0}) + \cos (El) \sin (Az) (Y - Y_{0})) \begin{matrix} \frac{x}{A z} = \frac{x}{s^{}} \frac{s^{}}{A z} & \frac{y}{A z} = \frac{y}{s^{}} \frac{s^{}}{A z} \end{matrix} \begin{matrix} \frac{x}{A z} = \frac{x_{l}}{A z} - \frac{x}{A z} & \frac{y}{A z} = \frac{y_{l}}{A z} - \frac{y}{A z} \end{matrix} & (4.16) \\ \frac{s^{}}{E l} = - (- \sin (El) \sin (Az) (X - X_{0}) + \sin (El) \cos (Az) (Y - Y_{0}) - \cos (El)) \begin{matrix} \frac{x}{El} = \frac{x}{s^{}} \frac{s^{}}{El} & \frac{y}{El} = \frac{y}{s^{}} \frac{s^{}}{El} \end{matrix} \begin{matrix} \frac{x}{El} = \frac{x_{l}}{El} - \frac{x}{El} & \frac{y}{El} = \frac{y_{l}}{El} - \frac{y}{El} \end{matrix} & (4.17) \\ \begin{matrix} \frac{s^{}}{Ro} = 0 & \frac{x}{Ro} = 0 & \frac{y}{Ro} = 0 \end{matrix} \begin{matrix} \frac{x}{Ro} = \frac{x_{l}}{Ro} & \frac{y}{Ro} = \frac{y_{l}}{Ro} \end{matrix} & (4.18) \end{matrix}$

(104) The present technology uses one-sided self-calibration bundle adjustment, and interior parameters participating in the adjustment are x.sub.p, y.sub.p, f, P.sub.1, P.sub.2, B.sub.1, B.sub.2, and one-sided radial distortion parameters k.sub.1ss.sub.1, k.sub.2ss.sub.2, and k.sub.3ss.sub.3 of S.sub.2.

(105) $\begin{matrix} \frac{x}{x_{p}} = - (K_{1 {ss}^{}} r^{2} + K_{2 {ss}^{}} r^{4} + K_{3 {ss}^{}} r^{6}) + \overline{x} (- 2 \overline{x}) (K_{1 {ss}^{}} + 2 K_{2 {ss}^{}} r^{2} + 3 K_{3 {ss}^{}} r^{4}) + P_{1} (- 6 \overline{x}) + 2 P_{2} (- \overline{y}) - B_{1} \frac{y}{x_{p}} = \overline{y} (- 2 \overline{x}) (K_{1 {ss}^{}} + 2 K_{2 {ss}^{}} r^{2} + 3 K_{3 {ss}^{}} r^{4}) + P_{2} (- 2 \overline{x}) + 2 P_{1} (- \overline{y}) \frac{x}{y_{p}} = \overline{x} (- 2 \overline{y}) (K_{1 {ss}^{}} + 2 K_{2 {ss}^{}} r^{2} + 3 K_{3 {ss}^{}} r^{4}) + P_{1} (- 2 \overline{y}) + 2 P_{2} (- \overline{x}) - B_{2} \frac{y}{y_{p}} = - (K_{1 {ss}^{}} r^{2} + K_{2 {ss}^{}} r^{4} + K_{3 {ss}^{}} r^{6}) + \overline{y} (- 2 \overline{y}) (K_{1 {ss}^{}} + 2 K_{2 {ss}^{}} r^{2} + 3 K_{3 {ss}^{}} r^{4}) + P_{2} (- 6 \overline{y}) + 2 P_{1} (- \overline{x}) \begin{matrix} \frac{x}{x_{p}} = 1 - \frac{x}{x_{p}} & \frac{y}{x_{p}} = - \frac{y}{x_{p}} \end{matrix} & (4.19) \\ \begin{matrix} \frac{x}{y_{p}} = - \frac{x}{y_{p}} & \frac{y}{y_{p}} = 1 - \frac{y}{y_{p}} \end{matrix} & (4.20) \\ \begin{matrix} \frac{x}{f} = - \frac{\overline{X}}{\overline{Z}} & \frac{v}{f} = - \frac{\overline{Y}}{\overline{Z}} \end{matrix} & (4.21) \\ \begin{matrix} \frac{K_{1 {ss}^{}}}{K_{1 {ss}_{2}}} = (1 -_{s^{}}) \frac{C_{s^{}}^{2}}{C_{s_{2}}^{2}} & \frac{K_{2 {ss}^{}}}{K_{2 {ss}_{2}}} = (1 -_{s^{}}) \frac{C_{s^{}}^{4}}{C_{s_{2}}^{4}} \end{matrix} \begin{matrix} \frac{K_{3 {ss}^{}}}{K_{3 {ss}_{2}}} = (1 -_{s^{}}) \frac{C_{s^{}}^{6}}{C_{s_{2}}^{6}} & \frac{x}{K_{1 {ss}_{2}}} = - \overline{x} r^{2} \frac{K_{1 {ss}^{}}}{K_{1 {ss}_{2}}} \end{matrix} \begin{matrix} \frac{x}{K_{2 {ss}_{2}}} = - \overline{x} r^{4} \frac{K_{2 {ss}^{}}}{K_{2 {ss}_{2}}} & \frac{x}{K_{3 {ss}_{2}}} = - \overline{x} r^{6} \frac{K_{3 {ss}^{}}}{K_{3 {ss}_{2}}} \end{matrix} & (4.22) \\ \begin{matrix} \frac{y}{K_{1 {ss}_{2}}} = - \overline{y} r^{2} \frac{K_{1 {ss}^{}}}{K_{1 {ss}_{2}}} & \frac{y}{K_{2 {ss}_{2}}} = - \overline{y} r^{4} \frac{K_{2 {ss}^{}}}{K_{2 {ss}_{2}}} \end{matrix} \frac{y}{K_{3 {ss}_{2}}} = - \overline{y} r^{6} \frac{K_{3 {ss}^{}}}{K_{3 {ss}_{2}}} & (4.23) \\ \begin{matrix} \frac{x}{P_{1}} = - (2 {\overline{x}}^{2} + r^{2}) & \frac{y}{P_{1}} = - 2 \overline{x y} \end{matrix} & (4.24) \\ \begin{matrix} \frac{y}{P_{2}} = - 2 \overline{x y} & \frac{y}{P_{2}} = - (2 {\overline{y}}^{2} + r^{2}) \end{matrix} & (4.25) \\ \begin{matrix} \frac{x}{B_{1}} = - \overline{x} & \frac{y}{B_{1}} = 0 \end{matrix} & (4.26) \\ \begin{matrix} \frac{x}{B_{2}} = - \overline{y} & \frac{y}{B_{2}} = 0 \end{matrix} & (4.27) \end{matrix}$

(106) Each term in B.sub.ij is solved through the following process:

(107) $\begin{matrix} \begin{matrix} \frac{x}{X} = - \frac{x}{X_{o}} & \begin{matrix} \frac{x}{Y} = - \frac{x}{Y_{o}} & \frac{x}{Z} = - \frac{x}{Z_{o}} \end{matrix} \end{matrix} & (4.28) \\ \begin{matrix} \begin{matrix} \frac{y}{X} = - \frac{y}{X_{o}} & \frac{y}{Y} = - \frac{y}{Y_{o}} & \frac{z}{Z} = - \frac{z}{Z_{o}} \end{matrix} \end{matrix} & (4.29) \end{matrix}$

(108) d3) performing a partitioning operation on the normal equation;

(109) Assuming that a dual-camera photogrammetry system performs photographing and measurement on n scale bar positions, the error equation set is:

(110) 0 $\begin{matrix} [\begin{matrix} v_{1, 11} \\ v_{1, 12} \\ v_{1, 21} \\ v_{1, 22} \\ .Math. \\ v_{1, i 1} \\ v_{1, i 2} \\ .Math. \\ v_{1, n 1} \\ v_{1, n 2} \\ v_{2, 11} \\ v_{2, 12} \\ v_{2, 21} \\ v_{2, 22} \\ .Math. \\ v_{2, j 1} \\ v_{2, j 2} \\ .Math. \\ v_{2, n 1} \\ v_{2, n 2} \end{matrix}] + [\begin{matrix} l_{1, 11} \\ l_{1, 12} \\ l_{1, 21} \\ l_{1, 22} \\ .Math. \\ l_{1, i 1} \\ l_{1, i 2} \\ .Math. \\ l_{1, n 1} \\ l_{1, n 2} \\ l_{2, 11} \\ l_{2, 12} \\ l_{2, 21} \\ l_{2, 22} \\ .Math. \\ l_{2, j 1} \\ l_{2, j 2} \\ .Math. \\ l_{2, n 1} \\ l_{2, n 2} \end{matrix}] = [\begin{matrix} A_{1, 11} & 0 & B_{1, 11} & 0 & 0 & 0 & .Math. & 0 & 0 & .Math. & 0 & 0 \\ A_{1, 12} & 0 & 0 & B_{1, 12} & 0 & 0 & .Math. & 0 & 0 & .Math. & 0 & 0 \\ A_{1, 21} & 0 & 0 & 0 & B_{1, 21} & 0 & .Math. & 0 & 0 & .Math. & 0 & 0 \\ A_{1, 22} & 0 & 0 & 0 & 0 & B_{1, 22} & .Math. & 0 & 0 & .Math. & 0 & 0 \\ .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. \\ A_{1, j 1} & 0 & 0 & 0 & 0 & 0 & .Math. & B_{1, j 1} & 0 & .Math. & 0 & 0 \\ A_{1, j 2} & 0 & 0 & 0 & 0 & 0 & .Math. & 0 & B_{1, j 2} & .Math. & 0 & 0 \\ .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. \\ A_{1, n 1} & 0 & 0 & 0 & 0 & 0 & .Math. & 0 & 0 & .Math. & B_{1, n 1} & 0 \\ A_{1, n 2} & 0 & 0 & 0 & 0 & 0 & .Math. & 0 & 0 & .Math. & 0 & B_{1, n 2} \\ 0 & A_{2, 11} & B_{2, 11} & 0 & 0 & 0 & .Math. & 0 & 0 & .Math. & 0 & 0 \\ 0 & A_{2, 12} & 0 & B_{2, 12} & 0 & 0 & .Math. & 0 & 0 & .Math. & 0 & 0 \\ 0 & A_{2, 21} & 0 & 0 & B_{2, 21} & 0 & .Math. & 0 & 0 & .Math. & 0 & 0 \\ 0 & A_{2, 22} & 0 & 0 & 0 & B_{2, 22} & .Math. & 0 & 0 & .Math. & 0 & 0 \\ .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. \\ 0 & A_{2, j 1} & 0 & 0 & 0 & 0 & .Math. & B_{2, j 1} & 0 & .Math. & 0 & 0 \\ 0 & A_{2, j 2} & 0 & 0 & 0 & 0 & .Math. & 0 & B_{2, j 2} & .Math. & 0 & 0 \\ .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. \\ 0 & A_{2, n 1} & 0 & 0 & 0 & 0 & .Math. & 0 & 0 & .Math. & B_{2, n 1} & 0 \\ 0 & A_{2, n 2} & 0 & 0 & 0 & 0 & .Math. & 0 & 0 & .Math. & 0 & B_{2, n 2} \end{matrix}] [\begin{matrix} _{1} \\ _{2} \\ {\dot{}}_{11} \\ {\dot{}}_{12} \\ {\dot{}}_{21} \\ {\dot{}}_{22} \\ .Math. \\ {\dot{}}_{j 1} \\ {\dot{}}_{j 2} \\ .Math. \\ {\dot{}}_{n 1} \\ {\dot{}}_{n 2} \end{matrix}] .Math. v + l = A \overset{.fwdarw.}{} & (4.30) \end{matrix}$

(111) where the symbol subscript i,jk represents the k.sup.th endpoint on (i=1,2, j=1,2, . . . n, k=1,2) the j.sup.th reference length for the i.sup.th camera.

(112) Assuming v.sup.Tv.fwdarw.min, the to-be-solved parameter increment is solved by using the following formula:
(A.sup.TA){right arrow over ()}=A.sup.Tl custom character N{right arrow over ()}=W(4.31)

(113) In view of the sparsity of A, a regular partitioned description may be made for A.sup.T A and A.sup.T l according to indexes of the image and the target point:

(114) $\begin{matrix} A^{T} A = N = [\begin{matrix} \underset{(32 32)}{N_{11}} & \underset{(32 6 n)}{N_{12}} \\ \underset{(6 n 32)}{N_{21}} & \underset{(6 n 6 n)}{N_{22}} \end{matrix}] & A^{T} l = W = [\begin{matrix} \underset{(32 1)}{w_{1}} \\ \underset{(6 n 1)}{w_{2}} \end{matrix}] \end{matrix}$ $N_{11} = [\begin{matrix} {.Math.}_{j = 1}^{n} {.Math.}_{k = 1}^{2} A_{1, jk}^{T} A_{1, jk} & 0 \\ 0 & {.Math.}_{j = 1}^{n} {.Math.}_{k = 1}^{2} A_{2, jk}^{T} A_{2, jk} \end{matrix}]$ $N_{22} = [\begin{matrix} {.Math.}_{i = 1}^{2} B_{i, 11}^{T} B_{i, 11} & 0 & 0 & 0 & .Math. & 0 & 0 & .Math. & 0 & 0 \\ 0 & {.Math.}_{i = 1}^{2} B_{i, 12}^{T} B_{i, 12} & 0 & 0 & .Math. & 0 & 0 & .Math. & 0 & 0 \\ 0 & 0 & {.Math.}_{i = 1}^{2} B_{i, 21}^{T} B_{i, 21} & 0 & .Math. & 0 & 0 & .Math. & 0 & 0 \\ 0 & 0 & 0 & {.Math.}_{i = 1}^{2} B_{i, 22}^{T} B_{i, 22} & .Math. & 0 & 0 & .Math. & 0 & 0 \\ .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. \\ 0 & 0 & 0 & 0 & .Math. & {.Math.}_{i = 1}^{2} B_{i, j 1}^{T} B_{i, j 1} & 0 & .Math. & 0 & 0 \\ 0 & 0 & 0 & 0 & .Math. & 0 & {.Math.}_{i = 1}^{2} B_{i, j 2}^{T} B_{i, j 2} & .Math. & 0 & 0 \\ .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. \\ 0 & 0 & 0 & 0 & .Math. & 0 & 0 & .Math. & {.Math.}_{i = 1}^{2} B_{i, n 1}^{T} B_{i, n 1} & 0 \\ 0 & 0 & 0 & 0 & .Math. & 0 & 0 & .Math. & 0 & {.Math.}_{i = 1}^{2} B_{i, n 2}^{T} B_{i, n 2} \end{matrix}]$ $N_{12} = [\begin{matrix} A_{1, 11}^{T} B_{1, 11} & A_{1, 12}^{T} B_{1, 12} & A_{1, 21}^{T} B_{1, 21} & A_{1, 22}^{T} B_{1, 22} & .Math. & A_{1, j 1}^{T} B_{1, j 1} & A_{1, j 2}^{T} B_{1, j 2} & .Math. & A_{1, n 1}^{T} B_{1, n 1} & A_{1, n 2}^{T} B_{1, n 2} \\ A_{2, 11}^{T} B_{2, 11} & A_{2, 12}^{T} B_{2, 12} & A_{2, 21}^{T} B_{2, 21} & A_{2, 22}^{T} B_{2, 22} & .Math. & A_{2, j 1}^{T} B_{2, j 1} & A_{2, j 2}^{T} B_{2, j 2} & .Math. & A_{2, n 1}^{T} B_{2, n 1} & A_{2, n 2}^{T} B_{2, n 2} \end{matrix}]$ $N_{21} = N_{12}^{T}$ $w_{1} = [\begin{matrix} {.Math.}_{j = 1}^{n} {.Math.}_{k = 1}^{2} A_{1, jk}^{T} l_{1, jk} \\ {.Math.}_{j = 1}^{n} {.Math.}_{k = 1}^{2} A_{2, jk}^{T} l_{2, jk} \end{matrix}]$ $w_{2} = [\begin{matrix} {.Math.}_{i = 1}^{2} B_{i, 11}^{T} l_{i, 11} \\ {.Math.}_{i = 1}^{2} B_{i, 12}^{T} l_{i, 12} \\ {.Math.}_{i = 1}^{2} B_{i, 21}^{T} l_{i, 21} \\ {.Math.}_{i = 1}^{2} B_{i, 22}^{T} l_{i, 22} \\ {.Math.}_{i = 1}^{2} B_{i, j 1}^{T} l_{i, j 1} \\ {.Math.}_{i = 1}^{2} B_{i, j 2}^{T} l_{i, j 2} \\ {.Math.}_{i = 1}^{2} B_{i, n 1}^{T} l_{i, n 1} \\ {.Math.}_{i = 1}^{2} B_{i, n 2}^{T} l_{i, n 2} \end{matrix}]$

(115) d4) obtaining self-calibration bundle adjustment with length constraints;

(116) Constraint conditions are described by using the following formula:

(117) $\begin{matrix} \underset{(n (32 + 6 n))}{C} \underset{(32 + 6 n) 1}{\overset{.fwdarw.}{}} + \underset{(n 1)}{W_{x}} = 0 & (4.32) \\ C = [\begin{matrix} 0 & 0 & \frac{L_{1}}{c_{11}} & \frac{L_{1}}{c_{12}} & 0 & 0 & .Math. & 0 & 0 & .Math. & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{L_{2}}{c_{21}} & \frac{L_{2}}{c_{22}} & .Math. & 0 & 0 & .Math. & 0 & 0 \\ .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. \\ 0 & 0 & 0 & 0 & 0 & 0 & .Math. & \frac{L_{j}}{c_{j 1}} & \frac{L_{j}}{c_{j 2}} & .Math. & 0 & 0 \\ .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. & .Math. \\ 0 & 0 & 0 & 0 & 0 & 0 & .Math. & 0 & 0 & .Math. & \frac{L_{n}}{c_{n 1}} & \frac{L_{n}}{c_{n 2}} \end{matrix}] \\ W_{x} = [\begin{matrix} L - L_{1} \\ L - L_{2} \\ .Math. \\ L - L_{j} \\ .Math. \\ L - L_{n} \end{matrix}] \end{matrix}$

(118) where L is a measured length value, L.sub.j is the length of the scale bar at the j.sup.th position reconstructed by the photogrammetry system, and a method for solving a partial derivative of the length relative to spatial coordinates of the corresponding endpoint is:

(119) $\frac{L_{j}}{c_{j 1}} = [\begin{matrix} \frac{x_{j 1} - x_{j 2}}{L_{j}} & \frac{y_{j 1} - y_{j 2}}{L_{j}} & \frac{z_{j 1} - z_{j 2}}{L_{j}} \end{matrix}]$ $\frac{L_{j}}{c_{j 2}} = [\begin{matrix} - \frac{x_{j 1} - x_{j 2}}{L_{j}} & - \frac{y_{j 1} - y_{j 2}}{L_{j}} & - \frac{z_{j 1} - z_{j 2}}{L_{j}} \end{matrix}]$

(120) Then, a least-squares solution and corresponding contact vector of the formula (4.31) under the constraint conditions described by the formula (4.32) are:

(121) $\begin{matrix} \overset{.fwdarw.}{} = (N^{- 1} - N^{- 1} C^{T} N_{cc}^{- 1} {CN}^{- 1}) W - N^{- 1} C^{T} N_{cc}^{- 1} Wx N_{cc} = {CN}^{- 1} C^{T} K_{s} = N_{cc}^{- 1} ({CN}^{- 1} W + W_{x}) & (4.33) \end{matrix}$

(122) d5) evaluating precision of the photogrammetry system;

(123) The mean error of weight unit for the photogrammetry system is:

(124) $\begin{matrix} {\hat{s}}_{0} = \sqrt{\frac{V^{T} V}{3 n - 32}} & (4.34) \end{matrix}$

(125) To ensure the orientation precision, it is recommended that scale bars be placed at 60 or more positions.

(126) An error covariance matrix of interior and exterior orientation parameters of the two cameras and spatial coordinates of feature points of the scale bar is:

(127) $\begin{matrix} \underset{((32 + 6 n) (32 + 6 n))}{C} = {\hat{s}}_{0}^{2} (\underset{((32 + 6 n) (32 + 6 n))}{N^{- 1} - N^{- 1} C^{T} N_{cc}^{- 1} {CN}^{- 1}}) & (4.35) \end{matrix}$

(128) d6) performing adaptive proportional adjustment of error equation terms;

(129) Unknown variables to be solved vary considerably in order of magnitude, especially for interior parameters. For example, generally terms K1, P1, P2, B1 and B2 have 10-5, term K2 has 10-7, and term K3 has 10-11. Such a considerable difference in order of magnitude leads to the ill-conditioned problem of the matrix A in the error equation (4.31). Due to the common problems in numerical computation such as machine precision, loss of low-order digits, and division of a large number by a small number, the computational result of the normal equation may be greatly affected. This often leads to problems such as matrix rank deficiency. To unify the orders of magnitude of all the terms in the error equation, an adaptive proportional adjustment technique is designed herein.

(130) Each time the adjustment iteration process begins, the order of magnitude Mj and an adjustment coefficient kj of each column in the error equation are calculated first, using the following method:
M.sub.j=average(round(ln(A(:,j)))) custom character k.sub.j=10.sup.M.sup.j(4.36)

(131) Then, the computational result of each column is multiplied by a corresponding proportionality coefficient, to obtain an adjusted error equation:

(132) $A^{} = A [\begin{matrix} k_{1} & 0 & .Math. & 0 \\ 0 & k_{2} & .Math. & 0 \\ .Math. & .Math. & .Math. & .Math. \\ 0 & 0 & .Math. & k_{u} \end{matrix}]$

(133) A relationship between an unknown variable {right arrow over ()} and {right arrow over ()} to be solved is calculated by using the formula (4.31):

(134) $\begin{matrix} \overset{.fwdarw.}{} = [\begin{matrix} k_{1} & 0 & .Math. & 0 \\ 0 & k_{2} & .Math. & 0 \\ .Math. & .Math. & .Math. & .Math. \\ 0 & 0 & .Math. & k_{u} \end{matrix}] \overset{.fwdarw.}{}, & (4.37) \end{matrix}$

(135) As proved by both simulation and experiments, such adaptive proportional adjustment eliminates the ill-conditioned problem of the normal equation in numerical computation, and provides higher adjustment stability and adjustment precision.

(136) Step e

(137) After the orientation process is completed, interior and exterior orientation parameters of the two cameras are fixed, and then spatial coordinates of points to be measured are solved by using a least-squares technique. In the process of solving the spatial coordinates, because there is no correlation between the measured points, the process may be performed point by point, thereby avoiding the time consumption caused by large-scale matrix operations and improving the performance of dynamic measurement. This is in essence a direct linear transformation method.

(138) According to an embodiment of the present invention, a specific method for performing traceable evaluation of orientation precision of the photogrammetry system is:

(139) e1) describing an error equation of a single imaging point;

(140) $\begin{matrix} {\hat{l}}_{ij} = l_{ij} + v_{ij} = K_{ij} {\hat{x}}_{ij} .Math. [\begin{matrix} x_{ij} -_{x} (X_{ij}^{0}) \\ y_{ij} -_{y} (X_{ij}^{0}) \end{matrix}] + [\begin{matrix} v_{xij} \\ v_{yij} \end{matrix}] = B_{ij} [\begin{matrix} X_{j} \\ Y_{j} \\ Z_{j} \end{matrix}] = B_{ij} {\dot{}}_{j} i = 1, 2; j = 1, 2, .Math. n & (5.1) \end{matrix}$

(141) The sparse matrix B.sub.ij is a partial derivative matrix of image plane coordinates relative to spatial coordinates, and is solved by using the formulas (4.30) and (4.31).

(142) e2) when the photogrammetry system includes two cameras, an error equation set corresponding to a spatial point [X Y Z] is:

(143) 0 $\begin{matrix} \underset{(4 1)}{[\begin{matrix} l_{1} \\ l_{2} \end{matrix}]} + \underset{(4 1)}{[\begin{matrix} v_{1} \\ v_{2} \end{matrix}]} = \underset{(4 3)}{[\begin{matrix} B_{1} \\ B_{2} \end{matrix}]} [\begin{matrix} X \\ Y \\ Z \end{matrix}] .Math. l + v = B \dot{} & (5.2) \end{matrix}$

(144) e3) a corresponding normal equation and coordinate correction amount are:
(B.sup.TB){dot over ()}=B.sup.Tl {dot over ()}=(B.sup.TB).sup.1B.sup.Tl(5.3)

(145) e4) the orientation precision evaluation is evaluation of the length measurement, and is performed by using data in the orientation process. Coordinates of endpoints of 2n scale bars are reconstructed by using the formula (5.3), and then a measured value of the length of the corresponding scale bar is solved:

(146) L.sub.1 L.sub.2 . . . L.sub.n

(147) e5) an average value of measurement results is:

(148) $\begin{matrix} \overline{L} = \frac{{.Math.}_{j = 1}^{n} L_{j}}{n} & (5.4) \end{matrix}$

(149) e6) performing scaling by using a ratio of an average length to a gage length as an orientation result (translation amount of exterior orientation parameter):

(150) $\begin{matrix} \begin{matrix} s = \frac{L}{\overline{L}} & T_{1}^{} = {sT}_{1} & T_{2}^{} = {sT}_{2} \end{matrix} & (5.5) \\ \begin{matrix} L_{1}^{} = {sL}_{1} & L_{2}^{} = {sL}_{2} & .Math. & L_{n}^{} = {sL}_{n} \end{matrix} & (5.6) \end{matrix}$

(151) and performing an uncertainty analysis on the length measurement result, to obtain an orientation precision evaluation:

(152) $\begin{matrix} u (L_{i}^{}) = \sqrt{\frac{{.Math.}_{i = 1}^{n} {(L_{i}^{} - L)}^{2}}{n - 1}} = \sqrt{\frac{{.Math.}_{i = 1}^{n} L_{i}^{2}}{n - 1}} = su (L_{i}) & (5.7) \end{matrix}$

(153) e7) a quality parameter of a measurement instrument is described as a multiple of the uncertainty of the length:
B=ku(L.sub.i)(5.8)

(154) e8) a relative error of each length is:

(155) $\begin{matrix} r_{i}^{} = \frac{L_{i}^{} - L}{L} = \frac{L_{i} - \overline{L}}{\overline{L}} = r_{i} & (5.9) \end{matrix}$

(156) e9) the relative error in length measurement is equal to the relative error in length measurement before scaling is performed, and therefore, an uncertainty of the relative error in length measurement is also equal to an uncertainty of the original relative error in length measurement:

(157) $\begin{matrix} u (r_{i}^{}) = u (r_{i}) = \frac{u (L_{i})}{\overline{L}} & (5.10) \end{matrix}$

(158) According to the method for implementing high-precision orientation and evaluating orientation precision of a large-scale dynamic photogrammetry system of the present invention, the following verification experiment was carried out:

(159) A photogrammetry system used in the experiment includes two industrial cameras. As shown in FIG. 3, a scale bar having a length of about 950 mm is used for orientation, and a reference length is indicated by backlight reflection code points at two ends of the scale bar. The volume of the measurement environment is about 3.6 m2.4 m2.4 m. The scale bar is moved in the measurement environment and imaged, self-calibration positioning and orientation of the photogrammetry system is implemented by using the method of the present invention, and the orientation precision of the system is evaluated by using the evaluation method according to claim 11. Results of three orientation experiments are shown in Table 1 and Table 2. Table 1 shows an orientation result of self-calibration bundle adjustment aimed at the smallest image plane error. Table 2 shows an orientation result of self-calibration bundle adjustment with multiple length constraints according to the present technology.

(160) TABLE-US-00001 TABLE 1 orientation result of conventional self-calibration bundle adjustment length measurement error standard deviation of stand- relative image plane error ard error direction direction devi- uncer- x y ation tainty orientation camera 2.85E05 4.24E04 3.800 1/250 experiment 1 1 camera 2.32E05 3.43E04 2 orientation the 1.60E05 4.72E04 6.271 1/150 experiment camera 2 1 camera 1.53E05 4.92E04 2 orientation camera 1.71E05 3.02E04 3.474 1/270 experiment 1 3 camera 1.43E05 2.88E04 2

(161) TABLE-US-00002 TABLE 2 orientation result of self-calibration bundle adjustment with multiple length constraints standard deviation of image plane error length direction measurement error x direction x orientation camera 1.80E04 4.60E04 0.060 1/15800 experiment 1 1 camera 1.85E04 4.18E04 2 orientation the 2.15E04 5.43E04 0.064 1/14800 experiment 2 camera 1 camera 2.28E04 5.32E04 2 orientation camera 1.47E04 3.68E04 0.051 1/18700 experiment 3 1 camera 1.89E04 3.34E04 2

(162) It can be seen from comparison of standard deviations of image plane errors of the three orientation experiments in Table 1 and Table 2 that in the conventional self-calibration method, the error in the direction x of the image plane is reduced to a ignorable degree, and is smaller than the image plane error in the other direction by one order of magnitude. This is obviously incorrect, because error levels in two directions of imaging should approximate. Such an adjustment method leads to an incorrect self-calibration result and incorrect positioning and orientation result, and the corresponding length error is large. For the self-calibration adjustment technique with multiple length constraints according to the present invention, due to the constraint of the space length, the image plane error is not the only factor that determines the adjustment result, the image plane error is more reasonable, and the corresponding length error is significantly reduced.

(163) The experimental results indicate that the use of the relative error in length measurement as an indicator of evaluation can ignore the analytical complexity caused by the reference length measurement error. If the length measurement is traceable, the precision of the photogrammetry system can also be evaluated according to an absolute measurement error, and the result of evaluation is also traceable.

(164) Based on the above, in the present invention, on-site positioning and orientation of a dual-camera large-scale dynamic photogrammetry system is completed by using a single scale bar. The measurement space is divided into regions so that the distance between the left and right cameras is equal to the region length of the divided space, and the gauge attitude and the camera position are adjusted at the same time during the photographing process, making the length orientation field more reliable. Preliminary orientation is performed by combining a five-point method with the space length ratio, so that the solution loss problem in special cases can be overcome, and a feasible solution can be accurately determined in multiple solutions. In addition, the correlation between parameters can be overcome by using self-calibration bundle adjustment with multiple space length constraints, thereby improving the overall orientation precision of the photogrammetry system. The statistical result of the length measurement error provides a traceable objective evaluation indicator. The experimental result shows that the relative error in length measurement of the present technology reaches 1/15000.

(165) The above descriptions are merely preferred embodiments of the present invention, but are not intended to limit the scope of implementation of the present invention. Hence, any equivalent variations or modifications made in accordance with the structures, features and principles described in the claims of the present invention shall fall within the scope of the claims of the present invention.

Method for implementing high-precision orientation and evaluating orientation precision of large-scale dynamic photogrammetry system

Assignee

Inventors

Cpc classification

Classification Explorer

G06T2200/21

PHYSICS

Classification Explorer

G01C11/14

PHYSICS

Classification Explorer

G06T2207/30244

PHYSICS

Classification Explorer

G06T7/77

PHYSICS

Classification Explorer

G06T3/40

PHYSICS

Classification Explorer

G06F17/18

PHYSICS

Classification Explorer

G01C11/10

PHYSICS

Classification Explorer

G06T7/85

PHYSICS

Classification Explorer

G06T7/75

PHYSICS

International classification

Classification Explorer

G06K9/36

PHYSICS

Classification Explorer

G06T7/77

PHYSICS

Classification Explorer

G06T3/40

PHYSICS

Classification Explorer

G01C11/14

PHYSICS

Classification Explorer

G01C11/10

PHYSICS

Classification Explorer

G06T7/80

PHYSICS

Classification Explorer

G06T7/73

PHYSICS

Abstract

Claims

Description