阅读说明：本技术 针对国产面阵摆扫成像系统获取的影像生产正射影像的方法 (Method for producing orthoimage aiming at image acquired by domestic area array swinging imaging system ) 是由 段延松 赵新博 张祖勋 柯涛 于 2021-08-24 设计创作，主要内容包括：本发明涉及一种针对国产面阵摆扫成像系统获取的影像生产正射影像的方法。首先选择合理的局部坐标系作为基准坐标系,将导航数据信息转换到基准坐标系中,然后在DEM辅助下对所有影像进行同名点自动提取,形成影像连接点网,之后使用连接点开展基于DEM的区域网平差,解算每张影像的外方位参数,最后基于最优化理论生产目标区域的整体最优正射影像。本发明引入GPS/IMU参数作为初值,使用DEM数据辅助平差,解决了传统摄影测量无法直接进行平差解算的问题。本发明可以更快速、更准确地完成影像拼接,生产合格的正射影像产品,填补国产面阵摆扫成像系统数据处理的空白。(The invention relates to a method for producing an orthoimage aiming at an image acquired by a domestic area array sweep imaging system. Firstly, a reasonable local coordinate system is selected as a reference coordinate system, navigation data information is converted into the reference coordinate system, then homonymous points are automatically extracted from all images under the assistance of a DEM (digital elevation model) to form an image connection point network, then, the connection points are used for carrying out adjustment of the regional network based on the DEM, the external orientation parameter of each image is solved, and finally, the overall optimal orthoscopic image of a target region is produced based on an optimization theory. The invention introduces GPS/IMU parameters as initial values, uses DEM data to assist adjustment, and solves the problem that the adjustment can not be directly calculated in the traditional photogrammetry. The method can complete image splicing more quickly and accurately, produce qualified orthoscopic image products, and fill the blank of data processing of a domestic area array scanning imaging system.)

1. A method for producing an orthoscopic image aiming at an image acquired by a domestic area array scanning imaging system is characterized by comprising the following steps:

step 1, establishing a local coordinate system as a unified reference coordinate system, and converting navigation data information into the reference coordinate system;

step 2, automatically extracting the same-name points of the image under the assistance of the rough DEM;

step 2.1, under the assistance of the general DEM, acquiring the adjacent relation between the images according to the scanning imaging periodicity and the sequence;

2.2, according to the adjacent relation of the images obtained in the step 2.1, carrying out homonymy point matching by adopting a method of carrying out dense matching on the stereo images based on two types of expansion;

step 3, adjusting the local area network based on DEM assistance, taking the navigation data information converted in the step 1 as an initial adjustment value, and calculating the external orientation element parameter of each image;

step 4, producing an overall optimal orthoimage;

step 4.1, calculating the ground coordinates corresponding to the four corner points of each original image;

step 4.2, selecting the optimal original image corresponding to the orthoimage to perform back projection;

and 4.3, producing the orthoimage by a back projection method.

2. The orthoscopic image production method for the domestic area array scanning imaging system as claimed in claim 1The image method is characterized in that: in the step 1, a local coordinate system adopts a fixed point tangent plane coordinate system, namely a TPC coordinate system, and the coordinates are defined as follows: selecting a geometric center point longitude L in a flight area under a WGS84 coordinate frame₀Latitude B₀As a plane origin, taking the height of an earth ellipsoid as 0 as an ellipsoid tangent plane, defining the east direction passing through an origin latitude line as an X coordinate axis, the north direction passing through an origin longitude line as a Y coordinate axis, and a vertical line passing through the origin and back to the geocenter as a Z coordinate axis; the TPC coordinate system and the image space coordinate system have a conversion relation as follows:

namely, the image space coordinate system → the load body coordinate system → the IMU body coordinate system → the navigation coordinate system → the geocentric geostationary coordinate system → the TPC coordinate system;

the calculation formula of the rotation matrix between the coordinate systems is as follows:

in the formula，Is an image attitude rotation matrix under the TPC coordinate system,three external orientation angle elements of the image are taken as parameters to be solved;is a transformation matrix between geocentric geostationary coordinate system and TPC coordinate system, from longitude L of TPC coordinate origin₀Latitude B₀Defining;is a conversion matrix between a navigation coordinate system and a geocentric geostationary coordinate system, which is determined by the longitude L of the position of the load_zLatitude B_zDefining;is a rotation matrix between a navigation coordinate system and an IMU body coordinate system, and (phi, theta and psi) are pitching, rolling and yawing numerical values recorded by the IMU;is a transformation matrix between the load body coordinate system and the IMU body coordinate system, also called a device placement matrix, placed on the load by the IMU_x、e_y、e_zDefining three angles;is a transformation matrix between the spatial coordinate system and the load body coordinate system;

