CN109903352B - Method for making large-area seamless orthoimage of satellite remote sensing image - Google Patents

Method for making large-area seamless orthoimage of satellite remote sensing image Download PDF

Info

Publication number
CN109903352B
CN109903352B CN201910057228.7A CN201910057228A CN109903352B CN 109903352 B CN109903352 B CN 109903352B CN 201910057228 A CN201910057228 A CN 201910057228A CN 109903352 B CN109903352 B CN 109903352B
Authority
CN
China
Prior art keywords
image
remote sensing
dom
point
sensing image
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910057228.7A
Other languages
Chinese (zh)
Other versions
CN109903352A (en
Inventor
刘斌
贾萌娜
辛鑫
邸凯昌
刘召芹
岳宗玉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Remote Sensing and Digital Earth of CAS
Original Assignee
Institute of Remote Sensing and Digital Earth of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Remote Sensing and Digital Earth of CAS filed Critical Institute of Remote Sensing and Digital Earth of CAS
Publication of CN109903352A publication Critical patent/CN109903352A/en
Application granted granted Critical
Publication of CN109903352B publication Critical patent/CN109903352B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Processing (AREA)

Abstract

The invention relates to a method for making a large-area seamless orthoimage of a satellite remote sensing image, which comprises the following steps: 1) constructing a remote sensing image imaging geometric model; 2) reading the remote sensing image to be matched and DEM data, and matching the DEM and the remote sensing image data to be registered; 3) orthorectifying the remote sensing image to obtain a corresponding DOM product; 4) obtaining the position relation between the orthoimages according to the step 3), determining the overlapping area between the images, matching in the overlapping area to obtain homonymous points, establishing an error equation based on a TPS (Perfect transform) model for each group of homonymous points, and obtaining the correction quantity of the geometric deviation between the images through iterative solution; 5) and (4) carrying out light evening and embedding on all corrected orthographic images to obtain an orthographic image map product of a large area. According to the method, the thin-plate spline model is introduced in the process of the orthoimage adjustment, so that the problem of low precision of the traditional geometric correction model caused by a small overlapping range between images is solved.

Description

