Self-adaptive layered and blocked true orthophoto orthorectification method

Abstract

In order to solve the problems of low exterior orientation element parameters and low correction efficiency which are solved by using ground control points in urban areas, the invention provides a self-adaptive layering and blocking true orthophoto orthorectification method which comprises the following steps: firstly, layering a building by adopting a self-adaptive radial scanning algorithm, and extracting a building corner point of each layer as a control point; calculating external orientation element parameters according to the control points of each floor, and orthogonally correcting the building; secondly, partitioning the original image according to a quadtree principle, determining the size of a generated orthoimage, then calculating geodetic coordinates of partitioned image blocks, inversely calculating image point coordinates of the original image, and endowing gray values of the image point coordinates to corresponding positions on the orthoimage; and finally, fusing the two parts of ortho-rectification images to generate a true orthophoto.

Description

Technical Field

[0001] The invention relates to the field of remote sensing image processing, and in particular to true orthogonal image generation. Background Art

[0002] Remote sensing images have the advantages of large monitoring area, fast data acquisition speed and high processing efficiency. Among them, remote sensing orthophotos have been widely used in mapping, land surveying, planning and design, etc. The traditional method uses digital differential correction method to solve and correct pixel by pixel to generate remote sensing orthophotos.

[0003] In the process of traditional remote sensing image orthorectification, the control points are all on the ground. However, in urban areas, high-rise buildings are relatively dense, and it is difficult to solve the ground control points to accurately calculate the exterior orientation elements (EOPs, Exterior Orientation Parameters). However, it is not realistic to collect control points on the roof of the building. At the same time, the traditional digital differential correction method operates on a single pixel, which is time-consuming and labor-intensive and cannot generate orthophotos in real time. The above problems seriously affect the accuracy and efficiency of orthophoto generation. Therefore, scholars at home and abroad have carried out a lot of research on the generation of true orthophotos.

[0004] In 1998, Amhar et al. from the National Research and Innovation Agency published a paper titled "Using 3D building models and traditional DTM to generate true orthophotos" in the journal "International Archives of Photogrammetry and Remote Sensing", proposing to generate orthophotos of building areas based on DBM and to generate orthophotos of terrain areas using DTM.

[0005] In 2021, Marseti et al. of SPACE SI published a paper titled "Automatic generation of true orthophotos from very high-resolution panchromatic satellite images based on occlusion detection based on image incidence angle" in the journal "IEEE Journal of Selected Topics inApplied Earth Observations and Remote Sensing", proposing a method for automatically generating true orthophotos from optical super-resolution satellite images. The automatic workflow includes five modules, starting with ground control point extraction, followed by geometric processing of image patches, occlusion detection, orthophotos, and finally generating a true orthophoto.

[0006] In 2003, Sheng et al. from Northeastern University published a paper titled "Generating true orthophotos of forest areas from large-scale aerial photographs" in the journal Photogrammetric Engineering & Remote Sensing. They used a model-based approach to reconstruct a true orthophoto model (CSM) instead of DEM to generate real photos of forest scenes. This method mainly focused on the occlusion and distortion caused by trees in forest areas.

[0007] The invention of Liu Jianming and others from Beijing Daoda Tianji Technology Co., Ltd., "Method and device for generating true orthogonal images based on oblique images of unmanned aerial vehicles" (patent number CN113566793A), generates true orthogonal images by obtaining the three-dimensional grid of the scene corresponding to the oblique images of unmanned aerial vehicles, mapping the texture, and orthographic projection.

[0008] The "True Orthogonal Image Generation Method and Device" (patent number CN112419443A) invented by Zhao Shuai and others from China Coal Survey and Remote Sensing Group Co., Ltd. is applied to the image generation server in the pixel factory. It generates DSM data based on the original aerial survey images, eliminates noise, and generates true orthogonal image data.

[0009] Shenzhen Pegasus Robotics Technology Co., Ltd., National Marine Environment Monitoring Center Wang Bangsong and others invented "A true radiographic image generation method" (patent number CN106875364A), which uses fully automated true radiographic processing of optical images from micro-UAVs to effectively solve the problems of collapsed and distorted high-resolution urban buildings and color inconsistency between images, and can realize true radiographic processing of images with large overlap.

