CN104931022A  Satellite image threedimensional area network adjustment method based on satelliteborne laser height measurement data  Google Patents
Satellite image threedimensional area network adjustment method based on satelliteborne laser height measurement data Download PDFInfo
 Publication number
 CN104931022A CN104931022A CN201510191096.9A CN201510191096A CN104931022A CN 104931022 A CN104931022 A CN 104931022A CN 201510191096 A CN201510191096 A CN 201510191096A CN 104931022 A CN104931022 A CN 104931022A
 Authority
 CN
 China
 Prior art keywords
 described
 control data
 image
 data
 satellite
 Prior art date
Links
 238000003384 imaging method Methods 0 abstract 2
 230000000875 corresponding Effects 0 abstract 1
 238000009740 moulding (composite fabrication) Methods 0 abstract 1
 230000036961 partial Effects 0 abstract 1
Classifications

 G—PHYSICS
 G01—MEASURING; TESTING
 G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
 G01C11/00—Photogrammetry or videogrammetry, e.g. stereogrammetry; Photographic surveying
Abstract
Description
Technical field
The present invention relates to a kind of satellite remotesensing 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, threedimensional 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 threedimensional 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 geolocation 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 onthespot 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 GPSRTK, 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 remotesensing 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 selfevident.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 multisource heterogeneous remote sensing image reference mark automatic acquiring method, can automatically extract remote sensing from multisource 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 remotesensing 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 threedimensional 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+BtL, the unification in adjustment processing by rigorous geometry model and rational function model, wherein
(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 _{0}, a _{1}, a _{2}, b _{0}, b _{1}, b _{2}) 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+BtL
Wherein,
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+BtL
Wherein,
Solve carry out combining threedimensional 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
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 _{m}be 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, σ _{m}be the mean square deviation of m Gaussian function, A _{m}the 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 lasermeasured 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 threedimensional 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 threedimensional localization precision of satellite, traditional mode generally adopts flat high point as reference mark, namely field GPS measurement point or highprecision topomap measures the threedimensional 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 threedimensional 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 lasermeasured 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 threedimensional 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 Automatching 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 multisource 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, threedimensional 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 expressionform is respectively (3) and (4),
Δc＝f _{1}(r,c)＝a _{0}(3)
Δr＝g _{1}(r,c)＝b _{0}
Δc＝f _{2}(r,c)＝a _{0}+a _{1}c+a _{2}r (4)
Δr＝g _{2}(r,c)＝b _{0}+b _{1}c+b _{2}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:
Abovementioned (4) formula is substituted into, is:
Then abovementioned (7) formula can be rewritten as V=Ax+BtL (8), wherein,
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, highprecision 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 _{1}X _{1}y _{1}z _{1}: be determined by the attitude control system of satellite (ACS), wherein initial point is at centroid of satellite, X1 axle, Y _{1}axle, Z _{1}axle gets three principal moments axles of satellite respectively.X _{1}axle along satellite transverse axis, Y _{1}axle points to satellite flight direction (not exclusively overlapping) along the longitudinal axis, Z _{1}axle is determined according to righthanded 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 _{1}the rotation of axle, roll (roll) is around body coordinate system Y _{1}the rotation of axle, course deviation (yaw) is around body coordinate system Z _{1}the rotation of axle.
(ψ _{x}) _{p}(ψ _{y}) _{p}the inclination angle in direction along ng a path and vertical track direction when representing chargecoupled 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 _{2}it is reverse that axle points to the earth's core, Y _{2}axle points to the direction of satellite motion on satellite orbit face, X _{2}axle is determined according to righthanded system rule.
O _{2}X _{2}y _{2}z _{2}the direction of each axle is by the position of satellite in t and speed two direction vectors are determined, as follows:
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 bodyfixed 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:
The conversion parameter table of ITRF90 to WGS84
Conversion formula is as follows:
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 Xaxis points to IERS guide meridian with mistake initial point and perpendicular to the intersection point of Z axis plane, Yaxis is determined by righthanded 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×(ll _{c})
Wherein, t _{c}for the photography time of central scan row, lsp is the sweep time of often going, l _{c}for the line number of central scan row
II, interpolation go out the position of t satellite with
III, the imaging inclination angle (ψ arranged by p _{x}) _{p}(ψ _{y}) _{p}calculate (l, p) imaging direction vector under navigational coordinate system
Imaging direction vector under IV, calculating orbital coordinate system
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}(t _{i}), a _{r}(t _{i}), a _{y}(t _{i}) be respectively t _{i}the 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
Imaging direction vector under V, calculating earth axes
Wherein, X _{2}, Y _{2}, Z _{2}for orbital coordinate system OX _{2}y _{2}z _{2}three change in coordinate axis direction vectors
VI, topocentric location:
As shown in Figure 4, measurement vector u _{3}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:
Wherein,
In conjunction with several transformation relations above, write as the form of similar collinearity condition equation:
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:
Can collinearity equation be obtained:
Namely
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 threedimensional adjustment, sight line vector corrects equation and is:
Wherein, (i, j) refers to the ranks number of picpointed coordinate.
Therefore formula (12), can be written as:
Be abbreviated as:
Wherein,
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 _{1}, v _{2}) be residual error, be the Three Dimensional Ground coordinate corrective value of nth, (a _{0m}, a _{1m}, a _{2m}, b _{0m}, b _{1m}, b _{2m}) be the sight line vector correction parameter that m opens image, k _{mn}be at nth 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+BtL, rigorous geometry model and rational function model are united in adjustment processing.
Wherein
Open on image for ith point in jth, then bring corresponding subscript respectively.Such as, open on image for ith point in jth, then V _{ij}=A _{ij}x _{i}+ B _{ij}t _{j}L _{ij}, weight matrix is P _{ij}.
Next, on the basis of the error equation derived at abovementioned 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 abovementioned 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
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+BtL, weight matrix is P; Can list normal equation by least square is:
NX＝Q (22)
Wherein,
Abovementioned 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: σ _{x}millimeter, then weigh P _{x}=1;
The precision of ground control point obserred coordinate value is: σ _{c}rice, then weigh wherein for photographic scale.Basic law is that precision is higher, i.e. σ _{c}less, 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 abovementioned 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 usecase, 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 abovementioned 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 threedimensional 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 threedimensional 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 lasermeasured 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 highprecision altitude figures, its laser footmark plane precision reaches 10m magnitude, and vertical accuracy is about 15cm.Within 2017, be about to launch ICESat2, 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'e2 " satellite of independent research have all carried highprecision 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 lasermeasured height system of earth observation, but list followup earth observation lasermeasured height satellite development plan in, as the followup star task of resource No. three high resolving power stereoplotting satellites, satellite surveying and mapping " 12 " planning also claims to the work of lasermeasured 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 highprecision 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 ISIPS (ICESAT ScienceInvestigatorled 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 GLA01GLA15.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 GLA01GLA04.L1B DBMS is corrected through device by L1A DBMS, data processing, based on attitude and position data carry out laser footpoint location formed, comprise GLA05GLA07.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 GLA08GLA15.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 threedimensional 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 threedimensional 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 echopeak, 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 bestfit effect.
In addition because laser firing pulses is through propagation in atmosphere, multiple processes such as target reflection, echo reception optoelectronic conversion are mixed into much noise, therefore need to estimate through noise, smothing filtering, Gaussian parameter preestimation 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 _{m}be 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, σ _{m}be the mean square deviation of m Gaussian function, A _{m}the 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 followup 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:
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:
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 longtail portion, right side.Larger negative value shows there is comparatively longtail 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 lasermeasured height data and threedimensional 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 lasermeasured height data is TOPEX/Poseidon, under needing first GLAS data to be transformed into WGS84 reference ellipsoid, be convenient to followuply 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).
In formula (25), B represents latitude, the difference of what da represented is WGS84 ellipsoid and TOPEX/Poseidon ellipsoid major axis, db represents the difference of WGS84 ellipsoid and TOPEX/Poseidon ellipsoid minor axis, dh is elevation correction value, the elevation of Hg relative TOPEX ellipsoid before being conversion, H _{g}for 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 lasermeasured 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 than20m 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 lasermeasured 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 threedimensional coordinate and the footmark size of lasermeasured height data, same place is found out from image in restriction planar range, finally obtain lasermeasured 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 warehousein, 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 lasermeasured height data 20072009, 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 WGS84, and elevation is WGS84 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
Table 4 GLAS point and LiDAR point elevation average discrepancy in elevation statistical form
Efficiency in Buildings in Tianjin Area test figure is the DEM result map of 1:2000 engineer's scale, and this DEM achievement is WGS84 through coordinate conversion back plane coordinate system, and height datum is WGS84 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 followup 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
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 lasermeasured 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 lasermeasured 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 pointtopoint 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 equalsized 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 _{1}＝g _{c,r}g _{c+1,r} (26)
dg _{2}＝g _{c,r}g _{c,r+1} (27)
dg _{3}＝g _{c,r}g _{c1,r} (28)
dg _{4}＝g _{c,r}g _{c,r1} (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 _{1},dg _{2},dg _{3},dg _{4}} (30)
For given threshold value T, if M>T, then (c, r) is the just reconnaissance of this intrazone 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 semiaxis is not more than 3.2 ~ 2.4, then can try to achieve Tq=0.32.
Wherein g _{u}and g _{v}the partial differential along u and v direction respectively
Comprehensively abovementioned 4 formulas.Can ask
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:
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 Ikonos2 and ZY3 two kinds of satellites in Australian Hobart area, Taiyuan, Shanxi, weinan altogether trizonal stereopsis carry out experiment.What wherein Hobart area adopted is the monoscape stereopsis of Ikonos2, and Taiyuan and Weinan Prefecture have employed the monoscape of ZY3 and region many scapes stereopsis respectively.
The Ikonos2 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 threeline 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 threeline imageries photographed November 23, 40 highprecision control points (followup 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 threedimensional 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 Ikonos2 image in Hobart area has added an elevation control data at zone line, and the ZY3 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 inorbit 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 lasermeasured 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 threedimensional 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 3501800 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 abovementioned experiment, for three regions, the satellite image of two kinds of different resolutions carries out single scape and the threedimensional 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 Ikonos2, and ZY3 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 ZY3 gets the interval being less than 50km) along rail by interval and lays multiple row vertical control point, can meet threedimensional adjustment to elevation domination number object demand.
Satellite photogrammetry and lasermeasured height are two important technical that satellite earth observation obtains threedimensional 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 lasermeasured height data assemblies carries out the target of resource No. three threedimensional 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)
Priority Applications (1)
Application Number  Priority Date  Filing Date  Title 

