CN104931022A - Satellite image three-dimensional area network adjustment method based on satellite-borne laser height measurement data - Google Patents

Abstract

The invention relates to a satellite image three-dimensional area network adjustment method based on satellite-borne laser height measurement data. The method comprises the following steps: respectively acquiring first control data and second control data corresponding to a connection point from a first control data source and a second control source; carrying out combined three-dimensional adjustment treatment for the first control data and second control data on the basis of a remote-sensing three-dimensional imaging model, acquiring a parameter correction value of the remote-sensing three-dimensional imaging model, wherein acquiring the first control data and the second control data includes acquiring a plane coordinate of the satellite-borne laser height measurement data; mapping the plane coordinate of the satellite-borne laser height measurement data to a satellite remote-sensing digital orthographic image, and acquiring a footprint image; forming second control data by the satellite-borne laser height measurement data and the footprint image or partial plane control data of the footprint image; acquiring homonymy points of a given quantity in the range of the footprint image or part of the footprint image as connection points, and using the control data of the homonymy points as first control data.

Description

Based on the satellite image stereoblock adjustment method of spaceborne laser altimeter system data

Technical field

The present invention relates to a kind of satellite remote-sensing image error compensation method, more specifically relate to a kind of satellite image stereoblock adjustment method based on spaceborne laser altimeter system data.

Background technology

Along with expanding economy, three-dimensional geographic information is used widely in each side such as digital earth, city planning, environmental protection.And the fast development of remote sensing satellite technology, space photogrammetry is made to have become the means of another quick obtaining three-dimensional geographic information after photogrammetric measurement, particularly in recent years along with the development of the technology such as three line scanner stereoscopic camera, the acquiring technology of 3 D Remote Sensing information achieves significant progress.

Such as, No. three, resource is the civilian stereoplotting satellite of China's first high precision, and succeeding in sending up and effectively applying of it has broken China for a long time to the dependence of external high precision satellite image, creates huge Social and economic benef@.But due to aspects such as optical stereo satellite stereoplotting mode and Satellite Attitude rail measuring accuracy, camera distortions, cause it under the condition of Pillarless caving, its measurement of higher degree precision is also difficult to the demand meeting high precision mapping.

Ground control point GCP (Ground Control Points) be satellite remote sensing image geometry correct and geo-location time important reference data sources.Correct in processing procedure at remotely sensing image geometric, for reaching certain correction precision, the ground control point of some is absolutely necessary, build video imaging model by the object coordinates at reference mark and corresponding picpointed coordinate and solving model parameter or compensation is optimized to existing imaging model solves compensating parameter, finally set up the correct transformational relation of object space and image space in imaging process.

Reference mark in traditional work process generally adopts the pattern of full field survey, need the operation through series of complexes such as " collect existing control data, survey that district is made an on-the-spot survey, reconnaissance is laying of markstone, field operation is measured, the arrangement of interior industry ", although along with the development of the advanced measuring techniques such as GPS-RTK, the workload that field operation is measured and complexity greatly reduce, but the reference mark field survey work of necessity is still inevitable.

In addition, because satellite remote-sensing image coverage is large (for resource No. three satellites, coverage is 50 kilometers × 50 kilometers), want to obtain equally distributed ground control point, the testing in the scope of hundreds of kilometer and even thousands of square kilometres of usual needs, its field process amount is big, labor intensive and material resources are self-evident.In addition, in the uninhabited areas such as disaster hotspots or virgin forest, marsh, desert such as earthquake, flood, rubble flow, survey crew usually cannot enter and measure on the spot.

Even if reference mark being by manually choosing remote sensing image and topomap same place obtains, also there is inefficiency, precision is difficult to the problems such as guarantee.And the different phase of areal or different sensors image are corrected, there will be again the situation of repetition reconnaissance.

Ground control point is generally divided into: flat high reference mark, planimetric control point and vertical control point.Chinese patent application 201310143369.3 discloses a kind of multi-source heterogeneous remote sensing image reference mark automatic acquiring method, can automatically extract remote sensing from multi-source heterogeneous image and control imaging point (reference mark image film), improve efficiency and the precision of control data acquisition, but it is in fact still a kind of automatic acquiring method of planimetric control point.

When ground control point data deficiencies, the area adjustment technology of satellite image can as a kind of important means of precise geometrical location.Such as, document " resource No. three cartographic satellite image planes and stereoblock adjustment compare " (" mapping journal " in April, 2014,43rd volume the 4th phase), for the problem that satellite remote-sensing image area adjustment under weak encounter conditions cannot correctly solve, the plane domain net adjusted data method that proposition utilizes digital elevation model (Digital Elevation Mode, DEM) to retrain as elevation is to improve its object positioning precision.But, because the method is only for image plane area adjustment, the problem of resource No. three cartographic satellite measurement of higher degree precision deficiencies therefore can not be solved.

Summary of the invention

In order to solve the problems of the technologies described above one of at least, according to embodiments of the invention, provide a kind of satellite image stereoblock adjustment method based on spaceborne laser altimeter system data, it comprises: obtain the first control data corresponding to tie point and the second control data respectively from the first control data source and the second control data source; Combine three-dimensional adjustment processing based on remote sensing stereopsis imaging model to the first control data and the second control data, obtain the parameter correction values of remote sensing stereopsis imaging model, wherein, the first control data source is storehouse, image reference mark; Second control data source is broad sense vertical control point storehouse; First control data is the control data of reference mark or reference mark image film; Second control data is the control data of broad sense vertical control point; The control data of this broad sense vertical control point comprises spaceborne laser altimeter system data, obtain the first control data corresponding to tie point and the second control data from the first control data source and the second control data source respectively, comprising: the planimetric coordinates obtaining spaceborne laser altimeter system data; The planimetric coordinates of spaceborne laser altimeter system data is mapped to satellite remote sensing digital orthophoto map and gathers footmark image; The plane control data of a part for spaceborne laser altimeter system data and footmark image or footmark image form the second control data jointly; The same place of the some in the scope of the part of footmark image or footmark image is obtained as tie point from storehouse, image reference mark, and using the control data of same place as the first control data.

According to the satellite image stereoblock adjustment method of the embodiment of the present invention, alternatively, remote sensing stereopsis imaging model is rational function model or rigorous geometry model, by expression formula V=AX+Bt-L, the unification in adjustment processing by rigorous geometry model and rational function model, wherein

$V=\left(\begin{array}{c}{v}_c\\ v_r\end{array}\right),$ $B=\left(\begin{array}{cccccc}1& c& r& 0& 0& 0\\ 0& 0& 0& 1& c& r\end{array}\right),L=\left(\begin{array}{c}{l}_c\\ l_r\end{array}\right)$

$t=\left(\begin{array}{c}d{a}_0\\ d{a}_1\\ d{a}_2\\ d{b}_0\\ d{b}_1\\ d{b}_2\end{array}\right),$

(c, r) is picpointed coordinate, for topocentric coordinates, (v _c, v _r) be picture point residual error, for Three Dimensional Ground coordinate corrective value, (a ₀, a ₁, a ₂, b ₀, b ₁, b ₂) correct parameter for the rigorous geometry model sight line vector of image corrects parameter or rational function model image space.

According to the satellite image stereoblock adjustment method of the embodiment of the present invention, alternatively, for the first control data, there is the first error equation:

V＝AX+Bt-L

Wherein, $A=\left(\begin{array}{c}\frac{\∂c}{\∂H}\\ \frac{\∂r}{\∂H}\end{array}\right),$ X=dH, or,

Wherein, a, b and d are the coefficients of plane equation,

For the second control data, then there is the second error equation:

V＝AX+Bt-L

Wherein,

Solve carry out combining three-dimensional adjustment processing by simultaneous first error equation and the second error equation.

According to the satellite image stereoblock adjustment method of the embodiment of the present invention, alternatively, adopt least square method to the first error equation and the second error equation simultaneous solution, comprise: by solving the normal equation corresponding with the first error equation and the second error equation: Nx=B obtains the model correction parameter of every scape image, wherein $N=\left(\begin{array}{cc}{A}^T\mathrm{PA}& A^T\mathrm{PB}\\ B^T\mathrm{PA}& B^{T}P{B}^T\end{array}\right),x=\left(\begin{array}{c}X\\ t\end{array}\right),B=\left(\begin{array}{c}{A}^T\mathrm{PL}\\ B^T\mathrm{PL}\end{array}\right),$ P is weight matrix.

According to the satellite image stereoblock adjustment method of the embodiment of the present invention, alternatively, for the second error equation, the mode adopting pointwise to cut unit forms the change normal equation relevant with the image space correction factor often opening image.

According to the satellite image stereoblock adjustment method of the embodiment of the present invention, alternatively, spaceborne laser altimeter system data are GLAS spaceborne laser altimeter system data.

According to the satellite image stereoblock adjustment method of the embodiment of the present invention, alternatively, after obtaining spaceborne laser altimeter system data from the second data source, through just slightly terrain data excluding gross error and echo waveform obtain the second control data being used for adjustment processing after screening.

According to the satellite image stereoblock adjustment method of the embodiment of the present invention, alternatively, comprised by echo waveform screening acquisition second control data: the background noise level of estimation echo waveform; Utilize the noise signal of smothing filtering algorithm filtering echo waveform; The Gauss curve fitting of filtered signal being carried out to initial parameter is estimated; Utilize the result that Gauss curve fitting is estimated, in conjunction with initial parameter, utilize nonlinear least square method to iterate and extract gaussian component parameter accurately.

According to the satellite image stereoblock adjustment method of the embodiment of the present invention, alternatively, adopt carry out Gauss curve fitting estimation, wherein w (t) be t backward energy value, W _mbe the backward energy value of m Gaussian function in t, Np is the number of Gaussian function, and ε is noise, and tm is the average of m Gaussian function, σ _mbe the mean square deviation of m Gaussian function, A _mthe peak value of m Gaussian function.

According to the satellite image stereoblock adjustment method of the embodiment of the present invention, alternatively, get off the imaging model mutually unification of GLAS laser-measured height data with stereoblock adjustment in coordinate system aspect by GLAS data being transformed into WGS84 coordinate system.

According to the satellite image stereoblock adjustment method based on spaceborne laser altimeter system data of the embodiment of the present invention, traditional area adjustment method basis is improved, adopt the method that plane is separated with elevation control data, improve adjustment Algorithm to the compatibility of all kinds of control data, unified general adjustment model is proposed, utilize spaceborne laser altimeter system data to overcome the problem that current Precision Elevation control data gathers difficulty, improve stereoplotting precision and the treatment effeciency of No. three, resource and subsequent satellites, and provide a kind of feasible solution route for improving No. three, resource in the area in hardship and the stereoplotting precision that cannot obtain control data area abroad.

