CN105374009A - Remote sensing image splicing method and apparatus - Google Patents

Abstract

The invention discloses a remote sensing image mosaic method and device, wherein the method includes: establishing a first forward calculation model and a first inverse calculation model between image points on each CCD and ground points according to the imaging parameters of each CCD , and generate the imaging parameters of the virtual CCD according to the imaging parameters of each slice of CCD, and establish the second positive calculation model and the second inverse calculation model between the image point on the virtual CCD and the ground point according to the imaging parameters of the virtual CCD; Calculation model, the first reverse calculation model, the second positive calculation model and the second reverse calculation model, establish the third normal calculation model and the third reverse calculation model between the image points on each CCD and the virtual CCD image points; read The image data of each piece of CCD, according to the 3rd forward calculation model and the 3rd inverse calculation model, determine the CCD corresponding to each image point on the virtual CCD and the image point coordinates on the CCD; Sampling processing to obtain the spliced image. Through the invention, the quality of remote sensing image mosaic is improved.

Description

Remote sensing image stitching method and device

technical field

The invention relates to the field of image processing, in particular to a remote sensing image splicing method and device.

Background technique

With the continuous improvement of the observation resolution of optical remote sensing satellites, due to the limitation of the number of pixels of a single CCD, the width of a single linear array CCD can no longer meet the observation requirements. In order to improve the efficiency of earth observation and ensure the acquisition of images with a certain width, it is an important trend in the development of spaceborne high-resolution optical cameras to use multiple CCDs to achieve a larger width through optical splicing or field of view splicing.

In related technologies, the splicing technology is to construct an image square transformation model based on connection points, so as to realize the registration of images in overlapping areas of adjacent CCDs. There are at least the following deficiencies in the image square splicing technology in the related art:

First, stitching algorithms all rely on the connection point information between slices. If the horizontal overlapping area of adjacent CCD images lacks texture features, this type of algorithm has application limitations.

Second, this type of algorithm lacks a rigorous theoretical basis as a guide, and cannot eliminate or reduce the internal image distortion caused by the geometric deformation of the sensor, etc., and cannot fundamentally improve the geometric quality of the virtual scanning scene product.

Third, the accuracy of the stitching process cannot be guaranteed for severe terrain fluctuations in the horizontal overlapping area on the original image.

Contents of the invention

Aiming at the problem of low image stitching quality in the related art, the present invention provides a remote sensing image stitching method and device to at least solve the above problems.

According to one aspect of the present invention, a remote sensing image stitching method is provided, including:

According to the imaging parameters of each sheet of CCD, the first forward calculation model and the first inverse calculation model between the image point and the ground point on the each sheet of CCD are established, and the imaging parameters of virtual CCD are generated according to the imaging parameters of each sheet of CCD. , establishing a second forward calculation model and a second inverse calculation model between the image point on the virtual CCD and the ground point according to the imaging parameters of the virtual CCD;

According to the first forward calculation model, the first inverse calculation model, the second forward calculation model and the second inverse calculation model, set up image points on each sheet of CCD and image points on the virtual CCD between the third forward model and the third inverse model;

Read the image data of each of the CCDs, determine the CCD corresponding to each image point on the virtual CCD and the image point coordinates on the CCD according to the third forward model and the third inverse model;

Resampling is performed on each determined image point to obtain a spliced image.

According to another aspect of the present invention, an image stitching device is provided, including:

The first building module is used to establish the first forward calculation model and the first inverse calculation model between the image point on the each piece of CCD and the ground point according to the imaging parameters of each piece of CCD, and according to the imaging parameters of each piece of CCD The parameters generate the imaging parameters of the virtual CCD, and set up the second forward calculation model and the second inverse calculation model between the image point and the ground point on the virtual CCD according to the imaging parameters of the virtual CCD;

The second building module is used to establish the relationship between the image points on each sheet of CCD according to the first forward calculation model, the first reverse calculation model, the second forward calculation model and the second reverse calculation model. The third forward calculation model and the third inverse calculation model between the image points on the virtual CCD;

A determining module, configured to read the image data of each of the CCDs, and determine the CCD corresponding to each image point on the virtual CCD and the CCD on the CCD according to the third forward model and the third inverse model. The image point coordinates;

The processing module is configured to perform resampling processing on each determined image point to obtain a spliced image.

Through the present invention, the image point coordinate conversion relationship between the virtual scanning scene and the original image is established, and then the seamless inner field of view splicing processing of the imaging data is realized, and the image splicing quality is improved.

Description of drawings

The accompanying drawings described here are used to provide a further understanding of the present invention and constitute a part of the application. The schematic embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute improper limitations to the present invention. In the attached picture:

Fig. 1 is a flowchart of a remote sensing image mosaic method according to an embodiment of the present invention;

Fig. 2 is a structural block diagram of a remote sensing image stitching device according to an embodiment of the present invention; and

Fig. 3 is a schematic flowchart of an optional implementation manner according to an embodiment of the present invention.

detailed description

Hereinafter, the present invention will be described in detail with reference to the drawings and examples. It should be noted that, in the case of no conflict, the embodiments in the present application and the features in the embodiments can be combined with each other.

An embodiment of the present invention provides a remote sensing image stitching method.

