CN104931022B - Satellite image stereoblock adjustment method based on spaceborne laser altimeter system data - Google Patents

Abstract

Satellite image stereoblock adjustment method based on spaceborne laser altimeter system data, including：The first control data and the second control data corresponding to tie point are obtained respectively from the first control data source and the second control data source；The first control data and the second control data are carried out based on remote sensing stereopsis imaging model to combine three-dimensional adjustment processing, obtain the parameter correction values of remote sensing stereopsis imaging model, wherein, obtaining the first control data and the second control data includes：Obtain the plane coordinates of spaceborne laser altimeter system data；The plane coordinates of spaceborne laser altimeter system data is mapped to satellite remote sensing digital orthophoto map and gathers footmark image；Spaceborne laser altimeter system data and the plane control data of footmark image or a part for footmark image collectively form the second control data；The a number of same place obtained from image control point storehouse in the range of the part of footmark image or footmark image is used as the first control data as tie point, and using the control data of same place.

Description

Satellite image stereoblock adjustment method based on spaceborne laser altimeter system data

Technical field

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

Background technology

With expanding economy, three-dimensional geographic information obtains in each side such as digital earth, urban planning, environmental protection To extensive use.And the fast development of remote sensing satellite technology, space photogrammetry is had become after aerophotogrammetry again The means of one quick obtaining three-dimensional geographic information, the particularly development recently as technologies such as three line scanner stereoscopic cameras, three The acquiring technology of dimension remote sensing information achieves significant progress.

For example, resource three is the civilian stereoplotting satellite of first of China high accuracy, it succeeding in sending up and effectively applying Break dependence of the China for a long time to external high-precision satellite image, generate huge social and economic benefit.But by In the reason for the optical stereo satellite stereoplotting mode and Satellite Attitude rail measurement accuracy, camera distortion etc., cause its Under conditions of Pillarless caving, its measurement of higher degree precision is also difficult to meet the needs of surveying and drawing in high precision.

Ground control point GCP (Ground Control Points) is satellite remote sensing image geometry correction and geo-location When important reference data sources.In remotely sensing image geometric corrects processing procedure, to reach certain correction precision, certain amount Ground control point be essential, pass through the object coordinates at control point and corresponding picpointed coordinate and build video imaging model And solving model parameter or to existing imaging model optimize compensation solve compensating parameter, be finally created as during object space With the correct transformational relation of image space.

Control point during traditional work is typically using the pattern of full field measurement, it is necessary to by " collection has controlled Data, survey area is made an on-the-spot survey, reconnaissance is laying of markstone, field operation measure, interior industry arrangement " etc. a series of complex process, although with GPS-RTK etc. The development of advanced e measurement technology, the workload and complexity of field operation measurement greatly reduce, but necessary control point field measurement work Make still inevitable.

Further, since satellite remote-sensing image coverage it is big (by taking No. three satellites of resource as an example, coverage is 50 kilometers × 50 kilometers), want to obtain equally distributed ground control point, it usually needs in hundreds of kilometer or even thousands of square kilometres of scope Interior testing, its field process amount is big, labor intensive and material resources are self-evident.In addition, in earthquake, flood, mudstone The uninhabited area such as natural calamity hotspots or virgin forest, marsh, desert, the survey crews such as stream usually cannot be introduced into reality Ground measures.

Even if control point is obtained by manually choosing remote sensing image and topographic map same place, there is also efficiency it is low, The problems such as precision is difficult to ensure.And different phase or the different sensors images of areal are corrected, weight occurs again The situation of final election point.

Ground control point is generally divided into：Flat high control point, horizontal control point and vertical control point.Chinese patent application 201310143369.3 disclose a kind of multi-source heterogeneous remote sensing image control point automatic acquiring method, can be from multi-source heterogeneous image Remote sensing control imaging point (control point image film) is automatically extracted, improves the efficiency and precision of control data acquisition, but its essence On be still a kind of automatic acquiring method of horizontal control point.

In the case of ground control point data deficiencies, the block adjustment technology of satellite image can be used as precise geometrical A kind of important means of positioning.For example, document《No. three cartographic satellite image planes of resource and stereoblock adjustment compare》 (《Survey and draw journal》In April, 2014, the 4th phase of volume 43), can not for satellite remote-sensing image block adjustment under weak encounter conditions Correct the problem of solving, propose what is constrained using digital elevation model (Digital Elevation Mode, DEM) as elevation Plane domain net adjusted data method improves its object positioning precision.But because this method is only for image plane region Net adjusted data, therefore can not solve the problems, such as No. three cartographic satellite measurement of higher degree precision deficiencies of resource.

The content of the invention

At least one in order to solve the above-mentioned technical problem, according to an embodiment of the invention, there is provided one kind is swashed based on spaceborne The satellite image stereoblock adjustment method of the high data of flash ranging, it includes：From the first control data source and the second control data Source obtains the first control data and the second control data corresponding to tie point respectively；Based on remote sensing stereopsis imaging model pair First control data and the second control data carry out combining three-dimensional adjustment processing, obtain the parameter of remote sensing stereopsis imaging model Correction value, wherein, the first control data source is image control point storehouse；Second control data source is broad sense vertical control point storehouse；The One control data is the control data of control point or control point image film；Second control data is the control of broad sense vertical control point Data processed；The control data of the broad sense vertical control point includes spaceborne laser altimeter system data, from the first control data source and second Control data source obtains the first control data and the second control data corresponding to tie point respectively, including：Obtain satellite borne laser Survey the plane coordinates of high data；The plane coordinates of spaceborne laser altimeter system data is mapped to satellite remote sensing digital orthophoto map simultaneously Gather footmark image；The plane control data of spaceborne laser altimeter system data and footmark image or a part for footmark image are common Form the second control data；Obtained from image control point storehouse certain in the range of the part of footmark image or footmark image The same place of quantity is used as the first control data as tie point, and using the control data of same place.

Satellite image stereoblock adjustment method according to embodiments of the present invention, alternatively, the imaging of remote sensing stereopsis Model is rational function model or rigorous geometry model, by expression formula V=AX+Bt-L, by rigorous geometry model and reasonable letter Exponential model is unified in adjustment processing, wherein,

(c, r) is picpointed coordinate,For topocentric coordinates, (v_c,v_r) it is picture point residual error,For Three Dimensional Ground coordinate corrective value, (a₀,a₁,a₂,b₀,b₁,b₂) correct and parameter or have for the rigorous geometry model sight line vector of image Manage function model image space correction parameter.

Satellite image stereoblock adjustment method according to embodiments of the present invention, alternatively, for the first control data, There is the first error equation：

V=AX+Bt-L

Wherein,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,

By the error equation of simultaneous first and the second error equation and solve to carry out combining three-dimensional adjustment processing.

Satellite image stereoblock adjustment method according to embodiments of the present invention, alternatively, using least square method To the first error equation and the second error equation simultaneous solution, including：By solving and the first error equation and the second error side Normal equation corresponding to journey：Nx=B obtains the model correction parameter of every scape image, wherein,P is weight matrix.

Satellite image stereoblock adjustment method according to embodiments of the present invention, alternatively, for the second error equation, The change normal equation relevant with the image space correction factor of every image is formed by the way of member is cut point by point.

Satellite image stereoblock adjustment method according to embodiments of the present invention, alternatively, spaceborne laser altimeter system data It is GLAS spaceborne laser altimeter system data.

Satellite image stereoblock adjustment method according to embodiments of the present invention, alternatively, from the second control data source After obtaining spaceborne laser altimeter system data, obtained after just omiting terrain data excluding gross error and echo waveform screening at adjustment Second control data of reason.

Satellite image stereoblock adjustment method according to embodiments of the present invention, alternatively, is screened by echo waveform Obtaining the second control data includes：Estimate the background noise level of echo waveform；Echo waveform is filtered out using smothing filtering algorithm Noise signal；The Gauss curve fitting that initial parameter is carried out to filtered signal is estimated；The result estimated using Gauss curve fitting, knot Initial parameter is closed, is iterated using nonlinear least square method and extracts accurate Gaussian component parameter.