CN201510191096.9A CN104931022B (en)  20150421  20150421  Satellite image stereoblock adjustment method based on spaceborne laser altimeter system data 
Applications Claiming Priority (1)
Application Number  Priority Date  Filing Date  Title 

CN201510191096.9A CN104931022B (en)  20150421  20150421  Satellite image stereoblock adjustment method based on spaceborne laser altimeter system data 
Publications (2)
Publication Number  Publication Date 

CN104931022A true CN104931022A (en)  20150923 
CN104931022B CN104931022B (en)  20180316 
Family
ID=54118300
Family Applications (1)
Application Number  Title  Priority Date  Filing Date 

CN201510191096.9A CN104931022B (en)  20150421  20150421  Satellite image stereoblock adjustment method based on spaceborne laser altimeter system data 
Country Status (1)
Country  Link 

CN (1)  CN104931022B (en) 
Cited By (14)
Publication number  Priority date  Publication date  Assignee  Title 

CN105527620A (en) *  20151201  20160427  中国科学院上海技术物理研究所  Aerosol thickness and laser radar ranging delay automatic calibration method 
CN105547254A (en) *  20160130  20160504  武汉大学  Normal equationbased quickbuilt superlarge scale block adjustment method and system 
CN105716580A (en) *  20160130  20160629  武汉大学  Superlargescale block adjustment method and system based on adjustment parameter fast resolution 
CN106248053A (en) *  20160831  20161221  铁道第三勘察设计院集团有限公司  Existing Railway Line plane center line manufacture method based on aerial survey mode 
CN106646430A (en) *  20161226  20170510  国家测绘地理信息局卫星测绘应用中心  Laser footprint center determining method based on ground detector 
CN106643804A (en) *  20161230  20170510  国家测绘地理信息局卫星测绘应用中心  Method of predetermining footprint position of satellite laser altimeter 
CN106885585A (en) *  20161230  20170623  国家测绘地理信息局卫星测绘应用中心  A kind of satellite borne photography measuring system integration calibration method based on bundle adjustment 
CN106908784A (en) *  20170216  20170630  中国电子科技集团公司第二十八研究所  A kind of multisource foreign peoples's radar collaborative detection method based on threedimensional grid 
CN107014399A (en) *  20170425  20170804  国家测绘地理信息局卫星测绘应用中心  A kind of spaceborne optical camera laser range finder combined system joint calibration method 
CN107421504A (en) *  20170802  20171201  中国科学院遥感与数字地球研究所  The shooting time computational methods of month base earth observation electrooptical photo 
CN107504981A (en) *  20170725  20171222  国家测绘地理信息局卫星测绘应用中心  A kind of attitude of satellite error correcting method and equipment based on lasermeasured height data 
WO2018223950A1 (en) *  20170605  20181213  中国测绘科学研究院  Method for extracting elevation control point with assistance of satellite laser altimetry data 
CN109741381A (en) *  20190123  20190510  张过  Spaceborne pushbroom type optical sensor high frequency error removing method based on parallel observation 
CN106885585B (en) *  20161230  20200121  自然资源部国土卫星遥感应用中心  Integrated calibration method of satelliteborne photogrammetry system based on light beam adjustment 
Citations (3)
Publication number  Priority date  Publication date  Assignee  Title 