FIG. 1 is a flowchart of a remote sensing image stitching method according to an embodiment of the present invention. As shown in FIG. 1 , the method includes steps 101 to 104 .

Step 101, establishing the first forward calculation model and the first inverse calculation model between the image point and the ground point on each piece of CCD according to the imaging parameters of each piece of CCD, and generating the imaging parameters of virtual CCD according to the imaging parameters of each piece of CCD, Establishing a second forward calculation model and a second inverse calculation model between the image point on the virtual CCD and the ground point according to the imaging parameters of the virtual CCD;

Step 102, according to the above-mentioned first forward calculation model, the first reverse calculation model, the second forward calculation model and the second reverse calculation model, establish the third forward calculation model and The third inverse model;

Step 103, read the image data of each CCD, determine the CCD corresponding to each image point on the virtual CCD and the image point coordinates on the CCD according to the above-mentioned third forward model and the third inverse model;

Step 104, performing resampling processing on each determined image point to obtain a spliced image.

In an optional implementation of the embodiment of the present invention, the same method can be used to determine the forward calculation model and the reverse calculation model for each piece of CCD and virtual CCD. In the above step 101, respective forward calculation models can be established according to their respective imaging parameters. Models and inverse models, including:

Establish a positive calculation model between each image point on the CCD and the ground point according to the imaging parameters of each image point;

Calculate the control point pair parameters according to the forward calculation model of each image point, establish the affine transformation model between each image point on the CCD and the ground point according to the control point pair parameters, and establish the relationship between each image point on the CCD and the ground point through multiple iterative calculations. Inverse model between points.

Optionally, for each image point, a positive calculation model between each image point on the CCD and the ground point is established according to the imaging parameters of each image point, including:

Obtain the position [X _S , Y _S , Z _S ] of the satellite in the agreed geocentric coordinate system at the time when the image point is shot;

Obtain the angle psiX between the main optical axis unit vector of the camera image point corresponding to the image and the X-axis of the satellite body coordinate system, and the angle psiY between the Y-axis;

Obtain the installation matrix M ₀ of the satellite body coordinate system relative to the camera;

Obtain the satellite-to-orbit coordinate system rotation matrix M ₁ at the shooting moment of the image point;

Obtain the rotation matrix M ₂ from the lower orbit to the J2000.0 coordinate system at the moment of image point shooting;

Obtain the J2000.0 to WGS84 coordinate system rotation matrix M ₃ at the moment of image point shooting;

The positive calculation model to determine the relationship between the image point and the ground point is:

$\left[\begin{array}{c}{x}_G\\ Y_G\\ Z_G\end{array}\right]=\left[\begin{array}{c}{x}_{\mathrm{the\; s}}\\ Y_{\mathrm{the\; s}}\\ Z_{\mathrm{the\; s}}\end{array}\right]+{\mathrm{uM}}_3{m}_2{m}_1{m}_0\left[\begin{array}{c}\mathrm{the\; tan}\mathrm{psiX}\\ -\mathrm{the\; tan}\mathrm{psiY}\\ 1\end{array}\right],$ Among them: [X _G , Y _G , Z _G ] is the coordinate of the ground target point corresponding to the image point in the agreed geocentric coordinate system; u is the scaling factor.

It should be noted that there is no sequence between the above steps of obtaining each parameter, and it can be adjusted according to the actual situation.

Optionally, calculate the control point pair parameters according to the forward calculation model of each image point, establish the affine transformation model between each image point on the CCD and the ground point according to the control point pair parameters, and establish each image point on the CCD through multiple iterative calculations. Inverse calculation model between image points and ground points, including:

Select four corner points within a predetermined radius around the coordinates of the center point of the image, and use the strict imaging model to calculate the latitude and longitude of the ground points corresponding to the four corner points;

Establish the affine transformation model Model1 between the image point and the ground point based on the four corner points and their corresponding ground point latitude and longitude;

Calculate the image point (i1, j1) corresponding to the ground point (lat, lon) according to the affine transformation model Model1;

Select four corner points within the predetermined radius around the image point (i1, j1), respectively calculate the latitude and longitude of the ground point corresponding to the four corner points, and establish the relationship between the image point and the ground by the four corner points and the corresponding ground point latitude and longitude Affine transformation model Model2 between points;

Calculate the image point (i2, j2) corresponding to the ground point (lat, lon) according to the affine transformation model Model2;

Determine the distance L between the image point (i2, j2) and the image point (i1, j1), and determine whether the distance L is greater than a preset threshold;

If the distance L is greater than the preset threshold, continue to select the affine transformation model near the point (i2, j2) to continue iterative calculation until the distance L is less than or equal to the preset threshold; if the distance L is less than the preset threshold, Then determine the ground point (lat, lon) as the ground point corresponding to the image point (i1, j1).

In an optional implementation manner of the embodiment of the present invention, for each image point, the above-mentioned step 102 establishes The third forward calculation model and the third inverse calculation model between the image points on each CCD and the virtual CCD, including:

1. Normal model

Calculate the ground point (lat, lon) at the image point (i, j) on the CCD according to the first forward calculation model of the CCD, and calculate the image corresponding to the virtual CCD at the ground point (lat, lon) according to the second inverse calculation model of the virtual CCD Point (i1, j1), get the third positive model of image point (i, j); and