[0010] The above technologies collect control points on the ground. It is difficult to calculate accurate EOPs using ground control points, and it is impossible to collect control points on the roof of a building. At the same time, the existing technology operates on a single pixel, which is time-consuming and labor-intensive, and cannot generate orthophotos in real time. Summary of the invention

[0011] In view of the above problems, the present invention proposes an adaptive layered and blocked true orthorectification method for orthorectification of true orthogonal images, which proposes an adaptive radial scanning algorithm to layer the buildings based on the building height difference of the digital building model (DBM) and the image plane distance of the buildings from the photography center, extracts the building corner points of each layer as control points, and orthorectifies the buildings based on the EOPs calculated from the control points of each layer.

[0012] The technical solution of the present invention is realized as follows:

[0013] The adaptive layered and block orthorectification algorithm is divided into two parts. One part performs layered orthorectification on buildings, and the other part performs block orthorectification on terrain. Finally, the two orthorectified images are fused to generate true orthorectified images.

[0014] S1: Building layered orthorectification

[0015] According to the law of error propagation, the EOPs error calculated by the control points with large elevation span will be correspondingly large. First, the buildings are layered based on the radial scanning algorithm, and then the corner points of each building are extracted as control points to calculate the EOPs. The buildings in each layer are orthorectified according to the EOPs calculated for each layer.

[0016] The image point displacement caused by traditional terrain undulation is shown in formula (1):

[0017] (1)

[0018] Where: is the image point displacement distance, is the radial distance from the image principal point to the displacement image point, is the elevation of the displaced object point, For the flight altitude.

[0019] From formula (1), we can see that when the flight altitude is constant, the image point displacement distance The size of the radial distance And object point elevation Decide when When the image point displacement is the maximum, To this end, the present invention proposes a building layering method: using a radial scanning algorithm to find Get the building with the maximum value and calculate the elevation interval based on it , to achieve building stratification.

[0020] First, a virtual orthographic projection plane is established, and the top corner points of the building are projected onto the virtual plane, and the bottom points are projected onto the virtual plane and set as , project the corner points of the building top surface onto the virtual plane The distance is set to , the image point displacement distance is required The maximum value of The maximum value of . However, direct calculation The maximum value of is very difficult, but according to and Proportional relationship, calculation The maximum value problem can be transformed into The maximum value of .

[0021] The radial scanning algorithm is used to solve The maximum value process is as follows:

[0022] (1) Calculate the value in a certain scanning direction The maximum value D

[0023] (2)

[0024] (2) Calculate the The maximum value D(k×△θ):

[0025] (3)

[0026] Where k is the number of increments and θ is the angle of each increment.

[0027] (3) Find the value of After the maximum value is obtained, the maximum value D among the maximum values ​​is obtained again. max :

[0028] (4)

[0029] (4) Get all the buildings The maximum value, the corresponding building elevation is , let the radial distance from the image point of the building to the image principal point be , maximum pixel displacement Calculate using formula (5)

[0030] (5)

[0031] (5) After obtaining the maximum image point displacement, calculate the maximum projection error δd of the orthophoto max :

[0032] (6)

[0033] (6) According to the photogrammetry principle and formula (6), the maximum elevation difference that makes the image point displacement appear less than one pixel on the orthophoto is: The calculation formula (7)

[0034] (7)

[0035] (7) The maximum elevation difference As a layered interval for building layers , for all buildings in DBM, see formula (8):

[0036] (8)

[0037] In the formula, is the minimum building height in DBM, and T is the number of building floors.

[0038] S2: Terrain Block Orthorectification

[0039] The principle of terrain block orthorectification is: first, the image is divided into blocks according to the quadtree principle, then the size of the generated orthophoto is determined, and then each image block is orthorectified, and the collinear condition equation (CCE) is used to orthorectify the image block.

[0040] First, the complexity of the features such as type, density, and size of each block is calculated based on the quadtree principle, and the blocks with scores higher than the threshold are divided into smaller sub-blocks. It can divide the image into multiple non-overlapping image blocks of different sizes. The complexity of the ground feature is calculated using formula (9): .

[0041] (9)

[0042] In the formula, Represents the grayscale standard deviation within the block, r u is the block number (u∈[1,4]), Represents the average grayscale value within the block. and The calculation methods of are shown in formulas (10) and (11) respectively.

[0043] (10)

[0044] Where g is the grayscale value and n is the size of the block.

[0045] (11)

