A method for generating orthophoto images for multi-machine cooperation

Abstract

The application discloses a kind of orthographic image generation methods for multi-machine cooperation, comprising: multi-machine cooperation gathers image and carries camera ID;Multiple-source image GPS coordinates are converted to ENU local coordinate system, and coordinates are unified;Grouping is constructed according to camera ID multiple parallel processing pipelines, and height model generation is carried out respectively by space triangulation solution;Based on digital elevation model, orthographic correction is carried out frame by frame, and single-frame orthographic image is generated;Each frame orthographic image is directly divided into tiles according to Mercator projection, and tiles are stored using hash mapping, directly fall into disk when Key does not exist, and then update and fall into disk after reading into old tile fusion when Key exists;For overlapping area, the optimal joint line is searched using energy function based on minimum cut, and multi-band fusion is carried out using Laplacian pyramid.The application does not need to generate global DOM large map, significantly reduces memory occupation, supports large-scale survey area and multi-source asynchronous input, effectively improves the system stability and processing efficiency of multi-machine cooperation orthographic image generation.

Description

Technical Field

[0001] This invention belongs to the field of UAV remote sensing mapping and image processing technology, specifically involving a multi-UAV collaborative orthophoto generation method. It is applicable to scenarios such as multi-UAV collaborative aerial surveying and large-scale survey area orthophoto (DOM) generation, and can be widely used in related fields such as UAV swarm collaborative mapping, large-scale terrain and geomorphology monitoring, and urban 3D modeling. Background Technology

[0002] In the field of UAV remote sensing and mapping, orthophotos are a core output, and their generation efficiency, processing scale, and system stability directly determine the quality and efficiency of surveying and mapping tasks. Traditional orthophoto generation methods are usually based on single-UAV aerial photography, and their standard processing flow is as follows: image acquisition, aerial triangulation / pose calculation, image projection, DOM generation, stitching and fusion, and tile segmentation. However, as surveying and mapping scenarios develop towards larger areas and higher resolutions, the traditional single-UAV mode is gradually failing to meet actual surveying and mapping needs. Therefore, multi-UAV collaborative acquisition and generation approaches are being gradually introduced, that is, acquiring data from large-scale scenes and generating orthophotos through multi-UAV collaboration to improve overall processing capabilities.

[0003] However, multi-drone collaborative orthophoto generation is not a simple transplant of the single-drone mode. In actual production operations, it faces the following problems: 1. Inconsistent coordinate system: The traditional single-drone mode has only a single data source, and the data coordinates are naturally unified. In the multi-drone collaborative mode, data from different sources must be unified under a single coordinate system; 2. Risk of memory overflow: In the traditional single-drone mode, the flight path is continuous and there is no geographical jump between two adjacent frames. However, in the multi-drone collaborative mode, two adjacent frames may come from different acquisition sources and have a large geographical gap. If memory management is still carried out in accordance with the traditional single-drone mode, it is easy to cause memory overflow; 3. Overlapping areas in images: In the multi-drone collaborative mode, the same ground feature may be observed multiple times by different drones (different internal parameters, different shooting angles, different time points), which can easily lead to problems such as geometric misalignment, texture distortion, or color banding in the resulting map. Summary of the Invention