US20080089610A1 (en) *  20040116  20080417  Geotango International Corp.  System, computer program and method for 3d object measurement, modeling and mapping from single imagery 
CN101609149A (en) *  20090713  20091223  北京航空航天大学  Method for improving attitude determination precision of airborne laser radar 
CN101846512A (en) *  20100603  20100929  东南大学  Method for implementing land patrol 

2015
 20150421 CN CN201510191096.9A patent/CN104931022B/en active IP Right Grant
Patent Citations (4)
Publication number  Priority date  Publication date  Assignee  Title 

US20080089610A1 (en) *  20040116  20080417  Geotango International Corp.  System, computer program and method for 3d object measurement, modeling and mapping from single imagery 
US7995862B2 (en) *  20040116  20110809  Microsoft Corporation  System, computer program and method for 3D object measurement, modeling and mapping from single imagery 
CN101609149A (en) *  20090713  20091223  北京航空航天大学  Method for improving attitude determination precision of airborne laser radar 
CN101846512A (en) *  20100603  20100929  东南大学  Method for implementing land patrol 
NonPatent Citations (2)
Title 

何钰: "基于月面CCD影像和激光测高数据的月球形貌测绘技术研究", 《信息科技辑》 * 
李国元等: "GLAS激光测高数据辅助的资源三号三线阵区域网平差研究", 《第三届高分辩率对地观测学术年会优秀论文集》 * 
Cited By (22)
Publication number  Priority date  Publication date  Assignee  Title 