Accompanying drawing explanation

In order to be illustrated more clearly in the technical scheme of the embodiment of the present invention, be briefly described below by the accompanying drawing of embodiment, apparently, the accompanying drawing in the following describes only relates to some embodiments of the present invention, but not limitation of the present invention.

Fig. 1 is according to an embodiment of the invention based on the indicative flowchart of the satellite image stereoblock adjustment method of spaceborne laser altimeter system data;

Fig. 2 shows picture point under rigorous geometry model and topocentric mathematical relation;

Fig. 3 shows the coordinate system that rigorous geometry model may be used;

Fig. 4 is the schematic diagram of the intersection of imaging direction and earth ellipsoid under rigorous geometry model;

Fig. 5 is the schematic diagram of spaceborne laser altimeter system reception and transmitted waveform;

Fig. 6 A ~ Fig. 6 C shows actual landform corresponding to GLAS footmark and echo waveform, and wherein Fig. 6 A is the actual image that footmark is corresponding; Fig. 6 B is the actual landform that footmark is corresponding; Fig. 6 C is GLAS actual ghosts waveform;

Fig. 7 schematically shows the process excavating the terrain information contained in return laser beam waveform further;

Fig. 8 is the schematic diagram that in database, GLAS surveys high data footprints image (left side) and echo waveform (right side); Fig. 9 is single scape stereopsis adjustment Scheme point position distribution schematic diagram;

Figure 10 is regional network three-dimensional adjustment Scheme point position distribution schematic diagram.

Embodiment

For making the object of the embodiment of the present invention, technical scheme and advantage clearly, below in conjunction with the accompanying drawing of the embodiment of the present invention, the technical scheme of the embodiment of the present invention is clearly and completely described.Obviously, described embodiment is a part of embodiment of the present invention, instead of whole embodiments.Based on described embodiments of the invention, the every other embodiment that those of ordinary skill in the art obtain under without the need to the prerequisite of creative work, all belongs to the scope of protection of the invention.

Unless otherwise defined, technical term used herein or scientific terminology should be in field belonging to the present invention the ordinary meaning that the personage with general technical ability understands." first ", " second " that use in patent application specification of the present invention and claims and similar word do not represent any order, quantity or importance, and are only used to distinguish different ingredients.Equally, the similar word such as " " or " " does not represent restricted number yet, but represents to there is at least one.

Because the track of optical satellite image, attitude and elements of interior orientation etc. all exist certain error, the Three Dimensional Ground coordinate therefore calculated by rigorous geometry model and actual coordinate value can there is some difference.In order to eliminate this species diversity, improve the geometric positioning accuracy of satellite, General Requirements introduces the ground control point of some.In order to improve the three-dimensional localization precision of satellite, traditional mode generally adopts flat high point as reference mark, namely field GPS measurement point or high-precision topomap measures the three-dimensional coordinate point of acquisition.

The control data at flat high reference mark generally includes panel data and altitude figures.Wherein, panel data can comprise longitude and the latitude of this ground control point; Altitude figures be this ground control point along pedal line direction to the distance of absolute datum, at survey field, be generally the height of ground control point on terrestrial gravitation direction, namely sea level elevation.

For traditional flat high reference mark, due to combination and this reference mark one_to_one corresponding of longitude, latitude and elevation, the control data including panel data and altitude figures also uniquely corresponds to concrete flat high reference mark.

And according to embodiments of the invention, propose a kind of three-dimensional error compensation method plane be separated with high process control.Specifically, when carrying out satellite image stereoblock adjustment, the control data adopted no longer is confined to the panel data corresponding with the flat high reference mark that a certain control data is originated and altitude figures, but can adopt the combination of the control data of originating from different control data.Described reference mark data of originating from different control data can have higher plane precision and higher vertical accuracy respectively, by carrying out local elevation constraint and carry out part plan constraint to the reference mark data with higher vertical accuracy carrying out simultaneous adjustment process to the reference mark data with higher level precision respectively, plane precision and/or the vertical accuracy of reference mark control data can be improved.

Further, the present invention proposes a kind of concept of new broad sense vertical control point.For this broad sense vertical control point, both traditional vertical control point can be comprised, also can be a planar region with approximate height value, wherein, the topomap mapping of different scale has different accuracy requirements to the topographic relief in this planar region or the discrepancy in elevation.So its altitude figures can be the altitude figures of concrete vertical control point, also can be by the elevation control data of laser-measured height satellite as various broad sense such as the altitude figures of ICESat/GLAS after waveform screens.

On the other hand, the sight line vector correcting method in tight imaging model and the image space correcting method in rational function model are unified by the satellite image stereoblock adjustment method according to the embodiment of the present invention, define unified three-dimensional adjustment Models.

Fig. 1 is according to an embodiment of the invention based on the indicative flowchart of the satellite image stereoblock adjustment method of spaceborne laser altimeter system data.

As shown in Figure 1, in step S100, obtain the control data of tie point (same place) respectively from different control data sources.

Specifically, by image Auto-matching from acquisition reference mark, storehouse, image reference mark or reference mark image film as tie point, and the control data (comprising plane control data and elevation control data) obtaining this tie point is as the first control data, wherein storehouse, image reference mark can be set up based on digital orthophoto map (Digital Orthophoto Map, DOM) by multi-source heterogeneous remote sensing image reference mark automatic acquiring method disclosed in Chinese patent application 201310143369.3.Such as, from broad sense vertical control point storehouse, spaceborne laser altimeter system database, obtain the control data (comprise plane control data and elevation control data) corresponding with this tie point as the second control data.Alternatively, according to embodiments of the invention, also first can obtain the second control data from broad sense vertical control point storehouse, then obtain the control data corresponding with this broad sense vertical control point as the first control data from storehouse, image reference mark.

Next, in step S300, based on remote sensing stereopsis imaging model, three-dimensional adjustment processing is combined to the first control data and the second control data, obtain the parameter correction values of imaging model.

The imaging model of remote sensing image can be divided into two large classes: physical model and universal model.

Physical model refer to be considered to as time cause the physical significance of deformation of image as factors such as surface relief, earth curvature, atmospheric refraction, camera distortions, then utilize these physical conditions to be construed as picture geometric model, the most representational is rigorous geometry model based on collinearity condition equation.

Universal model does not consider imaging mechanism, but direct mathematical function describes the geometric relationship of picture point and object point, and it has the feature of generality, confidentiality, high efficiency.Conventional universal imaging model has polynomial expression, direct linear transformation, affined transformation, rational function model (Rational FunctionModel is called for short RFM) etc.

Can unify to be applicable to rational function model and rigorous geometry model according to the stereoblock adjustment method that the plane of the embodiment of the present invention is separated with high process control, based on two kinds of models, the error compensation method according to the embodiment of the present invention is described respectively below.

RFM model is expressed as by picpointed coordinate (c, r) with correspondingly millet cake volume coordinate (wherein for latitude, λ is longitude, and H is elevation) be the polynomial ratio of independent variable:

Wherein, for standardized ground coordinate, (c _n, r _n) be standardized image coordinate.When there being reference mark, generally adopting translation, mode that the translation+image space such as rotation, affined transformation compensates revises rational function model, namely

Wherein (Δ c, Δ r) represents that image space compensates numerical value, and according to translation and an affined transformation, its expression-form is respectively (3) and (4),

Δc＝f ₁(r,c)＝a ₀(3)

Δr＝g ₁(r,c)＝b ₀

Δc＝f ₂(r,c)＝a ₀+a ₁c+a ₂r (4)

Δr＝g ₂(r,c)＝b ₀+b ₁c+b ₂r

By Taylor formula, by the compensating form of (4) formula affined transformation, once item is expanded to formula (2):

Being write as error equation is:

Above-mentioned (4) formula is substituted into, is:

Then above-mentioned (7) formula can be rewritten as V=Ax+Bt-L (8), wherein,

$V=\left(\begin{array}{c}{v}_c\\ v_r\end{array}\right),$ $B=\left(\begin{array}{cccccc}1& c& r& 0& 0& 0\\ 0& 0& 0& 1& c& r\end{array}\right),L=\left(\begin{array}{c}{l}_c\\ l_r\end{array}\right),$ $t=\left(\begin{array}{c}d{a}_0\\ d{a}_1\\ d{a}_2\\ d{b}_0\\ d{b}_1\\ d{b}_2\end{array}\right).$

For the plane of rule, can simultaneous equations:

wherein, a, b, d are the coefficient of plane equation.

Arranging into error equation is:

Being write as matrix form is:

Picture point and object space point strict geometric process one to one when what rigorous geometry model described is imaging, with conventional collinearity condition equation for representative.According to satellite almanac data, attitude data and time data, by imaging process, set up picture point (l, p) and ground point (X, Y, Z) or between mathematical relation.Fig. 2 shows picture point and topocentric mathematical relation.

Alternatively, high-precision satellite almanac data, attitude data and time data is provided by SPOT5 satellite image.

Coordinate system used by rigorous geometry model may comprise following coordinate system, as shown in Figure 3:

Navigational coordinate system (also claiming body coordinate system) O ₁-X ₁y ₁z ₁: be determined by the attitude control system of satellite (ACS), wherein initial point is at centroid of satellite, X1 axle, Y ₁axle, Z ₁axle gets three principal moments axles of satellite respectively.X ₁axle along satellite transverse axis, Y ₁axle points to satellite flight direction (not exclusively overlapping) along the longitudinal axis, Z ₁axle is determined according to right-handed system rule.Satellite attitude measurement carries out in body coordinate system, and three parameters describing its spatial attitude are: pitching (pitch) is around body coordinate system X ₁the rotation of axle, roll (roll) is around body coordinate system Y ₁the rotation of axle, course deviation (yaw) is around body coordinate system Z ₁the rotation of axle.

(ψ _x) _p(ψ _y) _pthe inclination angle in direction along ng a path and vertical track direction when representing charge-coupled image sensor (CCD) P components and parts imaging respectively.

Orbital coordinate system: be the interim coordinate system that navigational coordinate system and earth axes carry out transforming, the initial point of orbital coordinate system at centroid of satellite, Z ₂it is reverse that axle points to the earth's core, Y ₂axle points to the direction of satellite motion on satellite orbit face, X ₂axle is determined according to right-handed system rule.