[0004] To address the inherent problems of existing single-machine-based methods and simple porting of single-machine methods to multi-machine systems, this invention proposes a multi-machine collaborative orthophoto generation method. Its core objectives are: to adapt to the DOM generation requirements of ultra-large-scale survey areas, support asynchronous input and collaborative processing of multi-source data (multi-UAV), unify the coordinate system of multi-source data, eliminate global DOM memory overflow issues, and improve system stability through image boundary stitching and fusion. This method includes the following steps: Step 1: Multiple drones collaboratively collect images and generate single-frame orthophotos respectively; Step 2: For each frame of orthophoto, tile segmentation is performed directly according to Mercator projection. The segmented tiles are then stored using a hash map. Specifically, tile coordinates are converted into unique hash keys (Key), and the tile image data is used as the Value. For each tile in the current frame, if its Key does not exist in the hash map, the tile data is directly written to disk, and the correspondence between the Key and the tile's storage location is stored in the hash map. If the Key already exists, the corresponding old tile data is read from disk, merged with the current tile, and then re-written to disk. Step 3: After the finest layer tile of each frame image is generated, other lower-level tiles are generated layer by layer using the pyramid downsampling method; Step 4: For tiles with the same key from different pipes, perform full coverage or boundary overlap fusion based on minimum cut and multi-band fusion, depending on the coverage type. The energy function based on minimum cut is defined as: In the formula, The set of pixels in the overlapping region; It is a set of adjacent pixel pairs; For pixels The labels (0 indicates that the tile was taken from an old tile, and 1 indicates that it was taken from a new tile). Represents pixels Adjacent pixels; This represents the color data item, which measures the cost of color stitching for a single pixel. Represents the gradient data term, measuring the gradient stitching cost of a single pixel; This represents the smoothing term, which measures the cost of label transitions between adjacent pixels.

[0005] Furthermore, the multi-band fusion using the Laplace pyramid specifically includes: Laplacian pyramids LA and LB are constructed for the two images A and B to be fused, respectively, and a fusion weight mask is applied. Building the Gauss Pyramid ; At each level of the pyramid Above, the merged Laplace layer Calculate using the following formula: The merged layers of the Laplace Pyramid By restoring layer by layer, a final seamlessly blended image is obtained.

[0006] Further, step 1 includes: Step S1: Multiple drones collaboratively collect images. Each image frame carries at least GPS information, attitude information, camera intrinsic parameters, and a camera ID used to distinguish the source of the data collection. Step S2: Transform the GPS coordinates of all images to the local East-North-Sky ENU coordinate system with one of the images as the origin, thereby unifying the coordinate systems of the multi-source images; Step S3: Group the images according to camera ID, construct an independent parallel processing pipeline for each group of images, perform aerial triangulation calculation in each pipeline, and use the obtained 3D point cloud to generate a digital elevation model (DEM) of the survey area. Step S4: Based on the digital elevation model, perform frame-by-frame orthorectification of the original images in each pipeline to generate a single-frame orthorectified image.

[0007] Furthermore, in step S2, the coordinate unification process specifically includes: converting the WGS84 geographic coordinates of all images to the geocentric geofixed coordinate system ECEF as a transition; using one of the images as a reference point, calculating its ECEF coordinates; and converting the ECEF coordinates of each image point to ENU coordinates through a rotation matrix, thereby achieving coordinate unification of multi-source data.

[0008] Furthermore, in step S3, the aerial triangulation calculation performed in each parallel processing pipeline includes: extracting SIFT or ORB feature points from each group of images and performing feature matching between images; constructing an image connectivity graph, using bundle adjustment to optimize exterior orientation elements, and outputting the precise pose parameters of each frame of image; and generating a digital elevation model (DEM) of the survey area based on the matched 3D point cloud as a terrain constraint for orthorectification.