Method for making large-area seamless orthoimage of satellite remote sensing image
Technical Field
The invention relates to a method for making a large-area seamless orthoimage of a satellite remote sensing image, and relates to the technical field of remote sensing images in photogrammetry.
Background
Digital ortho-shadowgraph (DOM) is an important component of traditional surveying and mapping 4D products, is an important carrier for reflecting topography and feature, and is also a basic product for serving remote sensing data in various industries. The process of generating the Digital ortho-image generally includes performing pixel-by-pixel projection difference correction on the remote sensing image by using a Digital Elevation Model (DEM), eliminating various distortion and displacement errors, and finally obtaining a satellite remote sensing image product containing geographic information and various specialties. The originally obtained remote sensing image is affected by distortion of an optical system inside a sensor, nonlinearity of a scanning system, external attitude change, earth curvature, topographic relief and the like during imaging, and all the remote sensing image has distortion and distortion of different degrees, so that a projection image product generated by directly splicing images with certain overlapping degree is generated, and the situation of inconsistent geometric positions appears at the position of a joint edge. In addition, various errors introduced during the imaging process may cause positional deviation between the image and the existing elevation data, thereby affecting the accuracy of the projection difference and the terrain configuration. In the photogrammetry of lunar and mars orbiter, the lack of high-precision reference and ground control causes the problem of inconsistency between images and between the images and elevation data to be more prominent, and finally influences the engineering and scientific research application research of map products.
In the existing method, an existing DOM product is used as a medium, a control point is obtained between an image and a DEM, the control point and connection points between all images are simultaneously brought into adjustment of a local area network, and an imaging geometric model of the image is refined by adjustment calculation, so that the geometric consistency between the images and between the image and the DEM is finally improved; and then orthorectifying the image by using the refined geometric model and the DEM. However, when processing large-area image data, the problems of excessive calculation amount of the adjustment of the area network, unstable calculation, difficult convergence of iteration, unsatisfactory adjustment accuracy and the like are caused due to the fact that the related drawing range is wider, the number of included images is larger, the accuracy of control data does not meet the requirements and the like. In addition, if the overlapping range between images is too small with respect to the overall area of the images, it is difficult for the conventional block adjustment method based on the polynomial model to obtain the required adjustment accuracy.
Disclosure of Invention
In view of the above problems, an object of the present invention is to provide a method for making a large-area seamless ortho-image of a satellite remote sensing image, which facilitates convergence of adjustment iteration and is beneficial to improving adjustment accuracy.
In order to achieve the purpose, the invention adopts the following technical scheme: a method for making a large-area seamless orthoimage of a satellite remote sensing image comprises the following contents:
1) constructing a remote sensing image imaging geometric model;
2) reading the remote sensing image to be matched and DEM data, and matching the DEM and the remote sensing image data to be registered;
3) orthorectifying the remote sensing image to obtain a corresponding DOM product;
4) obtaining the position relation between the orthoimages according to the step 3), determining the overlapping area between the images, matching in the overlapping area to obtain homonymous points, establishing an error equation based on a TPS (Perfect transform) model for each group of homonymous points, and obtaining the correction quantity of the geometric deviation between the images through iterative solution;
5) and (4) carrying out light evening and embedding on all corrected orthographic images to obtain an orthographic image map product of a large area.
Further, the DEM and the remote sensing image data to be registered are matched, and the DEM and the image which are constrained by geometric and radiation information together are matched, and the specific process is as follows:
2.1) generating a simulation image by the DEM to be used as a reference image;
2.2) obtaining initial matching points of the reference image and the remote sensing image by using an affine-scale invariant feature transformation algorithm;
2.3) calculating initial refinement parameters of the rational function model by using the initial matching points;
2.4) obtaining an initial value of least square fine matching;
and 2.5) outputting the radiation deformation parameters and the geometric deformation parameters of the remote sensing images to be matched based on the least square matching of the global constraint.
Further, in the step 2.1), the analog image generated by the DEM is used as a reference image, and the specific process is as follows:
2.1.1) calculating gradients p and q of each pixel in x and y directions pixel by pixel on the DEM;
2.1.2) calculating the gradient and the slope direction value pixel by pixel according to the gradient value:
Figure BDA0001952915370000021
A=p/q
wherein S is a gradient value and A is a gradient value;
2.1.3) obtaining a solar azimuth angle and a solar altitude angle when the remote sensing image is obtained according to the head file information of the remote sensing image data, creating a null image with the same size as the DEM, calculating the corresponding simulated radiation value information pixel by pixel, storing the simulated radiation value information in the null image, and constructing to obtain a reference image.
Further, in the step 2.2), an affine-scale invariant feature transformation algorithm is used to obtain initial matching points of the reference image and the remote sensing image, and the specific process is as follows:
firstly, acquiring a sampling sequence of a longitude angle and a latitude angle of a camera, and simulating affine deformation possibly caused by an optical axis of the camera to realize image transformation;
then, carrying out oblique rotation transformation on the remote sensing image to be matched to generate a simulated image;
and finally, performing feature point detection and reference image matching on the generated simulated image according to a scale-invariant feature transformation algorithm to obtain an initial matching point between the remote sensing image and the reference image.
Further, the specific process of obtaining the initial value of the least square fine matching in the step 2.4) is as follows:
2.4.1) extracting the feature points on the reference image by using a feature point extraction operator;
2.4.2) calculating the characteristic points on the reference image to the remote sensing image by using the refined rational function model to be used as the initial value of least square matching.
Further, the step 2.5) outputs the radiation deformation parameters and the geometric deformation parameters of the remote sensing image to be matched based on the least square matching of the global constraint, and the specific process is as follows:
2.5.1) solving registration parameters of the remote sensing image and the reference image, simultaneously solving refinement parameters of a rational function model, and introducing a global constraint least square matching model as follows:
g1(x,y)+n1(x,y)=h0+h1g2(r,c)+n2(x,y)
r=ar 0+ar 1r′+ar 2c′+δr
c=bc 0+bc 1r′+bc 2c′+δc
in the formula, g1And g2Is a function of the radiation value of the corresponding coordinate position of the reference image and the remote sensing image connecting point to be registered, n1And n2Is a function of the noise of the reference image and the remote sensing image to be registered, (x, y) is the coordinate of the connecting point on the reference image, h0,h1Is a radiation deformation parameter between the reference image and the remote sensing image, ar 0……ar 2,bc 0……bc 2The geometric deformation parameters (refinement parameters), (r ', c') are projection coordinates of ground points back projected to the image side by RFM, (r, c) are coordinates of image points obtained by refinement of the projection coordinates back projected to the image side, and deltarcIs the model error;
2.5.2) establishing an error equation according to the global constraint least square matching model:
vg=h0+h1g2(r,c)-g1(x,y)
vr=ar 0+ar 1r′+ar 2c′+δr-r
vc=bc 0+bc 1r′+bc 2c′+δc-c
in the formula, vgError of radiation value between registered images, vrError of equation in column direction, vcIs the equation in the row directionAn error;
and performing first-order Taylor expansion on the error equation to obtain a linearized error equation:
Figure BDA0001952915370000031
Figure BDA0001952915370000034
Figure BDA0001952915370000035
in the formula (I), the compound is shown in the specification,
Figure BDA0001952915370000032
is the 0-order term of the Taylor series,
Figure BDA0001952915370000033
respectively representing the first partial derivatives of the coordinates in the column direction and the row direction;
2.5.3) calculating a normal equation point by point according to an error equation for each input matching initial value point;
2.5.4) solving a normal equation, correcting the unknown number, judging the correlation coefficient obtained by calculation after adjustment, if the sum of the correlation coefficients of all the points meets the requirement, outputting an adjustment result, and otherwise, entering the step 2.5.2 to carry out iterative solution;
2.5.5) outputting the radiation deformation parameters and the geometric deformation parameters of the remote sensing image, and completing the global least square matching between the remote sensing image and the reference image.
Further, the remote sensing image ortho-rectification in the step 3) adopts a forward solution method and a reverse solution method.
Further, the specific process of the remote sensing image orthorectification by adopting the inverse solution method comprises the following steps:
3.1) setting the X-direction and Y-direction resolutions (R) for generating an orthoimagexdom,Rydom) Using the original image size and imageThe geometric model of (2) calculating the coverage of the image, and setting the coordinates (X) of the upper left corner point of the image according to the range of the image0,Y0) The coordinate of the pixel point of the pixel center P of any point on the preset generated DOM is set as (x)dom,ydom) Calculating the ground coordinates (X, Y) corresponding to the P point;
X=X0+Rxdom*xdom
Y=Y0+Rydom*ydom
3.2) interpolating the elevation Z of the point by using (X, Y) and DEM, and calculating the corresponding image point coordinate p (X, Y) on the original image by using the established and refined image imaging geometric model;
x=f1(X,Y,Z)
y=f2(X,Y,Z)
wherein f1 and f2 represent the refined geometric model of the imaged image, and Z represents the elevation of the point P;
3.3) gray scale interpolation, obtaining the gray scale value g (x, y) of the image point p;
3.4) assigning the gray value of the image point P to the corrected image element, namely the P point on the DOM;
and 3.5) completing the operation on each pixel on the DOM in sequence to obtain the DOM product subjected to terrain correction.
Further, in the step 4), according to the position relationship between the orthoimages obtained in the step 3), an overlapping region between the images is determined, homonymous points are obtained by matching in the overlapping region, then an error method based on a TPS model is established for each group of homonymous points, and a correction amount of geometric deviation between the images is obtained by iterative solution, which specifically includes:
4.1) determining an overlapping area between the images according to the range information of the orthoimages obtained in the step 3), wherein the overlapping area of the images is an intersection of rectangular ranges of geographic coordinates, and homonymous points are obtained in the overlapping area;
4.2) reading the coordinates and the resolution of the upper left corner point of the DOM header file, solving the coordinate solution of the same name point obtained by matching to obtain the ground coordinates:
4.3) establishing an orthonormal adjustment error equation based on the TPS model, and solving by adopting a least square principle to obtain geometric refinement parameters among orthonormal images.
Further, in the step 4.3), an orthonormal adjustment error equation based on the TPS model is established, and a least square principle is used for solving to obtain geometric refinement parameters between orthonormal images, and the specific process is as follows:
4.3.1) modeling TPS:
Figure BDA0001952915370000051
Figure BDA0001952915370000052
in the formula, ax 0,ax 1,ax 2,by 0,by 1,by 2i(i ═ 1,2, … n) are TPS model parameters, x, y are the original coordinates, x ', y' are the corrected coordinates, n is the number of homologous points, ψ (r)i) For radial basis functions:
Figure BDA0001952915370000053
Figure BDA0001952915370000054
in the formula, r is the Euclidean distance from any point in the image to the point matched with the same name;
4.3.2) the error equation based on the TPS model finds the homonym to the object for each solution, wherein the error equation based on the TPS model:
Figure BDA0001952915370000055
Figure BDA0001952915370000056
in the formula (I), the compound is shown in the specification,
Figure BDA0001952915370000057
are respectively the coordinate residuals in the row and column directions,
Figure BDA0001952915370000058
is the original object space coordinate, and the object space coordinate,
Figure BDA0001952915370000059
i represents a point label, j represents an image number, and n represents the total number of points with the same name on the image;
expressed in matrix form as follows:
v=Ax-L
Figure BDA00019529153700000510
Figure BDA00019529153700000511
L=[Dx1 Dy1]T
4.3.3) setting the error equation shown above for each homonymous point on each image, and solving according to the least square principle to obtain the coefficient of the TPS model;
4.4) correcting and resampling each image according to the TPS model obtained by solving in the 4.3) to obtain a corrected orthoimage.
Further, the step 4.4) corrects and resamples each image by adopting a positive solution method or a reverse solution method; the specific process of adopting a positive solution method for correcting and resampling each image comprises the following steps:
4.4.1) determining the ground coordinates and the resolution of the upper left corner point of the corrected orthoimage:
Figure BDA0001952915370000061
Figure BDA0001952915370000062
in the formula (X)0,Y0) Is the ground coordinate of the upper left corner point of the original ortho-image, (X)0',Y0') is the top left corner point ground coordinates of the corrected orthographic image;
4.4.2) for each point (x) on the original ortho-imagedom,ydom) It is converted into ground coordinates (X, Y).
4.4.3) carrying out coordinate correction on the X and the Y to obtain corrected coordinates X 'and Y';
4.4.4) converting the corrected ground coordinates to the image coordinates (x) of the corrected DOMdom',ydom'):
xdom'=(X'-X0')/mx'
ydom'=(Y'-Y0')/my'
In the formula (m)x',my') resolution of corrected DOM, (X)0’,Y0') the top left corner point ground coordinates of the corrected DOM;
4.4.5) mixing (x)dom,ydom) Is given by (x)dom',ydom') finishing the correction process pixel by pixel, and finally performing gray resampling on the DOM to obtain the corrected orthoimage.
Due to the adoption of the technical scheme, the invention has the following advantages:
1. according to the large-area seamless ortho-image system, the original remote sensing image and the DEM are registered, and then the problem that the geometric accuracy between the ortho-images cannot meet the requirement due to insufficient control data accuracy is solved, so that the resolution of unknowns is decomposed, the number of unknowns solved at a time is reduced, the adjustment iteration is easy to converge, and the adjustment accuracy is improved;
2. the invention provides an orthoscopic image-based adjustment method, an original image is projected by using an imaging geometric model to achieve the effect of dimension reduction, and the problems of unstable adjustment and ill-condition solution caused by small intersection angle and insufficient geometric strength in the traditional adjustment method are solved.
3. According to the method, the thin-plate spline model is introduced in the process of the orthoimage adjustment, so that the problem of low precision of the traditional geometric correction model caused by a small overlapping range between images is solved;
4. the invention provides a method for directly matching a DEM (digital elevation model) and an image, which solves the problem that the prior DOM (document object model) assistance is required in the process of matching the DEM and the image, and simultaneously obtains the refinement parameters of an image imaging geometric model in the process of matching the DEM and the image, thereby improving the precision and quality of the image orthorectification.
Drawings
FIG. 1 is a schematic flow chart of a method for making a large-area seamless orthographic image of a satellite remote sensing image according to the invention;
fig. 2 is a schematic view of a 3 × 3 window according to an embodiment of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some, but not all, embodiments of the present invention. All other embodiments, which can be obtained by a person skilled in the art without any inventive step based on the embodiments of the present invention, are within the scope of the present invention.
As shown in fig. 1, the method for making a large-area seamless orthoimage of a satellite remote sensing image provided by the invention specifically comprises the following steps:
1. construction of remote sensing image imaging geometric model
The imaging geometric models of the remote sensing images commonly used at present are mainly divided into strict imaging geometric models and general imaging geometric models. The strict imaging geometric model is a mathematical model with strict theoretical basis, is mainly based on collinear equation, and can accurately express the strict geometric relation between image coordinates and ground point space coordinates. The general imaging geometric Model avoids the complex relation of the imaging process, and adopts a mathematical Model to fit the mutual relation between the image point coordinates and the three-dimensional coordinates of the object point, and the common fitting models comprise a general polynomial Model, a direct linear transformation Model, a Rational Function Model (RFM) and the like, wherein the Rational Function Model becomes the most widely applied mathematical Model in the general geometric Model of the remote sensing image due to the advantages of high fitting precision, good universality, convenient application and the like. The present invention may employ any rigorous imaging geometry model and general imaging geometry model.
The invention takes the imaging geometric model construction of the satellite image with the highest resolution ratio, namely the LRO NAC image, in the moon as an example, and the construction process of the strict imaging geometric model and the universal geometric model is explained in detail. The construction of the rigid imaging geometric model generally comprises two processes of inner orientation and outer orientation, while the construction of the universal geometric model needs to be based on the constructed rigid imaging geometric model, and the specific processes are as follows:
1.1) construction of LRO NAC rigid imaging geometry model
1.1.1) orientation within LRO NAC
Acquiring inner orientation parameters of the NAC camera from the IK auxiliary file of the LRO, for example: focal length, rank direction central coordinate, pixel size and distortion parameter etc. then according to LRO NAC's distortion model, carry out interior orientation to the NAC camera, wherein, LRO NAC's distortion model is:
xd=(sample-BORESIGHT_SAMPLE)*PIXEL_PITCH
r=xd
xc=xd/(1+k1*r2) (1)
in the formula, SAMPLE is a column coordinate of an image point on NAC EDR raw data, constraint _ SAMPLE is a center coordinate of a column direction, PIXEL _ PITCH is a PIXEL size of the column direction, xd is a coordinate (measurement coordinate) including distortion, k1 is a radial distortion parameter, r is a distance from the image point to a main point, and xc is a coordinate of the corrected image point on a focal plane, and a unit is mm. Since NAC is a CCD line scan camera, the same kind of parameter yd in the row direction is 0 and yc is 0.
1.1.2) LRO NAC Exo-orientation
(a) Establishing collinearity equations
After the internal orientation is finished, the coordinates of each pixel after distortion correction on a focal plane can be obtained, the relationship between a focal plane coordinate system and a star-fixed coordinate system is established by the external orientation, and the strict imaging geometric model can be expressed by a collinear equation:
Figure BDA0001952915370000081
wherein (xc, yc) is focal plane coordinates of an image point, f is focal length, (X, Y, Z) is coordinates of a corresponding object point in a star-solid coordinate system, (Xs, Ys, Zs) is coordinates of a photographing center in the star-solid coordinate system, and is a line element called an exterior orientation element, λ is a scale factor, and R is a rotation matrix from an image space coordinate system to the star-solid coordinate system, and is composed of three exterior orientation angle elements
Figure BDA0001952915370000084
And (4) forming.
(b) Initial exterior orientation element reading
To orient the image externally, the external orientation element at the time of imaging needs to be acquired first. The foreign orientation elements are obtained from aircraft position and attitude data obtained by orbit measurement, and the measured data is stored as auxiliary data in a SPICE kernel file of LRO NAC images, so that the foreign orientation elements of each image can be read from the corresponding SPICE kernel.
(c) Interpolating exterior orientation elements for each scan line
For a push-broom imaged orbiter image, each scan line has a corresponding exterior orientation element. However, the satellite orbit measurement time interval is larger than the scanning and imaging time interval of each line of image, and an interpolation mode is needed to obtain the external orientation element of each scanning line. An external orientation element function relative to the imaging time t is generally established by adopting a third-order polynomial, and the external orientation element of each scanning line can be obtained by interpolation according to the recorded CCD imaging time of each line.
Figure BDA0001952915370000082
In the formula, Xs(t),Ys(t),Zs(t) coordinates of the photographing center in the star-solid coordinate system at time t, i.e., elements of an external orientation line;
Figure BDA0001952915370000083
representing the attitude angle of the focal plane in the star-solid coordinate system at the moment t, namely an external orientation angle element; a isX 0……aκ 3Polynomial coefficients representing the corresponding parameters may be solved from the orbit measurement data using a least squares method.
(d) And the focal plane coordinate subjected to distortion correction can be converted into an object coordinate through a collinear equation and the obtained external orientation element, so that the establishment of a strict imaging geometric model of the sensor is completed.
1.2) establishment of LRO NAC rational function model
The establishment of the LRO NAC rational function model requires that a virtual control grid is established first, and then rational function model parameters are solved according to the generated virtual control points.
1.2.1) establishment of virtual control grids
When a virtual control grid is established, the elevation of an image area needs to be divided into a plurality of elevation surfaces, grid point coordinates of an image are generated at a certain distance in an image space and serve as image space virtual control points, and then grid point image coordinates are projected onto each elevation surface according to a strict geometric model to obtain object space virtual control points.
1.2.2) solution of rational function model parameters
The rational polynomial model establishes a one-to-one correspondence relationship between arbitrary ground point coordinates (lat, lon, h) and corresponding image coordinates (sample, line) by a ratio polynomial, and the expression form is as follows:
Figure BDA0001952915370000091
Figure BDA0001952915370000092
in the formula (I), the compound is shown in the specification,
NumL(P,L,H)=a1+a2L+a3P+a4H+a5LP+a6LH+a7PH+a8L2+a9P2
+a10H2+a11PLH+a12L3+a13LP2+a14LH2+a15L2P+a16P3+a17PH2
+a18L2H+a19P2H+a20H3
DenL(P,L,H)=b1+b2L+b3P+b4H+b5LP+b6LH+b7PH+b8L2+b9P2
+b10H2+b11PLH+b12L3+b13LP2+b14LH2+b15L2P+b16P3+b17PH2
+b18L2H+b19P2H+b20H3
Nums(P,L,H)=c1+c2L+c3P+c4H+c5LP+c6LH+c7PH+c8L2+c9P2
+c10H2+c11PLH+c12L3+c13LP2+c14LH2+c15L2P+c16P3+c17PH2
+c18L2H+c19P2H+c20H3
Dens(P,L,H)=d1+d2L+d3P+d4H+d5LP+d6LH+d7PH+d8L2+d9P2
+d10H2+d11PLH+d12L3+d13LP2+d14LH2+d15L2P+d16P3+d17PH2
+d18L2H+d19P2H+d20H3
wherein, ai,bi,ci,di(i is 1-20) is a rational function model parameter, b1And d1Typically 1, (P, L, H) normalized ground coordinates, (X, Y) normalized image coordinates,
the normalization is as follows:
Figure BDA0001952915370000101
Figure BDA0001952915370000102
in the formula, LINE _ SCALE, same _ OFF and LINE _ OFF are normalization parameters of image side coordinates; LAT _ OFF, LON _ OFF, HEIGHT _ OFF, LAT _ SCALE, LON _ SCALE, HEIGHT _ SCALE are object coordinate normalization parameters, LAT is latitude, LON is longitude, and h is elevation.
1.2.3) solving 78 rational function model parameters by least squares according to the virtual control points obtained in 1.2.1), and establishing a rational function model of each image.
2. Reading the remote sensing image to be matched and DEM data, and matching the DEM and the remote sensing image data to be registered. The matching method of the DEM and the image is explained by adopting a DEM and image matching method constrained by geometric and radiation information together, but the method is not limited to the method, and the specific process of the method is as follows:
2.1) generating a simulation image by the DEM to be used as a reference image.
2.1.1) on DEM, the gradient of each pixel in x and y directions is calculated pixel by pixel:
Figure BDA0001952915370000103
Figure BDA0001952915370000104
in the formula, e1,e2…e9In order to calculate the DEM pixel as the center, each pixel value corresponding to a 3 × 3 window is illustrated as a 3 × 3 window diagram in fig. 2.
2.1.2) calculating the gradient and the slope direction value pixel by pixel according to the gradient value:
Figure BDA0001952915370000105
A=p/q (8)
wherein S is a gradient value and A is a gradient value.
2.1.3) obtaining a solar azimuth angle (Azi) and a solar altitude angle (Alt) when the remote sensing image is obtained according to the header file information of the remote sensing image data, creating a null image with the same size as the DEM, calculating the corresponding simulated radiation value information pixel by pixel, storing the simulated radiation value information in the null image, and constructing to obtain a reference image.
I=255×(cosAlt×cosS+sinAlt×sinS×cosAzi-A) (9)
In the formula, I is a simulated radiation value of each pixel of the simulated reference image obtained by calculation, Alt is a solar altitude angle when the remote sensing image is acquired, and Azi is a solar azimuth angle when the remote sensing image is acquired.
And 2.2) obtaining initial matching points of the reference image and the remote sensing image by using an affine-scale invariant feature transform (ASIFT) algorithm.
An affine-scale invariant feature transform (ASIFT) algorithm is a common algorithm in the field of remote sensing technology, and the algorithm can obtain an initial matching point between a remote sensing image and a reference image, and specifically comprises the following steps:
firstly, acquiring a sampling sequence of a longitude angle and a latitude angle of a camera, and simulating affine deformation possibly caused by an optical axis of the camera to realize image transformation; then, carrying out oblique rotation transformation on the remote sensing image to be matched to generate a simulated image; and finally, performing feature point detection and reference image matching on the generated simulated image according to a Scale Invariant Feature Transform (SIFT) algorithm, and finally obtaining an initial matching point between the remote sensing image and the reference image.
And 2.3) calculating initial refinement parameters of the rational function model by using the initial matching points.
The invention adopts affine transformation of image space to correct the front and back projection error, and introduces refinement parameters (column parameters and row parameters) to correct the error:
Figure BDA0001952915370000111
Figure BDA0001952915370000112
wherein (r, c) is the coordinates of image points, (r ', c') is the projection coordinates of ground points back projected to the image side by RFM,
Figure BDA0001952915370000113
the parameters are refined for the RFM.
2.4) obtaining an initial value of the least square fine matching.
2.4.1) feature points are extracted on the reference image using feature point extraction operators, such as Harris, Forstner, etc.
2.4.2) calculating the characteristic points on the reference image to the remote sensing image by using the refined rational function model to be used as the initial value of least square matching.
Back projecting the feature points extracted in the step 2.4.1) to an image space through the rational function model obtained in the step 1 and the refined parameters calculated in the step 2.3), and obtaining the coordinates of the feature points on the image as initial values of global least square matching:
(a) calculating the reference image feature points extracted in the step 2.4.1) to geodetic coordinates according to header file information of the reference image:
Y=Y0+fy×Ry
X=X0+fx×Rx (11)
wherein Y and X are respectively latitude and longitude coordinates of the feature point, Y0And X0Latitude and longitude coordinates, f, of the upper left corner of the reference imagey,fxIs the row number, R, of the feature point on the reference imagey,RxThe resolution of the reference image in the y-direction and the x-direction is shown.
And according to the calculated longitude and latitude coordinates, extracting an elevation value Z corresponding to the longitude and latitude from the digital elevation model.
(b) Calculating image point coordinates (r, c) of the feature points on the remote sensing image to be registered through (X, Y, Z) according to the RPC and the refinement parameters calculated in the step 1 and the step 2.3) and the formula (4) and the formula (10);
2.5) based on the least square matching of global constraint, the radiation deformation parameters and the geometric deformation parameters of the remote sensing images to be matched are output, and the specific process is as follows:
2.5.1) solving registration parameters of the remote sensing image and the reference image, simultaneously solving refinement parameters of a rational function model, and introducing a global constraint least square matching model as follows:
Figure BDA0001952915370000121
in the formula, g1And g2Is a function of the radiation value of the corresponding coordinate position of the reference image and the remote sensing image connecting point to be registered, n1And n2Is a function of the noise of the reference image and the remote sensing image to be registered, (x, y) is the coordinate of the connecting point on the reference image, h0,h1Is a radiation deformation parameter between the reference image and the remote sensing image, ar 0……ar 2,bc 0……bc 2The geometric deformation parameters (refinement parameters), (r ', c') are projection coordinates of ground points back projected to the image side by RFM, (r, c) are coordinates of image points obtained by refinement of the projection coordinates back projected to the image side, and deltarcIs the model error.
2.5.2) establishing an error equation according to the global constraint least square matching model as follows:
Figure BDA0001952915370000122
in the formula, vgError of radiation value between registered images, vrError of equation in column direction, vcIs the error of the equation in the row direction.
By performing a first order taylor expansion on equation (13), a linearized error equation can be obtained:
Figure BDA0001952915370000123
Figure BDA0001952915370000124
in the formula (I), the compound is shown in the specification,
Figure BDA0001952915370000125
is the 0-order term of the Taylor series,
Figure BDA0001952915370000126
representing the first partial derivatives of the column-wise and row-wise coordinates, respectively.
Inputting the least square matching initial value obtained by calculation in the step 2.2) and all points in a window with the size of n multiplied by n by taking the initial value point as the center in the process of solving the matching parameters and the refinement parameters, and taking the points as input points, wherein the value of n is 9. And (3) establishing an error equation by using a linearized error equation (14) for all input points, establishing an error equation by using a linearized error equation (15) for the initial value points calculated in the step 2.2), and finally uniformly solving all established error equations to obtain matching parameters and rational function model refinement parameters.
2.5.3) calculating a normal equation point by point according to an error equation in the formula (15) for each input matching initial value point.
2.5.4) solving a normal equation, correcting the unknown number, judging correlation coefficients obtained by calculation after adjustment, for example, if the sum of the correlation coefficients of all points meets the requirement (the sum of the correlation coefficients of the current time is less than the sum of the correlation coefficients of the last iteration), outputting an adjustment result, and if the adjustment precision does not meet the requirement, repeating the steps 2.5.2),2.5.3) and 2.5.4), and performing iterative solution, wherein the correlation coefficients rho2
Figure BDA0001952915370000131
2.5.5) outputting the radiation deformation parameters and the geometric deformation parameters of the remote sensing image, and completing the global least square matching between the remote sensing image and the reference image.
3. Remote sensing image ortho-rectification
By using the geometric deformation parameters obtained by the solution in the step 2 and the image imaging geometric model established in the step 1, orthorectification can be performed on each image to obtain a corresponding DOM product. The image orthorectification is generally divided into a positive solution method and a reverse solution method, and the invention describes the specific implementation process by taking the reverse solution method as an example:
3.1) calculating the coordinates of the ground points.
First, X-direction and Y-direction resolutions (R) for an ortho-image to be generated are setxdom,Rydom) The resolution can be set according to the requirement, and can also be set according to the original resolution of the remote sensing image. Meanwhile, the coverage area of the image is calculated by using the size of the original image and the geometric model of the image. Setting up images according to their extentCoordinates of point at upper left corner (X)0,Y0). The coordinate of an image point of any point pixel center P on the preset generated DOM is set as (x)dom,ydom) Calculating the ground coordinates (X, Y) corresponding to the P point;
X=X0+Rxdom*xdom
Y=Y0+Rydom*ydom
3.2) calculating the coordinates of the image points
Interpolating the elevation Z of the point by using (X, Y) and DEM, and calculating the corresponding image point coordinate p (X, Y) on the original image by using the established and refined image imaging geometric model;
x=f1(X,Y,Z)
y=f2(X,Y,Z)
in the formula, f1 and f2 represent the refined geometric model of the imaged image, and Z represents the elevation of the point P.
3.3) Gray interpolation
Because the obtained pixel coordinate p (x, y) does not necessarily fall at the center of the pixel of the original image, gray scale interpolation is necessary for this purpose, and generally, bilinear interpolation can be adopted to obtain the gray scale value g (x, y) of the pixel p;
3.4) Gray scale assignment
Assigning the gray value of the image point P to the corrected pixel, namely the P point on the DOM;
and 3.5) completing the operation on each pixel on the DOM in sequence to obtain the DOM product subjected to terrain correction.
4. Orthonormal adjustment based on Thin Plate Spline (TPS) model: determining the overlapping area between the images according to the position relation between the orthoimages obtained in the step 3, matching the overlapping area to obtain the homonymous points, then establishing an error method based on a TPS (Perfect spline) model for each group of homonymous points, obtaining the correction quantity of the geometric deviation between the images through iterative solution, and further carrying out geometric position refinement on the orthoimages to ensure that the images achieve a seamless effect, wherein the specific process comprises the following steps:
4.1) determining the overlapping area between the images according to the range information of the orthoimages obtained in the step 3, wherein the overlapping area of the images is the intersection of the geographic coordinate rectangular ranges. And obtaining the homonymy points in the overlapping area by using matching methods such as least square matching, SIFT, SURF, A-SIFT and the like.
4.2) reading the coordinates and the resolution of the upper left corner point of the DOM header file, and solving the coordinate solution of the matched point with the same name according to the formula (11) to obtain the ground coordinates:
4.3) establishing an orthonormal adjustment error equation based on the TPS model, and solving by adopting a least square principle to obtain geometric refinement parameters among orthonormal images, wherein the geometric refinement parameters specifically comprise:
4.3.1) modeling TPS:
Figure BDA0001952915370000141
Figure BDA0001952915370000142
in the formula, ax 0,ax 1,ax 2,by 0,by 1,by 2i(i ═ 1,2, … n) are TPS model parameters, x, y are the original coordinates, x ', y' are the corrected coordinates, n is the number of homologous points, ψ (r)i) For radial basis functions:
Figure BDA0001952915370000143
Figure BDA0001952915370000144
in the formula, r is the Euclidean distance from any point in the image to the point matching the same name.
4.3.2) the TPS model based error equation solves for each solution to the homonym of the object, where the TPS model based error equation is as follows:
Figure BDA0001952915370000145
Figure BDA0001952915370000146
in the formula (I), the compound is shown in the specification,
Figure BDA0001952915370000147
are respectively the coordinate residuals in the row and column directions,
Figure BDA0001952915370000148
is the original object space coordinate, and the object space coordinate,
Figure BDA0001952915370000149
in the refined object coordinates, i represents a point number, j represents an image number, and n represents the total number of points of the same name on the image.
Expressed in matrix form as follows:
v=Ax-L
Figure BDA0001952915370000151
Figure BDA0001952915370000152
L=[Dx1 Dy1]T
4.3.3) for each same-name point on each image, an error equation shown above is listed, and the coefficient of the TPS model is obtained by solving according to the least square principle.
4.4) correcting and resampling each image according to the TPS model obtained by the solution in 4.3), wherein the solution process is similar to the process of generating the orthoimage in step 3, and can be divided into a positive solution method and a reverse connection method, and the positive solution method is taken as an example for explanation:
4.4.1) determining the ground coordinates and the resolution of the upper left corner point of the corrected orthoimage, wherein the default is generally the same as the resolution of the image before correction, and the coordinates of the upper left corner point are obtained by transforming the coordinates of the upper left corner point of the original DOM according to TPS model parameters:
Figure BDA0001952915370000153
Figure BDA0001952915370000154
in the formula (X)0,Y0) The ground coordinates of the upper left corner point of the original ortho-image can be read from the DOM header file (X)0',Y0') is the top left corner point ground coordinates of the corrected orthographic image.
4.4.2) for each point (x) on the original ortho-imagedom,ydom) It is converted into ground coordinates (X, Y) according to equation (11) using its header file information.
4.4.3) carrying out coordinate correction on the X and the Y to obtain corrected coordinates X 'and Y';
4.4.4) converting the corrected ground coordinates to the image coordinates (x) of the corrected DOMdom',ydom'):
xdom'=(X'-X0')/mx'
ydom'=(Y'-Y0')/my'
In the formula (m)x',my') resolution of corrected DOM, (X)0’,Y0') is the top left corner point ground coordinates of the corrected DOM.
4.4.5) mixing (x)dom,ydom) Is given by (x)dom',ydom') finishing the correction process pixel by pixel, and finally performing gray resampling on the DOM to obtain a corrected orthoimage.
5. And (4) carrying out light evening and inlaying on all corrected orthographic images to obtain an orthographic image map product of a large area. The map product can be used for topographic and geomorphic analysis, surface feature extraction, control reference of low-resolution images and the like.
Appropriate changes and modifications to the embodiments described above will become apparent to those skilled in the art from the disclosure and teachings of the foregoing description. Therefore, the present invention is not limited to the specific embodiments disclosed and described above, and some modifications and variations of the present invention should fall within the scope of the claims of the present invention. Furthermore, although specific terms are employed herein, they are used in a generic and descriptive sense only and not for purposes of limitation.