2. Inverse calculation model

A, calculate the ground point (lat, lon) corresponding to the image point (i1, j1) on the virtual CCD according to the second positive calculation model of the virtual CCD;

B, calculate the CCD image point (i2, j2) corresponding to the ground point (lat, lon) according to the first inverse calculation model of CCDM, wherein M is a preset piece number identification;

C, judge whether the ground point (lat, lon) is on the CCDM according to the value of i2, where,

If i2 is less than overlap/2, it is considered that the ground point (lat, lon) is on the left of CCDM, set the slice M to M-1, set the slice number identifier M to M-1, and return A, where overlap is two adjacent slices Overlapping image points between CCDs;

If i2 is greater than the width of CCDM minus the difference of overlap/2, then the ground point (lat, lon) is on the right side of the CCD, set the piece number identification M to M+1, and return A;

If i2 is between (overlap/2, the difference between the width of CCDM and overlap/2), then the ground point (lat, lon) is in CCDM, and the image point of CCDM corresponding to the image point (i1, j1) on the virtual CCD (i2, j2), to obtain the third inverse model of the image point (i1, j1).

In an optional implementation of the embodiment of the present invention, in the above step 101, generating the imaging parameters of the virtual CCD according to the imaging parameters of each piece of CCD includes: smoothing and fitting the imaging parameters of each piece of CCD to obtain the virtual CCD imaging parameters. Optionally, linear fitting is performed on the internal orientation parameters of each CCD by the least square method to establish internal orientation elements of the virtual CCD.

The embodiment of the present invention also provides a remote sensing image splicing device.

Fig. 2 is a structural block diagram of a remote sensing image stitching device according to the implementation strength of the present invention. As shown in Fig. 2, the device includes:

The first building module 201 is used to establish the first forward calculation model and the first inverse calculation model between the image point and the ground point on each piece of CCD according to the imaging parameters of each piece of CCD, and generate a virtual model according to the imaging parameters of each piece of CCD. The imaging parameters of the CCD, according to the imaging parameters of the virtual CCD, the second positive calculation model and the second inverse calculation model between the image point on the virtual CCD and the ground point are established;

The second building module 202 is connected with the first building module 201, and is used to set up image points on each sheet of CCD according to the above-mentioned first positive calculation model, the first inverse calculation model, the second positive calculation model and the second inverse calculation model The third forward calculation model and the third inverse calculation model between the image point on the virtual CCD;

Determining module 203, connected with the second building module 202, is used to read the image data of each slice CCD, determines the CCD corresponding to each image point on the described virtual CCD according to the 3rd forward model and the 3rd inverse model and Image point coordinates on the CCD;

The processing module 204 is connected with the determining module 203, and is configured to perform resampling processing on each determined image point to obtain a spliced image.

Optionally, the first building module 201 is used to set up a forward calculation model and an inverse calculation model according to the imaging parameters of each slice of CCD or the imaging parameters of the virtual CCD; the first building module 201 includes:

A positive calculation model determination unit is used to establish a positive calculation model between each image point on the CCD and the ground point according to the imaging parameters of each image point;

The inverse calculation model determination unit is used to calculate the control point pair parameters according to the forward calculation model of each image point, establish the affine transformation model between each image point on the CCD and the ground point according to the control point pair parameters, and calculate through multiple iterations Establish the inverse calculation model between each image point on the CCD and the ground point.

Optionally, the positive calculation model determines the unit, specifically for

Obtain the position [X _S , Y _S , Z _S ] of the satellite in the agreed geocentric coordinate system at the time when the image point is shot;

Obtain the angle psiX between the main optical axis unit vector of the camera image point corresponding to the image and the X-axis of the satellite body coordinate system, and the angle psiY between the Y-axis;

Obtain the installation matrix M ₀ of the satellite body coordinate system relative to the camera;

Obtain the satellite-to-orbit coordinate system rotation matrix M ₁ at the shooting moment of the image point;

Obtain the rotation matrix M ₂ from the lower orbit to the J2000.0 coordinate system at the moment of image point shooting;

Obtain the J2000.0 to WGS84 coordinate system rotation matrix M ₃ at the moment of image point shooting;

The positive calculation model to determine the relationship between the image point and the ground point is:

$\left[\begin{array}{c}{x}_G\\ Y_G\\ Z_G\end{array}\right]=\left[\begin{array}{c}{x}_{\mathrm{the\; s}}\\ Y_{\mathrm{the\; s}}\\ Z_{\mathrm{the\; s}}\end{array}\right]+{\mathrm{uM}}_3{m}_2{m}_1{m}_0\left[\begin{array}{c}\mathrm{the\; tan}\mathrm{psiX}\\ -\mathrm{the\; tan}\mathrm{psiY}\\ 1\end{array}\right],$ Among them: [X _G , Y _G , Z _G ] is the coordinate of the ground target point corresponding to the image point in the agreed geocentric coordinate system; u is the scaling factor.

Optionally, the above inverse calculation model determination unit is specifically used for:

Select four corner points within a predetermined radius around the coordinates of the center point of the image, and use the strict imaging model to calculate the latitude and longitude of the ground points corresponding to the four corner points;