[0046] In order to make the blocks as large as possible, a threshold of the complexity of the ground features is set The image block complexity C is compared with the adaptive threshold T to determine whether it needs to be divided into blocks again.

[0047] The block principle calculation is shown in formula (12).

[0048] (12)

[0049] In the formula represents the image block, x represents the xth image block, T represents the threshold, Represents an image block The complexity of The smaller the image block, the fewer features and the smoother it is. On the contrary, the image block is complex and needs to be divided again. , then it represents the image block The complexity is very small and no further division is required. Otherwise, it will continue to be divided into four smaller image blocks. And calculate one by one in the same way until no more division is required.

[0050] Depending on the resolution of the orthophoto, the extent of the DTM-based DOM is determined and blanks are created for it. The orthorectification process registers the original image into the selected map-based coordinate system, and the size of the original image always changes. In order to correctly set the storage space requirements when programming, the size of the space occupied by the resulting image must be determined in advance (upper left, lower left, upper right and lower right corners). These procedures are as follows:

[0051] Determination of the four corner coordinates: For a given ground resolution in the x and y directions of the original image and , assuming that the plane coordinates of any ground control point (GCP) are , whose corresponding position in the original image plane is Then, the coordinates of the four corner points can be determined conventionally. For example, for corner 1, its coordinates can be calculated by the following formula

[0052] (13)

[0053] Other angles can be calculated accordingly.

[0054] Determine the minimum and maximum coordinates from the four corners. For example, for the minimum x coordinate, you can calculate it by the following formula

[0055] (14)

[0056] This method can be used to find , , .

[0057] The size of the generated orthophoto is determined by the following formula:

[0058] (15)

[0059] After determining the range of the generated orthophoto, the CCE indirect method is used to orthorectify the above divided image blocks. The CCE indirect method is a method of obtaining scanning coordinates using geodetic coordinates. The CCE-based indirect method process can be divided into the following four parts:

[0060] (1) Calculate the geodetic coordinates of orthophoto pixels.

[0061] (16)

[0062] in is the geodetic coordinate of the pixel in the orthophoto; It is the geodetic coordinate of the lower left corner edge point of the orthophoto; is the scale denominator; are the column and row coordinates of the pixel in the orthophoto, is the sampling interval for columns and rows.

[0063] (2) Calculate the image plane coordinates.

[0064] (17)

[0065] in represents the camera's internal orientation parameter (IOP), is the spatial coordinate of the imaging position, Represents the rotation matrix, called the Exterior Orientation Parameters (EOP). The Z coordinate can be obtained by interpolating existing DTM data.

[0066] (3) Calculate the scanning coordinates.

[0067] (18)

[0068] in is the scanning coordinate, are the affine transformation parameters.

[0069] (4) Assign the grayscale value at the coordinates of the image point on the original image to the corresponding position on the orthophoto. Since the image point position of the pixel block corresponding to the original image may not fall at the center of the pixel, bilinear grayscale interpolation is required. Finally, the grayscale value on the original image is assigned to the corrected image. The above correction is performed on each pixel block in turn to complete the pixel block differential correction to generate a digital orthophoto.

[0070] Compared with the prior art, the present invention has significant improvements in that:

[0071] (1) The present invention divides high-rise buildings into layers, extracts the building corner points of each layer as control points, and performs orthorectification on the buildings of each layer respectively, so as to avoid inaccurate calculation of EOPs due to excessive difference in elevation of control points.

[0072] (2) Compared with the traditional pixel-by-pixel differential operation, the present invention uses the image block as a correction unit, which can improve the speed of orthorectification. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 The adaptive layered and block orthorectification method of the present invention

[0074] Figure 2 The building layered radial scanning algorithm of the present invention

[0075] Figure 3 The quadtree image block of the present invention

[0076] Figure 4 It is the pixel block correction of terrain orthorectification of the present invention DETAILED DESCRIPTION

[0077] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, preferred embodiments are given below and the specific implementation of the present invention is further described in detail in conjunction with the accompanying drawings.

[0078] Example:

[0079] Combination Figure 1 , the implementation process of the entire invention is described as follows:

[0080] (1) Building layer orthorectification