solving a rotation matrixThen, the inverse solution can obtain the attitude angle of the image in the TPC coordinate system, that is:

ω＝arcsin(-b₃) (9)

the GPS point position information can be losslessly transformed into any desired coordinate system using the earth coordinate definition model, and the TPC coordinates transformed from the original GPS coordinates (L, B, H) into the local tangent plane coordinate system are transformed as follows:

the original GPS data is first converted into the geocentric coordinate system:

X_c＝(N+H)cos B cos L (11)

Y_c＝(N+H)cos B sin L (12)

Z_c＝[N(1-e²)+H]sin B (13)

in the formula, N is the radius of the unitary mortise ring, a is the long semi-axis of the earth ellipsoid, b is the short semi-axis of the earth ellipsoid, (X)_c,Y_c,Z_c) Coordinates of a geocentric coordinate system;

then the origin coordinate (L) of the local tangent plane coordinate system₀,B₀And 0) carrying out the following steps (11) to (13) to find the origin of the local tangent plane coordinate system at the geocentric positionCoordinates in the system, noteThe calculation formula from the geocentric system coordinates to the TPC coordinate system is:

in the formula (I), the compound is shown in the specification,is a transformation matrix between the geocentric coordinate system and the TPC coordinate system, and can be obtained by calculation of the formula (2); (X)_TPC,Y_TPC,Z_TPC) Coordinates of a TPC coordinate system;

i.e. the navigation information converted into the reference coordinate system.

3. The method for producing an ortho image from an image acquired by a domestic area array sweep imaging system as claimed in claim 1, wherein: in the step 2, the adjacent relation among the images is obtained according to the scanning imaging periodicity and the sequence, namely N images obtained in one scanning period are divided into a group, each image is numbered from 1 to N according to the sequence of the image in the scanning period, and the images which are continuous in sequence are necessarily adjacent images; and the adjacent relation among the groups projects the acquired images onto the ground represented by the global DEM according to the GPS and IMU information, and then the image overlapping condition is obtained according to the ground range of each image, or the overlapping area of the original image is inversely calculated according to the overlapping area, so that the adjacent relation of the images among the groups is confirmed.

4. The method for producing an ortho image from an image acquired by a domestic area array sweep imaging system as claimed in claim 2, wherein: the step 3 of performing block adjustment based on DEM assistance comprises the following steps:

step 3.1, calculating the area network range by using the POS information of all the images, and acquiring corresponding DEM data from the global SRTM data according to the area network range;

step 3.2, calculating a rotation matrix of each image by using the external orientation angle element, and using the angle element obtained in the step 1 as an initial value of calculation;

step 3.3, respectively calculating the intersection points of the light rays connected with the image points and the DEM by taking the single-beam same-name light rays as an observation unit to obtain corresponding object space coordinates;

step 3.4, calculating the difference between object space points corresponding to the image points with the same name, and if the difference is greater than the theoretical object space positioning precision alpha derived from the angle element precision, marking the image points with the same name as rough difference points without participating in adjustment; otherwise, taking the average value of the object space coordinates corresponding to the image points with the same name as a true value to participate in adjustment;

step 3.5, constructing an error equation point by point, and resolving an external orientation element;

step 3.6, repeating steps 3.2 to 3.5 until the change of each corner element of the image does not exceed 10^-6Radian and the sum of the squares of all the image point corrections is less than a set threshold lambda.

5. The method for producing an ortho image from an image acquired by a domestic area array sweep imaging system as set forth in claim 4, wherein: in the step 3.3, the single beam of light is used as an observation unit, the intersection points of each connecting point and the DEM are respectively calculated, and a corresponding object space coordinate calculation formula is obtained as follows:

in the formula (X)_g,Y_g) The object space coordinate of the intersection point is the TCP coordinate; z is a radical of_demIs a crossDEM elevation of points; x_s、Y_s、Z_sIs the exterior orientation line element of the image, obtained in step 1 (X)_TPC,Y_TPC,Z_TPC) Can be used as an initial value of the exterior orientation line element; (x, y) are image point coordinates; f is the camera focal length.