Establish the affine transformation model Model1 between the image point and the ground point based on the four corner points and their corresponding ground point latitude and longitude;

Calculate the image point (i1, j1) corresponding to the ground point (lat, lon) according to the affine transformation model Model1;

Select four corner points within the predetermined radius around the image point (i1, j1), respectively calculate the latitude and longitude of the ground point corresponding to the four corner points, and establish the relationship between the image point and the ground by the four corner points and the corresponding ground point latitude and longitude Affine transformation model Model2 between points;

Calculate the image point (i2, j2) corresponding to the ground point (lat, lon) according to the affine transformation model Model2;

Determine the distance L between the image point (i2, j2) and the image point (i1, j1), and determine whether the distance L is greater than a preset threshold;

If the distance L is greater than the preset threshold, continue to select the affine transformation model near the point (i2, j2) to continue iterative calculation until the distance L is less than or equal to the preset threshold; if the distance L is less than the preset threshold, Then determine the ground point (lat, lon) as the ground point corresponding to the image point (i1, j1).

Optionally, the above-mentioned second building module 202 is specifically used for

Calculate the ground point (lat, lon) at the image point (i, j) on the CCD according to the first forward calculation model of the CCD, and calculate the image corresponding to the virtual CCD at the ground point (lat, lon) according to the second inverse calculation model of the virtual CCD Point (i1, j1), get the third positive model of image point (i, j); and

A, calculate the ground point (lat, lon) corresponding to the image point (i1, j1) on the virtual CCD according to the second positive calculation model of the virtual CCD;

B, calculate the CCD image point (i2, j2) corresponding to the ground point (lat, lon) according to the first inverse calculation model of CCDM, wherein M is a preset piece number identification;

C, judge whether the ground point (lat, lon) is on the CCDM according to the value of i2, where,

If i2 is less than overlap/2, it is considered that the ground point (lat, lon) is on the left of CCDM, set the slice M to M-1, set the slice number identifier M to M-1, and return A, where overlap is two adjacent slices Overlapping image points between CCDs;

If i2 is greater than the width of CCDM minus the difference of overlap/2, then the ground point (lat, lon) is on the right side of the CCD, set the piece number identification M to M+1, and return A;

If i2 is between (overlap/2, the difference between the width of CCDM and overlap/2), then the ground point (lat, lon) is in CCDM, and the image point of CCDM corresponding to the image point (i1, j1) on the virtual CCD (i2, j2), to obtain the third inverse model of the image point (i1, j1).

The parts that are the same as those of the above method will not be repeated here.

An optional implementation manner of the embodiment of the present invention will be described in detail below with reference to FIG. 3 .

Fig. 3 is a schematic flow chart of an optional implementation manner according to an embodiment of the present invention. In combination with what is shown in Fig. 3, this optional implementation manner may include the following processes:

1. Obtain auxiliary data such as orientation elements, orbits, attitudes, and travel times at each image point, and establish a positive calculation model between each image point and ground point, and a strict imaging model at each image point. Specifically include the following steps:

(1) Analyze the orbit, attitude, travel time and other data within the imaging time range of the image to be checked from the original data downloaded by the satellite;

(2) For any image point Pointsi{sample, line, lat, lon, height}, according to the image coordinates (sample, line) of the image point, obtain its corresponding photography time scanTime;

The imaging time corresponding to the image point can be directly analyzed from the auxiliary data downloaded from the satellite, and the imaging time corresponding to the auxiliary data in the line line is the imaging time corresponding to the image point.

(3) Use the Lagrange interpolation algorithm to interpolate and calculate the satellite orbit position (PX, PY, PZ, VX, VY, VZ) of the photography time scanTime; satellites download orbit data according to a certain frequency, therefore, the satellite corresponding to the photography time scanTime The orbital position needs to be interpolated using several sets of orbital data before and after the shooting time. The invention adopts the Lagrangian interpolation algorithm, and uses three sets of orbit data before and after the photographing time to calculate the satellite orbit position at the photographing time.

a) Starting from the first group of orbital data, judge the relationship between the generation time of this group of orbital data, the generation time of the next group of orbital data and scanTime, if scanTime is greater than the generation time of the i-th group of orbital data, and less than the i-th +1 generation time of track data, then record i is the sequence number of the track data closest to the shooting time.

b) Using the i-1, i, and i+1 groups of orbit data, calculate the orbit position and velocity at the time of photography based on the Lagrangian algorithm. The Lagrangian interpolation algorithm is expressed as follows:

For the function table (xi,f(xi))(i=0,1,…,n) of known y=f(x), for any x in the range of [xo,xn], there are:

${\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; *</span>\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; \&NotEqual;</span>\mathrm{<span\; class="notranslate">\; i</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{\mathrm{<span\; class="notranslate">\; \&Pi;</span>}}}\frac{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}}{{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}}<span\; class="notranslate">\; )</span>$

(4) Use the Lagrange algorithm to interpolate to calculate the three-axis attitude angle (Roll, Pitch, Yaw) of the scanTime camera relative to the orbital coordinate system;