[0081] The present invention uses aerial image data, digital ground model, and building model data of Denver. Figure 1 , the DBM data is layered using the radial scanning algorithm. According to the three main steps of the radial scanning algorithm (see formulas (2-4)), the building elevation corresponding to the maximum image point displacement and the radial distance from the corresponding image point to the image principal point are calculated. According to the resolution relationship between the original image and the orthophoto image, the maximum elevation difference can be calculated (see formula (7)). The DBM building data is layered according to the maximum elevation difference (see formula (8)). The top corner points of each building are extracted as control points to solve the EOPs parameters and perform orthorectification on the building.

[0082] (2) Terrain orthorectification

[0083] First, the complexity of the features such as the type, density, and size of each block is calculated based on the quadtree principle, and the blocks with scores higher than the threshold are divided into smaller sub-blocks. It can divide the image into multiple non-overlapping image blocks of different sizes. Secondly, according to the resolution of the orthophoto, the range of the DTM-based DOM is determined and blanks are created for it. After determining the range of the generated orthophoto, the CCE indirect method is used to orthorectify the above-mentioned divided image blocks. Orthorectify each image block until the orthorectification of the terrain part of the entire image is completed.

[0084] Finally, the orthorectification results of building layers and terrain are integrated to generate an orthophoto.

[0085] Comparison of experimental results with traditional methods:

[0086] 1. Taking images with a resolution of 1 meter as experimental data, the proposed algorithm takes 10.02 seconds to correct the image using a 5*5 pixel block size, while the traditional pixel-by-pixel differential correction method takes 14.81 seconds to correct the image.

[0087] 2. 18 check points were collected in the experimental area to verify the geometric accuracy of the orthophoto. The orthophoto results generated by the commonly used digital differential correction method were used as a control. The pixel mean error of the proposed algorithm in the X direction was 1.27 pixel, and the pixel mean error in the Y direction was 0.76 pixel. The pixel mean error of the traditional digital differential correction method in the X direction was 4.50 pixel, and the pixel mean error in the Y direction was 2.12 pixel.

[0088] Although the embodiments of the present invention have been disclosed as above, they are not limited to the applications listed in the specification and the implementation modes, and they can be fully applied to various fields suitable for the present invention. For those familiar with the art, additional modifications can be easily implemented. Therefore, without departing from the general concept defined by the claims and the scope of equivalents, the present invention is not limited to the specific details and the illustrations shown and described herein.

Claims

1. Adaptive hierarchical and block-based true orthorectification method, characterized by It is divided into two parts, one part is to perform layered orthorectification on buildings, and the other part is to perform block-by-block orthorectification on terrain. Finally, the two parts of the orthorectified images are fused to generate a true orthorectified image.

2. The method for orthorectification of self-adaptive hierarchical and block-based true orthorectification images according to claim 1, further comprising: Building layered orthorectification process: According to the law of error propagation, the EOPs error calculated by the control points with large elevation span will be correspondingly large. First, the buildings are layered based on the radial scanning algorithm, and then the corner points of each building are extracted as control points to calculate the EOPs. The orthorectification of each layered building is performed based on the EOPs calculated for each layer. The image point displacement caused by traditional terrain undulation is shown in formula (1): (1) Where: is the image point displacement distance, is the radial distance from the image principal point to the displacement image point, is the elevation of the displaced object point, is the flight altitude; From formula (1), we can see that when the flight altitude is constant, the image point displacement distance The size of the radial distance And object point elevation Decide when When the image point displacement is the maximum, The maximum value of; To this end, the present invention proposes a building layering method: using a radial scanning algorithm to find Get the building with the maximum value and calculate the elevation interval based on it , realize building stratification; First, a virtual orthographic projection plane is established, and the top corner points of the building are projected onto the virtual plane, and the bottom points are projected onto the virtual plane and set as , project the corner points of the building top surface onto the virtual plane The distance is set to , the image point displacement distance is required The maximum value of The maximum value of The maximum value of is very difficult, but according to and Proportional relationship, calculation The maximum value problem can be transformed into The maximum value of The radial scanning algorithm is used to solve The maximum value process is as follows: (1) Calculate the value in a certain scanning direction The maximum value D (2) (2) Calculate the The maximum value D(k×△θ): (3) In the formula, k is the number of increments, and θ is the angle of each increment; (3) Find the value of each scanning direction After the maximum value is obtained, the maximum value D among the maximum values ​​is obtained again. max : (4) (4) Get all the buildings The maximum value, the corresponding building elevation is , the radial distance from the image point to the image principal point of the building is , maximum pixel displacement Calculate using formula (5) (5) (5) After obtaining the maximum image point displacement, calculate the maximum projection error δd of the orthophoto max : (6) (6) According to the photogrammetry principle and formula (6), the maximum elevation difference that makes the image point displacement appear less than one pixel on the orthophoto is: The calculation formula (7) (7) (7) The maximum elevation difference As a layered interval for building layers , for all buildings in DBM, see formula (8): (8) In the formula, is the minimum building height in DBM, and T is the number of building floors.