[0009] Furthermore, in step S4, orthorectifying each frame of the original image specifically includes: Generate a regular grid on the orthographic projection plane according to the set DOM resolution; For any ground point coordinates The corresponding elevation value Z is read from the DEM, and the pixel coordinates of the ground point on the original image are calculated using the collinearity equation. : In the formula, ( The spatial coordinates of the camera center The rotation matrix components are composed of exterior orientation angle elements. The focal length of the camera; The color value of the non-integer pixel coordinates obtained by inverse calculation is calculated using bilinear interpolation and then filled into the DOM grid points.

[0010] Furthermore, before substituting into the collinearity equation, we first determine the image point coordinates of the original image. Lens distortion correction is performed using the following formula: In the formula, r is the distance from the image point to the principal point. The radial distortion coefficient is... Tangential distortion coefficient.

[0011] Further, in step S5, the tile segmentation uses Mercator projection to divide the tile layers, and the tile size is 512×512 pixels; the finest layer determines the maximum value of the Mercator projection zoom; the Mercator projection coordinates corresponding to the tile are... Convert to hash key value Key using the following formula: The offset M must satisfy the following condition. To ensure the uniqueness of the Key, M takes the following values: Furthermore, in step 2, the pyramid downsampling method is as follows: for each input frame of image, after all the corresponding finest layer tiles have been generated, the resolution is immediately reduced layer by layer and Gaussian filtering is applied to generate other lower-level tiles; during the generation process, the required fine-layer tiles are reread from the disk, and only the hash map and the orthophoto data of the current frame image are maintained in memory.

[0012] Furthermore, in step 3, when the new tile completely covers the effective area of ​​the old tile in space, a full coverage operation is performed, that is, the new tile is directly used to replace the old tile in the storage device without performing seam line search and multi-band fusion.

[0013] Compared with the prior art, the beneficial effects of the present invention are: 1. Significantly reduced memory usage: By eliminating the need to maintain a large global DOM image, the problem of memory overflow in large-scale test area processing is fundamentally solved, improving system stability.

[0014] 2. Supports large-scale survey areas: Since memory usage is not strongly coupled with the size of the survey area, and is only constrained by disk storage capacity, there is theoretically no DOM size limit, which can adapt to the orthophoto generation needs of large-scale survey areas.

[0015] 3. Naturally supports multi-drone collaboration: Adapts to multi-drone collaborative data acquisition mode, enabling asynchronous input and collaborative processing of multi-source data. It also supports multi-source data fusion and distributed processing systems, fully leveraging the operational advantages of multi-drone clusters to improve data acquisition and processing efficiency. Attached Figure Description

[0016] Figure 1A flowchart of the orthophoto generation and processing workflow for multi-machine collaboration; Figure 2 For multi-aircraft collaborative data acquisition route maps; Figure 3 The results of multi-machine collaborative orthophoto generation. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of the embodiments of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0018] Example This embodiment provides a method for generating orthophotos in a multi-drone collaborative manner. Taking a collaborative acquisition system consisting of two drones (drone A and drone B) as an example, it generates orthophotos of an area of ​​approximately 25 km². 2 Aerial survey missions were conducted in the survey area. The flight altitude was set to 120 m, the ground resolution (GSD) was approximately 5 cm, the forward overlap was 80%, and the lateral overlap was 70%. The mission included the following steps: Step 1: Multi-machine collaborative image acquisition UAV A and UAV B are used simultaneously to conduct image acquisition operations, with each UAV flying along different preset routes. Each acquired image frame contains standard GPS information (longitude λ, latitude φ, altitude h), attitude information (yaw angle, pitch angle, roll angle), and camera intrinsic parameters, plus a "camera ID" field to distinguish the image source. In this embodiment, camera ID=1 for UAV A and camera ID=2 for UAV B. All image data is transmitted to the ground station processing system in real time or in task batches.

[0019] Step 2: Establish a unified coordinate system To eliminate coordinate discrepancies between multi-source data, the WGS84 geographic coordinates of all imagery were first transformed to the Local Northeast-Sky (ENU) coordinate system, using the centroid coordinate system (ECEF) as an intermediary during the transformation. The first imagery with complete metadata acquired by UAV A was taken as the reference point, and its corresponding latitude, longitude, and altitude were set as follows: Then its corresponding ECEF coordinates Represented as: Where N is the radius of curvature of the zonal loop, the calculation formula is: In the formula, This is the semi-major axis (equatorial radius) of the WGS84 ellipsoid. This is the first eccentricity of the WGS84 ellipsoid.

[0020] Let any image point in the WGS84 geographic coordinate system The corresponding ECEF coordinates are Then its position in the ENU coordinate system can be expressed as: Through the above transformation, multi-source image data are unified into the same local spatial coordinate system, laying the geometric foundation for subsequent processing. Note that the ENU coordinates of any image point are used to unify multi-source images into the same local coordinate system, enabling them to participate in subsequent aerial triangulation, orthorectification, and other calculations.

