CN103217145B - A kind of Mars DEM makes and strip aerial triangulation method - Google Patents

Abstract

The present invention proposes a kind of Mars DEM to make and strip aerial triangulation method, described method is by data source based on Mars probes laser-measured height data, conversion between Mars geographic coordinate and Mars tangent plane coordinate is carried out to it, then interpolation generates high-precision Mars DEM, solves the blank of current China Mars DEM; Dual quaterion type is utilized to carry out strip aerial triangulation, by the attitude of time series analog image and the relative position unified process of dual quaterion, there is complete geometry implication, and interative computation process does not rely on initial value, the element of orientation and the closed points coordinate of high precision Mars sequential images can be obtained.

Description

A kind of Mars DEM makes and strip aerial triangulation method

Technical field

The invention belongs to deep space technical field of mapping, especially relate to a kind of Mars DEM and make and strip aerial triangulation method.

Background technology

Mars is first planet be positioned at outside Earth's orbit, has a kind of special attractive force for the mankind, and it is one of celestial body of the closest earth in the solar system.The angle of cut of mars equatorial plane and revolution orbit plane is in close proximity to the earth, and this makes it also have the four seasons of the similar earth alternately, is most interested, the detected event planet the most frequently of the mankind beyond the earth.Mars exploration is one of 5 priority research areass of 21 century survey of deep space, the history that mankind's usage space detector carries out mars exploration almost runs through whole mankind's space flight history, the U.S., Europe, Japan and Russia have had and had arrived Mars more than 30 detectors, and have beamed back a large amount of scientific exploration data comprising martian surface high resolution image, air and geologic data to the earth.Carry Russia's " Forbes-soil " detector of China's first Mars probes " No. one, the light of firefly " in lift-off smoothly on November 9, but this detector fails in flight course to realize according to plan becoming rail, detection mission failure.At present, China is actively preparing mars exploration project next time, strives the rocket and the TTC & DT Systems that within 2013, utilize oneself, launches independently, runs the Mars probes of China.Mapping martian surface landform is one of top priority of mars exploration, and it provides indispensable Back ground Information for subsequent probe Mars.Optical image is the general data source of current martian surface topography mapping, and the making of DEM and the hi-Fix of image remain key and the core of Mars optical image mapping.

Digital elevation model (DigitalElevationModel, DEM) be a kind of actual ground model representing ground elevation by one group of orderly array of values form, Mars DEM is the set of reaction Mars earth's surface elevation, is the indispensable Back ground Information such as Mars observation and prospecting Mars resource.The current U.S. has the most comprehensive Mars Controling network of covering, but United States Geological office not yet discloses its control point source, therefore must find suitable method and process available data, to obtaining high-precision Mars DEM.Before the Mars probes (such as " No. one, the light of firefly (YH-1) ") of China are launched, need to carry out a large amount of research work, according to prior art, data source required for research work, the such as acquisition of Mars DEM and analog image data is a great problems, due to the restriction by data source factor, the experimental work in later stage cannot be carried out, thus make the correctness of the principle of proposition can not get effective checking.

Aerotriangulation can be stereometric map, make image plane figure and orthophotoquad provides oriented control point and the element of orientation (determining the parameter of light shafts in photography moment time space position of photographing), is one of main path studying Mars landform further.Traditional aerotriangulation adopts rotation and the translation of Eulerian angle denotation coordination, and it is lower that precision often encrypted by process high spud angle image, and strong to the dependence of initial value.

Summary of the invention

Technical matters to be solved by this invention is to overcome the deficiencies in the prior art, proposes a kind of Mars DEM and makes and strip aerial triangulation method.

For solving the problems of the technologies described above, the technical solution adopted in the present invention is:

A kind of Mars DEM makes and strip aerial triangulation method, comprises the steps:

Steps A, carries out bilinear interpolation to laser-measured height data, makes Mars DEM;

Step B, the subpoint position of compute ray bundle on Mars DEM, generates Mars analog image data;

Step C, adopts dual quaterion to carry out strip aerial triangulation, obtains Mars image direction element and closed points coordinate.

Described steps A carries out bilinear interpolation to laser-measured height data, makes Mars DEM; Its process is as follows:

Steps A-1, the conversion between Mars geographic coordinate and Mars tangent plane coordinate is carried out to laser-measured height data:

$\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}{ccc}\sin B\cos L& \sin B\sin L& \cos B\\ -\sin L& \cos L& 0\\ \cos B\cos L& \cos B\sin L& \sin B\end{array}\right)\left(\begin{array}{c}(v+h)\cos L\cos B-X_0\\ (v+h)\cos L\sin B-Y_0\\ ((1-e^2)v+h)\sin L-Z_0\end{array}\right)$

In formula, (x, y, z) is Mars tangent plane coordinate, and (L, B, h) is the geographic coordinate of Mars ground point A, (X ₀, Y ₀, Z ₀) be the origin of tangent plane coordinate system, $v=a/\sqrt{(1-e^2{\sin }^{2}L)},e^2=(a^2-b^2)/a^2,$ A, b are respectively the length semiaxis of Mars ellipsoid;