Satellite image stereoblock adjustment method according to embodiments of the present invention, alternatively, useCarry out Gauss curve fitting estimation, wherein w (t) be t backward energy value, W_mFor m Individual Gaussian function is in the backward energy value of t, and Np is the number of Gaussian function, and ε is noise, t_mFor m-th Gaussian function Average, σ_mFor the mean square deviation of m-th of Gaussian function, A_mThe peak value of m-th of Gaussian function.

Satellite image stereoblock adjustment method according to embodiments of the present invention, alternatively, by the way that GLAS data are turned WGS84 coordinate systems are changed to get off GLAS laser-measured heights data and the imaging model phase of stereoblock adjustment in coordinate system aspect It is unified.

Satellite image stereoblock adjustment method based on spaceborne laser altimeter system data according to embodiments of the present invention, It is improved on the basis of traditional block adjustment method, the method separated using plane and elevation control data, improves adjustment Algorithm is proposed unified general adjustment model, overcome using spaceborne laser altimeter system data to the compatibility of all kinds of control datas Current Precision Elevation control data gathers the problem of difficult, improves resource three and the stereoplotting precision of subsequent satellites and place Efficiency is managed, and to improve resource three in the area in hardship and the external stereoplotting precision that can not obtain control data area A kind of feasible solution route is provided.

Brief description of the drawings

In order to illustrate the technical solution of the embodiments of the present invention more clearly, the accompanying drawing of embodiment will be simply situated between below Continue, it should be apparent that, drawings in the following description merely relate to some embodiments of the present invention, rather than limitation of the present invention.

Fig. 1 is that the satellite image solid region net according to an embodiment of the invention based on spaceborne laser altimeter system data is put down The indicative flowchart of difference method；

Fig. 2 shows the mathematical relationship of the picture point and ground point under rigorous geometry model；

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

Fig. 4 is the schematic diagram intersected 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 show actual landform and echo waveform corresponding to GLAS footmarks, and wherein Fig. 6 A are that footmark is corresponding Actual image；Fig. 6 B are actual landforms corresponding to footmark；Fig. 6 C are GLAS actual ghosts waveforms；

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

Fig. 8 is that GLAS surveys high data footprints image (left side) and the schematic diagram of echo waveform (right side) in database；Fig. 9 is single scape Stereopsis adjustment Scheme point bit distribution schematic diagram；

Figure 10 is regional network solid adjustment Scheme point bit distribution schematic diagram.

Embodiment

To make the purpose, technical scheme and advantage of the embodiment of the present invention clearer, below in conjunction with the embodiment of the present invention Accompanying drawing, the technical scheme of the embodiment of the present invention is clearly and completely described.Obviously, described embodiment is this hair Bright part of the embodiment, rather than whole embodiments.Based on described embodiments of the invention, ordinary skill The every other embodiment that personnel are obtained on the premise of without creative work, belongs to the scope of protection of the invention.

Unless otherwise defined, technical term or scientific terminology used herein should be in art of the present invention and had The ordinary meaning that the personage of general technical ability is understood.Used in present patent application specification and claims " the One ", " second " and similar word are not offered as any order, quantity or importance, and are used only to distinguish different Part.Equally, the similar word such as "one" or " one " does not indicate that quantity limits yet, but represents to exist at least one.

Because certain error all be present in the track of optical satellite image, posture and elements of interior orientation etc., therefore by strict There is some difference with the meeting of actual coordinate value for the Three Dimensional Ground coordinate that imaging model calculates.In order to eliminate this species diversity, improve The geometric positioning accuracy of satellite is general to require to introduce a number of ground control point.In order to improve the three-dimensional localization of satellite essence Degree, traditional mode as control point, i.e., are measured in field GPS measurement points or high-precision topographic map and obtained typically using flat high point The three-dimensional coordinate point obtained.

The control data for putting down high control point generally includes panel data and altitude data.Wherein, panel data can include The longitude and latitude of the ground control point；Altitude data is distance of the ground control point along plumb line direction to absolute datum, In survey field, usually height of the ground control point on terrestrial gravitation direction, that is, height above sea level.

For traditional flat high control point, because longitude, latitude and the combination of elevation and the control point correspond, The control data for including panel data and altitude data is also uniquely corresponding to specifically put down high control point.

And according to an embodiment of the invention, it is proposed that a kind of three-dimensional error compensation method for separating plane and high process control.Tool For body, when carrying out satellite image stereoblock adjustment, used control data is no longer limited to and a certain control number According to panel data and altitude data corresponding to the flat high control point in source, but can use from different control data sources The combination of control data.The control point data from different control data sources can have higher plane precision respectively Higher height accuracy, by carrying out local elevation constraint and to tool to the control point data with higher level precision respectively The control point data for having higher height accuracy carries out part plan constraint to carry out simultaneous adjustment processing, can improve control point control The plane precision and/or height accuracy of data processed.

Further, the present invention proposes a kind of concept of new broad sense vertical control point.For the high process control of this broad sense For point, it can both include traditional vertical control point or a planar region with approximate height value, wherein, The topographic map mapping of different scale has different required precisions to the hypsography or the discrepancy in elevation in this planar region.Then, its Altitude data can be the altitude data of specific vertical control point or by laser-measured height satellite such as ICESat/GLAS The elevation control data of the various broad sense such as the altitude data after waveform screens.

On the other hand, satellite image stereoblock adjustment method according to embodiments of the present invention is by tight imaging model Sight line vector correcting method and rational function model in image space correcting method unified, form unified solid and put down Differential mode type.

Fig. 1 is that the satellite image solid region net according to an embodiment of the invention based on spaceborne laser altimeter system data is put down The indicative flowchart of difference method.

As shown in figure 1, in step S100, the control number of tie point (same place) is obtained respectively from different control data sources According to.

Specifically, control point is obtained from image control point storehouse by image Auto-matching or control point image film is used as Tie point, and the control data (including plane control data and elevation control data) of the tie point is obtained as the first control number According to wherein image control point storehouse can the multi-source heterogeneous remote sensing image control according to disclosed in Chinese patent application 201310143369.3 System point automatic acquiring method is established based on digital orthophoto map (Digital Orthophoto Map, DOM).From broad sense elevation Control point storehouse, such as spaceborne laser altimeter system database, acquisition control data corresponding with the tie point (including plane control data With elevation control data) it is used as the second control data.Alternatively, according to an embodiment of the invention, can also be first from broad sense elevation Control point storehouse obtains the second control data, then obtains that broad sense vertical control point is corresponding controls number with this from image control point storehouse According to as the first control data.

Next, in step S300, number is controlled to the first control data and second based on remote sensing stereopsis imaging model According to carrying out combining three-dimensional adjustment processing, the parameter correction values of imaging model are obtained.

The imaging model of remote sensing image can be divided into two major classes：Physical model and universal model.

Physical model refers to consider physical significance such as surface relief, earth curvature, the air for causing deformation of image during imaging The factors such as refraction, camera distortion, then it is construed as using these physical conditions as geometrical model, the most representational is conllinear bar Rigorous geometry model based on part equation.

Universal model does not consider imaging mechanism, but the geometrical relationship of picture point and object point is directly described with mathematical function, It has the characteristics of generality, confidentiality, high efficiency.Conventional universal imaging model has multinomial, direct linear transformation, affine Conversion, rational function model (Rational Function Model, abbreviation RFM) etc..

Plane according to embodiments of the present invention and the stereoblock adjustment method of high process control separation can uniformly be applied to Rational function model and rigorous geometry model, illustrate adjustment side according to embodiments of the present invention based on two kinds of models individually below Method.

RFM models are that picpointed coordinate (c, r) is expressed as with corresponding ground point space coordinates(wherein For latitude, λ is longitude, and H is elevation) be independent variable polynomial ratio：