(5) According to the camera laboratory measurement data, read the optical axis pointing angle (psiX, psiY) of the first sample CCD detector corresponding to the image coordinates (sample, line) of the image point pair.

(6) Establish a strict imaging model at the image point, and use its internal orientation, orbit, attitude, travel time and other data to construct a strict geometric imaging model Model of remote sensing images for any image point

The strict imaging geometry model of the linear push-broom camera is described below.

$\left[\begin{array}{c}{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; G</span>}}\\ {\mathrm{<span\; class="notranslate">\; Y</span>}}_{\mathrm{<span\; class="notranslate">\; G</span>}}\\ {\mathrm{<span\; class="notranslate">\; Z</span>}}_{\mathrm{<span\; class="notranslate">\; G</span>}}\end{array}\right]<span\; class="notranslate">\; =</span>\left[\begin{array}{c}{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}\\ {\mathrm{<span\; class="notranslate">\; Y</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}\\ {\mathrm{<span\; class="notranslate">\; Z</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}\end{array}\right]<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; uM</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}{\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}\left[\begin{array}{c}\mathrm{<span\; class="notranslate">\; the\; tan</span>}\mathrm{<span\; class="notranslate">\; psiX</span>}\\ <span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; the\; tan</span>}\mathrm{<span\; class="notranslate">\; psiY</span>}\\ \mathrm{<span\; class="notranslate">\; 1</span>}\end{array}\right]$

in:

[X _S , Y _S , Z _S ] is the position of the satellite in the agreed geocentric coordinate system at this moment, that is, the orbital position at the image point (PX, PY, PZ);

[X _G , Y _G , Z _G ] are the coordinates of the ground target point corresponding to the pixel in the agreed geocentric coordinate system.

psiX and psiY are the angles between the main optical axis unit vector of the camera pixel corresponding to the image and the X-axis and Y-axis of the satellite body coordinate system;

u is a scaling factor.

M ₀ is the installation matrix of the satellite body coordinate system relative to the camera, which is obtained by ground measurement before the satellite is launched;

M ₁ is the rotation matrix from the satellite to the orbit coordinate system at this moment, which is composed of the attitude angle measured on the satellite.

M ₂ is the rotation matrix from the lower orbit to the J2000.0 coordinate system at this moment, which is composed of right ascension of ascending node, orbital inclination, argument, etc.

M ₃ is the rotation matrix of the J2000.0 to WGS84 coordinate system at this moment, which needs to be corrected for precession, nutation, Greenwich mean sidereal time and pole shift.

2. Calculate the control point pair parameters from the strict imaging model, establish the affine transformation model between the image point and the ground point, and establish the inverse calculation model through multiple iterative calculations.

(1) Select four corner points within a radius of 5*5 around the coordinates of the image center point, use the strict imaging model to calculate the latitude and longitude of the ground points corresponding to the four corner points, and build an image from the four corner points and their corresponding ground point latitude and longitude Affine transformation model Model1 between points and ground coordinates;

(2) Use the affine transformation model Model1 to calculate the image point coordinates (i1, j1) corresponding to the ground point coordinates (lat, lon);

(3) Select four corner points within a radius of 5*5 around the image point (i1, j1), calculate the latitude and longitude of the ground point corresponding to the four corner points, and build an image from the four corner points and their corresponding ground point latitude and longitude Affine transformation model Model2 between points and ground coordinates;

(4) Use the affine transformation model Model2 to calculate the image point coordinates (i2, j2) corresponding to the ground point coordinates (lat, lon);

(5) Calculate the distance L between (i2, j2) and (i1, j1). If L is greater than the set threshold, continue to select the affine transformation model near the (i2, j2) point to continue iterative calculation; if the value of L is less than the threshold, then end the iterative calculation, and the coordinates of the ground point calculated in step (4) are The corresponding coordinates of the image point.

3. Perform smoothing and fitting calculations based on the input internal orientation parameters of each CCD of the actual camera, as well as the actual orbit, attitude, and travel time data, and establish the imaging parameters of the ideal CCD.

The characteristic of this group of parameters is that each parameter presents smooth linear or low-order polynomial characteristics. The establishment goals of each parameter are as follows:

(1) Virtual CCD internal orientation elements: the real internal orientation parameters of each CCD of the camera are calibrated, and the least square method is used for linear fitting to establish an ideal undistorted internal orientation element;

(2) Time series of virtual imaging lines: the start and end times of imaging are equally divided according to the number of imaging times to ensure that the integration intervals between each line of images are equal;

(3) Virtual imaging track: smoothing the downlinked GPS measurement data (position and velocity) and fitting them into a straight line or low-order polynomial curve;

(4) Virtual imaging attitude: smooth the downloaded or calculated attitude angle data, and fit it into a straight line or a low-order polynomial curve;

4. According to the front and back calculation model of each CCD of the camera and the front and back calculation model of the virtual scanning CCD, the same ground coordinates are used for mapping, and the direct forward and reverse calculation solution of the image points on each CCD of the camera and the image points on the virtual CCD is established Model.

Forward calculation model: Let N represent the chip number of CCD. For a certain pixel (i, j) on a certain CCDN, i represents the row number and j represents the column number. The calculation process of the corresponding virtual CCD image point coordinates (i1, j1) is as follows:

a) Calculate the ground coordinates (lat, lon) at the image point (i, j) from the image point of the CCD and the ground point forward calculation model;

b) Calculate the image point coordinates (i1, j1) of the virtual CCD corresponding to the ground coordinates (lat, lon) from the image points of the virtual CCD and the inverse calculation model of the ground points.

Inverse calculation model: let N represent the chip number mark of CCD, for the virtual CCD image point coordinates (i1, j1), its corresponding CCD mark and the image point coordinates (i, j) on the CCD are calculated as follows:

a) calculate the ground point coordinates (lat, lon) corresponding to the image point coordinates (i1, j1) of the virtual CCD by the image point of the virtual CCD and the ground point forward calculation model;

b) Suppose the CCD chip number ccdID is N/2:

c) Calculate the CCD image point coordinates (i2, j2) corresponding to the ground point coordinates (lat, lon) from the image points of the CCD and the ground point inversion model;

d) judging whether the feature point is on the CCD according to the value of i2;

e) If i2 is less than overlap/2 (overlap is the overlapping pixel between two adjacent CCDs), then it is considered that the feature point is on the left side of the CCD, reset ccdID to ccdID-1, and repeat a);

f) If i2 is greater than (the width of the CCD-overlap/2), it is considered that the feature point is on the right side of the CCD, reset the ccdID to ccdID+1, and repeat a);

g) If i2 is between (overlap/2, (the width of the CCD-overlap/2)), it is considered that the feature point is in the CCD, and the virtual image point coordinates (i1, j1) correspond to the CCD The coordinates of the image point on are (i2, j2).

5. According to the mapping relationship between each CCD of the camera and the virtual CCD image points, grayscale resampling is performed on each image point of the virtual scanning scene image to obtain the spliced image. Specifically, read the original image data obtained by each piece of CCD, and for each pixel of the virtual CCD, obtain the piece number and pixel coordinates of its corresponding original CCD, perform resampling processing (multithreading) pixel by pixel, and output resampling After the virtual stitching image.

From the above description, it can be seen that the present invention achieves the following technical effects: using the inter-film photographic geometric constraint relationship of the CCD, based on the continuity of the object space, the image point coordinate conversion relationship between the virtual scanning scene and the original image is established, Furthermore, the seamless inner field of view splicing processing of the imaging data is realized; at the same time, the ideal CCD inner orientation parameters are obtained by fitting the least square method, which effectively improves the geometric distortion inside the image. It breaks through the limitations of traditional splicing methods, and provides a more rigorous and ideal technical solution for producing high-quality virtual scanning scene products.

Obviously, those skilled in the art should understand that each module or each step of the above-mentioned present invention can be realized by a general-purpose computing device, and they can be concentrated on a single computing device, or distributed in a network formed by multiple computing devices Alternatively, they may be implemented in program code executable by a computing device so that they may be stored in a storage device to be executed by a computing device, and in some cases in an order different from that shown here The steps shown or described are carried out, or they are separately fabricated into individual integrated circuit modules, or multiple modules or steps among them are fabricated into a single integrated circuit module for implementation. As such, the present invention is not limited to any specific combination of hardware and software.

The above descriptions are only preferred embodiments of the present invention, and are not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and changes. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention shall be included within the protection scope of the present invention.

Claims (10)

1. a remote sensing image joining method, is characterized in that, comprising:

Set up first between picture point and ground point on described each CCD according to the imaging parameters of each CCD and just calculate model and the first inverse model, and the imaging parameters of imaging parameters generating virtual CCD according to described each CCD, set up second between picture point and ground point on described virtual CCD according to the imaging parameters of described virtual CCD and just calculate model and the second inverse model;

According to described first just calculating model, model just calculated by described first inverse model, described second and described second inverse model, to set up on described each CCD in picture point and described virtual CCD between picture point the 3rd and just calculating model and the 3rd inverse model;

Read the image data of described each CCD, just calculate CCD that model and described 3rd inverse model determine that on described virtual CCD, each picture point is corresponding and the picpointed coordinate on this CCD according to the described 3rd;

Resampling process is carried out to determined each picture point, obtains spliced image.

2. method according to claim 1, is characterized in that, for described each CCD and described virtual CCD, sets up respective just calculation model and inverse model, comprising according to described respective imaging parameters:

The just calculation model on CCD between each picture point and ground point is set up according to the imaging parameters of each picture point;

Just calculation model according to each picture point described calculates dominating pair of vertices parameter, set up the affine Transform Model on CCD between each picture point and ground point according to dominating pair of vertices parameter, and calculated the inverse model set up on CCD between each picture point and ground point by successive ignition.

3. method according to claim 2, is characterized in that, for each picture point, sets up the just calculation model on CCD between each picture point and ground point, comprising according to the imaging parameters of each picture point:

Position [the X of satellite in agreement geocentric coordinate system is inscribed when obtaining picture point shooting _s, Y _s, Z _s];

Obtain camera picture point primary optical axis unit vector corresponding to image and the angle psiX of satellite body coordinate system X-axis, the angle psiY of Y-axis;

Obtain the installation matrix M of satellite body coordinate system relative to camera ₀;