Claims (8)

1. A method for making a large-area seamless orthoimage of a satellite remote sensing image is characterized by comprising the following steps:
1) constructing a remote sensing image imaging geometric model;
2) reading the remote sensing image to be matched and DEM data, and matching the DEM and the remote sensing image data to be registered;
3) orthorectifying the remote sensing image to obtain a DOM product corresponding to the primary rectification;
4) obtaining the position relation between the ortho images according to the step 3), determining the overlapping area between the images, obtaining the homonymy points in the overlapping area through matching, establishing an error equation based on a TPS model for each group of homonymy points, obtaining the correction quantity of the geometric deviation between the images through iterative solution, and further correcting the DOM preliminarily corrected in the step 3) by using the correction quantity of the geometric deviation between the images obtained through the solution to obtain the corrected DOM, wherein the specific process is as follows:
4.1) determining an overlapping area between the images according to the range information of the orthoimages obtained in the step 3), wherein the overlapping area of the images is an intersection of rectangular ranges of geographic coordinates, and homonymous points are obtained in the overlapping area;
4.2) reading the coordinates and the resolution of the upper left corner point of the DOM header file, solving the coordinate solution of the same name point obtained by matching to obtain the ground coordinates:
4.3) establishing an orthonormal adjustment error equation based on the TPS model, and solving by adopting a least square principle to obtain geometric refinement parameters among orthonormal images, wherein the specific process is as follows:
4.3.1) modeling TPS:
Figure FDA0002694248420000011
Figure FDA0002694248420000012
in the formula, ax 0,ax 1,ax 2,by 0,by 1,by 2i(i ═ 1,2, … n) are TPS model parameters, x, y are the original coordinates, x ', y' are the corrected coordinates, n is the number of homologous points, ψ (r)i) For radial basis functions:
Figure FDA0002694248420000013
Figure FDA0002694248420000014
in the formula, r is the Euclidean distance from any point in the image to the point matched with the same name;
4.3.2) the error equation based on the TPS model finds the homonym to the object for each solution, wherein the error equation based on the TPS model:
Figure FDA0002694248420000015
Figure FDA0002694248420000016
in the formula (I), the compound is shown in the specification,
Figure FDA0002694248420000017
are respectively the coordinate residuals in the row and column directions,
Figure FDA0002694248420000018
is the original object space coordinate, and the object space coordinate,
Figure FDA0002694248420000019
i represents a point label, j represents an image number, and n represents the total number of points with the same name on the image;
expressed in matrix form as follows:
v=Ax-L
Figure FDA0002694248420000021
Figure FDA0002694248420000022
L=[Dx1 Dy1]T
4.3.3) setting the error equation shown above for each homonymous point on each image, and solving according to the least square principle to obtain the coefficient of the TPS model;
4.4) correcting and resampling each image according to the TPS model obtained by solving in the step 4.3) to obtain a corrected orthoimage;
5) and (4) carrying out light evening and embedding on all corrected orthographic images to obtain an orthographic image map product of a large area.
2. The method for making the satellite remote sensing image large-area seamless orthoimage according to claim 1, wherein the DEM and the data of the remote sensing image to be registered are matched by adopting the DEM and the image which are constrained by geometric and radiation information together, and the specific process is as follows:
2.1) generating a simulation image by the DEM to be used as a reference image;
2.2) obtaining initial matching points of the reference image and the remote sensing image by using an affine-scale invariant feature transformation algorithm;
2.3) calculating initial refinement parameters of the rational function model by using the initial matching points;
2.4) obtaining an initial value of least square fine matching;
and 2.5) outputting the radiation deformation parameters and the geometric deformation parameters of the remote sensing images to be matched based on the least square matching of the global constraint.
3. The method for manufacturing the large-area seamless orthoimage of the satellite remote sensing image according to claim 2, wherein the step 2.2) of obtaining the initial matching point of the reference image and the remote sensing image by using an affine-scale invariant feature transformation algorithm comprises the following specific steps:
firstly, acquiring a sampling sequence of a longitude angle and a latitude angle of a camera, and simulating affine deformation possibly caused by an optical axis of the camera to realize image transformation;
then, carrying out oblique rotation transformation on the remote sensing image to be matched to generate a simulated image;
and finally, performing feature point detection and reference image matching on the generated simulated image according to a scale-invariant feature transformation algorithm to obtain an initial matching point between the remote sensing image and the reference image.
4. The method for making the satellite remote sensing image large-area seamless orthoimage according to claim 2, wherein the specific process of obtaining the initial value of least square fine matching in the step 2.4) is as follows:
2.4.1) extracting the feature points on the reference image by using a feature point extraction operator;
2.4.2) calculating the characteristic points on the reference image to the remote sensing image by using the refined rational function model to be used as the initial value of least square matching.
5. The method for making the satellite remote sensing image large-area seamless orthoimage according to claim 2, wherein the step 2.5) is based on global constraint least square matching, and outputs the radiation deformation parameters and the geometric deformation parameters of the remote sensing image to be matched, and the specific process is as follows:
2.5.1) solving registration parameters of the remote sensing image and the reference image, simultaneously solving refinement parameters of a rational function model, and introducing a global constraint least square matching model as follows:
g1(x,y)+n1(x,y)=h0+h1g2(r,c)+n2(x,y)
r=ar 0+ar 1r′+ar 2c′+δr
c=bc 0+bc 1r′+bc 2c′+δc
in the formula, g1And g2Is a function of the radiation value of the corresponding coordinate position of the reference image and the remote sensing image connecting point to be registered, n1And n2Is a function of the noise of the reference image and the remote sensing image to be registered, (x, y) is the coordinate of the connecting point on the reference image, h0,h1Is a radiation deformation parameter between the reference image and the remote sensing image, ar 0……ar 2,bc 0……bc 2The geometric deformation parameters (refinement parameters), (r ', c') are projection coordinates of ground points back projected to the image side by RFM, (r, c) are coordinates of image points obtained by refinement of the projection coordinates back projected to the image side, and deltarcIs the model error;
2.5.2) establishing an error equation according to the global constraint least square matching model:
vg=h0+h1g2(r,c)-g1(x,y)
vr=ar 0+ar 1r′+ar 2c′+δr-r
vc=bc 0+bc 1r′+bc 2c′+δc-c
in the formula, vgError of radiation value between registered images, vrIs a squareError of course in column direction, vcError of the equation in the row direction;
and performing first-order Taylor expansion on the error equation to obtain a linearized error equation:
Figure FDA0002694248420000031
vr=vr0+Δar 0+r′·Δar 1+c′·Δar 2+Δδr
vc=vc0+Δbc 0+r′·Δbc 1+c′·Δbc 2+Δδc
in the formula (I), the compound is shown in the specification,
Figure FDA0002694248420000033
is the 0-order term of the Taylor series,
Figure FDA0002694248420000032
respectively representing the first partial derivatives of the coordinates in the column direction and the row direction;
2.5.3) calculating a normal equation point by point according to an error equation for each input matching initial value point;
2.5.4) solving a normal equation, correcting the unknown number, judging the correlation coefficient obtained by calculation after adjustment, if the sum of the correlation coefficients of all the points meets the requirement, outputting an adjustment result, and otherwise, entering the step 2.5.2 to carry out iterative solution;
2.5.5) outputting the radiation deformation parameters and the geometric deformation parameters of the remote sensing image, and completing the global least square matching between the remote sensing image and the reference image.
6. The method for making a large-area seamless ortho image of a satellite remote sensing image according to claim 1, wherein the ortho correction of the remote sensing image in the step 3) adopts a forward solution method and a reverse solution method.
7. The method for making the satellite remote sensing image large-area seamless orthoimage according to claim 6, wherein the specific process of orthorectifying the remote sensing image by adopting the inverse solution method comprises the following steps:
3.1) setting the X-direction and Y-direction resolutions (R) for generating an orthoimagexdom,Rydom) Calculating the coverage of the image by using the size of the original image and the geometric model of the image, and setting the coordinates (X) of the upper left corner point of the image according to the range of the image0,Y0) The coordinate of the pixel point of the pixel center P of any point on the preset generated DOM is set as (x)dom,ydom) Calculating the ground coordinates (X, Y) corresponding to the P point;
X=X0+Rxdom*xdom
Y=Y0+Rydom*ydom
3.2) interpolating the elevation Z of the point by using (X, Y) and DEM, and calculating the corresponding image point coordinate p (X, Y) on the original image by using the established and refined image imaging geometric model;
x=f1(X,Y,Z)
y=f2(X,Y,Z)
wherein f1 and f2 represent the refined geometric model of the imaged image, and Z represents the elevation of the point P;
3.3) gray scale interpolation, obtaining the gray scale value g (x, y) of the image point p;
3.4) assigning the gray value of the image point P to the corrected image element, namely the P point on the DOM;
and 3.5) completing the operation on each pixel on the DOM in sequence to obtain the DOM product subjected to terrain correction.
8. The method for making the seamless orthoimage of the satellite remote sensing image large area according to claim 1, wherein the step 4.4) corrects and resamples each image by adopting a forward solution method or a reverse solution method; the specific process of adopting a positive solution method for correcting and resampling each image comprises the following steps:
4.4.1) determining the ground coordinates and the resolution of the upper left corner point of the corrected orthoimage:
Figure FDA0002694248420000041
Figure FDA0002694248420000042
in the formula (X)0,Y0) Is the ground coordinate of the upper left corner point of the original ortho-image, (X)0',Y0') is the top left corner point ground coordinates of the corrected orthographic image;
4.4.2) for each point (x) on the original ortho-imagedom,ydom) It is converted into ground coordinates (X, Y).
4.4.3) carrying out coordinate correction on the X and the Y to obtain corrected coordinates X 'and Y';
4.4.4) converting the corrected ground coordinates to the image coordinates (x) of the corrected DOMdom',ydom'):
xdom'=(X'-X0')/mx'
ydom'=(Y'-Y0')/my'
In the formula (m)x',my') resolution of corrected DOM, (X)0’,Y0') the top left corner point ground coordinates of the corrected DOM;
4.4.5) mixing (x)dom,ydom) Is given by (x)dom',ydom') finishing the correction process pixel by pixel, and finally performing gray resampling on the DOM to obtain the corrected orthoimage.
CN201910057228.7A 2018-12-24 2019-01-22 Method for making large-area seamless orthoimage of satellite remote sensing image Active CN109903352B (en)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN2018115805425 2018-12-24
CN201811580542 2018-12-24