Wherein,For the geographical coordinates of standardization, (c_n,r_n) it is the image coordinate standardized.There is control point In the case of, it is general that rational function model is repaiied by the way of the compensation of the image spaces such as translation, translation+rotation, affine transformation Just, i.e.,

Wherein (Δ c, Δ r) represent image space compensation numerical value, according to translation and an affine transformation, its expression-form difference For (3) and (4),

First order is expanded to by the compensating form of (4) formula affine transformation to formula (2) by Taylor formula：

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,

For rule plane,Can simultaneous equations：

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

Arranging into error equation is：

Being write as matrix form is：

V=AX-L

Rigorous geometry model describes during imaging picture point and object space point strict geometric process correspondingly, with conventional Collinearity condition equation is representative.According to satellite almanac data, attitude data and time data, by imaging process, picture is established Point (l, p) and ground point (X, Y, Z) orBetween mathematical relationship.Fig. 2 shows that the mathematics of picture point and ground point closes System.

Alternatively, high-precision satellite almanac data, attitude data and time data are provided by SPOT5 satellite images.

Coordinate system used in rigorous geometry model may include following coordinate system, as shown in Figure 3：

Navigational coordinate system (also referred to as body coordinate system) O₁-X₁Y₁Z₁：It is to be determined by the attitude control system (ACS) of satellite, Wherein origin is in centroid of satellite, X1 axles, Y₁Axle, Z₁Axle takes three principal moments axles of satellite respectively.X₁Axle is along satellite transverse axis, Y₁ Axle points to satellite flight direction (not exclusively overlapping), Z along the longitudinal axis₁Axle determines according to right-handed system rule.Satellite attitude measurement exists Carried out in body coordinate system, describing three parameters of its spatial attitude is：Pitching (pitch) is around body coordinate system X₁The rotation of axle Turn, it is around body coordinate system Y to roll (roll)₁The rotation of axle, course deviation (yaw) are around body coordinate system Z₁The rotation of axle.

(ψ_x)_p(ψ_y)_pDirection along ng a path and vertical when representing the imaging of the P component of charge coupling device (CCD) respectively The inclination angle of orbital direction.

Orbital coordinate system：The interim coordinate system that to be navigational coordinate system converted with earth axes, orbit coordinate The origin of system is in centroid of satellite, Z₂It is reverse that axle points to the earth's core, Y₂Axle points to the direction of satellite motion, X on satellite orbit face₂Axle Determined according to right-handed system rule.

O₂-X₂Y₂Z₂The direction of each axle is by satellite in the position of tAnd speedTwo direction vectors determinations, It is as follows：