Satellite is inscribed to orbital coordinate system rotation matrix M when obtaining the shooting of this picture point ₁;

Obtain picture point shooting moment lower railway to J2000.0 coordinate system rotation matrix M ₂;

J2000.0 to WGS84 coordinate system rotation matrix M is inscribed when obtaining picture point shooting ₃;

Determine that the just calculation model between picture point and ground point is:

$\left[\begin{array}{c}{X}_G\\ Y_G\\ Z_G\end{array}\right]=\left[\begin{array}{c}{X}_s\\ Y_s\\ Z_s\end{array}\right]+u{M}_3{M}_2{M}_1{M}_0\left[\begin{array}{c}\tan \mathrm{psiX}\\ -\tan \mathrm{psiY}\\ 1\end{array}\right],$ Wherein: [X _g, Y _g, Z _g] be the coordinate of terrain object point in agreement geocentric coordinate system corresponding to picture point; U is scale factor.

4. according to the method in claim 2 or 3, it is characterized in that, just calculation model according to each picture point described calculates dominating pair of vertices parameter, the affine Transform Model on CCD between each picture point and ground point is set up according to dominating pair of vertices parameter, and the inverse model set up on CCD between each picture point and ground point is calculated by successive ignition, comprising:

Choose four angle points in predetermined radii around heart point coordinate in the picture, utilize rigorous geometry model to calculate ground point longitude and latitude corresponding to four angle points respectively;

The affine Transform Model Model1 between picture point and ground point is set up according to the ground point longitude and latitude by four angle points and correspondence thereof;

Picture point (i1, j1) corresponding to ground point (lat, lon) is calculated according to affine Transform Model Model1;

At picture point (i1, j1) four angle points are chosen in described predetermined radii around, calculate ground point longitude and latitude corresponding to four angle points respectively, set up the affine Transform Model Model2 between picture point and ground point by the ground point longitude and latitude of these four angle points and correspondence thereof;

Picture point (i2, j2) corresponding to ground point (lat, lon) is calculated according to affine Transform Model Model2;

Determine the distance L between picture point (i2, j2) and picture point (i1, j1), and whether judging distance L is greater than predetermined threshold value;

If distance L is greater than described predetermined threshold value, continue to select affine Transform Model to continue iterative computation, until distance L is less than or equal to described predetermined threshold value near (i2, j2) point; If distance L is less than described predetermined threshold value, then millet cake (lat, lon) is the ground point of picture point (i1, j1) correspondence definitely.

5. method according to claim 1, for each picture point, according to described first just calculating model, model just calculated by described first inverse model, described second and described second inverse model, to set up on described each CCD in picture point and described virtual CCD between picture point the 3rd and just calculating model and the 3rd inverse model, comprising:

Just calculate model according to first of CCD and calculate picture point (i on CCD, j) ground point (lat at place, lon), the second inverse model according to virtual CCD calculates ground point (lat, lon) picture point (i1 of corresponding virtual CCD, j1), obtain picture point (i, j) the 3rd and just calculate model; And

A, is just calculating model according to second of virtual CCD and is calculating the ground point (lat, lon) that on virtual CCD, picture point (i1, j1) is corresponding;

B, calculates CCD picture point (i2, j2) corresponding to ground point (lat, lon) according to the first inverse model of CCDM, and wherein M is default sheet number mark;

C, judges ground point (lat, lon) whether on CCDM according to the value of i2, wherein,

If i2 is less than overlap/2, then think that ground point (lat, lon) is on the CCDM left side, is set to M-1 by sheet M, a sheet number mark M is set to M-1, and return A, wherein overlap is the overlapping picture point between adjacent two panels CCD;

If the width that i2 is greater than CCDM deducts the difference of overlap/2, then ground point (lat, lon) is on the right of this sheet CCD, a sheet number mark M is set to M+1, returns A;

If i2 is at (overlap/2, the width of CCDM and the difference of overlap/2) between, then ground point (lat, lon) in CCDM, the picture point (i2, j2) of the CCDM that picture point (i1, j1) is corresponding on virtual CCD, obtain the 3rd inverse model of picture point (i1, j1).

6. method according to claim 1, is characterized in that, according to the imaging parameters of the imaging parameters generating virtual CCD of described each CCD, comprising:

To the imaging parameters of described each CCD, smoothing and process of fitting treatment, obtains the imaging parameters of virtual CCD.

7. method according to claim 6, is characterized in that, the imaging parameters of described each CCD comprises the elements of interior orientation of described each CCD, and the imaging parameters of described virtual CCD comprises the elements of interior orientation of described virtual CCD; Smoothing and the process of fitting treatment to the imaging parameters of described each CCD, obtains the imaging parameters of virtual CCD, comprising:

To the interior orientation parameter of described each CCD, carry out linear fit by least square method, set up the elements of interior orientation of described virtual CCD.

8. a remote sensing image splicing apparatus, is characterized in that, comprising:

First sets up module, model and the first inverse model is just being calculated for setting up first between picture point and ground point on described each CCD according to the imaging parameters of each CCD, and the imaging parameters of imaging parameters generating virtual CCD according to described each CCD, set up second between picture point and ground point on described virtual CCD according to the imaging parameters of described virtual CCD and just calculate model and the second inverse model;