6. The method for producing an ortho image from an image acquired by a domestic area array sweep imaging system as set forth in claim 5, wherein: the calculation method of the positioning accuracy α of the theoretical object space in the step 3.4 is as follows:

in the formula, H is the height from the photographing position to the ground, epsilon is the IMU angle element precision of the imaging system, and theta is the included angle between the main optical axis and the gravity direction.

7. The method for producing an ortho image from an image acquired by a domestic area array sweep imaging system as set forth in claim 6, wherein: the point-by-point construction error equation in the step 3.5 is as follows:

wherein (x, y) is the coordinates of the image point, (x), (y) are the coordinates of the image point obtained by using the initial value of the unknown number, v_x、v_yIs the variable of the error, and is,and d represents error correction, and the external orientation element of each image can be updated after the correction number of the external orientation element is solved.

8. The method for producing an ortho image from an image acquired by a domestic area array sweep imaging system as claimed in claim 1, wherein: in said step 4.3The method for producing the orthoimage by the back projection method comprises the following steps: let the coordinate of any point P on the ortho-image be (X ', Y '), from the ground coordinate (X ') of the contour point at the lower left corner of the ortho-image₀,Y₀) Calculating the ground coordinates (X, Y) corresponding to the point P according to the orthographic image scale denominator M as follows:

X＝X₀+M×X′ (22)

Y＝Y₀+M×Y′ (23)

calculating the coordinate P (x, y) of the point P on the ortho image corresponding to the corresponding point on the original image as:

in the formula, X₀、Y₀The minimum value in the ground coordinates corresponding to the angular points of all the original images is taken as Z, the elevation of the P point is taken as Z, and the Z is obtained by interpolation of the DEM; a is_i,b_i,c_i(i ═ 1,2,3) is the rotation matrix element of the image, X_s、Y_s、Z_sIs the line element of the image, f is the focal length of the original image;

because the obtained coordinates of the image points do not necessarily exactly fall in the centers of the pixel elements, gray scale interpolation is needed for this purpose, a method of bilinear interpolation is adopted to obtain the gray scale value g (X, Y) of the image point P, and the gray scale value of the image point is assigned to the corresponding pixel point P (X ', Y') on the orthoimage, that is:

G(X′,Y′)＝g(x,y) (26)

and sequentially carrying out the operation on each pixel to obtain the orthoimage.

Technical Field

The invention belongs to the technical field of remote sensing mapping, and particularly relates to a method for producing an orthoimage aiming at an image acquired by a domestic area array sweep imaging system.

Background

The area array swinging imaging system is equipment for quickly acquiring ground image information, and compared with a traditional fixed frame imaging system and a multi-view inclined imaging system, the area array swinging imaging system is high in data acquisition efficiency, low in overall cost and high in practicability. Commercial systems have been introduced internationally as early as 2008, and are typically represented by the israel a3 area array swipe imaging system and the japanese SamplX area array swipe imaging system. These systems have been commercialized to provide a complete solution to the outside world, which involves a full flow production process operating specification from the aerial installation of the system to the final mapping product (including aerial triangulation, digital elevation model, orthophoto, digital line drawing), but does not provide relevant information on specific technical details, such as aerial triangulation processes and technical details. Manufacturers only provide device output result parameters such as camera resolution, camera focal length, image format and the like in device descriptions, do not provide more core parameters such as the working principle, parameters and the like of an imaging system POS system, and only provide software for developing and producing image information acquired by the system. Through the efforts of 10 years in scientific research institutions in related technical fields in China, a domestic area array scanning imaging system is successfully developed. The domestic area array swinging imaging system is composed of imaging hardware and processing software, wherein the imaging hardware is mainly composed of a lens, an imaging CCD panel, a GPS position measuring device, an IMU (inertial measurement unit) position measuring device, a swinging driving device, a storage device and other core components, and is shown in figure 1. The processing software mainly comprises an image signal decoding module, a GPS data processing module, an IMU data processing module and a mapping product production module.

In the imaging principle, the area array sweep imaging system has great difference with the traditional imaging system, and the maximum characteristic is continuous sweep imaging, in the sweep period, any two images can be approximately understood as concentric imaging, seamless splicing can be realized, the overlapped area is the same light imaging, and the three-dimensional positioning can not be carried out. In addition, the IMU used by the domestic area array swing scanning imaging system is limited by an autonomous technology, the precision can only be 5-10 arc seconds, and the accurate positioning and attitude determination of the image cannot be realized. Therefore, the processing software of the area array sweep imaging system needs to fully consider the above characteristics. In addition, each image of the area array scanning imaging system has GPS information and IMU information, and important information can be provided for processing software.

The characteristics of the swinging and scanning images, particularly the characteristic that the concentric imaging cannot be subjected to three-dimensional positioning, are not fully considered in the existing processing technology, so that the existing system cannot produce qualified products. Aiming at the characteristics of continuous image scanning, small image overlapping, low IMU precision and the like acquired by a domestic area array scanning imaging system, the method provided by the invention introduces DEM data auxiliary processing, finds balance between processing efficiency and achievement precision, can more quickly and accurately complete image splicing, produces qualified orthoscopic image products, and effectively fills the blank of data processing of the domestic area array scanning imaging system.

Disclosure of Invention

Aiming at the defects of the prior art, the invention provides a method for producing an orthoimage aiming at an image acquired by a domestic area array sweep imaging system. Firstly, a reasonable local coordinate system is selected as a reference coordinate system, navigation data information is converted into the reference coordinate system, then homonymous points are automatically extracted from all images under the assistance of a DEM (digital elevation model) to form an image connection point network, then, the connection points are used for carrying out adjustment of the regional network based on the DEM, the external orientation parameter of each image is solved, and finally, the overall optimal orthoscopic image of a target region is produced based on an optimization theory.

In order to achieve the above object, the present invention provides a method for producing an orthoscopic image from an image acquired by a domestic area array scanning imaging system, comprising the steps of:

step 1, establishing a local coordinate system as a unified reference coordinate system, and converting navigation data information into the reference coordinate system;

step 2, automatically extracting the same-name points of the image under the assistance of the rough DEM;

step 2.1, under the assistance of the general DEM, acquiring the adjacent relation between the images according to the scanning imaging periodicity and the sequence;

2.2, according to the adjacent relation of the images obtained in the step 2.1, carrying out homonymy point matching by adopting a method of carrying out dense matching on the stereo images based on two types of expansion;

step 3, adjusting the local area network based on DEM assistance, taking the navigation data information converted in the step 1 as an initial adjustment value, and calculating the external orientation element parameter of each image;

step 3.1, calculating the area network range by using the POS information of all the images, and acquiring corresponding DEM data from the global SRTM data according to the area network range;

step 3.2, calculating a rotation matrix of each image by using the external orientation angle element, and using the angle element obtained in the step 1 as an initial value of calculation;

step 3.3, respectively calculating the intersection points of the light rays connected with the image points and the DEM by taking the single-beam same-name light rays as an observation unit to obtain corresponding object space coordinates;

step 3.4, calculating the difference between object space points corresponding to the image points with the same name, and if the difference is greater than the theoretical object space positioning precision alpha derived from the angle element precision, marking the image points with the same name as rough difference points without participating in adjustment; otherwise, taking the average value of the object space coordinates corresponding to the image points with the same name as a true value to participate in adjustment;

step 3.5, constructing an error equation point by point, and resolving an external orientation element;

step 3.6, repeating steps 3.2 to 3.5 until the change of each corner element of the image does not exceed 10^-6Radian and the sum of squares of all image point correction numbers is smaller than a set threshold lambda;

step 4, producing an overall optimal orthoimage;

step 4.1, calculating the ground coordinates corresponding to the four corner points of each original image;

step 4.2, selecting the optimal original image corresponding to the orthoimage to perform back projection;

and 4.3, producing the orthoimage by a back projection method.

Moreover, in step 1, the local Coordinate System adopts a fixed point Tangent Plane Coordinate System, i.e. a TPC Coordinate System (referred to as "TPC Coordinate System"), and the coordinates are defined as follows: selecting a geometric center point longitude L in a flight area under a WGS84 coordinate frame₀Latitude B₀As a plane origin, the height of the earth ellipsoid is 0 to be taken as an ellipsoid tangent plane, the east direction passing through the latitude line of the origin is defined as an X coordinate axis, the north direction passing through the longitude line of the origin is defined as a Y coordinate axis, and the perpendicular line passing through the origin and back to the geocentric is defined as a Z coordinate axis, as shown in figure 4Shown in the figure.

The TPC coordinate system and the image space coordinate system have a conversion relation as follows:

namely, the image space coordinate system → the loading body coordinate system → the IMU body coordinate system → the navigation coordinate system → the geocentric geostationary coordinate system → the TPC coordinate system, and the definition and transformation of each coordinate system are shown in FIG. 5.

The calculation formula of the rotation matrix between the coordinate systems is as follows:

in the formula (I), the compound is shown in the specification,is an image attitude rotation matrix under the TPC coordinate system,three external orientation angle elements of the image are taken as parameters to be solved;is a transformation matrix between geocentric geostationary coordinate system and TPC coordinate system, from longitude L of TPC coordinate origin₀Latitude B₀Defining;is a conversion matrix between a navigation coordinate system and a geocentric geostationary coordinate system, which is determined by the longitude L of the position of the load_zLatitude B_zDefining;is a rotation matrix between a navigation coordinate system and an IMU body coordinate system, and (phi, theta and psi) are pitching, rolling and yawing numerical values recorded by the IMU;is a transformation matrix between the load body coordinate system and the IMU body coordinate system, also known as a device placement matrix, placed on the load by the IMU_x、e_y、e_zDefining three angles;is a transformation matrix between the spatial coordinate system and the load body coordinate system.

Solving a rotation matrixThen, the inverse solution can obtain the attitude angle of the image in the TPC coordinate system, that is:

ω＝arcsin(-b₃) (9)

the GPS point position information can be losslessly transformed into any desired coordinate system using the earth coordinate definition model, and the TPC coordinates transformed from the original GPS coordinates (L, B, H) into the local tangent plane coordinate system are transformed as follows:

the original GPS data is first converted into the geocentric coordinate system:

X_c＝(N+H)cos B cos L (11)

Y_c＝(N+H)cos B sin L (12)

Z_c＝[N(1-e²)+H]sin B (13)

in the formula, N is the radius of the unitary mortise ring, a is the long semi-axis of the earth ellipsoid, b is the short semi-axis of the earth ellipsoid, (X)_c,Y_c,Z_c) Is the coordinates of the geocentric coordinate system.