Steps A-2, the cosine law is utilized to calculate the size of the triangle interior angle on alternative 3rd summit, the point selecting the maximum corresponding is this leg-of-mutton 3rd summit, set up triangular grid, on the triangle at grid points place, the value of net point is asked for, interpolation create-rule Grid DEM by distance weighting counting backward technique; Circular is as follows:

If triangular apex coordinate is (X _i, Y _i, Z _i), i=1,2,3, interpolation point coordinate is (X ₀, Y ₀, Z ₀), and grid points is in this triangle;

Then the height value of this triangulation network lattice point is:

$Z_0=\frac{\underset{i=1}{\overset{3}{\mathrm{\Σ}}}\frac{1}{{(X_i-X_0)}^2+{(Y_i-Y_0)}^2}{Z}_i}{\underset{i=1}{\overset{3}{\mathrm{\Σ}}}\frac{1}{{(X_i-X_0)}^2+{(Y_i-Y_0)}^2}}$

In described step C, adopt dual quaterion to carry out strip aerial triangulation, obtain Mars image direction element and closed points coordinate; Its process is as follows:

Step C-1, dual quaterion $\hat{q}={\left[\begin{array}{cccc}{q}_0& q_x& q_y& q_z\end{array}\right]}^T+\mathrm{\ϵ}{\left[\begin{array}{cccc}{r}_0& r_x& r_y& r_z\end{array}\right]}^T$ The Mars image imaging equation represented is:

$\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]=\left[\begin{array}{c}2q_0{r}_x-2q_x{r}_0+{2q}_y{r}_z-2q_z{r}_y\\ 2q_0{r}_y-2q_x{r}_z-2q_y{r}_0+2q_z{r}_x\\ 2q_0{r}_z+2q_x{r}_y-2q_y{r}_x-2q_z{r}_0\end{array}\right]$

$\left[\begin{array}{c}{X}_A\\ Y_A\\ Z_A\\ 1\end{array}\right]=\left[\begin{array}{cccc}{a}_1& a_2& a_3& X_S\\ b_1& b_2& b_3& Y_S\\ c_1& c_2& c_3& Z_S\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]$

In formula,

$\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]=\left[\begin{array}{ccc}{q}_0^2+q_x^2-q_y^2-q_z^2& 2(q_x{q}_y-q_0{q}_z)& 2(q_x{q}_z+q_0{q}_y)\\ 2(q_x{q}_y+q_0{q}_Z)& q_0^2-q_x^2+q_y^2-q_z^2& 2(q_y{q}_z-q_0{q}_x)\\ 2(q_x{q}_z-q_0{q}_y)& 2(q_z{q}_y+q_0{q}_x)& q_0^2-q_x^2-q_y^2+q_z^2\end{array}\right];$

Q ₀, q _x, q _y, q _z, r ₀, r _x, r _y, r _zfor 8 elements of dual quaterion;

X _s, Y _s, Z _sfor element of orientation line element;

(X _a, Y _a, Z _a) be the coordinate of ground point A under Mars rectangular coordinate system;

(X, Y, Z) is the coordinate of ground point A under image space coordinate system;

Step C-2, input time, sequence Mars analog image data, were taken as 0 by the initial value of dual quaterion;

Step C-3, sets up dual quaterion error equation:

$V=\left[\begin{array}{cc}A& B\end{array}\right]\left[\begin{array}{c}t\\ X\end{array}\right]-L$

In formula,

$A=\left[\begin{array}{cccccccc}\frac{{\∂F}_x}{{\∂q}_0}& \frac{{\∂F}_x}{{\∂q}_x}& \frac{{\∂F}_x}{{\∂q}_y}& \frac{{\∂F}_x}{{\∂q}_z}& \frac{{\∂F}_x}{{\∂r}_0}& \frac{{\∂F}_x}{{\∂r}_x}& \frac{{\∂F}_x}{{\∂r}_y}& \frac{{\∂F}_x}{{\∂r}_z}\\ \frac{{\∂F}_y}{{\∂q}_0}& \frac{{\∂F}_y}{{\∂q}_x}& \frac{{\∂F}_y}{{\∂q}_y}& \frac{{\∂F}_y}{{\∂q}_z}& \frac{{\∂F}_y}{{\∂r}_0}& \frac{{\∂F}_y}{{\∂r}_x}& \frac{{\∂F}_y}{{\∂r}_y}& \frac{{\∂F}_y}{{\∂r}_z}\end{array}\right];$

$B=\left[\begin{array}{ccc}\frac{{\∂F}_x}{{\∂X}_A}& \frac{{\∂F}_x}{{\∂Y}_A}& \frac{{\∂F}_x}{{\∂Z}_A}\\ \frac{{\∂F}_y}{{\∂X}_A}& \frac{{\∂F}_y}{{\∂Y}_A}& \frac{{\∂F}_y}{{\∂Z}_A}\end{array}\right];$

$t={\left[\begin{array}{cccccccc}d{q}_0& {\mathrm{dq}}_x& {\mathrm{dq}}_y& {\mathrm{dq}}_z& {\mathrm{dr}}_0& {\mathrm{dr}}_x& {\mathrm{dr}}_y& {\mathrm{dr}}_z\end{array}\right]}^T,$ D is for differentiating;

$X={\left[\begin{array}{ccc}{\mathrm{dX}}_A& {\mathrm{dY}}_A& {\mathrm{dZ}}_A\end{array}\right]}^T;$

x, y are respectively the coordinate of picture point corresponding to ground point A under photo coordinate system;

$V={\left[\begin{array}{cc}{v}_x& v_y\end{array}\right]}^T;$

$F_x=-f\frac{a_1(X_A-X_S)+b_1(Y_A-Y_S)+c_1(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)};$

$F_y=-f\frac{a_2(X_A-X_S)+b_2(Y_A-Y_S)+c_2(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)};$

V _x, v _ybe respectively F _x, F _ycorrection;

F is the focal length of detector camera;

(F _x) ₀, (F _y) ₀by (X _a, Y _a, Z _a) and the element of orientation value substitute into formula F _xand F _yobtain;

In air strips, the element of orientation of all images meets restrictive condition:

B _Xt+W=0

In formula,

$W=\left[\begin{array}{c}{W}_1\\ W_2\\ \·\·\·\\ W_{(n-1)}\\ W_n\end{array}\right];$

$B_{\mathrm{Xi}}=\left[\begin{array}{cccccccc}{q}_{0i}& q_{\mathrm{xi}}& q_{\mathrm{yi}}& q_{\mathrm{zi}}& r_{0i}& r_{\mathrm{xi}}& r_{\mathrm{yi}}& r_{\mathrm{zi}}\\ 0& 0& 0& 0& q_{0i}& q_{\mathrm{xi}}& q_{\mathrm{yi}}& q_{\mathrm{zi}}\end{array}\right];$

$w_{1i}=\frac{1}{2}({q_{0i}}^2+{q_{\mathrm{xi}}}^2+{q_{\mathrm{yi}}}^2+{q_{\mathrm{zi}}}^2-1),w_{2i}=q_{0i}{r}_{0i}+q_{\mathrm{xi}}{r}_{\mathrm{xi}}+q_{\mathrm{yi}}{r}_{\mathrm{yi}}+q_{\mathrm{zi}}{r}_{\mathrm{zi}};$

$W_i={\left[\begin{array}{cc}{w}_{1i}& w_{2i}\end{array}\right]}^T;$

I represents the sequence number of Mars image;

Step C-4, sets up normal equation:

$\left[\begin{array}{cc}{A}^{T}A& A^{T}B\\ B^{T}A& B^{T}B\end{array}\right]\left[\begin{array}{c}t\\ X\end{array}\right]=\left[\begin{array}{c}{A}^{T}L\\ B^{T}L\end{array}\right]$

Step C-5, calculate the correction of Mars image direction element and closed points coordinate, its computing formula is:

$\left[\begin{array}{c}t\\ X\\ r\end{array}\right]-N^{-1}{W}_Y$

$N=\left[\begin{array}{cc}{C}^{T}C& {B_X}^T\\ B_X& 0\end{array}\right]{,W}_Y=\left[\begin{array}{c}{C}^{T}L\\ W\end{array}\right]$

In formula, $C=\left[\begin{array}{cc}{A}^{T}A& A^{T}B\\ B^{T}A& B^{T}B\end{array}\right],$ R is the redundant computation vector for trim both members matrix dimension in computation process;

Step C-6, upgrades the value of the element of orientation and pass point, checks whether the value of the correction of elements of exterior orientation and pass point is less than the threshold value of setting; If be less than threshold value, so iteration terminates, otherwise returns step C-3 continuation iteration, until correction is less than the threshold value of setting.

The invention has the beneficial effects as follows: the present invention proposes a kind of Mars DEM and make and strip aerial triangulation method, described method is by data source based on Mars whole world prospector's laser-measured height data, conversion between Mars geographic coordinate and Mars tangent plane coordinate is carried out to it, then interpolation generates high-precision Mars DEM, solves the blank of current China Mars DEM; Dual quaterion type is utilized to carry out strip aerial triangulation, by the attitude of time series analog image and the relative position unified process of dual quaterion, there is complete geometry implication, and interative computation process does not rely on initial value, the element of orientation and the closed points coordinate of high precision Mars sequential images can be obtained.

Accompanying drawing explanation

Fig. 1 is that a kind of Mars DEM makes and strip aerial triangulation method flow diagram.

Fig. 2 is Mars DEM schematic diagram.

Fig. 3 is dual quaterion strip aerial triangulation schematic diagram.

Fig. 4 is error equation structural drawing.

Fig. 5 is normal equation structural drawing.

Embodiment

Below in conjunction with accompanying drawing, a kind of Mars DEM making propose the present invention and strip aerial triangulation method are described in detail:

As shown in Figure 1, a kind of Mars DEM of the present invention makes and strip aerial triangulation method flow diagram, and its implementation process is as follows:

The first step: make Mars DEM, method is as follows:

1. with Mars north latitude 11.2 °-11.5 °, east longitude 41.2 ° of-42.5 ° of regions are example, and laser-measured height data discrete data point has 148.First the conversion between Mars geographic coordinate and Mars tangent plane coordinate is carried out to MOLA data:

$\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}{ccc}\sin B\cos L& \sin B\sin L& \cos B\\ -\sin L& \cos L& 0\\ \cos B\cos L& \cos B\sin L& \sin B\end{array}\right)\left(\begin{array}{c}(v+h)\cos L\cos B-X_0\\ (v+h)\cos L\sin B-Y_0\\ ((1-e^2)v+h)\sin L-Z_0\end{array}\right)---\left(1\right)$

In formula, (x, y, z) is Mars tangent plane coordinate, and (L, B, h) is the geographic coordinate of Mars ground point A, (X ₀, Y ₀, Z ₀) be the origin of tangent plane coordinate system $v=a/\sqrt{(1-e^2{\sin }^{2}L)},e^2=(a^2-b^2)/a^2,$ A, b are respectively the length semiaxis of Mars ellipsoid;

2. utilize the cosine law to calculate the size of the triangle interior angle on alternative 3rd summit, the point selecting the maximum corresponding is this leg-of-mutton 3rd summit, sets up triangular grid.On the triangle at grid points place, ask for the value of net point by distance weighting counting backward technique, interpolation create-rule Grid DEM.If triangular apex coordinate is (X _i, Y _i, Z _i), (i=1,2,3), interpolation point coordinate is (X ₀, Y ₀, Z ₀), and grid points is in this triangle.Then the height value of this triangulation network lattice point is:

$Z_0=\frac{\underset{i=1}{\overset{3}{\mathrm{\Σ}}}\frac{1}{{(X_i-X_0)}^2+{(Y_i-Y_0)}^2}{Z}_i}{\underset{i=1}{\overset{3}{\mathrm{\Σ}}}\frac{1}{{(X_i-X_0)}^2+{(Y_i-Y_0)}^2}}---\left(2\right)$

Fig. 2 is the schematic diagram of the DEM that Experimental Area generates, for follow-up image simulation data lay the foundation.

Second step: the subpoint position of compute ray bundle on DEM, method is as follows:

1. for Olympic Competition moral number (Odyssey) Mars probes, orbit altitude 400km, orbit inclination 93 °, positive circuit orbit, STK software emulation is utilized to generate Odyssey detector 2003-11-01 day 7:05:00.00 to 2003-11-02 day 7:05:00.00 position data according to above-mentioned nominal track parameter, sampling interval is 1 second, and Mars radius gets 3396.19km, then according to experimental formula emulation Odyssey detector attitude data below:

${\mathrm{\θ}}_{\mathrm{phi}}=a_{\mathrm{phi}}\cos \left(\frac{2\mathrm{\πt}}{u_{\mathrm{phi}}}\right)+b_{\mathrm{phi}}\sin \left(\frac{2\mathrm{\πt}}{u_{\mathrm{phi}}}\right)$

${\mathrm{\θ}}_{\mathrm{omega}}=a_{\mathrm{omega}}\cos \left(\frac{2\mathrm{\πt}}{u_{\mathrm{omega}}}\right)+b_{\mathrm{omega}}\sin \left(\frac{2\mathrm{\πt}}{u_{\mathrm{omega}}}\right)---\left(3\right)$

${\mathrm{\θ}}_{\mathrm{kapa}}=a_{\mathrm{kapa}}\cos \left(\frac{2\mathrm{\πt}}{u_{\mathrm{kapa}}}\right)+b_{\mathrm{kapa}}\sin \left(\frac{2\mathrm{\πt}}{u_{\mathrm{kapa}}}\right)$

In formula, θ _phi, θ _omega, θ _kapafor three attitude angle of detector, t is sampling instant, the coefficient a in formula _phi, a _omega, a _kapabe taken as respectively: 0.1,0,0.1, b _phi, b _omega, b _kapabe taken as 240,320,220, u respectively _phi, u _omega, u _kapabe taken as respectively :-1.0,0.5 ,-1.0, T _phi, T _omega, T _kapabe taken as 140 respectively, 220,120.

2. the detector locus generated with Odyssey sensor camera parameter and simulation, attitude data is for calculating parameter, utilize collinearity equation that DEM reference mark is projected to image, obtain corresponding picpointed coordinate, the image number that simulation generates is 3, amounts to 30 reference mark.

3rd step: carry out dual quaterion strip aerial triangulation, method is as follows:

As shown in Figure 3, be dual quaterion strip aerial triangulation schematic diagram.Point A is the intersection point of three light on image 1, image 2 and image 3, with O-S ₁light shafts are example, when detector motion is to position S ₁time, sensor shooting ground point forms photographic light flux, now remote sensor coordinate system S ₁-x ₁y ₁z ₁the angular range element rotated with Mars rectangular coordinate system O-XYZ can be write as 2 dual angles form, the analogy of line element of orientation vector is become S ₁-x ₁y ₁z ₁vow with the displacement of O-XYZ formula is utilized by dual quaterion can by S ₁-x ₁y ₁z ₁motion transform of spinning is in O-XYZ.

1. by the element of orientation of the time series analog image described in second step, use dual quaterion Unify legislation, carry out strip aerial triangulation; Its detailed process is as follows:

Dual quaterion $\hat{q}={\left[\begin{array}{cccc}{q}_0& q_x& q_y& q_z\end{array}\right]}^T+\mathrm{\ϵ}{\left[\begin{array}{cccc}{r}_0& r_x& r_y& r_z\end{array}\right]}^T$ The Mars image imaging equation represented is:

$\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]=\left[\begin{array}{c}2q_0{r}_x-2q_x{r}_0+{2q}_y{r}_z-2q_z{r}_y\\ 2q_0{r}_y-2q_x{r}_z-2q_y{r}_0+2q_z{r}_x\\ 2q_0{r}_z+2q_x{r}_y-2q_y{r}_x-2q_z{r}_0\end{array}\right]$

$\left[\begin{array}{c}{X}_A\\ Y_A\\ Z_A\\ 1\end{array}\right]=\left[\begin{array}{cccc}{a}_1& a_2& a_3& X_S\\ b_1& b_2& b_3& Y_S\\ c_1& c_2& c_3& Z_S\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]$

In formula,

$\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]=\left[\begin{array}{ccc}{q}_0^2+q_x^2-q_y^2-q_z^2& 2(q_x{q}_y-q_0{q}_z)& 2(q_x{q}_z+q_0{q}_y)\\ 2(q_x{q}_y+q_0{q}_Z)& q_0^2-q_x^2+q_y^2-q_z^2& 2(q_y{q}_z-q_0{q}_x)\\ 2(q_x{q}_z-q_0{q}_y)& 2(q_z{q}_y+q_0{q}_x)& q_0^2-q_x^2-q_y^2+q_z^2\end{array}\right];$

Q ₀, q _x, q _y, q _z, r ₀, r _x, r _y, r _zfor 8 elements of dual quaterion;

X _s, Y _s, Z _sfor element of orientation line element;

(X _a, Y _a, Z _a) be the coordinate of ground point A under Mars rectangular coordinate system;

(X, Y, Z) is the coordinate of ground point A under image space coordinate system;

2. input basic data, comprise time series Mars analog image data, the initial value of dual quaterion is taken as 0, then starts iteration.

3. dual quaterion error equation is set up:

$V=\left[\begin{array}{cc}A& B\end{array}\right]\left[\begin{array}{c}t\\ X\end{array}\right]-L$

In formula,

$A=\left[\begin{array}{cccccccc}\frac{{\∂F}_x}{{\∂q}_0}& \frac{{\∂F}_x}{{\∂q}_x}& \frac{{\∂F}_x}{{\∂q}_y}& \frac{{\∂F}_x}{{\∂q}_z}& \frac{{\∂F}_x}{{\∂r}_0}& \frac{{\∂F}_x}{{\∂r}_x}& \frac{{\∂F}_x}{{\∂r}_y}& \frac{{\∂F}_x}{{\∂r}_z}\\ \frac{{\∂F}_y}{{\∂q}_0}& \frac{{\∂F}_y}{{\∂q}_x}& \frac{{\∂F}_y}{{\∂q}_y}& \frac{{\∂F}_y}{{\∂q}_z}& \frac{{\∂F}_y}{{\∂r}_0}& \frac{{\∂F}_y}{{\∂r}_x}& \frac{{\∂F}_y}{{\∂r}_y}& \frac{{\∂F}_y}{{\∂r}_z}\end{array}\right];$

$B=\left[\begin{array}{ccc}\frac{{\∂F}_x}{{\∂X}_A}& \frac{{\∂F}_x}{{\∂Y}_A}& \frac{{\∂F}_x}{{\∂Z}_A}\\ \frac{{\∂F}_y}{{\∂X}_A}& \frac{{\∂F}_y}{{\∂Y}_A}& \frac{{\∂F}_y}{{\∂Z}_A}\end{array}\right];$

$t={\left[\begin{array}{cccccccc}d{q}_0& {\mathrm{dq}}_x& {\mathrm{dq}}_y& {\mathrm{dq}}_z& {\mathrm{dr}}_0& {\mathrm{dr}}_x& {\mathrm{dr}}_y& {\mathrm{dr}}_z\end{array}\right]}^T,$ D is for differentiating;