[0021] Step 3: Aerial triangulation and DEM generation Image data is grouped according to camera ID, and two independent processing pipelines are constructed: Pipeline 1 processes images with camera ID=1 (UAV A); Pipeline 2 processes images with camera ID=2 (UAV B). Aerial triangulation processing is then performed in parallel on each pipeline. With RTK, the image data from different pipelines have global coordinates, resulting in minimal pose misalignment in overlapping areas. Furthermore, graph-cut-based seam line finding and multi-band fusion strategies can mitigate visual differences caused by pose misalignment. The process involves extracting SIFT or ORB feature points, performing feature matching between images, constructing an image connectivity graph, optimizing exterior orientation elements using bundle adjustment, and outputting precise pose parameters for each frame. .

[0022] After aerial triangulation is completed, a digital elevation model (DEM) of the survey area is generated using the matched 3D point cloud. In this embodiment, the DEM resolution is set to 0.1m for terrain constraints in subsequent orthorectification.

[0023] Step 4: Orthophoto Generation (DOM) Based on the precise pose parameters and DEM data obtained from aerial triangulation, each frame of the original image is orthorectified to eliminate geometric deviations caused by terrain undulations, image tilt, and lens distortion, and to ensure color consistency of images collected by different UAVs. Finally, a single frame of standardized orthorectified image is generated as the image data source for subsequent tile segmentation.