Earth axes：The position of SPOT satellites and speed are all defined in earth Protocol Ref coordinate system (ITRS), ITRF is ITRS realization, is responsible for establishing and is safeguarded by the Terrestrial Reference Frame portion of IERS central offices, the reference frame Station coordinates and speed are the space geodesy technique such as to defend using very long baseline interferometry(VLBI (VLBI), Laser Measuring to be gathered Data calculate, be a continually changing reference frame, at present using ITRF90, calculate for convenience, one As using WGS84 coordinate systems come instead of ITRF frameworks determine coordinate system.Because WGS84 is ECEF coordinate system, ITRF is then inertial coodinate system, and final mapping is it is desirable that body-fixed coordinate system, and ITRS is differed with WGS84 only at 1 meter or so, from ITRF90 to WGS84 conversion parameter 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

ITRF90 to WGS84 conversion parameter table

Conversion formula is as follows：

WGS84 coordinate systems are one by global the earth's core coordinate frame of reference and one group of corresponding model and WGS84 the earth The measurement reference system that level surface is formed.The origin of the coordinate system is located at earth centroid, Z axis and IERS identical, the X with reference to pole sensing Axle points to IERS guide meridians with crossing origin and being determined perpendicular to the intersection point of Z axis plane, Y-axis by right-handed system rule；

It is as follows that imaging process derives (l, p) → (X, Y, Z)：

Imaging time corresponding to I, calculating l row p row pixels

T=t_c+lsp×(l-l_c)

Wherein, t_cFor the photography time of central scan row, lsp is often capable sweep time, l_cFor the line number of central scan row

II, interpolation go out the position of t satelliteWith

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

Imaging direction vector under IV, calculating orbital coordinate system

Wherein,WithIt is respectively the position of pixel (l, p) imaging moment t satellites and when velocity interpolation goes out t Carve the attitude angle of satellite：

a_p(t_i), a_r(t_i), a_y(t_i) it is respectively t_iThe attitude angle of moment satellite, it can be obtained from auxiliary data file.

By the imaging direction vector under navigational coordinate systemThe imaging direction vector being converted under orbital coordinate system

Imaging direction vector under V, calculating earth axes

Wherein, X₂,Y₂,Z₂For orbital coordinate system O-X₂Y₂Z₂Three change in coordinate axis direction vectors

The positioning of VI, ground point：

As shown in figure 4, measurement vector u₃Intersection point with ground is M, and the position vector of satellite isM is connected pair with the earth's core The vector answered isIt is again a+h in major axis because of M points, short axle is on b+h ellipsoid, so having after simultaneous equations：

Wherein,

With reference to several transformation relations above, write as the form similar to collinearity condition equation：

Obtain the relation of the observed direction of the imaging direction vector and earth axes under navigational coordinate system；Projection centre is (X_P,Y_P,Z_P), coordinate vector is (X, Y, Z), projection coefficient u, and spin matrix is:

Collinearity equation can be obtained：

I.e.

For stereopsis, two equations above can be listed on each image, if weighing twice It is folded, then four equations can be listed, then solve (X, Y, Z) three unknown numbers, obtain Three Dimensional Ground coordinate corresponding to image.

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

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

Therefore formula (12), can be written as：

It is abbreviated as：

Wherein,

More than for rigorous geometry model coordinate system analyze on the basis of, according to an embodiment of the invention, in order to it is preceding State the adjustment analysis based on rational function model mutually to unify, substitute (l, p) with picpointed coordinate (c, r), use topocentric coordinatesSubstitution (X, Y, Z) is stated.

Same image uses same sight line vector corrected parameter, then for stereoblock adjustment when, can be to (15) formula Listing error equation is:

In above formula, (v₁,v₂) it is residual error,For n-th point of Three Dimensional Ground coordinate corrective value, (a_0m, a_1m,a_2m,b_0m,b_1m,b_2m) for the sight line vector of m images correct parameter, k_mnFor n-th point of satellite corresponding with m images The distance between imaging center is worth.

(16) formula of contrast and foregoing (7) formula, remove lower target factor, it can be seen that two expression formulas are basically identical, are all Error equation on Three Dimensional Ground coordinate corrective value and imaging model correction.

Thus it is possible to by aforementioned expression (8), V=AX+Bt-L, rigorous geometry model and rational function model are existed United in adjustment processing.

Wherein

For i-th point on jth image, then corresponding subscript is taken respectively.For example, for i-th point in jth On image, then V_ij=A_ijx_i+B_ijt_j-L_ij, weight matrix P_ij。

Next, on the basis of error equation derived from above-mentioned imaging model, respectively for the first control data and Second control data drafts error equation.

, can be by high process control for the second control data as previously described from broad sense vertical control point data source Height value substitutes into above-mentioned expression formula (8) as actual value, then has

Or

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

Or

For the first control data, if control point image film, then using the plane equation of aforementioned expression (9), obtain To error equation (identical with previously described formula (11))：

V=Ax-L

Then, adjustment resolving is carried out using least square method.

For example, plane control and high process control object coordinates correction are similar into 0, virtual observation equation can be listed：

Cx=0 (21)

Wherein, C=(1 1 1)

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

NX=Q (22)

Wherein,

Above-mentioned normal equation is typical solving sparse linear equations., can be for above formula (18), using cutting member point by point for ease of calculating Mode form the change normal equation relevant with the image space correction factor of every image, for example, the exponent number of changing normal equation is 6m* 6m (wherein m is the image number for participating in adjustment), now changing normal equation is still sparse vectors.

Alternatively, certain weights P can also be assigned according to the measurement accuracy at control point.

The precision of picpointed coordinate observation：σ_xMillimeter, then weigh P_x=1；

The precision of ground control point obserred coordinate value is：σ_cRice, then weighWhereinFor photographic scale.Base This rule is that precision is higher, i.e. σ_cSmaller, weights are bigger, and the contribution to adjustment result is also bigger.

Finally solve the model correction parameter that sparse vectors obtain every scape image.For example, pass through the above-mentioned equation of simultaneous (18) solved with equation (19) and/or equation (20), obtain the model correction parameter of elevation coordinate and plane coordinates.Before It is formula to state equation (18), (19) and (20), can be represented for specific coordinate points with the mode of such as formula (17), to establish The corresponding relation of image of the coordinate points that need to be corrected with participating in adjustment, changes so as to obtain the model of every width image by equation solution Positive parameter.

It is presented above the adjustment processing of the control data of joint first according to embodiments of the present invention and the second control data Process.By the process, the correction to elevation coordinate and plane coordinates can be carried out simultaneously, while improves elevation coordinate peace The precision of areal coordinate.

It is alternatively possible to the correction of elevation coordinate is carried out only for the first control data, or only for the second control number According to the correction for carrying out its plane coordinates.That is, one in above-mentioned equation (18), (19) and (20) is resolved respectively When, only use control data corresponding to the equation.So, the correction of elevation coordinate or plane coordinates can be carried out respectively, respectively Improve precision.As it was previously stated, method according to embodiments of the present invention, the essence of elevation coordinate and plane coordinates can be improved simultaneously Degree.

As described above, satellite image stereoblock adjustment method according to embodiments of the present invention, it is combined from difference The control data in control data source, joint carry out three-dimensional adjustment resolving, can both carry out elevation coordinate or plane coordinates respectively Correction, the correction of elevation coordinate and plane coordinates can also be carried out simultaneously, this namely the present invention alleged by plane and height it is program control The processing of data separating processed.

According to an embodiment of the invention, spaceborne laser altimeter system data can be used as broad sense vertical control point database Control data.

Spaceborne laser altimeter system (Laser Altimeter) is the important technology that satellite earth observation obtains three-dimensional geographic information One of means.Spaceborne laser altimeter system is a kind of ground point High head ratio technology, it using satellite as platform, carry laser ceilometer from Space, time are observed to the earth, high accuracy, the distance between measure satellite to testee, and pass through data in real time Processing and analysis, obtain the information such as the landforms of earth surface, vegetative coverage situation, ocean surface topography.External spaceborne laser altimeter system skill Art is quickly grown, and world's Main Developed Countries are all carrying out the research of spaceborne laser altimeter system instrument.

For example, the U.S. transmitted ICESat satellites in 2003, the geoscience laser-measured height instrument system carried thereon (Geoscience Laser Altimeter System, GLAS) is that first, the whole world is used for the spaceborne sharp of the Continuous Observation earth Ligh-ranging system, is measurement polar ice sheet elevation and its change the purpose of its main science, distribution characteristics of cloud layer and aerosol etc.. ICESat satellites in orbit during, obtain a large amount of high-precision altitude datas, its laser footmark plane precision reaches 10m amounts Level, height accuracy about 15cm.ICESat-2 will be launched within 2017, the laser ceilometer of the satellite is by using 6 beams of 3*2 structures Laser, spot size will also narrow down to 10 meters, and the quantity and precision of the laser spots of acquisition greatly improve.Carried by U.S. NRC The LIST plans gone out will then launch in the year two thousand twenty or so, and the satellite will be entered using the photon detection array of 1000 beams, spot diameter One step narrows down to 5 meters, can quickly carry out the topographic survey of global seismic 5m spatial resolutions and 10cm vertical precisions.

Spaceborne laser altimeter system research is attached great importance in recent years in the country." goddess in the moon one for the moon exploration program series of satellites of independent research Number " and " Chang'e-2 " satellite be equipped with high-precision laser altitude gage, obtain moonscape 3 D stereoscopic image.Except Carry outside laser altimeter on " goddess in the moon " series of satellites, also without the laser-measured height system for earth observation, but be included in Follow-up earth observation laser-measured height satellite development plan, such as the follow-up star task of No. three high-resolution stereoplotting satellites of resource, Satellite surveying and mapping " 12 " planning also claims to the work of laser-measured height satellite technology parameter demonstration.

According to an embodiment of the invention, can be using GLAS spaceborne laser altimeter system data as broad sense vertical control point number According to the control data source in storehouse.

ICESat satellite orbital altitudes about 600km, orbit inclination angle are 94 °, latitude scope that can be ± 86 ° covering the whole world.System System working frequency is 40Hz, pulse width 4ns, each second while launch wavelength infrared laser and 532nm for 1064nm Green laser.The laser footmark that laser pulse is formed at the earth's surface is diameter 70m hot spot, between the serial hot spot of same track Every about 170m, adjacent orbit under the line on interval about 15km, 80 ° of high latitude areas at intervals of 2.5km.GLAS is used High-precision waveform collection equipment, efficient Echo Processing algorithm and precision satellite fix and determine appearance equipment, make it System has very high precision, and wherein plane precision is up to 10m magnitudes, height accuracy about 15cm.ICESat satellite data products By Ge Dade space flight centers I-SIPS (ICESAT Science Investigator-led Processing System) Generation, the echo information received according to GLAS extract every physical parameter, can obtain the other normal data product of 15 species, Naming method is GLA01-GLA15.System data file is by American National ice and snow data center (National Snow and Ice Data Center, Distributed Active Archive Center, NASA NSIDC DAAC) externally announce, point For L1 levels and the DBMS of L2 levels two, wherein L1 DBMSs are divided into two ranked datas of L1A and L1B again.L1A DBMSs are without number According to processing, the full resolution initial data based on time series, including GLA01-GLA04 are recorded.L1B DBMSs are by L1A DBMS corrects by device, data processing, and carrying out laser footpoint based on posture and position data positions what is formed, including GLA05-GLA07.L2 DBMSs are that L1B DBMSs are more accurate afterwards by various correcting process, and are divided into differently The elevation and position data of table type, and for distance, the various amendment data of elevation correction, including GLA08-GLA15.Its Middle GLA12 data files are ice sheet altitude data, and GLA13 is sea ice altitude data, and GLA14 is land altitude data, and GLA15 is Ocean altitude data.

As it was previously stated, the imaging mode of three-dimensional intersection cause satellite image without control precision particularly height accuracy from mapping Using there is gap, it is necessary to improve its precision by ground control point data.Spaceborne laser altimeter system is high without control height accuracy, about 0.15m, and satellite can obtain global altitude data, can be used as the elevation control data of satellite image, The work at artificial surface deployment measurement control point is reduced, and High Precision Stereo mapping regional overseas can be realized.

In the step S100 of foregoing three-dimensional adjustment processing, the second control data is obtained from broad sense vertical control point storehouse, so Control data corresponding with broad sense vertical control point is obtained as the first control data from image control point storehouse afterwards.If using GLAS spaceborne laser altimeter system data, because its corresponding ground footmark is relatively large in diameter, about 70m, thus plane precision is not special It is not high, about 10m, therefore the plane coordinates for being primarily based on spaceborne laser altimeter system data is mapped to known DOM striographs and gathers foot Image is printed, the footmark image makes reference use for carrying same place in image control point data base, it may not be necessary to by footmark image The image film in point data base is controlled to match with image, it is only necessary to be controlled from image in point data base and be chosen at footmark image model Interior same place is enclosed, also therefore, the known DOM striographs precision need not be very high, 20 meters or so cans of precision.Generally For including, be exactly first plane coordinates based on spaceborne laser altimeter system data obtain a footmark image (plane information) substantially, Then a number of same place in footmark image capturing range is obtained from image control point storehouse, and with the control number of the same place According to as the first control data.

Preferably, first pass around and just omit terrain data (such as SRTM) excluding gross error, waveform screens step to satellite borne laser Survey high data to be selected, can so improve the precision of subsequent treatment.

The spaceborne laser altimeter system shape information being exemplified below with reference to ICESat satellites is entered to the features of terrain in footmark Row analysis, filters out available data as elevation control data.

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

Because laser firing pulses are being Gaussian Profile along the direction of the launch and Vertical Launch direction, and terrain and its features is to returning Wave action is then mainly in position etc. corresponding to echo-peak, broadening, crest, so still can use Gauss for echo waveform Function is fitted analysis, only for the landform of complexity, it may be necessary to which multiple Gaussian functions, which are combined, can be only achieved most preferably Fitting effect.

Pass through multiple processes such as propagation in atmosphere, target reflection, echo reception opto-electronic conversion additionally, due to laser firing pulses It is mixed into much noise, it is therefore desirable to estimate by noise, smothing filtering, Gaussian parameter pre-estimation and fitting extracting parameter.

According to the characteristic distributions of echo-signal and the property of Gaussian function, determining the processing method of echo-signal includes：Estimate Calculate background noise level；Most noise signals are filtered out using smothing filtering algorithm；Based on the property of Gaussian function to filtering Signal afterwards carries out the estimation of initial parameter；With reference to initial parameter, extraction is iterated accurately using nonlinear least square method Gaussian component parameter.

For the echo waveform of reality, can useIt is fitted, wherein w (t) For the backward energy value of t, W_mFor m-th of Gaussian function t backward energy value, Np be Gaussian function number, ε For noise, t_mFor the average of m-th of Gaussian function, σ_mFor the mean square deviation of m-th of Gaussian function, A_mThe peak of m-th of Gaussian function Value.By the method for least square can the parameter of each Gaussian function obtained of solution, pair between quantitative analysis landform and waveform It should be related to.Alternatively, can be by calculating the skewness and kurtosis value of echo waveform, to returning before Gauss curve fitting estimation is carried out The waveform gaussian shape of ripple is estimated, in order to which follow-up Gaussian parameter calculates.

The Skewness degrees of bias are the measurements of statistics distribution skew direction and degree, are that statistics is distributed asymmetric journey The numerical characteristic of degree.For the degree of bias of n sample value, common degree of bias estimation formulas：

Wherein：N is sample number, and σ is sample standard deviation,For sample average.

Under general scenario, when statistics is that right avertence is distributed, Sk>0, and Sk values are bigger, right avertence degree is higher；Work as system Count when being distributed for left avertence, Sk<0, and Sk values are smaller, left avertence degree is higher；When statistics is symmetrical, it is clear that have Sk=0.

Kurtosis kurtosis is the index for reflecting the sharp high and steep or flat degree in curve of frequency distribution top, for measuring number According in center aggregation extent.Sometimes the arithmetic average of two groups of data, standard deviation and the coefficient of skew are all identical, but they are distributed song The towering degree on line top is different.Its estimation formulas is：

In the case of normal distribution, coefficient of kurtosis value is 0.Positive coefficient of kurtosis illustrates that observed quantity is more concentrated, and has and compares normal state It is distributed longer afterbody；Negative coefficient of kurtosis illustrates that observed quantity is less concentrated, and has the afterbody more shorter than normal distribution, is similar to Rectangle is uniformly distributed.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 For examining normality.If the ratio absolute value is more than 2, normality will be refused.The coefficient of skewness (Skewness) is used for measuring Whether distribution is symmetrical.Normal distribution or so is symmetrical, the coefficient of skewness 0.It is larger on the occasion of show the distribution have right side compared with Long afterbody.Larger negative value shows there is the longer afterbody in left side.The ratio of the coefficient of skewness and its standard error may also used to examine Normality.

According to the principle of laser ranging, it is 15cm that 1ns intervals, which correspond to actual range, on the time, therefore is passed through after waveform decomposition Time corresponding to crest can calculate corresponding discrepancy in elevation information, the optical image data with reference to corresponding to the footmark, can be achieved The decomposition that becomes more meticulous is carried out to the elevation in footmark.Fig. 6 A~6C show actual landform corresponding to GLAS footmarks and echo ripple Shape, wherein Fig. 6 A are actual images corresponding to footmark；Fig. 6 B are actual landforms corresponding to footmark；Fig. 6 C are GLAS actual ghosts ripples Shape.In figure 6 c, axis of abscissas is the relative time for receiving echo, and unit nanosecond (ns), axis of ordinates is echo strength, unit Volt (volts).

Alternatively, using the echo analysis model to become more meticulous, the landform contained in return laser beam waveform can further be excavated Information, its illustrative process are as shown in Figure 7.

Further, it is necessary in coordinate system aspect by GLAS laser-measured heights data with being imaged used by three-dimensional local net adjusted data Unified Model gets up.For example, GLAS data are transformed under WGS84 coordinate systems.

Specifically, because the reference ellipsoid of GLAS laser-measured height data is TOPEX/Poseidon, it is necessary to first GLAS Data are transformed under WGS-84 reference ellipsoids, are easy to subsequently carry out precision test and application.Because two reference ellipsoids are in parallel Direction difference is smaller, warp direction without correction, therefore when carrying out data conversion can only consider height value amendment.Typically adopt Converted with empirical equation, the calculating of data comparison between two different ellipsoids can be greatly simplified, such as formula (25).

In formula (25), B represents latitude, and what da was represented is the difference of WGS-84 ellipsoids and TOPEX/Poseidon ellipsoid major axis, Db represents the difference of WGS-84 ellipsoids and TOPEX/Poseidon ellipsoid short axles, and dh is elevation correction value, and Hg is relative before conversion The elevation of TOPEX ellipsoids, H_GFor the GLAS data elevations after conversion.

GLAS laser beams can not pass through cirrus and Bao Yun in communication process, in addition by mountain region, hills etc. intricately The influence of shape, part GLAS laser-measured height data values may be caused to be significantly larger than real top elevation.Alternatively, for This part rough error value, is rejected by being contrasted with dem data first.For example, using SRTM data and the high numbers of passes of GLAS According to being contrasted.SRTM data are using WGS84 as reference ellipsoid, using the positive height of EGM96 geoid.For the ease of right Than SRTM is transformed into EGM2008 from EGM96.It is 20m to set elimination of rough difference threshold value, i.e., by GLAS altitude datas and SRTM DEM depth displacement is less than -20m or the point more than 20m removes as rough error point.By the depth displacement of two kinds of data is counted Analysis, discovery meet normal distribution, and the standard deviation for further setting ± 2 times only retains the GLAS elevations within threshold value as threshold value Data.Scalping is carried out using rough terrain data such as SRTM to GLAS laser-measured heights data before waveform screening is carried out to pick Remove, the precision of waveform subsequent screening can be improved.

Available ICESat/GLAS spaceborne laser altimeter systems data are obtained by waveform screening, then according to its substantially flat Areal coordinate finds image blocks from known DOM figures, surveys high data and collectively forms broad sense elevation with image blocks plane control data Control the control data of point data base.For existing GLAS spaceborne laser altimeter systems data, the plane of 70 rice diameter hot spots About 10 meters of error, with reference to the image blocks obtained from DOM figures the error of itself, composition error is at 20 meters or so.According to of the invention real The method for applying example, the plane coordinates of the image blocks is without accurate, and it is only as the reference for choosing tie point (same place).

As it was previously stated, the plane coordinates of broad sense vertical control point is inaccurate for elevation, but can be with base In laser spots geometric coordinate corresponding with broad sense vertical control point and footmark spot size, it is aided with footmark image, satellite can be predicted The same place initial position of image, and pass through the same place of image matching technology acquisition satellite image.That is, pass through laser The three-dimensional coordinate and footmark size of high data are surveyed, same place is found out from image limiting in planar range, finally gives laser High data vertical control point pair is surveyed, this is to the object space plane coordinates and image pair including accurate object space height value, just slightly The image space ranks coordinate value answered.

It is exemplified below and how establishes broad sense vertical control point database.

The name of table, which should try one's best, in database uses directly perceived, succinct mark, and the classification and its name of table are included in name Claim.Among the title prefixing " TB " of table, classification and title plus " _ " is split.

For the ease of the operation and management of data, GLAS footmarks tables of data " TB_GLAS_PT " and GLAS footmark shadows are created As table " TB_IMG ".It is mainly used to the image block message that high data grid technology is surveyed in storage with GLAS in table " IMG_TABLE ".Table knot Structure is as shown in Table 1 and Table 2.

The GLAS_PT_TABLE table structures of table 1

The IMG_TABLE table structures of table 2

If developed using Oracle data platforms, existed by data and data base administration subsystem keeper GCPDB table spaces are created in oracle database, table space sets as follows：

● initial file size 5GB；

● AutoExtend1024KB is set；

● largest amount 20GB is 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 after Administrator by accusing Know user, and set GCPAdmin accounts that there is complete operation authority to GCPDB table spaces.

● database is named：GCPDB

● database user：GCPAdmin and with control point image storage, the note with management system operating right Volume user.

In data base design process, by distinguishing different visitors, different access types and different data pair As being treated respectively to carry out database security design.

For concrete example, broad sense vertical control point storehouse control point data source can use GLAS laser-measured height data 2007-2009 data, wherein, echo data is included in GLAS01 data, and altitude data is GLA14 data；Data model Enclose for 64 million figure amplitude ranges within Chinese territory, totally 1202338 points after preliminary screening.Footmark image data uses flat The national RGB orthography result map (DOM) of the resource that about 15 meters of face precision three, the result map is without elevation information, and plane Precision does not have image control point data base precision high.

Program is write using C# in the platforms of Microsoft Visual Studio 2010, reads the GLAS of binary format High data are surveyed, and stores and arrives database.

GLAS laser footmark image datas have chosen the image data of the regional extent for precision test and Preliminary Applications. Being produced and updated in basic surveying and mapping product based on the dominating pair of vertices of image has very important science and realistic meaning, while right Play the role of in the relation of analysis GLAS echo waveforms and features of terrain very big, therefore footmark image blocks data are as important number Be put in storage according to source.Program is write using C# in the platforms of Microsoft Visual Studio 2010, utilizes the grid increased income Conversion of Spatial Data storehouse GDAL (Geospatial Data Abstraction Library), surveying high data in using GLAS In the image blocks data Cun Chudao databases of the heart.Fig. 8 is that GLAS surveys high data footprints image (left side) and echo waveform in database The schematic diagram on (right side).Similar to Fig. 6 C, in the echo waveform figure of Fig. 8 right parts, axis of abscissas is the relative time for receiving echo, Unit nanosecond (ns), axis of ordinates are echo strengths, unit volt (volts).Footmark image does not join directly with echo waveform System, echo waveform is only relevant with the features of terrain in footmark and characters of ground object, as it was previously stated, footmark image is according to the flat of laser What areal coordinate obtained from DOM.

In order to verify the precision of GLAS broad sense vertical control points, a range of airborne LiDAR of Taiyuan region is collected into (Light Detection And Ranging, laser radar) data and Efficiency in Buildings in Tianjin Area dem data, to GLAS vertical control points Carry out precision test.

Taiyuan region LiDAR point cloud data break is about 3m, and coordinate system uses WGS-84, and elevation system is WGS-84 the earth It is high.Plane precision is about 0.5~0.8m, and absolute elevation precision is better than 0.15m.

For Taiyuan region LiDAR point data, precision evaluation is carried out to GLAS vertical control points using two methods.

If the GLAS by preliminary screening, which surveys high data, also has rough error value, footmark image and actual landform can be combined Excluding gross error value.

First method is to calculate the LiDAR that plan range is nearest in the range of the elevation and its footmark of GLAS vertical control points The difference of the elevation of data point.Statistical analysis, error is 0.603m in the discrepancy in elevation, worst error 1.617m, and minimal error is 0.017m.Statistical result is as shown in table 3.

Because the elevation of GLAS points is the elevation average in the range of a diameter of 70m hot spot, thus second method with Centered on GLAS vertical control points, the elevation average value of all LiDAR points in the range of radius 35m footmarks is calculated, obtains GLAS points Height value and LiDAR point elevation average value difference.Through statistics, error is 0.531m, worst error 1.645m in the discrepancy in elevation, Minimal error is 0.016m.Statistical result is as shown in table 4.

The GLAS points of table 3 and LiDAR closest approach discrepancy in elevation statistical forms

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

The GLAS points of table 4 and LiDAR point elevation average discrepancy in elevation statistical form

Discrepancy in elevation absolute value (m)	Quantity	Percentage
			<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 data is 1:The DEM result maps of 2000 engineer's scales, DEM achievements plane after Coordinate Conversion are sat Mark system is WGS-84, and height datum is WGS-84 geodetic heights.Through calculating DEM mesh spacing about 1.5m, absolute elevation precision expires The code requirement of the national principal scale DEM result maps of foot.

Be discrete vector point data because GLAS surveys high data, dem data is for raster data, it is necessary to by two kinds of data Form, which carries out unification, can just be analyzed.It is high using the Spatial Analyst Tool in ArcToolBox, extraction DEM Height value in DEM raster datas is extracted GLAS and surveys high number by journey value to point (extract values to points) In, as the attribute field of GLAS vector point datas, follow-up comparison and statistical analysis are carried out.

The difference of statistical analysis, the elevation of GLAS laser spots of the Efficiency in Buildings in Tianjin Area after screening and the height value extracted Middle error is 0.706m, worst error 1.767m, minimal error 0.0004m.Statistical result is as shown in table 5.

The GLAS points of table 5 and DEM elevational point discrepancy in elevation statistical forms

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

By carrying out GLAS laser footmarks image with reference to analysis, it is found that the larger laser spots of error are all to be located at mountain area, Or landform changes and has the region of significant small ditch, intuitively shows as the relation of elevation difference and landform, explanation The precision of GLAS laser-measured height data has close relationship with features of terrain.In the region that relief is larger, elevation difference Substantially, for example the larger GLAS points of Taiyuan region error are all located at mountain region substantially；In the area that landform is flat, elevation difference is universal It is smaller.Aggregate analysis, the GLAS laser-measured heights data after preliminary screening have higher precision, can be used as 1:50000 Or 1:2.5 ten thousand high process controls use.

Generally when using ground control point, control can be clearly recognized by the description of station at corresponding control point The point position of point and distribution characteristics.Because GLAS laser footmarks are diameter 70m hot spots, not corresponding point-to-point mark, Corresponding point position can not be accurately found in the range of hot spot, therefore is needed in test to the feature in GLAS laser faculas Point is matched.

According to an embodiment of the invention, alternatively, based on Forstner operator extraction characteristic points, GLAS broad sense elevations are realized Image feature point automatically extracting and matching in the range of the hot spot of control point, can not only so solve the tired of artificial visual differentiation Difficulty, while can also improve operating efficiency.

Its essence is a weighted operator for Forstner operators.It is the Robert's gradients and picture by calculating each pixel The gray scale covariance matrix of a window (such as 5 × 5) centered on plain (c, r), found in image it is minimum, close to round The point of error ellipse is as characteristic point.Because Forstner operators are more complicated, a substitution method is, first with simple poor Divide operator extraction just reconnaissance, for first reconnaissance using Forstner operators in 3 × 3 window calculation interest values, and select candidate point, Finally extraction extreme point is characterized a little.Its calculating process is as follows：

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

2) difference operator of pixel (c, r) in each piece of region is calculated.Pixel (c, r) is calculated in four direction up and down Grey scale difference absolute value dg1, 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) just reconnaissance and candidate point are determined, by dg1, dg2, dg3, dg4 according to being ranked up from low to high, and according to given Threshold value T extraction just reconnaissance.

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