O ₂-X ₂y ₂z ₂the direction of each axle is by the position of satellite in t and speed two direction vectors are determined, as follows:

${\stackrel{\→}{Z}}_2=\frac{\stackrel{\→}{P}\left(t\right)}{\left|\right|\stackrel{\→}{P}\left(t\right)\left|\right|}$

${\stackrel{\→}{X}}_2=\frac{\stackrel{\→}{V}\left(t\right)\mathrm{\Λ}{\stackrel{\→}{Z}}_2}{\left|\right|\stackrel{\→}{V}\left(t\right)\mathrm{\Λ}{\stackrel{\→}{Z}}_2\left|\right|}$

${\stackrel{\→}{Y}}_2={\stackrel{\→}{Z}}_2\mathrm{\Λ}{\stackrel{\→}{X}}_2$

Earth axes: the position of SPOT satellite and speed are all be defined in earth Protocol Ref coordinate system (ITRS), ITRF is the realization of ITRS, be responsible for setting up and safeguarding by the Terrestrial Reference Frame portion of IERS central office, the station coordinates of this reference frame and speed utilize such as very long baseline interferometry(VLBI (VLBI), Laser Measuring to defend etc. that data that space geodesy technique gathers calculate, it is a reference frame constantly changed, current employing be ITRF90, conveniently calculate, the coordinate system that general employing WGS84 coordinate system replaces ITRF framework to determine.This is because WGS84 is ECEF coordinate system, ITRF is then inertial coordinates system, and final mapping is it is desirable that body-fixed coordinate system, and ITRS and WGS84 difference is only at about 1 meter, and the conversion parameter from ITRF90 to WGS84 is as shown in the table:

T1(m)	T2(m)	T3(m)	D(ppm)	R1(arcsec.)	R2(arcsec.)	R3(arcsec.)
							0.060	-0.517	-0.223	-0.011	0.0183	-0.0003	0.0070

The conversion parameter table of ITRF90 to WGS84

Conversion formula is as follows:

$\left(\begin{array}{c}{X}_{\mathrm{WGS}84}\\ Y_{\mathrm{WGS}84}\\ Z_{\mathrm{WGS}84}\end{array}\right)=\left(\begin{array}{c}{T}_1\\ T_2\\ T_3\end{array}\right)+(1+D)\left(\begin{array}{ccc}0& -R_3& R_2\\ R_3& 0& -R_1\\ -R_2& R_1& 0\end{array}\right)\left(\begin{array}{c}{X}_{\mathrm{ITRF}}\\ Y_{\mathrm{ITRF}}\\ Z_{\mathrm{ITRF}}\end{array}\right)$

WGS84 coordinate system is a measurement reference system be made up of global the earth's core coordinate frame of reference and one group of corresponding model and WGS84 geoid surface.The initial point of this coordinate system is positioned at earth centroid, and Z axis and IERS point to identical with reference to pole, and X-axis points to IERS guide meridian with mistake initial point and perpendicular to the intersection point of Z axis plane, Y-axis is determined by right-handed system rule;

(l, p) → (X, Y, Z) is as follows in imaging process derivation:

I, calculate imaging time corresponding to l capable p row pixel

t＝t _c+lsp×(l-l _c)

Wherein, t _cfor the photography time of central scan row, lsp is the sweep time of often going, l _cfor the line number of central scan row

II, interpolation go out the position of t satellite with

$\stackrel{\→}{P}\left(t\right)=\underset{j=1}{\overset{8}{\mathrm{\Σ}}}\frac{\stackrel{\→}{P}\left(t_j\right)\×\underset{i\≠j}{\overset{8}{\underset{i=1}{\mathrm{\Π}}}}(t-t_i)}{\underset{i\≠j}{\overset{8}{\underset{i=1}{\mathrm{\Π}}}}(t-t_i)}$

$\stackrel{\→}{V}\left(t\right)=\underset{j=1}{\overset{8}{\mathrm{\Σ}}}\frac{\stackrel{\→}{V}\left(t_j\right)\×\underset{i\≠j}{\overset{8}{\underset{i=1}{\mathrm{\Π}}}}(t-t_i)}{\underset{i\≠j}{\overset{8}{\underset{i=1}{\mathrm{\Π}}}}(t-t_i)}$

III, the imaging inclination angle (ψ arranged by p _x) _p(ψ _y) _pcalculate (l, p) imaging direction vector under navigational coordinate system

${\stackrel{\→}{u}}_1^{\text{'}}=\left(\begin{array}{c}-\tan \left({\mathrm{\Ψ}}_Y\right)\\ +\tan \left({\mathrm{\Ψ}}_X\right)\\ -1\end{array}\right)$

${\stackrel{\→}{u}}_1=\frac{{\stackrel{\→}{u}}_1^{\text{'}}}{\left|\right|{\stackrel{\→}{u}}_1^{\text{'}}\left|\right|}$

Imaging direction vector under IV, calculating orbital coordinate system

${\stackrel{\→}{Z}}_2=\frac{\stackrel{\→}{P}\left(t\right)}{\left|\right|\stackrel{\→}{P}\left(t\right)\left|\right|}$

${\stackrel{\→}{X}}_2=\frac{\stackrel{\→}{V}\left(t\right)\mathrm{\Λ}{\stackrel{\→}{Z}}_2}{\left|\right|\stackrel{\→}{V}\left(t\right)\mathrm{\Λ}{\stackrel{\→}{Z}}_2\left|\right|}$ ${\stackrel{\→}{Y}}_2={\stackrel{\→}{Z}}_2\mathrm{\Λ}{\stackrel{\→}{X}}_2$

Wherein, with the attitude angle that the position of pixel (l, p) imaging moment t satellite and velocity interpolation go out t satellite respectively:

$a_p\left(t\right)=a_p\left(t_i\right)+(a_p\left(t_{i+1}\right)-a_p\left(t_i\right))\frac{t-t_i}{t_{i+1}-t_i}$

$a_r\left(t\right)=a_r\left(t_i\right)+(a_r\left(t_{i+1}\right)-a_r\left(t_i\right))\frac{t-t_i}{t_{i+1}-t_i}$

$a_y\left(t\right)=a_y\left(t_i\right)+(a_y\left(t_{i+1}\right)-a_y\left(t_i\right))\frac{t-t_i}{t_{i+1}-t_i}$

A _p(t _i), a _r(t _i), a _y(t _i) be respectively t _ithe attitude angle of moment satellite, can obtain from auxiliary data file.

By the imaging direction vector under navigational coordinate system be converted into the imaging direction vector under orbital coordinate system

${\stackrel{\→}{u}}_2=\frac{{\stackrel{\→}{u}}_2^{\text{'}}}{\left|\right|{\stackrel{\→}{u}}_2^{\text{'}}\left|\right|}$

${\stackrel{\→}{u}}_2^{\text{'}}=M_p\·M_r\·M_y\·{\stackrel{\→}{u}}_1$

$M_p=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(a_p\left(t\right)\right)& \sin \left(a_p\left(t\right)\right)\\ 0& -\sin \left(a_p\left(t\right)\right)& \cos \left(a_p\left(t\right)\right)\end{array}\right]$

$M_r=\left[\begin{array}{ccc}\cos \left(a_r\left(t\right)\right)& 0& -\sin \left(a_r\left(t\right)\right)\\ 0& 1& 0\\ \sin \left(a_r\left(t\right)\right)& 0& \cos \left(a_r\left(t\right)\right)\end{array}\right]$

$M_y=\left[\begin{array}{ccc}\cos \left(a_y\left(t\right)\right)& -\sin \left(a_y\left(t\right)\right)& 0\\ \sin \left(a_y\left(t\right)\right)& \cos \left(a_y\left(t\right)\right)& 0\\ 0& 0& 1\end{array}\right]$

Imaging direction vector under V, calculating earth axes

${\stackrel{\→}{u}}_3=\left[\begin{array}{ccc}{\left(X_2\right)}_X& {\left(Y_2\right)}_X& {\left(Z_2\right)}_X\\ {\left(X_2\right)}_Y& {\left(Y_2\right)}_Y& {\left(Z_2\right)}_Y\\ {\left(X_2\right)}_Z& {\left(Y_2\right)}_Z& {\left(Z_2\right)}_Z\end{array}\right]\·{\stackrel{\→}{u}}_2$

Wherein, X ₂, Y ₂, Z ₂for orbital coordinate system O-X ₂y ₂z ₂three change in coordinate axis direction vectors

VI, topocentric location:

As shown in Figure 4, measurement vector u ₃be M with the intersection point on ground, the position vector of satellite is m and the earth's core connect corresponding vector again because M point is a+h at major axis, minor axis is on the ellipsoid of b+h, so have after simultaneous equations:

$\stackrel{\→}{O_{3}M}=\stackrel{\→}{P}\left(t\right)+\mathrm{\μ}\×u_3\⇒\left\{\begin{array}{c}X=X_P+\mathrm{\μ}\×{\left(u_3\right)}_X\\ Y=Y_P+\mathrm{\μ}\×{\left(u_3\right)}_Y\\ Z=Z_P+\mathrm{\μ}\×{\left(u_3\right)}_Z\end{array}\right.$

$\frac{X^2+Y^2}{A^2}+\frac{Z^2}{B^2}=1$

Wherein,

$\left\{\begin{array}{c}A=a+h\\ B=b+h\end{array}\right.$

In conjunction with several transformation relations above, write as the form of similar collinearity condition equation:

$u_1=\left(\begin{array}{c}-\tan {\left({\mathrm{\Ψ}}_y\right)}_p\\ +\tan {\left({\mathrm{\Ψ}}_x\right)}_p\\ -1\end{array}\right)=\frac{1}{u}{R}_{21}^{-1}.R_{32}^1.\left[\begin{array}{c}X-X_P\\ Y-Y_P\\ Z-Z_P\end{array}\right]$

The relation of the imaging direction vector under acquisition navigational coordinate system and the observed ray of earth axes; Projection centre is (X _p, Y _p, Z _p), coordinate vector is (X, Y, Z), and projection coefficient is u, and rotation matrix is:

$R=R_{21}^{-1}[a_r\left(t\right),a_p\left(t\right),a_y\left(t\right)]R_{32}^{-1}[P\left(t\right),V\left(t\right)]=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right]$

Can collinearity equation be obtained:

$\frac{r_{11}(X-X_P)+r_{12}(Y-Y_p)+r_{13}(Z-Z_p)}{r_{31}(X-X_P)+r_{32}(Y-Y_p)+r_{33}(Z-Z_p)}-\tan {\left({\mathrm{\ψ}}_y\right)}_p=0$

Namely $\frac{r_{21}(X-X_P)+r_{22}(Y-Y_p)+r_{23}(Z-Z_p)}{r_{31}(X-X_P)+r_{32}(Y-Y_p)+r_{33}(Z-Z_p)}+\tan {\left({\mathrm{\ψ}}_x\right)}_p=0$

$\left\{\begin{array}{c}{f}_y=\arctan \left(\frac{r_{11}(X-X_P)+r_{12}(Y-Y_P)+r_{13}(Z-Z_P)}{r_{31}(X-X_P)+r_{32}(Y-Y_P)+r_{33}(Z-Z_P)}\right)-\left({\mathrm{\ψ}}_y\right)\≡0\\ f_x=\arctan \left(\frac{r_{21}(X-X_P)+r_{22}(Y-Y_P)+r_{23}(Z-Z_P)}{r_{31}(X-X_P)+r_{32}(Y-Y_P)+r_{33}(Z-Z_P)}\right)+\left({\mathrm{\ψ}}_x\right)\≡0\end{array}\right.---\left(12\right)$

For stereopsis, two equations above each image can be listed, if overlapping twice, then can list four equations, then solve (X, Y, Z) three unknown numbers, obtain the Three Dimensional Ground coordinate that image is corresponding.

Introduce sight line vector correction equation and carry out three-dimensional adjustment, sight line vector corrects equation and is:

$\begin{array}{c}{\left({\mathrm{\ψ}}_y\right)}_{\mathrm{ij}}^{\text{'}}={\left({\mathrm{\ψ}}_y\right)}_{\mathrm{ij}}+\mathrm{\Δ}{\left({\mathrm{\ψ}}_y\right)}_{\mathrm{ij}}={\left({\mathrm{\ψ}}_y\right)}_{\mathrm{ij}}-\frac{(a_0+a_{1}i+a_{2}j)}{k_{\mathrm{ij}}}\\ {\left({\mathrm{\ψ}}_x\right)}_{\mathrm{ij}}^{\text{'}}={\left({\mathrm{\ψ}}_x\right)}_{\mathrm{ij}}+\mathrm{\Δ}{\left({\mathrm{\ψ}}_x\right)}_{\mathrm{ij}}={\left({\mathrm{\ψ}}_x\right)}_{\mathrm{ij}}-\frac{(b_0+b_{1}i+b_{2}j)}{k_{\mathrm{ij}}}\end{array}---\left(13\right)$

Wherein, (i, j) refers to the ranks number of picpointed coordinate.

Therefore formula (12), can be written as:

$\begin{array}{c}{f}_{\mathrm{yij}}=\arctan \left(\frac{r_{11}(X-X_P)+r_{12}\left({Y-Y}_P\right)+r_{13}(Z-Z_P)}{r_{31}(X-X_P)+r_{32}(Y-Y_P)+r_{33}(Z-Z_P)}\right)-({\left({\mathrm{\ψ}}_y\right)}_{\mathrm{ij}}-\frac{a_0+a_{1}i+a_{2}j}{k_{\mathrm{ij}}})=0\\ f_{\mathrm{xij}}=\arctan \left(\frac{r_{21}(X-X_P)+r_{22}\left({Y-Y}_P\right)+r_{23}(Z-Z_P)}{r_{31}(X-X_P)+r_{32}(Y-Y_P)+r_{33}(Z-Z_P)}\right)+({\left({\mathrm{\ψ}}_x\right)}_{\mathrm{ij}}-\frac{b_0+b_{1}i+b_{2}j}{k_{\mathrm{ij}}})=0\end{array}---\left(14\right)$

Be abbreviated as:

$\left\{\begin{array}{c}{F}_1=k_{\mathrm{ij}}{f}_1+(a_0+a_{1}i+a_{2}j)=0\\ F_2=k_{\mathrm{ij}}{f}_2+(b_0+b_{1}i+b_{2}j)=0\end{array}\right.---\left(15\right)$

Wherein,

$\left\{\begin{array}{c}{f}_1=\arctan \left(\frac{h_1}{h_3}\right)-{\left({\mathrm{\ψ}}_y\right)}_{\mathrm{ij}}\\ f_2=\arctan \left(\frac{h_2}{h_3}\right)+{\left({\mathrm{\ψ}}_x\right)}_{\mathrm{ij}}\end{array}\right.,$

$\left\{\begin{array}{c}{h}_1=r_{11}(X-X_P)+r_{12}(Y-Y_P)+r_{13}(Z-Z_P)\\ h_2=r_{21}(X-X_P)+r_{22}(Y-Y_P)+r_{23}(Z-Z_P)\\ h_3=r_{31}(X-X_P)+r_{32}(Y-Y_P)+r_{33}(Z-Z_P)\end{array}\right.$

Above for rigorous geometry model coordinate system analyze basis on, according to embodiments of the invention, in order to the unification mutually of the aforementioned adjustment analysis based on rational function model, with picpointed coordinate (c, r) replace (l, p), use topocentric coordinates replace (X, Y, Z) to state.

Same image adopts same sight line vector corrected parameter, then, time for stereoblock adjustment, can list error equation to (15) formula is:

In above formula, (v ₁, v ₂) be residual error, be the Three Dimensional Ground coordinate corrective value of n-th, (a _0m, a _1m, a _2m, b _0m, b _1m, b _2m) be the sight line vector correction parameter that m opens image, k _mnbe at n-th and open the distance value between satellite imagery center corresponding to image with m.

Contrast (16) formula and aforementioned (7) formula, the lower target factor of removing, can find out that two expression formulas are basically identical, be all the error equation about Three Dimensional Ground coordinate corrective value and imaging model correction.

So, can aforementioned expression (8) be passed through, V=AX+Bt-L, rigorous geometry model and rational function model are united in adjustment processing.

Wherein

Open on image for i-th point in jth, then bring corresponding subscript respectively.Such as, open on image for i-th point in jth, then V _ij=A _ijx _i+ B _ijt _j-L _ij, weight matrix is P _ij.

Next, on the basis of the error equation derived at above-mentioned imaging model, error equation is drafted for the first control data and the second control data respectively.

For foregoing the second control data from broad sense vertical control point data source, the height value of high process control can be substituted into above-mentioned expression formula (8) as actual value, then have

or

For the first control data, if planimetric control point, can using the planimetric coordinates of plane control as actual value, then

$\begin{array}{cc}A=\left(\begin{array}{ccc}1& 1& \frac{\∂c}{\∂H}\\ 1& 1& \frac{\∂r}{\∂H}\end{array}\right)& x=\left(\begin{array}{c}0\\ 0\\ \mathrm{dH}\end{array}\right)\end{array}$ Or $\begin{array}{cc}A=\left(\begin{array}{c}\frac{\∂c}{\∂H}\\ \frac{\∂r}{\∂H}\end{array}\right)& x=\mathrm{dH}\end{array}---\left(19\right)$

For the first control data, if reference mark image film, then use the plane equation of aforementioned expression (9), obtain error equation (identical with previously described formula (11)):

Then, adopt least square method to carry out adjustment to resolve.

Such as, plane control and high process control object coordinates correction are similar to 0, can virtual observation equation be listed:

Cx＝0 (21)

Wherein, C=(1 1 1)

For V=Ax+Bt-L, weight matrix is P; Can list normal equation by least square is:

NX＝Q (22)

Wherein, $N=\left(\begin{array}{cc}{A}^T\mathrm{PA}& A^T\mathrm{PB}\\ B^T\mathrm{PA}& B^{T}P{B}^T\end{array}\right),X=\left(\begin{array}{c}x\\ t\end{array}\right),Q=\left(\begin{array}{c}{A}^T\mathrm{PL}\\ B^T\mathrm{PL}\end{array}\right)$

Above-mentioned normal equation is typical solving sparse linear equations.For ease of calculating, can for above formula (18), the mode adopting pointwise to cut unit forms the change normal equation relevant with the image space correction factor often opening image, such as, the exponent number of changing normal equation is 6m*6m (wherein m is the image number participating in adjustment), and now changing normal equation is still sparse vectors.

Alternatively, certain weights P can also be given according to the measuring accuracy at reference mark.

The precision of picpointed coordinate observed reading: σ _xmillimeter, then weigh P _x=1;

The precision of ground control point obserred coordinate value is: σ _crice, then weigh wherein for photographic scale.Basic law is that precision is higher, i.e. σ _cless, weights are larger, also larger to the contribution of adjustment result.

Finally separate the model correction parameter that sparse vectors obtains every scape image.Such as, solved by the above-mentioned equation of simultaneous (18) and equation (19) and/or equation (20), obtain the model correction parameter of elevation coordinate and planimetric coordinates.Aforementioned equation (18), (19) and (20) are general formulas, concrete coordinate points can be represented such as formula the mode of (17) by use-case, to set up the corresponding relation of coordinate points and the image of participation adjustment that need correct, thus obtained the model correction parameter of every width image by equation solution.

Be presented above the process according to associating first control data of the embodiment of the present invention and the adjustment processing of the second control data.By this process, the correction to elevation coordinate and planimetric coordinates can be carried out simultaneously, improve the precision of elevation coordinate and planimetric coordinates simultaneously.

Alternatively, only can carry out the correction of elevation coordinate for the first control data, or only the second control data be carried out to the correction of its planimetric coordinates.That is, when respectively in above-mentioned equation (18), (19) and (20) being resolved, the control data that this equation is corresponding is only used.Like this, the correction of elevation coordinate or planimetric coordinates can be carried out respectively, improve precision respectively.As previously mentioned, according to the method for the embodiment of the present invention, the precision of elevation coordinate and planimetric coordinates can be improved simultaneously.

As mentioned above, according to the satellite image stereoblock adjustment method of the embodiment of the present invention, it combines the control data from different control data source, combine and carry out three-dimensional adjustment and resolve, both the correction of elevation coordinate or planimetric coordinates can have been carried out respectively, also can carry out the correction of elevation coordinate and planimetric coordinates, the process that this plane namely alleged by the present invention is separated with elevation control data simultaneously.

According to embodiments of the invention, spaceborne laser altimeter system data can be adopted as the control data of broad sense vertical control point database.

Spaceborne laser altimeter system (Laser Altimeter) is that satellite earth observation obtains one of important technical of three-dimensional geographic information.Spaceborne laser altimeter system is a kind of ground point High head ratio technology, it take satellite as platform, carry laser ceilometer from space, the time observes the earth, high precision, measure distance between satellite to testee in real time, and by data process&analysis, obtain the information such as landforms, vegetative coverage situation, ocean surface topography of earth surface.External spaceborne laser altimeter system technical development is rapid, and world's Main Developed Countries is all in the research carrying out spaceborne laser altimeter system instrument.

Such as, the U.S. transmitted ICESat satellite in 2003, geoscience laser-measured height instrument system (the Geoscience Laser Altimeter System that it carries, GLAS) be first the satellite borne laser range measurement system for the Continuous Observation earth in the whole world, its main science object measures polar ice sheet elevation and change thereof, cloud layer and aerocolloidal distribution characteristics etc.ICESat satellite is period in orbit, and obtain a large amount of high-precision altitude figures, its laser footmark plane precision reaches 10m magnitude, and vertical accuracy is about 15cm.Within 2017, be about to launch ICESat-2, the laser ceilometer of this satellite will adopt 6 bundle laser of 3*2 structure, and spot size also will narrow down to 10 meters, and the quantity of the laser spots of acquisition and precision have and significantly improve.The LIST plan proposed by U.S. NRC then will be launched in about the year two thousand twenty, and the photon detection array that employing 1000 is restrainted by this satellite, spot diameter narrows down to 5 meters further, can carry out the topographical surveying of global seismic 5m spatial resolution and 10cm vertical precision fast.

Domestic attach great importance in recent years spaceborne laser altimeter system research.Moon exploration program series of satellites " Chang'e I " and " Chang'e-2 " satellite of independent research have all carried high-precision laser altitude gage, obtain moonscape 3 D stereoscopic image.Except " goddess in the moon " series of satellites carries except laser altimeter, also not used for the laser-measured height system of earth observation, but list follow-up earth observation laser-measured height satellite development plan in, as the follow-up star task of resource No. three high resolving power stereoplotting satellites, satellite surveying and mapping " 12 " planning also claims to the work of laser-measured height satellite technology parameter demonstration.

According to embodiments of the invention, the spaceborne laser altimeter system data of GLAS can be adopted to originate as the control data of broad sense vertical control point database.

ICESat satellite orbital altitude is about 600km, and orbit inclination is 94 °, can the latitude scope of covering the whole world ± 86 °.System operating frequency is 40Hz, and pulse width is 4ns, and p.s., emission wavelength was the infrared laser of 1064nm and the green laser of 532nm simultaneously.The laser footmark that laser pulse is formed at the earth's surface is the hot spot of diameter 70m, and the serial spot separation of same track is about 170m, and the interval that adjacent orbit is gone up under the line is about 15km, is spaced apart 2.5km the high latitude area of 80 °.GLAS have employed the satnav of high-precision waveform collection equipment, high efficiency Echo Processing algorithm and precision and determines appearance equipment, and make its system have very high precision, wherein plane precision can reach 10m magnitude, and vertical accuracy is about 15cm.ICESat satellite data product is generated by Ge Dade space flight center I-SIPS (ICESAT ScienceInvestigator-led Processing System), every physical parameter is extracted according to the echo information that GLAS receives, can obtain the normal data product of 15 kinds, naming method is GLA01-GLA15.System data file is by American National ice and snow data center (National Snow andIce Data Center, Distributed Active Archive Center, NASA NSIDC DAAC) externally announce, be divided into L1 level and L2 level two DBMS, wherein L1 DBMS is divided into again L1A and L1B two ranked datas.L1A DBMS is without data processing, records based on seasonal effect in time series full resolution raw data, comprises GLA01-GLA04.L1B DBMS is corrected through device by L1A DBMS, data processing, based on attitude and position data carry out laser footpoint location formed, comprise GLA05-GLA07.L2 DBMS is that L1B DBMS is more accurate after various correcting process, and is divided into elevation and the position data of different ground surface type, and for distance, elevation correction various correction data, comprise GLA08-GLA15.Wherein GLA12 data file is ice sheet altitude figures, and GLA13 is sea ice altitude figures, and GLA14 is land altitude figures, and GLA15 is ocean altitude figures.

As previously mentioned, the imaging mode of three-dimensional intersection cause satellite image without control precision particularly vertical accuracy have gap from mapping application, its precision must be improved by ground control point data.Spaceborne laser altimeter system is high without control vertical accuracy, about 0.15m, and satellite can obtain global altitude figures, and the elevation control data of satellite image can be it can be used as to use, reduce the work at artificial surface deployment Survey control point, and the High Precision Stereo mapping in area overseas can be realized.

In the step S100 of aforementioned three-dimensional adjustment processing, obtain the second control data from broad sense vertical control point storehouse, then obtain the control data corresponding with broad sense vertical control point as the first control data from storehouse, image reference mark.If adopt the spaceborne laser altimeter system data of GLAS, because the ground footmark diameter of its correspondence is larger, about 70m, thus plane precision is not high especially, about 10m, therefore be first mapped to known DOM striograph based on the planimetric coordinates of spaceborne laser altimeter system data and gather footmark image, this footmark image is used for carrying same place at image reference mark database and makes reference use, can not do with the image film in the database of image reference mark and mate by footmark image, only need the same place be chosen at from the database of image reference mark in footmark image capturing range, also therefore, described known DOM striograph precision does not need very high, precision about 20 meters is just passable.In short, first obtain a footmark image (plane information) roughly based on the planimetric coordinates of spaceborne laser altimeter system data exactly, then the same place of the some in footmark image capturing range is obtained from storehouse, image reference mark, and using the control data of this same place as the first control data.

Preferably, first through just omiting the steps such as terrain data (as SRTM) excluding gross error, waveform screening, spaceborne laser altimeter system data being selected, the precision of subsequent treatment can be improved like this.

Illustrate that the spaceborne laser altimeter system shape information in conjunction with ICESat satellite is analyzed the terrain feature in footmark below, filter out available data as elevation control data.

From principle, distance measurement value is transmitted waveform and the half receiving light velocity product in waveform center of gravity moment difference and vacuum, and Fig. 5 is the schematic diagram of spaceborne laser altimeter system reception and transmitted waveform, and wherein abscissa line is time shaft.

Because laser firing pulses is being Gaussian distribution along transmit direction and Vertical Launch direction, terrain and its features is then main in echo-peak, broadening, position etc. that crest is corresponding on echo impact, so Gaussian function still can be adopted to carry out Fitting Analysis for echo waveform, just for the landform of complexity, multiple Gaussian function may be needed to carry out combining and just can reach best-fit effect.

In addition because laser firing pulses is through propagation in atmosphere, multiple processes such as target reflection, echo reception opto-electronic conversion are mixed into much noise, therefore need to estimate through noise, smothing filtering, Gaussian parameter pre-estimation and matching extracting parameter.

According to the characteristic distributions of echoed signal and the character of Gaussian function, determine that the disposal route of echoed signal comprises: estimation background noise level; Utilize smothing filtering algorithm filtering overwhelming majority noise signal; Based on the character of Gaussian function, filtered signal is carried out to the estimation of initial parameter; In conjunction with initial parameter, utilize nonlinear least square method to iterate and extract gaussian component parameter accurately.

For the echo waveform of reality, can adopt carry out matching, wherein w (t) be t backward energy value, W _mbe the backward energy value of m Gaussian function in t, Np is the number of Gaussian function, and ε is noise, and tm is the average of m Gaussian function, σ _mbe the mean square deviation of m Gaussian function, A _mthe peak value of m Gaussian function.Can the parameter of each Gaussian function obtained of solution by the method for least square, the corresponding relation between quantitative test landform and waveform.Alternatively, before carrying out Gauss curve fitting estimation, by calculating the degree of bias and the kurtosis value of echo waveform, the waveform gaussian shape of echo can be estimated, so that follow-up Gaussian parameter calculates.

The Skewness degree of bias is the tolerance of statistics distribution skew direction and degree, is the numerical characteristic of statistics distribution degree of asymmetry.For the degree of bias of n sample value, common degree of bias estimation formulas:

$\mathrm{Sk}=\frac{1}{n-1}\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^3}{{\mathrm{\σ}}^3}---\left(23\right)$

Wherein: n is sample number, σ is sample standard deviation, for sample average.

Under general situation, when statistics is right avertence distribution, Sk>0, and Sk value is larger, right avertence degree is higher; When statistics is left avertence distribution, Sk<0, and Sk value is less, left avertence degree is higher; When statistics is symmetrical, obviously there is Sk=0.

Kurtosis kurtosis is used to the index reflecting the high and steep or flat degree of curve of frequency distribution top point, is used for metric data at center aggregation extent.Sometimes the arithmetical mean of two groups of data, standard deviation are all identical with the coefficient of skew, but their the towering degree on distribution curve top is different.Its estimation formulas is:

$K=\frac{1}{n-1}\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_i-\stackrel{\&OverBar;}{x})}^4}{{\mathrm{\&sigma;}}^4}-3---\left(24\right)$

In normal distribution situation, coefficient of kurtosis value is 0.Positive coefficient of kurtosis illustrates that observed quantity is more concentrated, has the afterbody longer than normal distribution; Negative coefficient of kurtosis illustrates that observed quantity is so unconcentrated, has the afterbody shorter than normal distribution, is similar to being uniformly distributed of rectangle.The standard error of coefficient of kurtosis is used for judging the normality of distribution.The ratio of coefficient of kurtosis and its standard error is used for checking normality.If this ratio absolute value is greater than 2, normality will be refused.Whether the coefficient of skewness (Skewness) is used for tolerance distribution symmetrical.Be symmetrical about normal distribution, the coefficient of skewness is 0.Larger on the occasion of showing that this distribution has comparatively long-tail portion, right side.Larger negative value shows there is comparatively long-tail portion, left side.Can be used for same with the ratio of its standard error of the coefficient of skewness checks normality.

According to the principle of laser ranging, on time, the corresponding actual range in 1ns interval is 15cm, therefore waveform can calculate corresponding discrepancy in elevation information by the time corresponding to crest after decomposing, and with reference to the optical image data that this footmark is corresponding, can realize becoming more meticulous decomposition to the elevation in footmark.Fig. 6 A ~ 6C shows actual landform corresponding to GLAS footmark and echo waveform, and wherein Fig. 6 A is the actual image that footmark is corresponding; Fig. 6 B is the actual landform that footmark is corresponding; Fig. 6 C is GLAS actual ghosts waveform.In figure 6 c, abscissa axis is the relative time receiving echo, and unit nanosecond (ns), axis of ordinates is echo strength, unit volt (volts).

Alternatively, adopt the echo analysis model become more meticulous, can excavate the terrain information contained in return laser beam waveform further, its illustrative process as shown in Figure 7.

Further, need in coordinate system aspect, the imaging model that GLAS laser-measured height data and three-dimensional LAN (Local Area Network) adjustment adopt to be united.Such as, under GLAS data being transformed into WGS84 coordinate system.

Specifically, because the reference ellipsoid of GLAS laser-measured height data is TOPEX/Poseidon, under needing first GLAS data to be transformed into WGS-84 reference ellipsoid, be convenient to follow-uply carry out precision test and application.Because two reference ellipsoids are less in weft direction difference, warp direction, without correction, therefore only can consider the correction of height value when carrying out data transformations.General employing experimental formula transforms, and greatly can simplify the calculating of Data Comparison between two different ellipsoids, such as formula (25).

$\left\{\begin{array}{c}\mathrm{dh}={\cos }^2\mathrm{Bda}+{\sin }^2\mathrm{Bdb}\\ H_G=H_g+\mathrm{dh}\end{array}\right.---\left(25\right)$

In formula (25), B represents latitude, the difference of what da represented is WGS-84 ellipsoid and TOPEX/Poseidon ellipsoid major axis, db represents the difference of WGS-84 ellipsoid and TOPEX/Poseidon ellipsoid minor axis, dh is elevation correction value, the elevation of Hg relative TOPEX ellipsoid before being conversion, H _gfor the GLAS data elevation after conversion.

GLAS laser beam cannot pass through cirrus and Bao Yun in communication process, is subject to mountain region in addition, the impact of the complicated landform such as hills, may cause part GLAS laser-measured height data value far away higher than real top elevation.Alternatively, for this part rough error value, first by carrying out contrasting being rejected with dem data.Such as, SRTM data and GLAS altitude figures is adopted to contrast.SRTM data take WGS84 as reference ellipsoid, adopts the positive height of EGM96 geoid surface.For the ease of contrast, SRTM is transformed into EGM2008 from EGM96.Arranging elimination of rough difference threshold value is 20m, and the point that the difference of elevation by GLAS altitude figures and SRTM DEM is less than-20m or is greater than 20m is removed as rough error point.Through the difference of elevation of two kinds of data is carried out statistical study, find to meet normal distribution, the standard deviation of setting ± 2 times is further threshold value, only retains the GLAS altitude figures within threshold value.Utilize rough terrain data before carrying out waveform screening, to carry out scalping rejecting to GLAS laser-measured height data as SRTM, the precision of waveform subsequent screening can be improved.

Available ICESat/GLAS spaceborne laser altimeter system data are obtained by waveform screening, then find image blocks according to its planimetric coordinates roughly from known DOM figure, survey the control data that high data and image blocks plane control data form broad sense vertical control point database jointly.For existing GLAS spaceborne laser altimeter system data, the plane error of 70 rice diameter hot spots about 10 meters, in conjunction with the error of the image blocks obtained from DOM figure self, composition error is at about 20 meters.According to the method for the embodiment of the present invention, the planimetric coordinates of described image blocks is without the need to accurately, and it is only as the reference choosing tie point (same place).

As previously mentioned, the planimetric coordinates of broad sense vertical control point is coarse relative to elevation, but can based on the laser spots geometric coordinate corresponding with broad sense vertical control point and footmark spot size, be aided with footmark image, the same place initial position of measurable satellite image, and the same place being obtained satellite image by image matching technology.That is, by three-dimensional coordinate and the footmark size of laser-measured height data, same place is found out from image in restriction planar range, finally obtain laser-measured height data vertical control point pair, this is to comprising image space ranks coordinate figure corresponding to accurate object space height value, the object space planimetric coordinates just omited and image.

Illustrate below and how to set up broad sense vertical control point database.

The name shown in database should adopt directly perceived, succinct mark as far as possible, name will comprise classification and the title thereof of table.The title prefixing " TB " of table, classification and title centre add " _ " and split.

For the ease of operation and the management of data, create GLAS footmark tables of data " TB_GLAS_PT " and GLAS footmark shadow table " TB_IMG ".Be mainly used in table " IMG_TABLE " store the image blocks information surveying high data grid technology with GLAS.List structure as shown in Table 1 and Table 2.

Table 1 GLAS_PT_TABLE list structure

Table 2 IMG_TABLE list structure

If use Oracle data platform to develop, then in oracle database, create GCPDB table space by data and data base administration subsystem keeper, table space arranges as follows:

● initial file size 5GB;

● AutoExtend1024KB is set;

● largest amount 20GB or unrestricted;

● the block size in table space can not be less than 8192 bytes;

And increasing Oracle user GCPAdmin newly by data and data base administration subsystem, password informs user by after Administrator, and arranges GCPAdmin account and have complete operation authority to GCPDB table space.

● database is named: GCPDB

● database user: GCPAdmin and there is reference mark image warehouse-in, there is management system

The registered user of system operating right.

In data base design process, by distinguishing different visitors, different access types and different data objects, treat respectively to carry out database security design.

Concrete example, reference mark, broad sense vertical control point storehouse data source can adopt the data of GLAS laser-measured height data 2007-2009, and wherein, echo data is included in GLAS01 data, and altitude figures is GLA14 data; Data area is 64 1,000,000 figure amplitude ranges within Chinese territory, totally 1202338 points after preliminary screening.The national true color orthography result map (DOM) of No. three, resource of the plane precision about 15 meters that footmark image data adopts, this result map is without elevation information, and plane precision does not have image reference mark database precision high.

Utilize C# coding at Microsoft Visual Studio 2010 platform, the GLAS reading binary format surveys high data, and is stored into database.

GLAS laser footmark image data have chosen the image data of the regional extent for precision test and Preliminary Applications.Dominating pair of vertices based on image has very important science and realistic meaning in basic surveying and mapping product production and renewal, have very large effect for the relation analyzing GLAS echo waveform and terrain feature, therefore footmark image blocks data are put in storage as important data source simultaneously.C# coding is utilized at Microsoft VisualStudio 2010 platform, utilize the raster spatial data transformation warehouse GDAL (Geospatial Data Abstraction Library) increased income, the image blocks data surveying high data grid technology with GLAS are stored in database.Fig. 8 is the schematic diagram that in database, GLAS surveys high data footprints image (left side) and echo waveform (right side).Be similar to Fig. 6 C, in the echo waveform figure of Fig. 8 right part, abscissa axis is the relative time receiving echo, and unit nanosecond (ns), axis of ordinates is echo strength, unit volt (volts).Footmark image and echo waveform are not contacted directly, and echo waveform is only relevant with characters of ground object with the terrain feature in footmark, and as previously mentioned, footmark image obtains from DOM according to the planimetric coordinates of laser.

In order to verify the precision of GLAS broad sense vertical control point, collect airborne LiDAR (the Light Detection And Ranging of Taiyuan region certain limit, laser radar) data and Efficiency in Buildings in Tianjin Area dem data, carry out precision test to GLAS vertical control point.

Taiyuan region LiDAR point cloud data break is about 3m, and coordinate system adopts WGS-84, and elevation is WGS-84 geodetic height.Plane precision is about 0.5 ~ 0.8m, and absolute elevation precision is better than 0.15m.

For Taiyuan region LiDAR point data, two kinds of methods are adopted to carry out precision evaluation to GLAS vertical control point.

If survey high data through the GLAS of preliminary screening also there is rough error value, can in conjunction with footmark image and actual landform excluding gross error value.

First method is the difference of the elevation calculating the elevation of GLAS vertical control point and the nearest LiDAR data point of its footmark scope inner plane.Statistical analysis, discrepancy in elevation medial error is 0.603m, and maximum error is 1.617m, and least error is 0.017m.Statistics is as shown in table 3.

Due to the elevation average within the scope of the hot spot of the elevation of GLAS point to be diameter be 70m, therefore second method is centered by GLAS vertical control point, calculate the elevation mean value of all LiDAR point within the scope of radius 35m footmark, obtain the height value of GLAS point and the difference of LiDAR point elevation mean value.Through statistics, discrepancy in elevation medial error is 0.531m, and maximum error is 1.645m, and least error is 0.016m.Statistics is as shown in table 4.

Table 3 GLAS point and LiDAR closest approach discrepancy in elevation statistical form

Discrepancy in elevation absolute value (m)	Quantity	Number percent
			<0.20	11	26.8％
<0.40	22	53.7％
			<0.50	30	73.2％
<0.80	34	82.9％

Table 4 GLAS point and LiDAR point elevation average discrepancy in elevation statistical form

Discrepancy in elevation absolute value (m)	Quantity	Number percent
			<0.20	11	26.8％
<0.40	28	68.3％
			<0.50	30	73.2％
<0.80	34	82.9％

Efficiency in Buildings in Tianjin Area test figure is the DEM result map of 1:2000 engineer's scale, and this DEM achievement is WGS-84 through coordinate conversion back plane coordinate system, and height datum is WGS-84 geodetic height.Be about 1.5m through measuring and calculating DEM mesh spacing, absolute elevation precision meets the code requirement of national primary scale DEM result map.

Be discrete vector point data because GLAS surveys high data, dem data is raster data, needs to be carried out by the form of two kinds of data unifying just can be analyzed.Utilize the Spatial Analyst Tool in ArcToolBox, extract DEM height value to point (extract values to points), extracting GLAS by the height value in DEM raster data surveys in high data, as the attribute field of GLAS vector point data, carry out follow-up comparison and statistical study.

Statistical analysis, the elevation of the GLAS laser spots of Efficiency in Buildings in Tianjin Area after screening is 0.706m with the difference medial error of the height value extracted, and maximum error is 1.767m, and least error is 0.0004m.Statistics is as shown in table 5.

Table 5 GLAS point and DEM spot elevation discrepancy in elevation statistical form

Discrepancy in elevation absolute value (m)	Quantity	Number percent
			<0.20	156	19.8％
<0.40	339	43.1％
			<0.50	397	50.5％
<0.80	563	71.6％

Analyze through carrying out reference to GLAS laser footmark image, find that the laser spots that error is larger is all be positioned at mountain area, or landform changes and has the region of significant little ditch, visualize is the relation of elevation difference and landform, illustrates that the precision of GLAS laser-measured height data and terrain feature have close relationship.In the region that relief is larger, elevation difference is obvious, and the GLAS point that such as Taiyuan region error is larger is substantially all positioned at mountain region; In the area that landform is smooth, elevation difference is general less.Aggregate analysis, the GLAS laser-measured height data after preliminary screening have higher precision, can use as the high process control of 1:5 ten thousand or 1:2.5 ten thousand.

Usually, when utilizing ground control point, some position and the distribution characteristics at reference mark can clearly be recognized by the description of station at corresponding reference mark.Due to the hot spot that GLAS laser footmark is diameter 70m, corresponding point-to-point mark, namely cannot not find accurately and put position accordingly, therefore need in test to mate the unique point in GLAS laser facula within the scope of hot spot.

According to embodiments of the invention, alternatively, based on Forstner operator extraction unique point, realize automatic extraction and the coupling of image feature point within the scope of GLAS broad sense vertical control point hot spot, the difficulty that artificial visual differentiates can not only be solved like this, also can improve operating efficiency simultaneously.

Forstner operator its essence is a weighted operator.It is the gray scale covariance matrix of a window (as 5 × 5) centered by Robert's gradient by calculating each pixel and pixel (c, r), finds minimum in image, close to the point of the error ellipse of circle as unique point.Because Forstner operator is more complicated, an alternative method is, first extracts just reconnaissance by simple difference operator, and adopt Forstner operator in 3 × 3 window calculation interest value for first reconnaissance, and select candidate point, finally extracting extreme point is unique point.Its computation process is as follows:

1) first image is divided into equal-sized region.

2) difference operator of pixel (c, r) in each block region is calculated.Calculate the grey scale difference absolute value dg1 of pixel (c, r) at four direction up and down, dg2, dg3, dg4:

dg ₁＝|g _c,r-g _c+1,r| (26)

dg ₂＝|g _c,r-g _c,r+1| (27)

dg ₃＝|g _c,r-g _c-1,r| (28)

dg ₄＝|g _c,r-g _c,r-1| (29)

3) determining just reconnaissance and candidate point, by dg1, dg2, dg3, dg4 according to sorting, and extracting just reconnaissance according to given threshold value T from low to high.

M＝mid{dg ₁,dg ₂,dg ₃,dg ₄} (30)

For given threshold value T, if M>T, then (c, r) is the just reconnaissance of this intra-zone one feature; Otherwise (c, r) is not unique point.Then, then in 3x3 window centered by (c, r), calculate the circularity qc of covariance matrix N and error ellipse according to Forstner operator, according to the roundness threshold Tq of error ellipse, r, then judges whether it is a feature candidate point.If restraining error ellipse is long, the ratio of minor semi-axis is not more than 3.2 ~ 2.4, then can try to achieve Tq=0.32.

$N=\left[\begin{array}{cc}\mathrm{\&Sigma;}{g}_u^2& \mathrm{\&Sigma;}{g}_u{g}_v\\ \mathrm{\&Sigma;}{g}_v{g}_u& \mathrm{\&Sigma;}{g}_v^2\end{array}\right]---\left(31\right)$

Wherein g _uand g _vthe partial differential along u and v direction respectively

$\mathrm{\&Sigma;}{g}_u^2=\underset{i=r-k}{\overset{r+k-1}{\mathrm{\&Sigma;}}}\underset{j=c-k}{\overset{c+k-1}{\mathrm{\&Sigma;}}}{(g_{i+1,j+1}-g_{i,j})}^2---\left(32\right)$

$\mathrm{\&Sigma;}{g}_v^2=\underset{i=r-k}{\overset{r+k-1}{\mathrm{\&Sigma;}}}\underset{j=c-k}{\overset{c+k-1}{\mathrm{\&Sigma;}}}{(g_{i+1,j}-g_{i,j-1})}^2---\left(33\right)$

$\mathrm{\&Sigma;}{g}_u{g}_v=\underset{i=r-k}{\overset{r+k-1}{\mathrm{\&Sigma;}}}\underset{j=c-k}{\overset{c+k-1}{\mathrm{\&Sigma;}}}(g_{i+1,j+1}-g_{i,j})(g_{i+1,j}-g_{i,j+1})---\left(34\right)$

$k=\mathrm{INT}\left(\frac{s}{2}\right)---\left(35\right)$

Comprehensively above-mentioned 4 formulas.Can ask

$q_{c,r}=\frac{4\mathrm{DetN}}{{\left(\mathrm{trN}\right)}^2}---\left(36\right)$

Wherein Det is the value of determinant, and tr is matrix trace.Because N is nonsingular matrix, i.e. DetN ≠ 0.If

Q _c,r< T _q=0.32, then this pixel is a feature candidate point.

4) determine that unique point take weights as foundation, the extreme point of the graticule mesh divided in extraction step (1) is as unique point, and the pixel namely extracted in this graticule mesh corresponding to maximum weight person is its unique point, determines that the foundation weighed is as follows:

$w=\left\{\begin{array}{c}0(q_{c,r}<T_q)\\ \frac{\mathrm{DetN}}{\mathrm{trN}}(q_{c,r}>T_q)\end{array}\right.---\left(37\right)$

5) foundation weights are foundation, and choosing the point weighing the maximum in a suitable window is unique point.

Therefore, in the process utilizing Forstner operator extraction unique point, need to determine two threshold values: one is the threshold value T weighing equal value coefficient.Another is the threshold value T of interest value _q.The determination of these two threshold values is on the registration accuracy impact that fixed output quota is raw very large really, and Threshold selection too young pathbreaker increases the calculated amount in unique point selection course greatly, and threshold value arranged senior general increases the probability that unique point leaks choosing, wrong choosing.

The implementation result according to the stereoblock adjustment method of the embodiment of the present invention is introduced in citing below.

As an example, have chosen Ikonos-2 and ZY-3 two kinds of satellites in Australian Hobart area, Taiyuan, Shanxi, weinan altogether trizonal stereopsis carry out experiment.What wherein Hobart area adopted is the mono-scape stereopsis of Ikonos-2, and Taiyuan and Weinan Prefecture have employed the mono-scape of ZY-3 and region many scapes stereopsis respectively.

The Ikonos-2 stereopsis acquisition time in Australia Hobart area is in February, 2003, image resolution is 1m, measured by International Photography and remote sensing association ISRPS website provides free, wherein reference mark data are obtained by field GPS field survey by University of Melbourne, vertical accuracy is better than 0.15m, have chosen 21 reference mark altogether.

Taiyuan region scope of experiment is: 37.5448 °--38.4856 ° of N, 112.3624 °--112.9165 ° of E, experimental data is resource No. three three-line imageries photographed on October 31st, 2012, wherein plane control is data from the DOM image of 0.5m resolution, as the GLAS of high process control data from 1 grade and 14 grades of products in October, 2009, elevation checks that data are the airborne LIDAR data of the 3 meters of graticule mesh in this area, and precision is about 15cm.All data all transfer the utm projection of central meridian 111 degree under WGS84 ellipsoid to.

Weinan Prefecture's scope of experiment is: 33.955 °--35.323 ° of N, 107.862 °--109.278 ° of E, area about 20,000 square kilometres, experimental data is on July 28th, 2012, November 13, the three rail resource three-line imageries photographed November 23, 40 high-precision control points (follow-up as planimetric control point and checkpoint use) adopting GPS field operation to measure acquisition are had in whole survey district, there is on March 31st, 2007, on March 4th, 2008, on March 28th, 2009 totally three GLAS track datas, by 25 kilometers, interval, 50 kilometers, 100 kilometers pick 14 GLAS points altogether as vertical control point.

the three-dimensional adjustment experiment of single scape

For single scape stereopsis, deploy to ensure effective monitoring and control of illegal activities scheme in the corner that have employed in photogrammetric measurement, and in order to verify that plane is separated with high process control, only used as plane control at the reference mark of corner, direction along ng a path selects two as vertical control point simultaneously.The Ikonos-2 image in Hobart area has added an elevation control data at zone line, and the ZY-3 satellite image of Taiyuan region adds cloth plane control data in centre.Adjustment Scheme point position distribution schematic diagram is shown in Fig. 9, and concrete scheme is described as follows:

adjustment Scheme P1:0 planimetric control point, 0 vertical control point;

adjustment Scheme P2:4 angle point cloth 4 planimetric control points, 0 vertical control point;

an adjustment Scheme H1:0 planimetric control point, 3 vertical control points along rail direction;

an adjustment Scheme H2:4 planimetric control point, 1 vertical control point along rail direction;

an adjustment Scheme H3:4 planimetric control point, 2 vertical control points along rail direction;

an adjustment Scheme H4:4 planimetric control point, 3 vertical control points along rail direction.

The Ikonos image in Hobart area is the stereogram product after certain geometric manipulations that ISRPS provides, therefore it can reach plane 3.561m without control positioning precision, elevation 3.612m, it also has certain gap without control vertical accuracy from 1:10000 stereoplotting precision.In experiment, adopt adjustment mode P2, adopt separately 4 planimetric control points plane precision can be brought up to 0.837m, and vertical accuracy is not improved; Within only adopting the adjustment mode H1 of 3 vertical control points vertical accuracy can be brought up to 1m, and plane precision does not improve.On the basis having four planimetric control points, direction along ng a path progressively adds 1,2 and 3 vertical control points, plane precision remains unchanged substantially, and 2 vertical control points are obviously better than only using 1 vertical control point, and the adjustment result of 3 vertical control points then with 2 is suitable.

Test block, Taiyuan height above sea level is from 700 meters to more than 1400 meters, a part belongs to Plain, a part belongs to little high mountain region, this test block possible is close geometry calibration district in-orbit just in time, therefore it is higher without control vertical accuracy as can be seen from Table 1, even mountain area still can reach 5.415m, after the plane control data adding 4 angle points and 1 central point, plane positioning precision significantly promotes, but vertical accuracy has no change; Select the GLAS laser-measured height data after screening as elevation control data, if add a vertical control point, vertical accuracy has certain lifting, but effect is not obvious especially, selecting to add two (selecting two points along rail direction by two) then to promote precision the good results are evident, can 2.932m be reached.

Table 6 Hobart and Taiyuan region list scape stereopsis adjustment result statistical form

the three-dimensional adjustment experiment of regional network

For 9 scape resource No. three satellite images in area, weinan in experimentation, having carried out a kind respectively has the adjustment of control to test without control and 7 kinds.These 8 kinds of adjustment cloth point modes are respectively:

Scheme P1: without the adjustment of Free Networks at reference mark;

Scheme P2: only lay 4 planimetric control points altogether at 4 angle points;

7 planimetric control points are laid altogether along rail direction in scheme P3:4 angle point and centre;

Scheme P4: lay 9 planimetric control points altogether by three row three row;

Scheme H1: only along GLAS trajectory direction equispaced about 25 kilometers altogether lay 14 vertical control points;

Scheme H2: on scheme P4 basis, about 100 kilometers of reconnaissances along GLAS trajectory direction interval, lay 9 planimetric control points and 6 vertical control points altogether;

Scheme H3: on scheme P4 basis, about 50 kilometers of reconnaissances along GLAS trajectory direction interval, lay 9 planimetric control points and 8 vertical control points altogether;

Scheme H4: scheme P4 and the combination of H1, lay 9 planimetric control points and 14 vertical control points altogether.

7 kinds next have the some position distribution of control adjustment as shown in Figure 10, and 8 kinds of adjustment result are in table 7, and the unit of medial error is rice.

Resource No. three area adjustment result statistical forms that table 7.GLAS altitude figures is auxiliary

Test block, Weinan belongs to alpine region substantially, and height above sea level is from 350--1800 rice, and the discrepancy in elevation is up to 1500 meters, and poor without control plane positioning precision, about 13.284m, elevation location precision is 13.532m, and has obvious Systematic Errors.Add 4 angle points as after planimetric control point, plane precision significantly promotes, and vertical accuracy is unchanged, is undertaken only having the mode of plane control to carry out adjustment more afterwards by 7 points, 9 points, and plane precision increases but significantly, vertical accuracy is change no; When not having plane control, adopt 14 GLAS vertical control points, plane precision does not improve, but vertical accuracy is significantly improved, and reaches 3.066m; Based on 9 plane control, adopt again 100 kilometers, interval, 50 kilometers, 25 kilometers add 6,8,14 GLAS vertical control points, plane precision with only adopt the adjustment result of 9 plane control suitable, vertical accuracy improves constantly along with the reduction at GLAS vertical control point interval, when wherein adopting 8 GLAS vertical control points during about 50 kilometers, interval, vertical accuracy can meet the requirement of 1:5 ten thousand stereoplotting, and interval 25 km choose 14 GLAS vertical control point adjustments after vertical accuracy can reach 2.903m.

In above-mentioned experiment, for three regions, the satellite image of two kinds of different resolutions carries out single scape and the three-dimensional adjustment experiment of many scene area net, its result can show: 1) for satellite stereo image, plane control can be separated with high process control, a kind of control data of independent employing and adopt the flat high result controlled in the direction of accuracy improvements and be numerically consistent, namely during the high process control of simple employing, plane precision does not improve, but vertical accuracy is significantly improved and suitable with flat high Hybrid mode result, like this for overseas region, only use spaceborne laser altimeter system data, without any other ground control datas, also the stereoplotting (plane precision relies on the plane precision with satellite image itself) of satellite image can be carried out, simple employing planimetric control point, vertical accuracy does not change substantially, and plane precision is significantly improved, 2) 4 angle points or 4 angle points is adopted to add the planar cloth prosecutor formula of 1 central point for single scape, many scapes adopt survey district four angle point and image overlap district certain intervals to carry out plane and deploy to ensure effective monitoring and control of illegal activities, can reach good plane control effect, plane precision can have greatly improved, 3) adopt the elevation control data of same accuracy level, elevation can be brought up to about 0.8m to Ikonos-2, and ZY-3 can only bring up to 2.9m.Control data of the same type is described, carries out adjustment to the satellite stereo image of different resolution, its result precision level is distinguishing, and with geometric resolution positive correlation, this can be following 1:10000 stereoplotting satellite of launching and provides reference; 4) plane to be separated with high process control carry out adjustment time, single scape adopts lays two vertical control points along rail direction at the two ends of image, the employing of many scene area net is less than every scape overlay length (as ZY-3 gets the interval being less than 50km) along rail by interval and lays multiple row vertical control point, can meet three-dimensional adjustment to elevation domination number object demand.

Satellite photogrammetry and laser-measured height are two important technical that satellite earth observation obtains three-dimensional geographic information, and different imaging modes determines and is better than vertical accuracy without the former plane precision under control condition, and the latter's vertical accuracy is better than plane precision.According to embodiments of the invention, propose a kind of plane and be separated the new method of carrying out adjustment with high process control, set up the part plan constraint of vertical control point and the local elevation constraint of planimetric control point, according to precision, different weights are invested to inhomogeneous control information, effectively can realize the image control strip set up at present and laser-measured height data assemblies carries out the target of resource No. three three-dimensional adjustments as control source, the final stereoplotting production efficiency improving No. three, resource.According to embodiments of the invention, give the distribution scheme of the control data of concrete elevation and planar separation.Such as, corner, region adds center and lays 5 planimetric control points, plane control demand can be met, select spaceborne laser altimeter system data can meet elevation demand for control as vertical control point being less than 50km along interval, rail direction, both combinationally use, and finally can meet stereoplotting demand.

The above is only exemplary embodiment of the present invention, but not for limiting the scope of the invention, protection scope of the present invention is determined by appended claim.

Claims (10)

1., based on a satellite image stereoblock adjustment method for spaceborne laser altimeter system data, it is characterized in that, comprising:

The first control data corresponding to tie point and the second control data is obtained respectively from the first control data source and the second control data source;