Then the origin coordinate (L) of the local tangent plane coordinate system₀,B₀And 0) taking the expressions (11) to (13), obtaining the coordinates of the origin of the local tangent plane coordinate system in the geocentric coordinate system, and recording asThe calculation formula from the geocentric system coordinates to the TPC coordinate system is:

in the formula (I), the compound is shown in the specification,between the earth-centered earth-fixed coordinate system and the TPC coordinate systemA conversion matrix calculated by the formula (2); (X)_TPC,Y_TPC,Z_TPC) Are the coordinates of the TPC coordinate system.

I.e. the navigation information converted into the reference coordinate system.

In step 2, the adjacent relationship between the images is obtained according to the sweep imaging periodicity and the sequence, that is, N images obtained by sweeping one period are divided into a group, each image is numbered from 1 to N according to the sequence of the image in the sweep period, and the images which are continuous in sequence are necessarily adjacent images. The adjacent relation between each group can project the acquired images to the ground represented by the global DEM according to the GPS and IMU information, then the image overlapping condition is obtained according to the ground range of each image, and the overlapping area of the original image can be inversely calculated according to the overlapping area, so that the adjacent relation of the images among the groups is confirmed.

In step 3.3, the single beam of light is used as an observation unit, and the intersection points of each connecting point and the DEM are respectively calculated to obtain corresponding object coordinates, wherein the calculation formula is as follows:

in the formula (X)_g,Y_g) Object coordinates (i.e., TCP coordinates) that are the intersection points; z is a radical of_demDEM elevation as intersection point; x_s、Y_s、Z_sIs the exterior orientation line element of the image, obtained in step 1 (X)_TPC,Y_TPC,Z_TPC) Can be used as an initial value of the exterior orientation line element; (x, y) are image point coordinates; f is the camera focal length.

In step 3.4, the theoretical object positioning accuracy α is calculated as follows:

in the formula, H is the height from the photographing position to the ground, epsilon is the IMU angle element precision of the imaging system, and theta is the included angle between the main optical axis and the gravity direction.

Moreover, the point-by-point construction error equation in the step 3.5 is as follows:

wherein (x, y) is the coordinates of the image point, (x), (y) are the coordinates of the image point obtained by using the initial value of the unknown number, v_x、v_yIs the variable of the error, and is,and d represents error correction, and the external orientation element of each image can be updated after the correction number of the external orientation element is solved.

Moreover, the method for producing an orthoimage by the back projection method in the step 4.3 is as follows: let the coordinate of any point P on the ortho-image be (X ', Y '), from the ground coordinate (X ') of the contour point at the lower left corner of the ortho-image₀,Y₀) Calculating the ground coordinates (X, Y) corresponding to the point P according to the orthographic image scale denominator M as follows:

X＝X₀+M×X′ (22)

Y＝Y₀+M×Y′ (23)

calculating the coordinate P (x, y) of the point P on the ortho image corresponding to the corresponding point on the original image as:

in the formula, X₀、Y₀The minimum value in the ground coordinates corresponding to the angular points of all the original images is taken as Z, the elevation of the P point is taken as Z, and the Z is obtained by interpolation of the DEM; a is_i,b_i,c_i(i ═ 1,2,3) is the rotation matrix element of the image, X_s、Y_s、Z_sIs the line element of the image, and f is the focal length of the original image.

Because the obtained image point coordinate does not exactly fall in the center of the pixel, and the gray interpolation is needed for the purpose, the invention adopts a bilinear interpolation method to obtain the gray value g (X, Y) of the image point P, and the gray value of the image point is assigned to the corresponding pixel point P (X ', Y') on the orthoimage, namely:

G(X′,Y′)＝g(x,y) (26)

and sequentially carrying out the operation on each pixel to obtain the orthoimage.

Compared with the prior art, the invention has the following advantages:

(1) the domestic area array scanning imaging system continuously scans and images, utilizes the characteristic that images in an imaging period are continuous, namely adjacent, and uses GPS/IMU and DEM data to obtain an image adjacent matrix and an overlapping area, so that image matching is pertinently carried out, and image connection points are quickly extracted.

(2) The POS precision of the domestic area array scanning imaging system is not good enough, the POS precision is limited when the POS system is used for direct geographic positioning, and an orthoimage cannot be directly produced, and adjustment processing is still needed; however, the system has small image overlapping degree and weak intersection condition, the GPS/IMU parameters are introduced as initial values, and the problems that the adjustment calculation cannot be directly carried out by the traditional photogrammetry by using DEM data to assist the adjustment are solved.

(3) The method can complete image splicing more quickly and accurately, produce qualified orthoscopic image products, and fill the blank of data processing of a domestic area array scanning imaging system.

Drawings

Fig. 1 is a schematic diagram of core components of a domestic area array sweep imaging system according to an embodiment of the present invention.

FIG. 2 is a flowchart of an embodiment of the present invention.

Fig. 3 is a diagram illustrating international standard definition of IMU rotation angle according to an embodiment of the present invention.

Fig. 4 is a schematic diagram illustrating the definition of a local tangent plane coordinate system (TPC coordinate system) according to an embodiment of the present invention.

Fig. 5 is a schematic diagram of converting the IMU coordinate system into the TPC coordinate system according to an embodiment of the present invention.

FIG. 6 is a diagram illustrating a GPS/IMU parameter conversion result of the sweep imaging of the domestic area array sweep imaging system according to an embodiment of the present invention.

FIG. 7 is a diagram illustrating an image adjacency matrix according to an embodiment of the invention.

FIG. 8 shows an embodiment of an ortho-image data result.

FIG. 9 is a partially enlarged view of the ortho-image data stitching region according to the embodiment of the present invention.

Detailed Description

The invention provides a method for producing an ortho-image aiming at an image acquired by a domestic area array sweep imaging system.

The technical solution of the present invention is further explained with reference to the drawings and the embodiments.

As shown in fig. 2, the process of the embodiment of the present invention includes the following steps:

step 1, establishing a local coordinate system as a unified reference coordinate system, and converting navigation data information into the reference coordinate system.

The domestic area array sweep imaging system is provided with GPS equipment and IMU equipment, the GPS equipment records aircraft position information including longitude, latitude and elevation, the IMU records aircraft attitude information including pitching, rolling and yawing 3 angle values, and an IMU coordinate system is defined according to international standards, as shown in figure 3.

The local Coordinate System adopts a fixed point Tangent Plane Coordinate System, i.e. a TPC Coordinate System (indexed plant Coordinate System)) The coordinates are defined as follows: selecting a geometric center point longitude L in a flight area under a WGS84 coordinate frame₀Latitude B₀As a plane origin, the ellipsoidal height of the earth is 0 to serve as an ellipsoidal tangential plane, the east direction passing through the latitude line of the origin is defined as an X coordinate axis, the north direction passing through the longitude line of the origin is defined as a Y coordinate axis, and the perpendicular line passing through the origin and back to the geocentric is defined as a Z coordinate axis, as shown in fig. 4.