For given threshold value T, if M>T, then (c, r) is a feature just reconnaissance inside the region；Otherwise, (c, r) is not Characteristic point.Then, in the 3x3 windows then centered on (c, r), covariance matrix N and error are calculated according to Forstner operators Oval circularity qc, r, then judges whether it is a feature candidate point according to the roundness threshold Tq of error ellipse.If limitation misses Poor oval long, the ratio between semi-minor axis is not more than 3.2~2.4, then can try to achieve Tq=0.32.

Wherein g_uAnd g_vIt is the partial differential along u and v directions respectively

The formula of summary 4.It can ask

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

q_c,r<T_q=0.32, then the pixel is a feature candidate point.

4) characteristic point is determined using weights as foundation, and the extreme point of the grid divided in extraction step (1) is as special Point is levied, that is, it is its characteristic point to extract the pixel in the grid corresponding to maximum weight person, it is determined that the foundation of power is as follows：

5) it is foundation according to weights, the point for choosing power the maximum in an appropriate window is characterized a little.

Therefore, during using Forstner operator extraction characteristic points, it is thus necessary to determine that two threshold values：One is that power is equal The threshold value T of value coefficient.Another is the threshold value T of interest value_q.To registration accuracy, fixed output quota life is very big really for the determination of the two threshold values Influence, threshold value selects the too small amount of calculation that will be greatly increased in characteristic point selection course, and threshold value set senior general to increase feature Point leakage choosing, the probability of wrong choosing.

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