Publications (2)

Publication Number Publication Date
CN109903352A CN109903352A (en) 2019-06-18
CN109903352B true CN109903352B (en) 2021-03-30

Family

ID=66943936

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910057228.7A Active CN109903352B (en) 2018-12-24 2019-01-22 Method for making large-area seamless orthoimage of satellite remote sensing image

Country Status (1)

Country Link
CN (1) CN109903352B (en)

Families Citing this family (35)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110503604B (en) * 2019-07-31 2022-04-29 武汉大学 Aviation area array image real-time orthotropic splicing method based on high-precision POS
CN110555817B (en) * 2019-09-10 2022-06-24 中国科学院遥感与数字地球研究所 Geometric normalization method and device for remote sensing image
CN110853140A (en) * 2019-10-11 2020-02-28 北京空间机电研究所 DEM (digital elevation model) -assisted optical video satellite image stabilization method
CN112750075A (en) * 2019-10-31 2021-05-04 中国科学院长春光学精密机械与物理研究所 Low-altitude remote sensing image splicing method and device
CN110763205B (en) * 2019-11-05 2021-07-02 新疆维吾尔自治区测绘科学研究院 Method for generating orthophoto map of border narrow and long area by digital photogrammetric system
CN110929782B (en) * 2019-11-20 2021-05-18 天津大学 River channel abnormity detection method based on orthophoto map comparison
CN111156969A (en) * 2020-02-07 2020-05-15 黄文超 Wide remote sensing image stereo mapping method and system
CN111611525B (en) * 2020-05-14 2022-07-29 中国科学院空天信息创新研究院 Remote sensing data elevation calculation method based on object space matching elevation deviation iterative correction
CN111754458B (en) * 2020-05-18 2023-09-15 北京吉威空间信息股份有限公司 Satellite image three-dimensional space reference frame construction method for geometric fine processing
CN111861934A (en) * 2020-07-29 2020-10-30 贵阳欧比特宇航科技有限公司 Hyperspectral satellite image data production, mosaic and metadata manufacturing method
CN112017108B (en) * 2020-08-07 2022-04-15 武汉大学 Satellite image color relative correction method based on adjustment of independent model method
CN111982015B (en) * 2020-08-18 2022-06-21 深圳大学 Method for monitoring structure geometry
CN112016178A (en) * 2020-09-02 2020-12-01 浙江清环智慧科技有限公司 Method for building pipe network model by aid of online monitoring network data
CN112132029B (en) * 2020-09-23 2023-07-11 中国地震局地震预测研究所 Unmanned aerial vehicle remote sensing image rapid positioning method for earthquake emergency response
CN112288650B (en) * 2020-10-28 2021-07-20 武汉大学 Multi-source remote sensing satellite image geometric and semantic integrated processing method and system
CN112561805A (en) * 2020-11-12 2021-03-26 中国资源卫星应用中心 Correction method and device for geometric rough correction image
CN112270320B (en) * 2020-11-16 2023-08-22 云南电网有限责任公司昆明供电局 Power transmission line tower coordinate calibration method based on satellite image correction
CN112422917A (en) * 2020-11-18 2021-02-26 南通市测绘院有限公司 Road accumulated water monitoring method and system based on combination of monitoring video and high-precision DEM
CN112465986A (en) * 2020-11-27 2021-03-09 航天恒星科技有限公司 Method and device for inlaying satellite remote sensing image
CN114152267A (en) * 2021-02-26 2022-03-08 武汉大学 Mars orbit camera image simulation method and system
CN112884342B (en) * 2021-03-10 2024-03-12 陕西九州遥感信息技术有限公司 Quality evaluation and cross calibration method for water-color satellite atmospheric roof radiation product
CN113160071B (en) * 2021-03-11 2023-11-07 北京师范大学 Satellite image automatic geometric correction method, system, medium and terminal equipment
CN113096047B (en) * 2021-04-25 2022-04-29 华中师范大学 Geometric fine correction method and system for generalized cloud driving and radiation cooperative remote sensing image
CN113624207A (en) * 2021-08-06 2021-11-09 苏州九宇遥感科技有限公司 Satellite image three-dimensional area network adjustment method based on orbit constraint
CN113762098A (en) * 2021-08-19 2021-12-07 北京和德宇航技术有限公司 Satellite remote sensing image matching method
CN113739767B (en) * 2021-08-24 2022-09-13 武汉大学 Method for producing orthoscopic image aiming at image acquired by domestic area array swinging imaging system
CN113724165B (en) * 2021-09-02 2022-05-24 中国科学院空天信息创新研究院 Synthetic aperture radar mosaic image color homogenizing processing method
CN114187179B (en) * 2021-12-14 2023-02-03 广州赋安数字科技有限公司 Remote sensing image simulation generation method and system based on video monitoring
CN114255051B (en) * 2021-12-21 2023-04-18 中科星通(廊坊)信息技术有限公司 Authenticity inspection method of orthometric product based on stereo mapping satellite
CN114331137B (en) * 2021-12-29 2023-04-07 中国人民解放军32021部队 Data processing method and device for equipment efficiency evaluation
CN114399541B (en) * 2021-12-29 2022-10-21 北京师范大学 Regional coordinate conversion method and device
CN114252060B (en) * 2021-12-31 2023-12-08 中铁第一勘察设计院集团有限公司 Large scene manufacturing method based on space satellite images
CN114964169B (en) * 2022-05-13 2023-05-30 中国科学院空天信息创新研究院 Remote sensing image adjustment method for image space object space cooperative correction
CN114998399B (en) * 2022-05-20 2023-04-11 中国人民解放军61540部队 Heterogeneous optical remote sensing satellite image stereopair preprocessing method
CN115423696B (en) * 2022-07-29 2024-06-18 上海海洋大学 Remote sensing orthographic image parallel generation method of self-adaptive thread parameters

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102750697B (en) * 2012-06-08 2014-08-20 华为技术有限公司 Parameter calibration method and device
CN103489176B (en) * 2012-06-13 2016-02-03 中国科学院电子学研究所 A kind of SAR image for serious geometric distortion carries out the method for same place extraction
CN106228521B (en) * 2016-07-25 2018-12-11 哈尔滨工业大学 A kind of barrier feature extracting method based on thin-plate spline interpolation
CN106548462B (en) * 2016-11-02 2019-02-19 西安电子科技大学 Non-linear SAR image geometric correction method based on thin-plate spline interpolation
CN107341778B (en) * 2017-07-10 2020-11-10 自然资源部国土卫星遥感应用中心 SAR image orthorectification method based on satellite control point library and DEM
CN108305237B (en) * 2018-01-23 2021-09-21 中国科学院遥感与数字地球研究所 Multi-stereo image fusion drawing method considering different illumination imaging conditions
CN108830889B (en) * 2018-05-24 2022-05-31 中国科学院遥感与数字地球研究所 Global geometric constraint-based remote sensing image and reference image matching method