Based on remote sensing stereopsis imaging model, three-dimensional adjustment processing is combined to described first control data and described second control data, obtains the parameter correction values of described remote sensing stereopsis imaging model,

Wherein,

Described first control data source is storehouse, image reference mark;

Described second control data source is broad sense vertical control point storehouse;

Described first control data is the control data of reference mark or reference mark image film;

Described second control data is the control data of broad sense vertical control point;

The control data of this broad sense vertical control point comprises spaceborne laser altimeter system data,

Describedly obtain the first control data corresponding to tie point and the second control data respectively from the first control data source and the second control data source, comprising:

Obtain the planimetric coordinates of spaceborne laser altimeter system data;

The planimetric coordinates of described spaceborne laser altimeter system data is mapped to satellite remote sensing digital orthophoto map and gathers footmark image;

The plane control data of a part for described spaceborne laser altimeter system data and described footmark image or described footmark image form described second control data jointly;

The same place of the some in the scope of the described part of described footmark image or described footmark image is obtained as tie point from storehouse, described image reference mark, and using the control data of described same place as described first control data.

2. satellite image stereoblock adjustment method according to claim 1,

It is characterized in that,

Described remote sensing stereopsis imaging model is rational function model or rigorous geometry model,

By expression formula V=AX+Bt-L, the unification in described adjustment processing by described rigorous geometry model and described rational function model,