As an example, two kinds of satellites of Ikonos-2 and ZY-3 are have chosen in Australian Hobart areas, Shanxi too Trizonal stereopsis carries out experiment altogether for former, weinan.Wherein Hobart areas are stood using the mono- scapes of Ikonos-2 The more scape stereopsis of the mono- scapes of ZY-3 and region have been respectively adopted in body image, Taiyuan and Weinan Prefecture.

The Ikonos-2 stereopsis acquisition time in Australian Hobart areas is 2 months 2003, and image resolution is 1m, by International Photography measurement with being provided free on remote sensing association ISRPS websites, wherein control point data is passed through by University of Melbourne Field GPS field surveys obtain, and height accuracy is better than 0.15m, have chosen 21 control points altogether.

Taiyuan region scope of experiment is：37.5448 ° -- 38.4856 ° of N, 112.3624 ° -- 112.9165 ° of E, experimental data For No. three three-line imageries of resource photographed on October 31st, 2012, wherein plane control data source is in 0.5m resolution ratio DOM images, for the GLAS data sources as high process control in 1 grade of in October, 2009 and 14 grades of products, elevation checks that data are The airborne LIDAR data of 3 meters of this area grid, precision about 15cm.All data switch to central meridian 111 under WGS84 ellipsoids The utm projection of degree.