CN105527620B (en) *  20151201  20171027  中国科学院上海技术物理研究所  The automatic calibration method that a kind of aerosol thickness postpones with laser radar range 
CN105527620A (en) *  20151201  20160427  中国科学院上海技术物理研究所  Aerosol thickness and laser radar ranging delay automatic calibration method 
CN105547254B (en) *  20160130  20180112  武汉大学  Ultralarge block adjustment method and system based on normal equation rapid build 
CN105547254A (en) *  20160130  20160504  武汉大学  Normal equationbased quickbuilt superlarge scale block adjustment method and system 
CN105716580A (en) *  20160130  20160629  武汉大学  Superlargescale block adjustment method and system based on adjustment parameter fast resolution 
CN105716580B (en) *  20160130  20180123  武汉大学  The ultralarge block adjustment method and system quickly resolved based on adjustment parameter 
CN106248053B (en) *  20160831  20181019  中国铁路设计集团有限公司  Existing Railway Line plane center line production method based on aerial survey mode 
CN106248053A (en) *  20160831  20161221  铁道第三勘察设计院集团有限公司  Existing Railway Line plane center line manufacture method based on aerial survey mode 
CN106646430A (en) *  20161226  20170510  国家测绘地理信息局卫星测绘应用中心  Laser footprint center determining method based on ground detector 
CN106885585A (en) *  20161230  20170623  国家测绘地理信息局卫星测绘应用中心  A kind of satellite borne photography measuring system integration calibration method based on bundle adjustment 
CN106643804B (en) *  20161230  20191122  自然资源部国土卫星遥感应用中心  A kind of method of the footmark position of predetermined spaceborne laser altimeter system instrument 
CN106643804A (en) *  20161230  20170510  国家测绘地理信息局卫星测绘应用中心  Method of predetermining footprint position of satellite laser altimeter 
CN106885585B (en) *  20161230  20200121  自然资源部国土卫星遥感应用中心  Integrated calibration method of satelliteborne photogrammetry system based on light beam adjustment 
CN106908784A (en) *  20170216  20170630  中国电子科技集团公司第二十八研究所  A kind of multisource foreign peoples's radar collaborative detection method based on threedimensional grid 
CN106908784B (en) *  20170216  20190326  中国电子科技集团公司第二十八研究所  A kind of multisource foreign peoples's radar collaborative detection method based on threedimensional grid 
CN107014399A (en) *  20170425  20170804  国家测绘地理信息局卫星测绘应用中心  A kind of spaceborne optical camera laser range finder combined system joint calibration method 
WO2018223950A1 (en) *  20170605  20181213  中国测绘科学研究院  Method for extracting elevation control point with assistance of satellite laser altimetry data 
CN107504981B (en) *  20170725  20200103  国家测绘地理信息局卫星测绘应用中心  Satellite attitude error correction method and device based on laser height measurement data 
CN107504981A (en) *  20170725  20171222  国家测绘地理信息局卫星测绘应用中心  A kind of attitude of satellite error correcting method and equipment based on lasermeasured height data 
CN107421504B (en) *  20170802  20191112  中国科学院遥感与数字地球研究所  The shooting time calculation method of month base earth observation electrooptical photo 
CN107421504A (en) *  20170802  20171201  中国科学院遥感与数字地球研究所  The shooting time computational methods of month base earth observation electrooptical photo 
CN109741381A (en) *  20190123  20190510  张过  Spaceborne pushbroom type optical sensor high frequency error removing method based on parallel observation 
Also Published As
Publication number  Publication date 