Wherein,

$V=\left(\begin{array}{c}{v}_c\\ v_r\end{array}\right),$ $B=\left(\begin{array}{cccccc}1& c& r& 0& 0& 0\\ 0& 0& 0& 1& c& r\end{array}\right),L=\left(\begin{array}{c}{l}_c\\ l_r\end{array}\right)$

$t=\left(\begin{array}{c}{\mathrm{da}}_0\\ {\mathrm{da}}_1\\ {\mathrm{da}}_2\\ {\mathrm{db}}_0\\ {\mathrm{db}}_1\\ {\mathrm{db}}_2\end{array}\right),$

(c, r) is picpointed coordinate, for topocentric coordinates, (v _c, v _r) be picture point residual error, for Three Dimensional Ground coordinate corrective value, (a ₀, a ₁, a ₂, b ₀, b ₁, b ₂) correct parameter for the rigorous geometry model sight line vector of image corrects parameter or rational function model image space.

3. satellite image stereoblock adjustment method according to claim 2,

It is characterized in that,

For described first control data, there is the first error equation:

V＝AX+Bt-L

Wherein, $A=\left(\begin{array}{c}\frac{\&PartialD;c}{\&PartialD;H}\\ \frac{\&PartialD;r}{\&PartialD;H}\end{array}\right),$ X=dH, or,