The TPC coordinate system and the image space coordinate system have a conversion relation as follows:

namely, the image space coordinate system → the loading body coordinate system → the IMU body coordinate system → the navigation coordinate system → the geocentric geostationary coordinate system → the TPC coordinate system, and the definition and transformation of each coordinate system are shown in FIG. 5.

The calculation formula of the rotation matrix between the coordinate systems is as follows:

in the formula (I), the compound is shown in the specification,is an image attitude rotation matrix under the TPC coordinate system,three external orientation angle elements of the image are taken as parameters to be solved;is a transformation matrix between geocentric geostationary coordinate system and TPC coordinate system, from longitude L of TPC coordinate origin₀Latitude B₀Defining;is a conversion matrix between a navigation coordinate system and a geocentric geostationary coordinate system, which is determined by the longitude L of the position of the load_zLatitude B_zDefining;is a rotation matrix between a navigation coordinate system and an IMU body coordinate system, and (phi, theta and psi) are pitching, rolling and yawing numerical values recorded by the IMU;is a transformation matrix between the load body coordinate system and the IMU body coordinate system, also known as a device placement matrix, placed on the load by the IMU_x、e_y、e_zDefining three angles;is a transformation matrix between the spatial coordinate system and the load body coordinate system.

Solving a rotation matrixThen, the inverse solution can obtain the attitude angle of the image in the TPC coordinate system, that is:

ω＝arcsin(-b₃) (9)

the GPS point position information can be losslessly transformed into any desired coordinate system using the earth coordinate definition model, and the TPC coordinates transformed from the original GPS coordinates (L, B, H) into the local tangent plane coordinate system are transformed as follows:

the original GPS data is first converted into the geocentric coordinate system:

X_c＝(N+H)cos B cos L (11)

Y_c＝(N+H)cos B sin L (12)

Z_c＝[N(1-e²)+H]sin B (13)

in the formula, N is the radius of the unitary mortise ring, a is the long semi-axis of the earth ellipsoid, b is the short semi-axis of the earth ellipsoid, (X)_c,Y_c,Z_c) Is the coordinates of the geocentric coordinate system.

Then the origin coordinate (L) of the local tangent plane coordinate system₀,B₀And 0) taking the expressions (11) to (13), obtaining the coordinates of the origin of the local tangent plane coordinate system in the geocentric coordinate system, and recording asThe calculation formula from the geocentric system coordinates to the TPC coordinate system is:

in the formula (I), the compound is shown in the specification,is a transformation matrix between the geocentric coordinate system and the TPC coordinate system, and can be obtained by calculation of the formula (2); (X)_TPC,Y_TPC,Z_TPC) Are the coordinates of the TPC coordinate system.

I.e. the navigation information converted into the reference coordinate system.

And 2, automatically extracting the image homonymous points under the assistance of the general DEM.

Aiming at the images acquired by a domestic area array scanning imaging system, the periodic scanning rule and GPS/IMU parameters can be fully utilized, and the global DEM data (mainly SRTM data) is matched to realize quick and effective image matching. The schematic diagram of the image of the area array sweep imaging system and the projection of the GPS and IMU parameters thereof to the global DEM data processing principle and effect is shown in FIG. 6.

And 2.1, acquiring the adjacent relation between the images according to the scanning imaging periodicity and the sequence under the assistance of the general DEM.

First, N (in this embodiment, N is 36) images obtained by sweeping one period are divided into a group, each image is numbered 1 to N according to the order in the sweeping period, and the sequentially consecutive images are necessarily adjacent images. The adjacent relation between each group can be obtained by projecting the acquired images onto the ground represented by the global DEM according to the GPS and IMU information, then obtaining the image overlapping condition according to the ground range of each image, or reversely calculating the overlapping area of the original image according to the overlapping area, so as to confirm the adjacent relation of the images between each group, and finally confirming the adjacency matrix as shown in FIG. 7.

And 2.2, carrying out homonymy point matching by adopting a method for carrying out dense matching on the stereo images based on two types of expansion according to the adjacent relation of the images obtained in the step 2.1.

And 3, carrying out adjustment of the area network based on DEM assistance, taking the navigation data information converted in the step 1 as an initial adjustment value, and calculating the external orientation element parameter of each image.

And 3.1, calculating the range of the local area network by using the POS information of all the images, and acquiring corresponding DEM data from the global SRTM data according to the range of the local area network.