Second sets up module, for according to described first just calculating model, model just calculated by described first inverse model, described second and described second inverse model, to set up on described each CCD in picture point and described virtual CCD between picture point the 3rd and just calculating model and the 3rd inverse model;

Determination module, for reading the image data of described each CCD, is just calculating CCD that model and described 3rd inverse model determine that on described virtual CCD, each picture point is corresponding and the picpointed coordinate on this CCD according to the described 3rd;

Processing module, for carrying out resampling process to determined each picture point, obtains spliced image.

9. device according to claim 8, is characterized in that, described first sets up module, is just calculating model and inverse model for setting up according to the imaging parameters of each CCD or the imaging parameters of virtual CCD; Described first sets up module, comprising:

Just calculating model determining unit, for setting up the just calculation model on CCD between each picture point and ground point according to the imaging parameters of each picture point;

Inverse model determining unit, dominating pair of vertices parameter is calculated for the just calculation model according to each picture point described, set up the affine Transform Model on CCD between each picture point and ground point according to dominating pair of vertices parameter, and calculated the inverse model set up on CCD between each picture point and ground point by successive ignition.

Or

Described second sets up module, specifically for

Just calculate model according to first of CCD and calculate picture point (i on CCD, j) ground point (lat at place, lon), the second inverse model according to virtual CCD calculates ground point (lat, lon) picture point (i1 of corresponding virtual CCD, j1), obtain picture point (i, j) the 3rd and just calculate model; And

A, is just calculating model according to second of virtual CCD and is calculating the ground point (lat, lon) that on virtual CCD, picture point (i1, j1) is corresponding;

B, calculates CCD picture point (i2, j2) corresponding to ground point (lat, lon) according to the first inverse model of CCDM, and wherein M is default sheet number mark;

C, judges ground point (lat, lon) whether on CCDM according to the value of i2, wherein,

If i2 is less than overlap/2, then think that ground point (lat, lon) is on the CCDM left side, is set to M-1 by sheet M, a sheet number mark M is set to M-1, and return A, wherein overlap is the overlapping picture point between adjacent two panels CCD;

If the width that i2 is greater than CCDM deducts the difference of overlap/2, then ground point (lat, lon) is on the right of this sheet CCD, a sheet number mark M is set to M+1, returns A;

If i2 is at (overlap/2, the width of CCDM and the difference of overlap/2) between, then ground point (lat, lon) in CCDM, the picture point (i2, j2) of the CCDM that picture point (i1, j1) is corresponding on virtual CCD, obtain the 3rd inverse model of picture point (i1, j1).

10. device according to claim 9, is characterized in that, is describedly just calculating model determining unit, for

Position [the X of satellite in agreement geocentric coordinate system is inscribed when obtaining picture point shooting _s, Y _s, Z _s];

Obtain camera picture point primary optical axis unit vector corresponding to image and the angle psiX of satellite body coordinate system X-axis, the angle psiY of Y-axis;

Obtain the installation matrix M of satellite body coordinate system relative to camera ₀;

Satellite is inscribed to orbital coordinate system rotation matrix M when obtaining the shooting of this picture point ₁;

Obtain picture point shooting moment lower railway to J2000.0 coordinate system rotation matrix M ₂;

J2000.0 to WGS84 coordinate system rotation matrix M is inscribed when obtaining picture point shooting ₃;

Determine that the just calculation model between picture point and ground point is:

$\left[\begin{array}{c}{X}_G\\ Y_G\\ Z_G\end{array}\right]=\left[\begin{array}{c}{X}_s\\ Y_s\\ Z_s\end{array}\right]+u{M}_3{M}_2{M}_1{M}_0\left[\begin{array}{c}\tan \mathrm{psiX}\\ -\tan \mathrm{psiY}\\ 1\end{array}\right],$ Wherein: [X _g, Y _g, Z _g] be the coordinate of terrain object point in agreement geocentric coordinate system corresponding to picture point; U is scale factor.

Or

Described inverse model determining unit, for:

Choose four angle points in predetermined radii around heart point coordinate in the picture, utilize rigorous geometry model to calculate ground point longitude and latitude corresponding to four angle points respectively;

The affine Transform Model Model1 between picture point and ground point is set up according to the ground point longitude and latitude by four angle points and correspondence thereof;

Picture point (i1, j1) corresponding to ground point (lat, lon) is calculated according to affine Transform Model Model1;

At picture point (i1, j1) four angle points are chosen in described predetermined radii around, calculate ground point longitude and latitude corresponding to four angle points respectively, set up the affine Transform Model Model2 between picture point and ground point by the ground point longitude and latitude of these four angle points and correspondence thereof;

Picture point (i2, j2) corresponding to ground point (lat, lon) is calculated according to affine Transform Model Model2;

Determine the distance L between picture point (i2, j2) and picture point (i1, j1), and whether judging distance L is greater than predetermined threshold value;

If distance L is greater than described predetermined threshold value, continue to select affine Transform Model to continue iterative computation, until distance L is less than or equal to described predetermined threshold value near (i2, j2) point; If distance L is less than described predetermined threshold value, then millet cake (lat, lon) is the ground point of picture point (i1, j1) correspondence definitely.