$X={\left[\begin{array}{ccc}{\mathrm{dX}}_A& {\mathrm{dY}}_A& {\mathrm{dZ}}_A\end{array}\right]}^T;$

x, y are respectively the coordinate of picture point corresponding to ground point A under photo coordinate system;

$V={\left[\begin{array}{cc}{v}_x& v_y\end{array}\right]}^T;$

$F_x=-f\frac{a_1(X_A-X_S)+b_1(Y_A-Y_S)+c_1(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)};$

$F_y=-f\frac{a_2(X_A-X_S)+b_2(Y_A-Y_S)+c_2(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)};$

V _x, v _ybe respectively F _x, F _ycorrection;

F is the focal length of detector camera;

(F _x) ₀, (F _y) ₀by (X _a, Y _a, Z _a) and the element of orientation value substitute into formula F _xand F _yobtain;

In air strips, the element of orientation of all images meets restrictive condition:

B _Xt+W=0

In formula,

$W=\left[\begin{array}{c}{W}_1\\ W_2\\ \·\·\·\\ W_{(n-1)}\\ W_n\end{array}\right];$

$B_{\mathrm{Xi}}=\left[\begin{array}{cccccccc}{q}_{0i}& q_{\mathrm{xi}}& q_{\mathrm{yi}}& q_{\mathrm{zi}}& r_{0i}& r_{\mathrm{xi}}& r_{\mathrm{yi}}& r_{\mathrm{zi}}\\ 0& 0& 0& 0& q_{0i}& q_{\mathrm{xi}}& q_{\mathrm{yi}}& q_{\mathrm{zi}}\end{array}\right];$

$w_{1i}=\frac{1}{2}({q_{0i}}^2+{q_{\mathrm{xi}}}^2+{q_{\mathrm{yi}}}^2+{q_{\mathrm{zi}}}^2-1),w_{2i}=q_{0i}{r}_{0i}+q_{\mathrm{xi}}{r}_{\mathrm{xi}}+q_{\mathrm{yi}}{r}_{\mathrm{yi}}+q_{\mathrm{zi}}{r}_{\mathrm{zi}};$

$W_i={\left[\begin{array}{cc}{w}_{1i}& w_{2i}\end{array}\right]}^T;$

I represents the sequence number of Mars image;

4. normal equation is set up:

$\left[\begin{array}{cc}{A}^{T}A& A^{T}B\\ B^{T}A& B^{T}B\end{array}\right]\left[\begin{array}{c}t\\ X\end{array}\right]=\left[\begin{array}{c}{A}^{T}L\\ B^{T}L\end{array}\right]$

5. calculate the correction of Mars image direction element and closed points coordinate, its computing formula is:

$\left[\begin{array}{c}t\\ X\\ r\end{array}\right]=-N^{-1}{W}_Y$

$N=\left[\begin{array}{cc}{C}^{T}C& {B_X}^T\\ B_X& 0\end{array}\right]{,W}_Y=\left[\begin{array}{c}{C}^{T}L\\ W\end{array}\right]$

In formula, $C=\left[\begin{array}{cc}{A}^{T}A& A^{T}B\\ B^{T}A& B^{T}B\end{array}\right],$ R is the redundant computation vector for trim both members matrix dimension in computation process;

6. upgrade the value of the element of orientation and pass point, check whether the value of the correction of elements of exterior orientation and pass point is less than the threshold value of setting; If be less than threshold value, so iteration terminates, otherwise returns step and 3. continue iteration, until correction is less than the threshold value of setting.

Error equation structure as shown in Figure 4, because the reference mark number in experiment is more, typical several reference mark is only listed in Fig. 4, point 1 is obtained by two light on image one and two, therefore this point can list two error equations often opening on image, and nonzero element occupies two little black patches, and one piece relevant with image one, another block is relevant with image two, for coordinate corrective value then two little black patches be linked to be a vertical rectangle.Point 5 is the intersection point of three light, and have 3 black patches in Fig. 4, they are relevant with three width images respectively.The matrix of coefficients structural drawing of corresponding normal equation unknown number as shown in Figure 5.It is the sparse matrix of symmetrical band shape.For often opening A image ^tthe rank of A are 8 × 8, concerning a point to be located, and B ^tthe rank of B are 3 × 3, A ^tthe rank of B are 8 × 3.Air strips number is 1, has 3 photos, so the A of the whole district ^tthe exponent number of A is 24 × 24.Whole district's point to be located has n ' individual, then corresponding B ^tthe exponent number of B is 3n ' × 3n '.In figure, on principal diagonal, the exponent number of generous black patch is 8 × 8, and the exponent number of little side's black patch is 3 × 3, non-principal diagonal reflects the content of often opening about point to be located on image.Table 1 is for utilizing the precision result of dual quaterion strip aerial triangulation:

Table 1 dual quaterion strip aerial triangulation result

Table 1 result shows, dual quaterion strip aerial triangulation physical plane precision can reach 0.79cm, vertical accuracy can reach 1.51cm, and realistic accuracy and theoretical precision basically identical, result proves the reliable dem data of the method by making, realize the aerotriangulation of high-precision Mars image dual quaterion, thus be that the three-dimensional mapping of Mars image high-precision provides technical support.

Claims (1)

1. Mars DEM makes and a strip aerial triangulation method, and it is characterized in that, described method comprises the steps:

Steps A, carries out bilinear interpolation to laser-measured height data, makes Mars DEM; Its detailed process is as follows:

Steps A-1, the conversion between Mars geographic coordinate and Mars tangent plane coordinate is carried out to laser-measured height data:

$\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}{ccc}\sin B\cos L& \sin B\sin L& \cos B\\ -\sin L& \cos L& 0\\ \cos B\cos L& \cos B\sin L& \sin B\end{array}\right)\left(\begin{array}{c}(v+h)\cos L\cos B-X_0\\ (v+h)\cos L\sin B-Y_0\\ ((1-e^2)v+h)\sin L-Z_0\end{array}\right)$

In formula, (x, y, z) is Mars tangent plane coordinate, and (L, B, h) is the geographic coordinate of Mars ground point A, (X ₀, Y ₀, Z ₀) be the origin of tangent plane coordinate system, e ²=(a ²-b ²)/a ², a, b are respectively the length semiaxis of Mars ellipsoid;

Steps A-2, the cosine law is utilized to calculate the size of the triangle interior angle on alternative 3rd summit, the point selecting the maximum corresponding is this leg-of-mutton 3rd summit, set up triangular grid, on the triangle at grid points place, the value of net point is asked for, interpolation create-rule Grid DEM by distance weighting counting backward technique; Circular is as follows:

If triangular apex coordinate is (X _i, Y _i, Z _i), i=1,2,3, interpolation point coordinate is (X ₀, Y ₀, Z ₀), and grid points is in this triangle;

Then the height value of this triangulation network lattice point is:

$Z_0=\frac{\underset{i=1}{\overset{3}{\mathrm{\Σ}}}\frac{1}{{(X_i-X_0)}^2+{(Y_i-Y_0)}^2}{Z}_i}{\underset{i=1}{\overset{3}{\mathrm{\Σ}}}\frac{1}{{(X_i-X_0)}^2+{(Y_i-Y_0)}^2}};$

Step B, the subpoint position of compute ray bundle on Mars DEM, generates Mars analog image data;

Step C, adopts dual quaterion to carry out strip aerial triangulation, obtains Mars image direction element and closed points coordinate; Its detailed process is as follows:

Step C-1, dual quaterion $\hat{q}={\left[\begin{array}{cccc}{q}_0& q_x& q_y& q_z\end{array}\right]}^T+\mathrm{\ϵ}{\left[\begin{array}{cccc}{r}_0& r_x& r_y& r_z\end{array}\right]}^T$ The Mars image imaging equation represented is:

$\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]=\left[\begin{array}{c}{2q}_0{r}_x-{2q}_x{r}_0+{2q}_y{r}_z-{2q}_z{r}_y\\ {2q}_0{r}_y-{2q}_x{r}_z-{2q}_y{r}_0+{2q}_z{r}_x\\ {2q}_0{r}_z+{2q}_x{r}_y-{2q}_y{r}_x-{2q}_z{r}_0\end{array}\right]$

$\left[\begin{array}{c}{X}_A\\ Y_A\\ Z_A\\ 1\end{array}\right]=\left[\begin{array}{cccc}{a}_1& a_2& a_3& X_S\\ b_1& b_2& b_3& Y_S\\ c_1& c_2& c_3& Z_S\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]$

In formula,

$\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]=\begin{array}{}\end{array}\left[\begin{array}{ccc}{q}_0^2+q_x^2-q_y^2-q_z^2& 2(q_x{q}_y-q_0{q}_z)& 2(q_x{q}_z+q_0{q}_y)\\ 2(q_x{q}_y+q_0{q}_z)& q_0^2-q_x^2+q_y^2-q_z^2& 2(q_y{q}_z-q_0{q}_x)\\ 2(q_x{q}_z-q_0{q}_y)& 2(q_z{q}_y+q_0{q}_x)& q_0^2-q_x^2-q_y^2+q_z^2\end{array}\right];$

Q ₀, q _x, q _y, q _z, r ₀, r _x, r _y, r _zfor 8 elements of dual quaterion;

X _s, Y _s, Z _sfor element of orientation line element;

(X _a, Y _a, Z _a) be the coordinate of ground point A under Mars rectangular coordinate system;

(X, Y, Z) is the coordinate of ground point A under image space coordinate system;

Step C-2, input time, sequence Mars analog image data, were taken as 0 by the initial value of dual quaterion;