Wherein, a, b and d are the coefficients of plane equation,

For described second control data, then there is the second error equation:

V＝AX+Bt-L

Wherein,

Solve carry out the three-dimensional adjustment processing of described associating by the first error equation described in simultaneous and described second error equation.

4. satellite image stereoblock adjustment method according to claim 3,

It is characterized in that,

Adopt least square method to described first error equation and described second error equation simultaneous solution, comprising:

By solving the normal equation corresponding with described first error equation and described second error equation: Nx=B obtains the model correction parameter of every scape image, wherein, $N=\left(\begin{array}{cc}{A}^T\mathrm{PA}& A^T\mathrm{PB}\\ B^T\mathrm{PA}& B^T{\mathrm{PB}}^T\end{array}\right),x=\left(\begin{array}{c}X\\ t\end{array}\right),$ $B=\left(\begin{array}{c}{A}^T\mathrm{PL}\\ B^T\mathrm{PL}\end{array}\right),$ P is weight matrix.

5. satellite image stereoblock adjustment method according to claim 3,

It is characterized in that,

For described second error equation, the mode adopting pointwise to cut unit forms the change normal equation relevant with the image space correction factor often opening image.

6. satellite image stereoblock adjustment method according to claim 1,

It is characterized in that,

Described spaceborne laser altimeter system data are GLAS spaceborne laser altimeter system data.