[0024] Specifically, a regular grid is generated on the orthographic projection plane according to the set DOM resolution (5 cm in this embodiment). For any ground point coordinates... The corresponding elevation value Z is read from the DEM, and the pixel coordinates of the ground point on the original image are calculated using the collinearity equation. : in, These are the coordinates of the image point in the image coordinate system. Let be the object coordinates of the ground point, ( The spatial coordinates of the camera center Let F be the rotation matrix component consisting of exterior orientation elements, and F be the camera focal length.

[0025] Before substituting into the collinearity equation, it is usually necessary to first determine the image point coordinates of the original image. The formula for lens distortion correction is as follows: In the formula, r is the distance from the image point to the principal point. The radial distortion coefficient is... Tangential distortion coefficient.

[0026] Because the pixel coordinates are obtained by inverse calculation using the collinearity equation Typically, these are non-integer positions, requiring bilinear interpolation to calculate their color values. The formula is as follows: The final pixel value is obtained by using the grayscale values ​​of the four surrounding integer pixels, weighted according to their distance, and then filled into the DOM grid points. yes The pixel value of the top left corner is represented as ; yes The pixel in the top right corner has the following pixel value: ; yes The pixel value of the bottom left corner is represented as , yes The pixel in the bottom right corner has the following pixel value: .

[0027] Step 5: Real-time image segmentation and hash mapping To avoid the memory overflow problem caused by maintaining a large global DOM image in traditional methods, this embodiment does not perform global fusion and stitching, but instead directly performs unified tile segmentation on each frame of orthophoto from different sources.

[0028] Mercator projection was used to divide the tile layers, with a tile size of 512×512 pixels. At a zoom level of 20, each pixel on the equator represents approximately 7.46 cm of geographic length. For each frame of orthophoto, its covered area was mapped to the Mercator coordinate system. Only the standard tiles of the finest layer were calculated and segmented. This finest layer determined the maximum value of the Mercator projection zoom (i.e., the highest resolution layer), thus obtaining the Mercator projection coordinates corresponding to each tile. ,in It is the row and column number of the tile in the Mercator projection coordinate system, and This represents the pixel row and column number (position) in the image coordinate system. The subscript 'm' here is only used as a coordinate identifier and does not represent which frame of the image. It is converted into a hash key value using the following formula: To ensure key uniqueness and prevent hash collisions, the offset M must satisfy the following conditions: In this embodiment, we take... The image data of each tile is used as the Value, and key-value pairs are stored in a hash map structure to enable on-demand indexing.

[0029] Step 6: Tiles Placement and Replacement After each frame of image is tiled and hash-mapped, immediately perform the following disk write and update operations: For each tile generated in the current frame, its key is searched in the hash map: if the key does not exist, it means the tile is new, and the tile data is directly written to disk, and the correspondence between the key and the tile location information is stored in the hash map; if the key already exists, it means the tile needs to be updated. At this time, the corresponding old tile data is read from the storage device into memory, and merged with the newly segmented tile. After the merging is completed, the updated tile is written back to disk.

[0030] After the current image frame is processed and all the finest layer tiles are generated, a pyramid downsampling method (i.e., reducing the resolution layer by layer and applying a slight Gaussian filter to avoid jagged edges) is used to generate other lower-level tile data sequentially. It's important to note that lower-level tiles are not generated uniformly after all images have been processed. Instead, for each input image frame, after all the corresponding 512×512 finest layer tiles have been generated, lower-level tiles are generated immediately based on this. The processing sequence is: input an image frame, generate the finest layer tiles, then generate lower-level tiles, and finally input the next image frame, and so on. When generating lower-level tiles, the finest layer tiles used are reread from disk. Throughout the entire process, only the hash map structure and the orthophoto data corresponding to the currently processed image frame are maintained in memory; the global DOM large image is not saved, thus significantly reducing memory usage and supporting ultra-large-scale survey areas. In other words, the size of the map is allocated on demand based on the actual number of tiles being processed, and only the key-value pairs of the tiles to be processed or updated are stored. There is no need to maintain a complete DOM map, which fundamentally reduces memory usage.

[0031] Step 7: Overlapping Region Coverage and Adaptive Fusion In multi-machine collaborative mode, tiles generated by different pipelines may have the same Key value (i.e., correspond to the same geographic spatial region). This embodiment adopts an adaptive processing strategy based on coverage type, specifically divided into two cases: 1. Full Coverage (Full Update): When the newly generated tiles completely cover the effective area of ​​the original tiles in space (i.e., the image fills the entire view without transparent or invalid edges; since the tile size for Mercator projection is fixed at 512*512, the system only needs to determine whether the size of the newly generated tiles can fill this area), the system determines that the new data has higher timeliness or completeness and performs a full coverage operation, that is, directly replacing the old tiles in the storage device with the new tiles. 2. Boundary Overlap (Merge and Mosaic): When the new tiles only partially cover the original area and there is overlap between the new and old tiles, that is, if the newly generated tiles cannot fill the geographical area covered by the original tiles, the system determines that there is overlap (i.e., at the boundary of the viewpoints of the two acquisition sources), and the system triggers the tile merging mechanism. The mechanism consists of two core steps: First, the optimal seam line is searched in the overlapping area based on the color features of the image to avoid the edges of significant features; then, a multi-band fusion algorithm is used to perform smooth transition processing on both sides of the seam line to eliminate color breaks and splicing traces caused by differences in lighting and angle, and finally generate seamless fused tiles and place them on the disk.

[0032] To address the aforementioned boundary overlap issues, this invention further employs the following method to achieve high-quality splicing: Based on the minimum cut theory, an energy function is defined to find the optimal seam line. The optimal seam line is obtained by optimizing this energy function to minimize it. Let the set of pixels in the overlapping area of ​​the old and new tiles be denoted as... The set of adjacent pixel pairs (where adjacent pixel pairs refer to overlapping areas, but the adjacent pixels come from different tiles) is denoted as . Using labels Mark the pixels in the overlapping area of ​​old and new tiles, if the pixel is taken from the old tile pixel. If the pixel is taken from a new tile pixel Calculate the total cost of the seam, i.e., the total energy function. for: in, This represents the total cost of the entire seam, in order to minimize the total energy function. Optimize the total energy function for the objective Find a pixel allocation method for the overlapping areas of old and new tiles, so that... Minimum. This represents the color data item, measuring the cost of color stitching for a single pixel. This represents the gradient data term, measuring the gradient concatenation cost of a single pixel. This represents the smoothing term, which measures the cost of label transitions between adjacent pixels.

[0033] The color data item reflects how pixels are... The simplified formula for calculating the color difference between two tile images at a given pixel, based on the absolute value of the color difference at that pixel, can be expressed as: The smaller the difference in this value, the lower the cost for that pixel to become part of the seam. Similar to the color data, the gradient data term should identify regions with small gradient changes as seam lines, represented as: The smoothing effect ensures the continuity of the seams and prevents jagged or extremely fine noise. If two adjacent pixels... and They were assigned the same label ( This indicates that the seam does not pass through the middle of them, in which case the smoothing cost is 0. If their labels are different ( This indicates that the seam has been cut. and Between these points, a larger penalty value will be introduced. This forces the algorithm to minimize unnecessary label jumps, making the seams smoother.

[0034] Once the optimal seam line in the sense of minimum cut is found, fusion and splicing can be performed at the seam line. A multi-band fusion algorithm is used here, utilizing Gaussian and Laplace pyramids for layering. The operation steps are as follows: First, a Laplacian pyramid is constructed to decompose the overlapping images to be fused into multiple sub-bands of different frequencies. The low-frequency layer contains the basic color and brightness of the image, while the high-frequency layer contains the edge and texture details of the image. Then, layered fusion is performed. A smooth gradient transition is used in the low-frequency band to eliminate exposure differences and color bands, while the high-frequency band is directly selected or finely fused to preserve the sharpness of object edges and avoid blurring the image. Finally, reconstruction is performed to recombine the fused pyramid layers and restore the final seamless panoramic image.

[0035] Let the original image be Then the first Image of the Gaussian pyramid layer By the upper layer Gaussian smoothing (convolution) and downsampling yielded: in, It is a normalized odd Gaussian convolution kernel.

[0036] After obtaining the Gaussian pyramid, the Laplace pyramid is obtained by difference calculus. Laplace Pyramid Image The calculation formula is: in, This indicates an upsampling operation, that is, first sampling... Upsampling is performed to increase the size, followed by smoothing with a Gaussian kernel for the top layer. Directly equal to the top-level Gaussian image .

[0037] Suppose there are two images to be merged. and The Pyramid of Laplace was constructed separately. and Simultaneously, the fusion weight mask Building the Gauss Pyramid On each level of the pyramid Above, the merged Laplace layer Calculated using the following weighted formula: In the formula, , For two images in The Pyramid of Laplace, A Gaussian pyramid with a weighted mask.

[0038] The physical meaning of this formula lies in: at the low-frequency level The changes are very gradual, achieving a smooth transition in brightness and color and eliminating exposure differences; the original sharp edges and texture details of each image are preserved in the high-frequency layers, avoiding the blurring effect brought by traditional feathering blending.

[0039] Finally, the merged layers of the Laplace Pyramid By restoring layer by layer, a final seamless panoramic image is obtained. : The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for generating orthophotos for multi-machine collaboration, characterized in that, Includes the following steps: Step 1: Multiple drones collaboratively collect images and generate single-frame orthophotos respectively; Step 2: For each frame of orthophoto, tile segmentation is performed directly according to Mercator projection. The segmented tiles are then stored using a hash map. Specifically, tile coordinates are converted into unique hash keys (Key), and the tile image data is used as the Value. For each tile in the current frame, if its Key does not exist in the hash map, the tile data is directly written to disk, and the correspondence between the Key and the tile's storage location is stored in the hash map. If the Key already exists, the corresponding old tile data is read from disk, merged with the current tile, and then re-written to disk. Step 3: After the finest layer tile of each frame image is generated, other lower-level tiles are generated layer by layer using the pyramid downsampling method; Step 4: For tiles with the same key from different pipes, perform full coverage or boundary overlap fusion based on minimum cut and multi-band fusion, depending on the coverage type. The energy function based on minimum cut is defined as: In the formula, The set of pixels in the overlapping region; It is a set of adjacent pixel pairs; For pixels The labels (0 indicates that the tile was taken from an old tile, and 1 indicates that it was taken from a new tile). Represents pixels Adjacent pixels; This represents the color data item, which measures the cost of color stitching for a single pixel. Represents the gradient data term, measuring the gradient stitching cost of a single pixel; This represents the smoothing term, which measures the cost of label transitions between adjacent pixels.

2. The method according to claim 1, characterized in that, The use of the Laplace pyramid for multi-band fusion specifically includes: Laplacian pyramids LA and LB are constructed for the two images A and B to be fused, respectively, and a fusion weight mask is applied. Building the Gauss Pyramid ; At each level of the pyramid Above, the merged Laplace layer Calculate using the following formula: The merged layers of the Laplace Pyramid By restoring layer by layer, a final seamlessly blended image is obtained.

3. The method according to claim 1, characterized in that, Step 1 includes: Step S1: Multiple drones collaboratively collect images. Each image frame carries at least GPS information, attitude information, camera intrinsic parameters, and a camera ID used to distinguish the source of the data collection. Step S2: Transform the GPS coordinates of all images to the local East-North-Sky ENU coordinate system with one of the images as the origin, thereby unifying the coordinate systems of the multi-source images; Step S3: Group the images according to camera ID, construct an independent parallel processing pipeline for each group of images, perform aerial triangulation calculation in each pipeline, and use the obtained 3D point cloud to generate a digital elevation model (DEM) of the survey area. Step S4: Based on the digital elevation model, perform frame-by-frame orthorectification of the original images in each pipeline to generate a single-frame orthorectified image.

4. The method according to claim 3, characterized in that, In step S2, the coordinate unification process specifically includes: converting the WGS84 geographic coordinates of all images to the geocentric-fixed coordinate system ECEF as a transition; using one frame of the image as a reference point to calculate its ECEF coordinates; and converting the ECEF coordinates of each image point to ENU coordinates through a rotation matrix, thereby achieving coordinate unification of multi-source data.

5. The method according to claim 3, characterized in that, In step S3, the aerial triangulation calculation performed in each parallel processing pipeline includes: extracting SIFT or ORB feature points from each group of images and performing feature matching between images; constructing an image connectivity graph, using bundle adjustment to optimize exterior orientation elements, and outputting the precise pose parameters of each frame of image; and generating a digital elevation model (DEM) of the survey area based on the matched 3D point cloud as a terrain constraint for orthorectification.

6. The method according to claim 3, characterized in that, In step S4, orthorectification of each original image frame specifically includes: Generate a regular grid on the orthographic projection plane according to the set DOM resolution; For any ground point coordinates The corresponding elevation value Z is read from the DEM, and the pixel coordinates of the ground point on the original image are calculated using the collinearity equation. : In the formula, ( The spatial coordinates of the camera center The rotation matrix components are composed of exterior orientation angle elements. The focal length of the camera; The color value of the non-integer pixel coordinates obtained by inverse calculation is calculated using bilinear interpolation and then filled into the DOM grid points.

7. The method according to claim 6, characterized in that, Before substituting into the collinearity equation, first determine the image point coordinates of the original image. Lens distortion correction is performed using the following formula: In the formula, r is the distance from the image point to the principal point. The radial distortion coefficient is... Tangential distortion coefficient.

8. The method according to claim 1, characterized in that, In step 2, tile segmentation uses Mercator projection to divide tile layers, with each tile measuring 512×512 pixels. The finest layer determines the maximum value of the Mercator projection zoom. The Mercator projection coordinates corresponding to the tiles are then... Convert to hash key value Key using the following formula: The offset M must satisfy the following condition. To ensure the uniqueness of the Key, M takes the following values: .

9. The method according to claim 1, characterized in that, In step 3, the pyramid downsampling method is as follows: for each input frame of image, after all the corresponding finest layer tiles have been generated, the resolution is immediately reduced layer by layer and Gaussian filtering is applied to generate other lower-level tiles; during the generation process, the required fine-layer tiles are reread from the disk, and only the hash map and the orthophoto data of the current frame image are maintained in memory.

10. The method according to claim 1, characterized in that, In step 4, when the new tile completely covers the effective area of ​​the old tile in space, a full coverage operation is performed, that is, the new tile is directly used to replace the old tile in the storage device without performing seam line search and multi-band fusion.

Images