Weinan Prefecture's scope of experiment is：33.955 ° -- 35.323 ° of N, 107.862 ° -- 109.278 ° of E, area about 20,000 are flat Fang Gongli, experimental data are the three rail resource three-line imageries photographed on July 28th, 2012, November 13, November 23, The whole area that surveys has 40 high-precision control points obtained using the measurement of GPS field operations (subsequently to make as horizontal control point and checkpoint With), there are on March 31st, 2007, on March 4th, 2008, on March 28th, 2009 totally three GLAS track datas, it is public by interval 25 In, 50 kilometers, 100 kilometers pick 14 GLAS points altogether as vertical control point.

Single scape solid adjustment experiment

For single scape stereopsis, employ the corner in aerophotogrammetry and deploy to ensure effective monitoring and control of illegal activities scheme, in order to verify plane and height Process control separates, and the control point of corner is only used as into plane control uses, while direction along ng a path selects two as high program control Point processed.The Ikonos-2 images in Hobart areas have added an elevation control data in intermediate region, and the ZY-3 of Taiyuan region is defended Star image adds one plane control data of cloth in centre.Adjustment Scheme point bit distribution schematic diagram is shown in Fig. 9, specific scheme explanation It is as follows：

Adjustment Scheme P1：0 horizontal control point, 0 vertical control point；

Adjustment Scheme P2：4 angle point cloth, 4 horizontal control points, 0 vertical control point；

Adjustment Scheme H1：0 horizontal control point, along 3, rail direction vertical control point；

Adjustment Scheme H2：4 horizontal control points, along 1, rail direction vertical control point；

Adjustment Scheme H3：4 horizontal control points, along 2, rail direction vertical control point；

Adjustment Scheme H4：4 horizontal control points, along 3, rail direction vertical control point.

The Ikonos images in Hobart areas are the stereogram products after certain geometric manipulations that ISRPS is provided, Therefore it can reach plane 3.561m, elevation 3.612m without control positioning precision, and it is without control height accuracy from 1:10000 three-dimensional surveys Figure precision also has certain gap.In experiment, using adjustment mode P2, plane precision can individually be carried using 4 horizontal control points Height arrives 0.837m, and height accuracy is not improved；Height accuracy can be improved only with the adjustment mode H1 of 3 vertical control points To within 1m, and plane precision does not improve.On the basis of having four horizontal control points, direction along ng a path be gradually added 1, 2 and 3 vertical control points, plane precision are held essentially constant, and 2 vertical control points are substantially better than and only use 1 elevation Control point, 3 vertical control points are then suitable with the adjustment result of 2.

For Taiyuan test block height above sea level from 700 meters to more than 1400 meters, a part belongs to Plain, and a part belongs to small high mountain region, The possible test block is just close to in-orbit geometry calibration area, therefore it is higher without control height accuracy as can be seen from Table 1, even if Be mountain area still up to 5.415m, after the plane control data of 4 angle points and 1 central point are added, plane positioning precision is notable Lifting, but height accuracy has no change；The GLAS laser-measured heights data after screening are selected as elevation control data, if A vertical control point is added, height accuracy has certain lifting, but effect is not particularly evident, and selection addition two ( Along rail direction two points are selected by two) then lifting precision, the good results are evident, up to 2.932m.

The Hobart of table 6 and Taiyuan region list scape stereopsis adjustment result statistical form

Regional network solid adjustment is tested

For No. three satellite images of 9 scape resource in weinan area in experimentation, a kind has been carried out respectively without control and 7 The adjustment that kind has control is tested.This 8 kinds of adjustment stationing modes are respectively：

Scheme P1：Adjustment of Free Networks without control point；

Scheme P2：Only 4 horizontal control points are laid altogether in 4 angle points；

Scheme P3：7 horizontal control points are laid in 4 angle points and centre altogether along rail direction；

Scheme P4：9 horizontal control points of laying altogether are arranged by three rows three；

Scheme H1：Only along GLAS course bearings equispaced 14 vertical control points are laid altogether for about 25 kilometers；

Scheme H2：On the basis of scheme P4, along the reconnaissance of about 100 kilometers of GLAS course bearings interval, 9 plane controls are laid altogether System point and 6 vertical control points；

Scheme H3：On the basis of scheme P4, along the reconnaissance of about 50 kilometers of GLAS course bearings interval, 9 plane controls are laid altogether System point and 8 vertical control points；

Scheme H4：Scheme P4 and H1 combination, 9 horizontal control points and 14 vertical control points are laid altogether.

7 kinds next have control adjustment point bit distribution it is as shown in Figure 10,8 kinds of adjustment results are shown in Table 7, and the unit of middle error is Rice.

No. three block adjustment result statistical forms of resource of table 7.GLAS altitude datas auxiliary

Weinan test block substantially belongs to alpine region, and height above sea level is from 350--1800 rice, and the discrepancy in elevation is up to 1500 meters, no control plane Positioning precision is poor, about 13.284m, and elevation location precision is 13.532m, and has obvious Systematic Errors.Add 4 angle points After horizontal control point, plane precision is obviously improved, and height accuracy is unchanged, is only had again by 7 points, 9 points afterwards The mode of plane control carries out adjustment, and plane precision increases but not notable, height accuracy no change；In no plane control In the case of, using 14 GLAS vertical control points, plane precision does not improve, but height accuracy is significantly improved, and reaches 3.066m；Based on 9 plane controls, then using 100 kilometers of interval, 50 kilometers, 25 kilometers of additions 6,8,14 GLAS vertical control points, plane precision is suitable with the adjustment result only with 9 plane controls, and height accuracy is with GLAS elevations The reduction at control point interval and improve constantly, height accuracy when using 8 GLAS vertical control points when about 50 kilometers among it 1 can be met:The requirement of 50000 stereoplottings, and be spaced 25 kms and choose height after 14 GLAS vertical control point adjustments Cheng Jingdu is up to 2.903m.