3. The method for orthorectification of self-adaptive hierarchical and block-based true orthorectification images according to claim 1, further comprising: Terrain block orthorectification process: The principle of terrain block orthorectification is: first, the image is divided into blocks according to the quadtree principle, then the size of the generated orthoimage is determined, and then each image block is orthorectified, and the collinear condition equation (CCE) is used to orthorectify the image block; The main steps are as follows: First, the complexity of the features such as type, density, and size of the objects in each block is calculated based on the quadtree principle, and the blocks with scores higher than the threshold are divided into smaller sub-blocks; it can divide the image into multiple non-overlapping image blocks of different sizes; the complexity of the features of the objects is calculated using formula (9): ; (9) In the formula, represents the grayscale standard deviation within the block, ru is the number of blocks (u∈[1,4]), Represents the average grayscale value within the block; and The calculation methods of are shown in formulas (10) and (11) respectively; (10) Where g is the gray value and n is the size of the block; (11) In order to make the blocks as large as possible, a threshold of the complexity of the ground features is set Perform block division; compare the image block complexity C with the adaptive threshold T to determine whether it needs to be divided into blocks again; The calculation of the block principle is shown in formula (12); (12) In the formula represents the image block, x represents the xth image block, T represents the threshold, Represents an image block The complexity of The smaller the image block, the fewer features and the smoother it is. On the contrary, the image block is complex and needs to be divided again. , then it represents the image block The complexity is very small and no further division is required; otherwise, it will continue to be divided into four smaller image blocks And calculate them one by one in the same way until no more division is needed; Depending on the resolution of the orthophoto, determine the extent of the DTM-based DOM and create blanks for it; the orthorectification process registers the original image to the selected map-based coordinate system, and the size of the original image always changes; in order to correctly set the storage space requirements when programming, the size of the space occupied by the generated image must be determined in advance (upper left corner, lower left corner, upper right corner and lower right corner); these procedures are as follows: Determination of the four corner coordinates: For a given ground resolution in the x and y directions in the original image and , assuming that the plane coordinates of any ground control point (GCP) are , whose corresponding position in the original image plane is ; Then, the coordinates of the four corner points can be determined conventionally; for example, for corner 1, its coordinates can be calculated by the following formula (13) The other angles can be calculated accordingly; Determine the minimum and maximum coordinates from the four corners; for example, the minimum x coordinate can be calculated by the following formula (14) This method can be used to find , , ; The size of the generated orthophoto can be calculated by the following formula: (15) After determining the range of the generated orthophoto, the CCE indirect method is used to orthorectify the above-mentioned image blocks; the CCE indirect method is a method of obtaining scanning coordinates using geodetic coordinates; the CCE-based indirect method process can be divided into the following four parts: (1) Calculate the geodetic coordinates of orthophoto pixels (16) in is the geodetic coordinate of the pixel in the orthophoto; It is the geodetic coordinate of the lower left corner edge point of the orthophoto; is the scale denominator; are the column and row coordinates of the pixel in the orthophoto, is the sampling interval of columns and rows; (2) Calculate image plane coordinates (17a in represents the camera's internal orientation parameter (IOP), is the spatial coordinate of the imaging position, Represents the rotation matrix, called the Exterior Orientation Parameters (EOP); the Z coordinate can be obtained by interpolating the existing DTM data; (3) Calculate scanning coordinates (18) in is the scanning coordinate, is the affine transformation parameter; (4) Assign the grayscale value of the image point coordinates on the original image to the corresponding position on the orthophoto; since the pixel block The image point position corresponding to the original image may not fall at the center of the pixel. Therefore, grayscale interpolation is performed using formula (18). Finally, the grayscale value of the original image is assigned to the corrected image. The above correction is performed on each pixel block in turn to complete the pixel block differential correction and generate a digital orthophoto.

Images