7. the satellite image stereoblock adjustment method according to claim 1 or 6,

It is characterized in that,

After obtaining spaceborne laser altimeter system data from described second data source, through just slightly terrain data excluding gross error and echo waveform obtain described second control data being used for described adjustment processing after screening.

8. satellite image stereoblock adjustment method according to claim 7,

It is characterized in that,

Obtain described second control data by described echo waveform screening to comprise:

The background noise level of estimation echo waveform;

Utilize the noise signal of echo waveform described in the filtering of smothing filtering algorithm;

The Gauss curve fitting of filtered signal being carried out to initial parameter is estimated;

Utilize the result that described Gauss curve fitting is estimated, in conjunction with initial parameter, utilize nonlinear least square method to iterate and extract gaussian component parameter accurately.

9. satellite image stereoblock adjustment method according to claim 8,

It is characterized in that,

Adopt carry out described Gauss curve fitting estimation, wherein w (t) be t backward energy value, W _mbe the backward energy value of m Gaussian function in t, Np is the number of Gaussian function, and ε is noise, and tm is the average of m Gaussian function, σ _mbe the mean square deviation of m Gaussian function, A _mthe peak value of m Gaussian function.

10. satellite image stereoblock adjustment method according to claim 6,

It is characterized in that, getting off the imaging model mutually unification of described GLAS laser-measured height data with stereoblock adjustment in coordinate system aspect by GLAS data being transformed into WGS84 coordinate system.