Step C-3, sets up dual quaterion error equation:

$V=\left[\begin{array}{cc}A& B\end{array}\right]\left[\begin{array}{c}t\\ X\end{array}\right]-L$

In formula,

$A=\left[\begin{array}{cccccccc}\frac{{\∂F}_x}{{\∂q}_0}& \frac{{\∂F}_x}{{\∂q}_x}& \frac{{\∂F}_x}{{\∂q}_y}& \frac{{\∂F}_x}{{\∂q}_z}& \frac{{\∂F}_x}{{\∂r}_0}& \frac{{\∂F}_x}{{\∂r}_x}& \frac{{\∂F}_x}{{\∂r}_y}& \frac{{\∂F}_x}{{\∂r}_z}\\ \frac{{\∂F}_y}{{\∂q}_0}& \frac{{\∂F}_y}{{\∂q}_x}& \frac{{\∂F}_y}{{\∂q}_y}& \frac{{\∂F}_y}{{\∂q}_z}& \frac{{\∂F}_y}{{\∂r}_0}& \frac{{\∂F}_y}{{\∂r}_x}& \frac{{\∂F}_y}{{\∂r}_y}& \frac{{\∂F}_y}{{\∂r}_z}\end{array}\right];$

$B=\left[\begin{array}{ccc}\frac{{\∂F}_x}{{\∂X}_A}& \frac{{\∂F}_x}{{\∂Y}_A}& \frac{{\∂F}_x}{{\∂Z}_A}\\ \frac{{\∂F}_y}{{\∂X}_A}& \frac{{\∂F}_y}{{\∂Y}_A}& \frac{{\∂F}_y}{{\∂Z}_A}\end{array}\right];$

T=[dq ₀dq _xdq _ydq _zdr ₀dr _xdr _ydr _z] ^t, d is for differentiating;

X＝[dX _AdY _AdZ _A] ^T；

L=[x-(F _x) ₀y-(F _y) ₀] ^t, x, y are respectively the coordinate of picture point corresponding to ground point A under photo coordinate system;

V＝[v _xv _y] ^T；

$F_x=-f\frac{a_1(X_A-X_S)+b_1(Y_A-Y_S)+c_1(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)};$

$F_y=-f\frac{a_2(X_A-X_S)+b_2(Y_A-Y_S)+c_2(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)};$

V _x, v _ybe respectively F _x, F _ycorrection;

F is the focal length of detector camera;

(F _x) ₀, (F _y) ₀by (X _a, Y _a, Z _a) and the value of the element of orientation substitute into F respectively _xand F _yobtain _;

In air strips, the element of orientation of all images meets restrictive condition:

B _Xt+W＝0

In formula,

$B_{\mathrm{Xi}}=\left[\begin{array}{cccccccc}{q}_{0i}& q_{\mathrm{xi}}& q_{\mathrm{yi}}& q_{\mathrm{zi}}& r_{0i}& r_{\mathrm{xi}}& r_{\mathrm{yi}}& r_{\mathrm{zi}}\\ 0& 0& 0& 0& q_{0i}& q_{\mathrm{xi}}& q_{\mathrm{yi}}& q_{\mathrm{zi}}\end{array}\right];$

$w_{1i}=\frac{1}{2}({q_{0i}}^2+{q_{\mathrm{xi}}}^2+{q_{\mathrm{yi}}}^2+{q_{\mathrm{zi}}}^2-1),$ w _2i＝q _0ir _0i+q _xir _xi+q _yir _yi+q _zir _zi；

W _i＝[w _1iw _2i] ^T；

I represents the sequence number of Mars image;

Step C-4, sets up normal equation:

$\left[\begin{array}{cc}{A}^{T}A& A^{T}B\\ B^{T}A& B^{T}B\end{array}\right]\left[\begin{array}{c}t\\ X\end{array}\right]=\left[\begin{array}{c}{A}^{T}L\\ B^{T}L\end{array}\right]$

Step C-5, calculate the correction of Mars image direction element and closed points coordinate, its computing formula is:

$\left[\begin{array}{c}t\\ X\\ r\end{array}\right]=-N^{-1}{W}_Y$

$\begin{array}{cc}N=\left[\begin{array}{cc}{C}^{T}C& {B_X}^T\\ B_X& 0\end{array}\right]& W_Y=\left[\begin{array}{c}{C}^{T}L\\ W\end{array}\right]\end{array}$

In formula, $C=\left[\begin{array}{cc}{A}^{T}A& A^{T}B\\ B^{T}A& B^{T}B\end{array}\right],$ R is the redundant computation vector for trim both members matrix dimension in computation process;

Step C-6, upgrades the value of the element of orientation and pass point, checks whether the value of the correction of elements of exterior orientation and pass point is less than the threshold value of setting; If be less than threshold value, so iteration terminates, otherwise returns step C-3 continuation iteration, until correction is less than the threshold value of setting.