CN104931022B (en)  20180316 
Similar Documents
Publication  Publication Date  Title 

Baltsavias et al.  Radiometric and geometric evaluation of Ikonos GEO images and their use for 3D building modelling  
Konecny  Geoinformation: remote sensing, photogrammetry and geographic information systems  
EP1709396B1 (en)  System,computer program and method for 3d object measurement, modeling and mapping from single imagery  
Muller et al.  MISR stereoscopic image matchers: Techniques and results  
Toutin  Comparison of stereoextracted DTM from different highresolution sensors: SPOT5, EROSA, IKONOSII, and QuickBird  
Mukherjee et al.  Evaluation of vertical accuracy of open source Digital Elevation Model (DEM)  
JP2012511697A (en)  How to generate a geodetic reference database  
Rexer et al.  Comparison of free high resolution digital elevation data sets (ASTER GDEM2, SRTM v2. 1/v4. 1) and validation against accurate heights from the Australian National Gravity Database  
Fernández‐Hernandez et al.  Image‐based modelling from unmanned aerial vehicle (UAV) photogrammetry: an effective, low‐cost tool for archaeological applications  
Li  Potential of highresolution satellite imagery for national mapping products  
AlRousan et al.  Automated DEM extraction and orthoimage generation from SPOT level 1B imagery  
Berthier et al.  Glacier topography and elevation changes derived from Pléiades submeter stereo images  
KR20090024146A (en)  Method and system for generating an imagetextured digital surface model (dsm) for a geographical area of interest  
US8427505B2 (en)  Geospatial modeling system for images and related methods  
US7194389B2 (en)  Fusion of data from differing mathematical models  
US7944547B2 (en)  Method and system of generating 3D images with airborne oblique/vertical imagery, GPS/IMU data, and LIDAR elevation data  
Altmaier et al.  Digital surface model generation from CORONA satellite images  
Toutin  Error tracking of radargrammetric DEM from RADARSAT images  
Niethammer et al.  Open source imageprocessing tools for lowcost UAVbased landslide investigations  
Zhang et al.  Photogrammetric processing of low‐altitude images acquired by unpiloted aerial vehicles  
Kaichang et al.  Geometric processing of Ikonos stereo imagery for coastal mapping applications  
Ridley et al.  Evaluating the potential of the forthcoming commercial US highresolution satellite sensor imagery at the ordnance survey (R)  
Raggam et al.  Assessment of the stereoradargrammetric mapping potential of TerraSARX multibeam spotlight data  
Habib et al.  Error budget of LiDAR systems and quality control of the derived data  
Takaku et al.  GENERATION OF HIGH RESOLUTION GLOBAL DSM FROM ALOS PRISM. 
Legal Events
Date  Code  Title  Description 

C06  Publication  
PB01  Publication  
C10  Entry into substantive examination  
SE01  Entry into force of request for substantive examination  
GR01  Patent grant  
GR01  Patent grant 