In above-mentioned experiment, single scape and more scene areas are carried out for the satellite image in three regions, two kinds of different resolutions Three-dimensional adjustment experiment is netted, its result can be shown that：1) for satellite stereo image, plane control and high process control are to separate , individually using a kind of control data and the result using flat high control the direction of accuracy improvements and be numerically it is consistent, When using high process control merely, plane precision does not improve, but height accuracy is significantly improved and control is mixed with flat height and ties Fruit is suitable, so for region overseas, only uses spaceborne laser altimeter system data, without any other ground control data, also may be used To carry out the stereoplotting of satellite image (plane precision relies on and the plane precision of satellite image in itself)；Use plane control merely Processed, height accuracy does not change substantially, and plane precision is significantly improved；2) 4 angle points or 4 angle points are used for single scape Adding the planar cloth prosecutor formula of 1 central point, more scapes carry out plane using the survey angle point of area four and image overlap area certain intervals and deployed to ensure effective monitoring and control of illegal activities, Good plane control effect can be reached, plane precision can have greatly improved；3) using the high program control of same accuracy level Elevation is data processed, 0.8m or so can be brought up to to Ikonos-2, and ZY-3 can only bring up to 2.9m.Illustrate the control of same type Data processed, adjustment is carried out to the satellite stereo image of different resolution, its result precision level is distinguishing, is differentiated with geometry Rate positive correlation, this can be following transmitting 1:10000 stereoplotting satellites provide reference；4) progress is separated with high process control in plane During adjustment, single scape uses along rail direction and lays two vertical control points at the both ends of image, between more scene area nets along rail using pressing Every less than every scape overlay length (such as ZY-3 takes the interval less than 50km) laying multiple row vertical control point, you can meet three-dimensional adjustment To the demand of high process control number.

Satellite photogrammetry and laser-measured height are two important technology hands that satellite earth observation obtains three-dimensional geographic information Section, the former plane precision is better than height accuracy under the conditions of different imaging modes determines no control, and the latter's height accuracy is better than Plane precision.According to an embodiment of the invention, it is proposed that the new method of a kind of plane and high process control separation progress adjustment, establish The part plan constraint of vertical control point and the local elevation of horizontal control point constrain, to inhomogeneous control information according to precision Different weights are invested, can effectively realize that the image control sheet established at present and laser-measured height data group cooperation are entered for voltage input The target of No. three three-dimensional adjustments of row resource, the final stereoplotting production efficiency for improving resource three.According to the implementation of the present invention Example, give the distribution scheme of the control data of specific elevation and planar separation.For example, region corner adds center cloth If 5 horizontal control points, can meet plane control demand, it is being less than 50km selection spaceborne laser altimeter system data along rail direction interval Elevation demand for control can be met as vertical control point, both are combined use, may finally meet stereoplotting demand.

Described above is only the exemplary embodiment of the present invention, not for limiting the scope of the invention, this hair Bright protection domain is determined by appended claim.

Claims (10)

1. A kind of 1. satellite image stereoblock adjustment method based on spaceborne laser altimeter system data, it is characterised in that including：

   Obtained respectively corresponding to the first control data of tie point and second from the first control data source and the second control data source Control data；

   It is three-dimensional that joint is carried out to first control data and second control data based on remote sensing stereopsis imaging model Adjustment processing, the parameter correction values of the remote sensing stereopsis imaging model are obtained,

   Wherein,

   The first control data source is image control point storehouse；

   The second control data source is broad sense vertical control point storehouse；

   First control data is the control data of control point or control point image film；

   Second control data is the control data of broad sense vertical control point；

   The control data of the broad sense vertical control point includes spaceborne laser altimeter system data,

   It is described from the first control data source and the second control data source obtain respectively corresponding to tie point the first control data and Second control data, including：

   Obtain the plane coordinates of spaceborne laser altimeter system data；

   The plane coordinates of the spaceborne laser altimeter system data is mapped to satellite remote sensing digital orthophoto map and gathers footmark shadow Picture；

   The spaceborne laser altimeter system data and the footmark image or the plane control data of a part for the footmark image Collectively form second control data；

   Obtained from image control point storehouse in the range of the part of the footmark image or the footmark image A number of same place is used as first control data as tie point, and using the control data of the same place；

   Wherein, the spaceborne laser altimeter system is a kind of ground point High head ratio technology, and it carries laser-measured height using satellite as platform Instrument is observed from space, time to the earth；

   The satellite image stereoblock adjustment method is by the sight line vector correcting method in tight imaging model and reasonable letter Image space correcting method in exponential model is unified, and forms unified three-dimensional adjustment Models.
2. 2. satellite image stereoblock adjustment method according to claim 1,

   Characterized in that,

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

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

   Wherein,

   (c, r) is picpointed coordinate,For topocentric coordinates, (v_c,v_r) it is picture point residual error,For ground three Dimension coordinate correction, (a₀,a₁,a₂,b₀,b₁,b₂) for the rigorous geometry model sight line vector of image correct parameter or rational function Model image space corrects parameter.
3. 3. satellite image stereoblock adjustment method according to claim 2,

   Characterized in that,

   For first control data, there is the first error equation：

   V=Ax+Bt-L

   Wherein,X=dH, or,

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

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

   V=Ax+Bt-L

   Wherein,

   By the first error equation described in simultaneous and second error equation and solve to carry out at the three-dimensional adjustment of joint Reason.
4. 4. satellite image stereoblock adjustment method according to claim 3,

   Characterized in that,

   Using least square method to first error equation and the second error equation simultaneous solution, including：

   By solving normal equation corresponding with first error equation and second error equation：Nx=B obtains every scape The model correction parameter of image, wherein, P is weight matrix.
5. 5. satellite image stereoblock adjustment method according to claim 3,

   Characterized in that,

   For second error equation, formed by the way of member is cut point by point relevant with the image space correction factor of every image Changing normal equation.
6. 6. satellite image stereoblock adjustment method according to claim 1,

   Characterized in that,

   The spaceborne laser altimeter system data are GLAS spaceborne laser altimeter system data.
7. 7. the satellite image stereoblock adjustment method according to claim 1 or 6,

   Characterized in that,

   After obtaining spaceborne laser altimeter system data from the second control data source, by just omiting terrain data excluding gross error and echo Second control data for the adjustment processing is obtained after waveform screening.
8. 8. satellite image stereoblock adjustment method according to claim 7,

   Characterized in that,

   Obtaining second control data by echo waveform screening includes：

   Estimate the background noise level of echo waveform；

   The noise signal of the echo waveform is filtered out using smothing filtering algorithm；

   The Gauss curve fitting that initial parameter is carried out to filtered signal is estimated；

   The result estimated using the Gauss curve fitting, with reference to initial parameter, iterated extraction using nonlinear least square method Accurate Gaussian component parameter.
9. 9. satellite image stereoblock adjustment method according to claim 8,

   Characterized in that,

   UsingThe Gauss curve fitting estimation is carried out, wherein w (t) is the echo energy of t Value, W_mIt is m-th of Gaussian function in the backward energy value of t, Np is the number of Gaussian function, and ε is noise, t_mFor m-th The average of Gaussian function, σ_mFor the mean square deviation of m-th of Gaussian function, A_mThe peak value of m-th of Gaussian function.
10. 10. satellite image stereoblock adjustment method according to claim 6,

    Characterized in that, got off the GLAS laser by the way that GLAS data are transformed into WGS84 coordinate systems in coordinate system aspect Survey high data and the imaging model of stereoblock adjustment is mutually unified.