And 3.2, calculating a rotation matrix of each image by using the external orientation angle elements, and using the angle elements obtained in the step 1 as initial values of calculation.

The rotation matrix of the image calculated by the external orientation angle element is shown as the following formula:

and 3.3, respectively calculating the intersection points of the light rays connected with the image points and the DEM by taking the single-beam same-name light rays as an observation unit to obtain corresponding object space coordinates, wherein the calculation formula is as follows:

in the formula (X)_g,Y_g) Object coordinates (i.e., TCP coordinates) that are the intersection points; z is a radical of_demDEM elevation as intersection point; x_s、Y_s、Z_sIs the exterior orientation line element of the image, obtained in step 1 (X)_TPC,Y_TPC,Z_TPC) Can be used as an initial value of the exterior orientation line element; (x, y) are image point coordinates; f is the camera focal length.

Step 3.4, calculating the difference between object space points corresponding to the image points with the same name, and if the difference is greater than the theoretical object space positioning precision alpha derived from the angle element precision, marking the image points with the same name as rough difference points without participating in adjustment; otherwise, taking the average value of the object coordinates corresponding to the image points with the same name as a true value to participate in adjustment.

If the IMU angle element precision of the imaging system is epsilon, and the included angle between the main optical axis and the gravity direction is theta, the theoretical object space positioning precision alpha is as follows:

in the formula, H is the height from the imaging position to the ground.

Step 3.5, constructing an error equation point by point, resolving an unknown number correction number, wherein the error equation is as follows:

wherein (x, y) is the coordinates of the image point, (x), (y) are the coordinates of the image point obtained by using the initial value of the unknown number, v_x、v_yIs the variable of the error, and is,and d represents error correction, and the external orientation element of each image can be updated after the correction number of the external orientation element is solved.

Step 3.6, repeating steps 3.2 to 3.5 until the change of each corner element of the image does not exceed 10^-6Radian and the sum of the squares of all the image point corrections is less than a set threshold lambda.

And 4, producing the integral optimal orthoimage.

And 4.1, calculating the ground coordinates corresponding to the four corner points of each original image.

And (3) calculating the ground coordinates corresponding to the four corner points of each original image to form a quadrangle by using the formulas (18) and (19), and solving the gravity center of each quadrangle.

And 4.2, selecting the optimal original image corresponding to the orthoimage to perform back projection.

And (3) setting the ground coordinates corresponding to any point P on the orthoimage as (X, Y), and selecting the original image corresponding to the center of gravity of the quadrangle closest to the point P (X, Y) as the optimal original image for back projection.

And 4.3, producing the orthoimage by a back projection method.

Let the coordinate of any point P on the ortho-image be (X ', Y '), from the ground coordinate (X ') of the contour point at the lower left corner of the ortho-image₀,Y₀) Calculating the ground coordinates (X, Y) corresponding to the point P according to the orthographic image scale denominator M as follows:

X＝X₀+M×X′ (22)

Y＝Y₀+M×Y′ (23)

calculating the coordinate P (x, y) of the point P on the ortho image corresponding to the corresponding point on the original image as:

in the formula, X₀、Y₀The minimum value in the ground coordinates corresponding to the angular points of all the original images is taken as Z, the elevation of the P point is taken as Z, and the Z is obtained by interpolation of the DEM; a is_i,b_i,c_i(i ═ 1,2,3) is the rotation matrix element of the image, X_s、Y_s、Z_sIs the line element of the image, and f is the focal length of the original image.

Because the obtained image point coordinate does not exactly fall in the center of the pixel, and the gray interpolation is needed for the purpose, the invention adopts a bilinear interpolation method to obtain the gray value g (X, Y) of the image point P, and the gray value of the image point is assigned to the corresponding pixel point P (X ', Y') on the orthoimage, namely:

G(X′,Y′)＝g(x,y) (26)

and sequentially carrying out the operation on each pixel to obtain the orthoimage.

The step 4.1 to the step 4.3 are adopted to produce the integral orthoimage, each pixel is only calculated once, redundant calculation is avoided, the processing efficiency is greatly improved, and the integral orthoimage finally obtained is also optimized because each pixel selects the optimal original image.

In specific implementation, the above process can adopt computer software technology to realize automatic operation process.

The specific embodiments described herein are merely illustrative of the spirit of the invention. Various modifications or additions may be made to the described embodiments or alternatives may be employed by those skilled in the art without departing from the spirit or ambit of the invention as defined in the appended claims.