Also Published As

Publication number Publication date
CN109903352A (en) 2019-06-18

Similar Documents

Publication Publication Date Title
CN109903352B (en) Method for making large-area seamless orthoimage of satellite remote sensing image
CN108830889B (en) Global geometric constraint-based remote sensing image and reference image matching method
CN110388898B (en) Multisource multiple coverage remote sensing image adjustment method for constructing virtual control point constraint
CN108305237B (en) Multi-stereo image fusion drawing method considering different illumination imaging conditions
US5606627A (en) Automated analytic stereo comparator
KR101965965B1 (en) A method of automatic geometric correction of digital elevation model made from satellite images and provided rpc
CN112017224B (en) SAR data area network adjustment processing method and system
CN107644435B (en) Attitude correction-considered agile optical satellite field-free geometric calibration method and system
CN110500995A (en) The method for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter
CN109709551B (en) Area network plane adjustment method for satellite-borne synthetic aperture radar image
Teo et al. DEM-aided block adjustment for satellite images with weak convergence geometry
CN103063200A (en) High-resolution optical satellite ortho-rectification image generation method
CN113900125B (en) Satellite-ground combined linear array imaging remote sensing satellite full-autonomous geometric calibration method and system
CN110006452B (en) Relative geometric calibration method and system for high-resolution six-size wide-view-field camera
CN113514829B (en) InSAR-oriented initial DSM area network adjustment method
CN113739767B (en) Method for producing orthoscopic image aiming at image acquired by domestic area array swinging imaging system
CN109029379B (en) High-precision small-base-height-ratio three-dimensional mapping method
CN114241064A (en) Real-time geometric calibration method for internal and external orientation elements of remote sensing satellite
CN110631555A (en) Historical image ortho-rectification method based on adjustment of second-order polynomial control-point-free area network
CN116753916B (en) Multi-view satellite image area network adjustment method and system
CN109579796B (en) Area network adjustment method for projected image
CN107705272A (en) A kind of high-precision geometric correction method of aerial image
CN110503604B (en) Aviation area array image real-time orthotropic splicing method based on high-precision POS
CN111275773A (en) Method and system for calibrating field-free geometry
Jacobsen Orientation of high resolution optical space images

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant