A satellite remote sensing orthographic image automatic mosaic method considering view angle and ground object height

Abstract

The present application relates to the field of remote sensing image processing and photogrammetry, and particularly relates to a satellite remote sensing orthographic image automatic mosaic method considering view angle and ground object height. The method obtains rational polynomial model parameters, a digital elevation model, a digital surface model, a digital orthographic image effective range and a resolution corresponding to two satellite images, constructs a grid coordinate set and determines an overlapping grid index; performs elevation interpolation and coordinate calculation on the grid points, calculates a displacement vector of a ground object vertex relative to a ground bottom point, and generates a displacement vector diagram of each scene; calculates a vector difference and a module length of the displacement vectors of the two scenes based on the overlapping grid index, and forms a reverse displacement vector module length diagram; determines overlapping units and start and end intersection points, and takes the module length diagram as a geometric cost to perform path iteration, and generates a mosaic line. The present application makes the geometric displacement formed by the ground object height and the view angle difference participate in the generation of the mosaic line, and reduces the geometric connection discontinuity at the edge.

Description

Technical Field

[0001] This invention relates to the fields of remote sensing image processing and photogrammetry, and in particular to an automatic mosaicking method for satellite remote sensing orthophotos that takes into account both viewing angle and ground feature height. Background Technology

[0002] In the fields of remote sensing image processing and photogrammetry, mosaicking of satellite remote sensing orthophotos typically involves processing multiple digital orthophotos, the effective image range, overlapping areas, and mosaicking lines. Existing solutions mostly determine mosaicking lines based on grayscale differences, color differences, texture continuity, the central axis of overlapping areas, or manual editing. These solutions can complete image edge stitching and mosaicking in typical scenarios where there are minimal variations in ground feature height and insignificant differences in imaging perspective. However, in densely built-up urban areas, forest areas, mountainous areas, or areas with dense structures, limitations arise such as inconsistent positions of identical ground features in two orthophotos, mosaicking lines crossing areas where tall ground features lean, and geometric discontinuities at the edge stitching points.

[0003] In current orthophoto production processes, digital orthophotos are typically orthorectified based on digital elevation models (DEMs). However, DEMs primarily represent ground elevation and struggle to reflect the surface elevation of vertices of features such as buildings, tree canopies, bridges, and structures. When different satellite imagery have different imaging perspectives, the direction and distance of the planar offset of the same feature vertex relative to the ground floor may differ. If mosaic lines are still primarily determined based on image grayscale, color, or the geometry of overlapping areas, the lack of a stable connection between the preceding orthorectification stage and the subsequent mosaic line generation stage, addressing differences in feature height and perspective, can easily lead to local geometric misalignments within overlapping areas that are difficult to identify during the mosaic line generation process.

[0004] For the joint processing of RPC model parameters, DEM, and DSM in the overlapping area of ​​two satellite remote sensing orthophotos, existing technologies still suffer from common shortcomings, such as difficulty in quantifying vertex tilt misalignment, difficulty in incorporating geometric misalignment information into the generation of mosaic lines for overlapping units, and insufficient connection between edge tilt misalignment limits and path iteration processes. Therefore, it is necessary to address the technical challenges of generating mosaic lines based on RPC model parameters, DEM, and DSM in the overlapping area of ​​two satellite remote sensing orthophotos. Summary of the Invention

[0005] To address the aforementioned technical problems, this invention provides an automatic mosaicking method for satellite remote sensing orthophotos that takes into account both viewing angle and ground feature height, comprising:

[0006] S100: Obtain the RPC model parameters, DEM, DSM, effective range of digital orthophoto DOM, and DOM resolution corresponding to the two satellite images; construct a set of DOM grid coordinates for each image based on the effective range of DOM and the DOM resolution, determine the overlapping grid index, and generate DOM grid mapping data.

[0007] S200. Based on the DOM grid mapping data, perform DEM elevation interpolation, DSM elevation interpolation, and RPC coordinate calculation on each scene's DOM grid points; calculate the offset vector of the feature vertex relative to the ground bottom point, and generate the DOM offset vector map for each scene.

[0008] S300. Based on the overlapping grid index and the DOM offset vector map of each scene, extract the offset vector of the corresponding grid point; perform vector difference calculation and magnitude calculation on the offset vectors of the two scenes to generate the reverse misalignment vector magnitude map.

[0009] S400. Based on the effective range of the two DOMs, perform topological intersection and connected component partitioning to determine overlapping units and start and end intersection points; associate the backward misalignment vector magnitude graph with the overlapping units to generate overlapping units with additional geometric costs.

[0010] S500: Based on the overlapping units with the additional geometric cost, the start and end intersections and the edge misalignment limit, perform path iteration according to the magnitude of the misalignment vector of the neighborhood grid to generate tessellation lines.

[0011] Furthermore, the generation of the DOM grid mapping data includes:

[0012] Based on the boundary coordinates of the effective range of each scene's DOM and the DOM resolution, generate a set of DOM grid coordinates for each scene;

[0013] Transform the coordinate sets of each scene's DOM grid to the same plane coordinate system;

[0014] Spatially overlay the effective ranges of the two DOM scenes, extract the grid point numbers located within the overlapping area of ​​the two DOM scenes, and generate an overlapping grid index.

[0015] Furthermore, the DEM elevation interpolation, DSM elevation interpolation, and RPC coordinate calculation include:

[0016] For each DOM grid point, the ground elevation is interpolated from the DEM, and the plane coordinates of the ground bottom point are calculated based on the ground elevation and RPC model parameters.

[0017] For each scene's DOM grid points, interpolate the feature vertex elevations from the DSM, and calculate the feature vertex planar coordinates based on the feature vertex elevations and RPC model parameters;

[0018] The difference between the planar coordinates of the vertex of the feature and the planar coordinates of the ground bottom point is used as the offset vector, and the offset vector is organized according to the DOM grid coordinate set to generate a DOM offset vector map.

[0019] Furthermore, the generation of the DOM offset vector map also includes:

[0020] The elevation difference of each scene's DOM grid points is marked based on the difference between the DSM elevation and the DEM elevation.

[0021] Set the offset vector of the DOM grid points marked with height difference as low feature grid to zero vector;

[0022] Mark DOM grid points with null DEM or DSM elevation values ​​as invalid grid points and remove them from the calculation range of the backward misalignment vector magnitude map.

[0023] Furthermore, the generation of the backward misalignment vector magnitude map includes:

[0024] Based on the overlapping grid index, the first offset vector is extracted from the first DOM offset vector map, and the second offset vector is extracted from the second DOM offset vector map;

[0025] The vector difference between the first offset vector and the second offset vector is used as the reverse misalignment vector;

[0026] Calculate the magnitude of the backward misalignment vector, and organize the magnitudes of the backward misalignment vector according to the overlapping grid positions to generate a backward misalignment vector magnitude map.

[0027] Furthermore, the determination of the overlapping units and their start and end points includes:

[0028] Perform topological intersection on the effective ranges of the two DOMs to obtain the boundary of the overlapping region;

[0029] The boundaries of the overlapping regions are divided into connected components to obtain multiple overlapping units;

[0030] Extract the intersection points of each overlapping unit boundary and the effective range boundaries of the two DOMs. When there are two intersection points corresponding to the overlapping unit boundary, take one of the intersection points as the starting intersection point and the other intersection point as the ending intersection point.

[0031] When there are more than two intersection points corresponding to the boundary of an overlapping unit, the intersection points are grouped into pairs according to the order in which they appear along the boundary of the effective range of the DOM, and each group of intersection points is used as the start and end intersection points of an overlapping unit.

[0032] Further, associating the magnitude map of the inverted misalignment vector with the overlapping unit includes:

[0033] The magnitude map of the inverted misalignment vector is cropped according to the boundary range of the overlapping unit to obtain the magnitude sub-map of the overlapping unit;

[0034] The magnitude of the backward misalignment vector of each grid point in the overlapping unit magnitude subgraph is used as the geometric cost of the neighboring grid within the overlapping unit.

[0035] The overlapping units, the start and end points, and the geometric cost are associated to generate overlapping units with additional geometric costs.

[0036] Further, the step of generating mosaic lines by iterating the path according to the magnitude of the backward misalignment vector of the neighborhood grid includes:

[0037] The pathfinding trend direction is determined based on the start and end points, and candidate neighborhood grids are determined based on the current pathfinding node and the pathfinding trend direction.

[0038] When there is a grid in the candidate neighborhood grid where the magnitude of the backward misalignment vector is not greater than the edge backward misalignment limit, the grid with the smallest vertical distance to the start and end intersection points is selected as the next pathfinding node.

[0039] When the magnitude of the backward misalignment vector of the candidate neighborhood grid is greater than the edge backward misalignment limit, the grid with the smallest magnitude of the backward misalignment vector is selected as the next pathfinding node.

[0040] Write the next pathfinding node into the path node sequence, and update the current pathfinding node based on the next pathfinding node, until the path node sequence reaches the termination point, generating a mosaic line.

[0041] Furthermore, the path iteration also includes:

[0042] When the path node sequence reaches the boundary of the overlapping unit and the boundary point is not the termination point, the boundary point is taken as the new starting point and the pathfinding trend direction is redefined.

[0043] When a candidate neighborhood grid is an invalid grid point, the candidate neighborhood grid is removed from the selection range of the next pathfinding node;

[0044] When duplicate nodes exist in the path node sequence, the candidate neighborhood grid corresponding to the duplicate node is deleted based on the visited node record.

[0045] Furthermore, after generating the mosaic lines, the process also includes:

[0046] Merge the mosaic lines corresponding to each overlapping unit to generate a DOM mosaic line set;

[0047] The set of DOM mosaic lines is used as the pixel selection boundary between the two DOM scenes;

[0048] Based on the pixel selection boundary, pixel selection, color balancing, and feathering are performed on the two DOM scenes to generate DOM mosaicking results.

[0049] The key innovations of this invention include:

[0050] (1) Couple the DOM grid coordinate set with the RPC model parameters, DEM and DSM, and form the coordinate calculation relationship between the ground bottom point and the ground feature vertex on the same DOM grid point, and organize the DOM offset vector map of each scene by the planar difference between the ground feature vertex and the ground bottom point.

[0051] (2) Establish the correspondence between the offset vector maps of two scenes using the overlapping grid index, calculate the vector difference and magnitude of the offset vectors of the two scenes at the same overlapping grid point, and form a offset vector magnitude map that represents the degree of offset of the vertices of the two scenes.

[0052] (3) Associate the backward misalignment vector magnitude graph with the overlapping units obtained by topological intersection and connected domain partitioning, and combine the start and end intersection points and the backward misalignment limit of the connecting edge to perform path iteration on the neighborhood grid to form a mosaic line constrained by geometric cost.

[0053] The following are its main beneficial effects:

[0054] (1) To address the problem that the elevation of the top of the ground feature in digital orthophotos is not fully expressed by ordinary ground elevation, the parameters of the DEM, DSM and RPC models are introduced into the same DOM grid coordinate relationship to form the DOM offset vector map of each scene, so that the planar offset of the top of the ground feature relative to the ground bottom point can be expressed in a gridded manner.

[0055] (2) To address the problem that it is difficult to quantify the inversion misalignment of the same ground features in the overlapping area of ​​two orthophoto images, the overlapping grid index is used to extract the corresponding offset vectors of the two images and calculate the vector difference to form an inversion misalignment vector magnitude map, so that the subsequent mosaic line generation can obtain geometric misalignment data that can be directly called.

[0056] (3) To address the problem of insufficient connection between geometric misalignment information and path selection in the existing mosaic line generation process, the magnitude of the backward misalignment vector is used as the geometric cost of the neighborhood grid within the overlapping unit, and path iteration is performed in combination with the edge backward misalignment limit to reduce the edge geometric discontinuity caused by the mosaic line crossing the high backward misalignment region. Attached Figure Description

[0057] Figure 1 This is a flowchart illustrating an automatic mosaicking method for satellite remote sensing orthophotos that takes into account both viewing angle and ground feature height, provided as an embodiment of this application. Detailed Implementation

[0058] Example 1: Refer to Figure 1 This is a flowchart illustrating an automatic mosaicking method for satellite remote sensing orthophotos that takes into account viewing angle and ground feature height, provided by an embodiment of the present invention. The process may include at least steps S100-S500:

[0059] S100: Obtain the RPC model parameters, DEM, DSM, effective range of digital orthophoto DOM, and DOM resolution corresponding to the two satellite images; construct a set of DOM grid coordinates for each image based on the effective range of DOM and the DOM resolution, determine the overlapping grid index, and generate DOM grid mapping data.

[0060] S200. Based on the DOM grid mapping data, perform DEM elevation interpolation, DSM elevation interpolation, and RPC coordinate calculation on each scene's DOM grid points; calculate the offset vector of the feature vertex relative to the ground bottom point, and generate the DOM offset vector map for each scene.

[0061] S300. Based on the overlapping grid index and the DOM offset vector map of each scene, extract the offset vector of the corresponding grid point; perform vector difference calculation and magnitude calculation on the offset vectors of the two scenes to generate the reverse misalignment vector magnitude map.

[0062] S400. Based on the effective range of the two DOMs, perform topological intersection and connected component partitioning to determine overlapping units and start and end intersection points; associate the backward misalignment vector magnitude graph with the overlapping units to generate overlapping units with additional geometric costs.

[0063] S500: Based on the overlapping units with the additional geometric cost, the start and end intersections and the edge misalignment limit, perform path iteration according to the magnitude of the misalignment vector of the neighborhood grid to generate tessellation lines.

[0064] S100: Obtain the RPC model parameters, DEM, DSM, effective range of the digital orthophoto DOM, and DOM resolution corresponding to the two satellite images; construct the DOM grid coordinate set for each image based on the effective range and resolution of the DOM, determine the overlapping grid index, and generate DOM grid mapping data.

[0065] In this embodiment, the automatic mosaicking method for satellite remote sensing orthophotos is executed by a remote sensing image processing platform, which includes an image data reading unit, a grid construction unit, a geometric calculation unit, an overlapping area processing unit, and a mosaicking line pathfinding unit. After receiving the processing task of two satellite images to be mosaicked, the image data reading unit reads the corresponding RPC model parameters, DEM, DSM, effective range of the digital orthophoto map (DOM), and DOM resolution for the two satellite images. Here, RPC stands for Rational Polynomial Coefficients; DEM stands for Digital Elevation Model; DSM stands for Digital Surface Model; and DOM stands for Digital Orthophoto Map. The RPC model parameters typically come from the satellite image metadata file, the DEM and DSM come from elevation data files under the same coordinate reference, the effective range of the DOM comes from the effective pixel boundary of the orthophoto or image range vector, and the DOM resolution comes from the ground resolution configuration of the DOM product to be generated.

[0066] When reading input objects, the image data reading unit first checks whether the effective extents of the two DOMs, DEM, and DSM are located in the same plane coordinate system. If the coordinate references of the effective extents of the two DOMs are different, the coordinate transformation parameters are called to transform the effective extents of the two DOMs to the same plane coordinate system; if the coordinate references of the DEM or DSM are inconsistent with the effective extents of the DOM, elevation data coordinate transformation is performed first, and then the grid construction process begins. If the RPC model parameters lack image row and column normalization parameters or ground coordinate normalization parameters, the image data reading unit marks the corresponding satellite image as an image unusable for geometric solution and stops the automatic generation process of mosaic lines for that image pair. This process avoids grid points in the subsequent S200 that cannot complete RPC coordinate solution.

[0067] After obtaining the effective ranges of two DOM scenes in the same planar coordinate system, the grid construction unit establishes regular grids for the first and second DOM effective ranges according to their respective DOM resolutions. Specifically, the grid construction unit uses the outer rectangle of the DOM effective range as the grid generation boundary and the DOM resolution as the grid spacing to generate grid points along the horizontal and vertical axes of the planar coordinate system. Then, it clips the outer rectangular grid based on the DOM effective range boundary, retaining grid points located within or on the boundary of the DOM effective range. For grid points near the boundary that fall outside the effective range but are adjacent to effective pixels, the grid construction unit does not consider them as grid points eligible for geometric calculations to prevent the subsequent offset vector map from forming pseudo-grid values ​​in invalid image areas.

[0068] In one implementation, the grid spacing is consistent with the DOM resolution; in another implementation, the grid spacing is an integer multiple of the DOM resolution, and the grid density required for mosaicking pathfinding is restored through neighborhood interpolation when generating the backward misalignment vector modulus map. In the latter implementation, the grid building unit writes the ratio between the grid spacing and the DOM resolution into the DOM grid mapping data, ensuring that subsequent steps S300 and S500 use the same grid scale for overlapping grid localization and path iteration.

[0069] After establishing the coordinate sets of the two DOM grids, the grid construction unit spatially overlays the effective ranges of the two DOMs, extracting the overlapping portion to form the overlapping region of the two DOMs. For grid points located within this overlapping region, the grid construction unit records their grid point numbers in the first DOM grid coordinate set, their grid point numbers in the second DOM grid coordinate set, and their corresponding planar coordinates. If the grid points of the first and second DOMs cannot coincide one-to-one in planar position, the target grid in a unified planar coordinate system is used as an intermediate index, and the proximity mapping relationship between the target grid point and the grid points of the two DOMs is recorded. The resulting overlapping grid index is not simply a range record, but a spatial correspondence for subsequent point-by-point comparisons of the offset vectors of the two scenes.

[0070] When there is no effective overlap between the overlapping areas of two DOM scenes, the grid construction unit does not generate an overlapping grid index and marks the processing task as completed with no overlapping area. When the overlapping area of ​​two DOM scenes consists of multiple disconnected regions, the grid construction unit still generates an overlapping grid index, but retains the connected component numbers of the overlapping areas in the index for S400 to divide the overlapping units. The final DOM grid mapping data generated by S100 includes the set of DOM grid coordinates for each scene, the range of the overlapping area between the two DOM scenes, the overlapping grid index, the grid point numbers, and the correspondence between the planar coordinates. When S200 performs DEM elevation interpolation, DSM elevation interpolation, and RPC coordinate calculation, it will directly read the grid point numbers and planar coordinates from the DOM grid mapping data without rebuilding the spatial index.

[0071] S200. Based on the DOM grid mapping data, perform DEM elevation interpolation, DSM elevation interpolation, and RPC coordinate calculation on each DOM grid point; calculate the offset vector of the feature vertex relative to the ground bottom point, and generate the DOM offset vector map for each scene:

[0072] S200 is executed by the geometry solving unit. After receiving the DOM grid mapping data output by S100, the geometry solving unit reads the planar coordinates of grid points in the first and second scene DOM grid coordinate sets, respectively, according to the scene. Since this method needs to express the planar offset at the same DOM grid point caused by the height of the ground object and the imaging viewpoint, the geometry solving unit does not directly use the DOM pixel grayscale or texture value as the mosaic basis. Instead, it first splits each grid point into two geometric objects: a ground bottom point and a ground object vertex. The ground bottom point is the ground corresponding point determined based on the DEM elevation, and the ground object vertex is the top of the surface object or the corresponding point on the surface determined based on the DSM elevation.

[0073] For each DOM grid point, the geometric solution unit first performs elevation interpolation in the DEM. During interpolation, if the DOM grid point is located between DEM grid nodes, the ground elevation is calculated based on the surrounding DEM elevation values; if there are missing values ​​in the surrounding DEM elevation, the elevation values ​​within the effective neighborhood are used to fill in the gaps; if the effective neighborhood is insufficient to complete the elevation interpolation, the DOM grid point is marked as an invalid grid point. After completing the DEM elevation interpolation, the geometric solution unit combines the planar coordinates of the DOM grid point with the ground elevation to obtain the ground bottom point coordinates, and performs coordinate calculation based on the RPC model parameters of the current foreground satellite image to obtain the corresponding position of the ground bottom point in the current foreground DOM plane.

[0074] After completing the ground baseline calculation, the geometric calculation unit interpolates the feature vertex elevations of the same DOM grid point from the DSM. The DSM elevation reflects the top or surface elevation of surface objects; therefore, in areas with buildings, tree canopies, bridges, and structures, there is usually a difference between the DSM elevation and the DEM elevation. The geometric calculation unit combines the planar coordinates of the DOM grid point with the DSM elevation to obtain the feature vertex coordinates, and uses the RPC model parameters of the same satellite imagery to calculate the coordinates, obtaining the corresponding position of the feature vertex in the current foreground DOM plane. If the DSM elevation is null but the DEM elevation is valid, the grid point is not included in the feature vertex offset calculation; if both the DSM and DEM elevations are valid but the difference between them is marked as a low-feature grid, the offset vector of the grid point is treated as a zero vector.

[0075] In one implementation, the geometric solution unit marks DOM grid points with height differences based on the difference between the DSM elevation and the DEM elevation. These height difference marks distinguish between low-feature grids and grids where feature height influences the selection. For low-feature grids, the offset between feature vertices and ground datum points in the orthographic plane has little impact on mosaic line selection; the geometric solution unit sets the offset vector of these grid points to zero. For grid points whose height difference marks do not belong to low-feature grids, the geometric solution unit retains the planar coordinates of feature vertices and ground datum points obtained through RPC coordinate calculation and calculates the difference between these two coordinates, using this difference as the offset vector. The offset vector includes a lateral offset component and a longitudinal offset component, both of which together characterize the tilt direction and distance of feature vertices relative to ground datum points in the DOM plane from the perspective of the current foreground image.

[0076] The DOM offset vector maps for each scene are organized according to the DOM grid coordinate set. Specifically, the geometric calculation unit uses the grid point number as an index to write the offset vector corresponding to each valid grid point into the current scene's DOM offset vector map; for invalid grid points, the invalid grid point mark is retained but no offset vector value is generated. If a grid point cannot obtain the planar coordinates of the feature vertex or the ground bottom point due to the non-convergence of the RPC model parameter calculation, the geometric calculation unit marks the grid point as a geometric calculation failure grid point and deletes it from the calculation range of the backward misalignment vector magnitude map in subsequent S300. After this processing, when S300 reads the DOM offset vector maps for each scene, it only needs to extract the offset vectors that are valid for both scenes based on the overlapping grid index, avoiding invalid grid points from participating in the vector difference calculation.

[0077] The DOM offset vector maps output by S200 maintain the same grid point numbering system as the DOM grid mapping data formed by S100. The first and second DOM offset vector maps not only record the offset vector values ​​but also retain the corresponding grid point numbers and valid status. After receiving this result, S300 will establish a one-to-one offset vector extraction relationship between the two DOM offset vector maps using the overlapping grid index as a channel.

[0078] S300. Based on the overlapping grid index and the DOM offset vector map of each scene, extract the offset vector of the corresponding grid point; perform vector difference calculation and magnitude calculation on the offset vectors of the two scenes to generate the magnitude map of the inverted misalignment vector:

[0079] S300 is performed by the inverted misalignment calculation unit. This unit does not recalculate the terrain or image geometry; instead, it reads the overlapping grid index formed by S100 and the DOM offset vector maps of each scene formed by S200, and compares the offset vectors of the two scenes point by point within the overlapping area of ​​the two DOMs. The comparison object here is not the grayscale value, color value, or texture feature of the two images, but the planar offset vector formed by the height of the ground features and the imaging viewpoint at the same overlapping grid point of the two images.

[0080] In practice, the backward misalignment calculation unit reads the overlapping grid point numbers item by item according to the overlapping grid index. For each overlapping grid point, the first offset vector is extracted from the first scene's DOM offset vector map, and the second offset vector is extracted from the second scene's DOM offset vector map. If the offset vector corresponding to any scene is a null value, an invalid grid point, or a grid point where the geometric solution failed, then that overlapping grid point does not participate in the backward misalignment vector magnitude calculation and remains in an invalid state in the backward misalignment vector magnitude map. If both scene offset vectors are valid, then the vector difference calculation begins.

[0081] The backward misalignment vector is formed by the vector difference between the first offset vector and the second offset vector. Because the satellite imaging perspectives of the first and second scenes are different, the backward direction and backward distance of the same feature vertex relative to the ground bottom point in the two scene DOMs may differ. Vector difference calculation unifies these direction and distance differences into a single plane coordinate system. The backward misalignment calculation unit calculates the magnitude of the vector difference to obtain the backward misalignment vector magnitude of the overlapping grid point. The backward misalignment vector magnitude is a non-negative number, its spatial location is determined by the overlapping grid index, and its value participates in the mosaic line iteration in S500 as the path cost of the neighboring grid.

[0082] In one implementation, the backward misalignment calculation unit performs zero-vector inheritance processing on overlapping grid points marked as low-feature grids by height difference. If the corresponding grid points of two scenes are both low-feature grids, then the offset vectors of both scenes are zero vectors, and the magnitude of the backward misalignment vector is zero; if one scene is a low-feature grid while the other scene has a valid offset vector, then the magnitude is calculated according to the vector difference between the offset vectors of the two scenes. This processing ensures that low-feature areas are not misclassified as high backward misalignment areas, while retaining the constraint on mosaic pathfinding when there are elevation differences within a single scene.

[0083] The backward misalignment vector magnitude map is organized according to the overlapping grid positions. Each valid overlapping grid point corresponds to a backward misalignment vector magnitude, and each invalid overlapping grid point corresponds to an invalid state. For regions with a few invalid grid points but continuous surrounding valid grid points, in one implementation, neighborhood completion can be performed based on the magnitudes of neighboring valid grid points; for regions with continuously distributed invalid grid points that cannot reflect the local geometric state, the backward misalignment calculation unit retains them as invalid regions that cannot participate in path selection. When performing overlapping unit association, S400 reads the distribution of valid and invalid grid points in the backward misalignment vector magnitude map to ensure that subsequent mosaicking pathfinding does not cross regions lacking geometric cost basis.

[0084] The backward misalignment vector magnitude map output by S300 lies in the same coordinate system as the overlapping area of ​​the two DOM scenes, and its grid position is consistent with the overlapping grid index of S100. This map is not the final mosaic line, but rather the geometric cost basis required by S400 and S500. S400 will determine the overlapping units based on the effective range of the two DOM scenes, and will crop or index the backward misalignment vector magnitude map into the range of each overlapping unit, so that each overlapping unit has independent local magnitude data.

[0085] S400. Based on the effective range of the two DOMs, perform topological intersection and connected component partitioning to determine overlapping units and their start and end points; associate the backward misalignment vector magnitude graph with the overlapping units to generate overlapping units with additional geometric costs.

[0086] S400 is executed by the overlapping region processing unit. The overlapping region processing unit receives the backward misalignment vector magnitude map generated by S300 and simultaneously reads the effective ranges of the two DOM scenes from S100. Since the effective ranges of the two DOM scenes may be irregular polygons in actual production, and the overlapping region may consist of multiple disconnected blocks, the overlapping region processing unit first performs topological intersection on the effective ranges of the two DOM scenes to obtain the overlapping region boundary, and then divides the overlapping region boundary into connected components to form one or more overlapping units. Each overlapping unit corresponds to a segment of effective space for subsequent independent pathfinding.

[0087] During the topological intersection process, the overlapping region processing unit uses the effective boundaries of two DOM scenes as input to calculate their overlapping regions. If the overlapping region is a single closed region, one overlapping unit is generated; if the overlapping region is divided into multiple closed regions by invalid extents, image holes, or discontinuous boundaries, multiple overlapping units are generated according to connectivity. For boundary self-intersections, repeated edges, or small fragmented regions, the overlapping region processing unit first performs boundary cleanup, deleting fragmented regions with areas smaller than the area corresponding to the grid resolution, and retaining regions that intersect with the overlapping grid index. The resulting overlapping units have a corresponding grid range to the backward misalignment vector magnitude map in S300.

[0088] After determining the overlapping units, the overlapping region processing unit extracts the intersection points of each overlapping unit boundary with the effective boundaries of the two DOMs. These intersection points define the start and end positions of the mosaic lines within the overlapping unit. When there are two intersection points corresponding to the overlapping unit boundary, one is designated as the start intersection point and the other as the end intersection point. When there are more than two intersection points corresponding to the overlapping unit boundary, the intersection points are paired and grouped according to their order of appearance along the DOM's effective boundaries, with each pair of intersection points serving as the start and end intersection points of an overlapping unit. If the intersection points form a line segment intersection due to boundary overlap, the overlapping region processing unit selects the endpoints of this line segment intersection as candidate intersection points and groups them according to their order of arrangement on the overlapping unit boundary.

[0089] After the overlapping units are determined, the overlapping region processing unit associates the backward misalignment vector magnitude map with each overlapping unit. Specifically, the overlapping region processing unit trims the backward misalignment vector magnitude map according to the boundary range of the overlapping unit to obtain the overlapping unit magnitude submap. For grid points that fall on the boundary of the overlapping unit, if the center coordinates of the grid point are located inside or on the boundary of the overlapping unit, it is included in the overlapping unit magnitude submap. If the effective magnitude grid points within the range of a certain overlapping unit are insufficient to form a continuous candidate path, the overlapping unit is marked as an overlapping unit to be reviewed and does not enter the automatic path iteration, or in another implementation, the candidate neighborhood is expanded before entering S500.

[0090] The overlapping region processing unit associates overlapping cells, start and end points, and geometric costs to generate overlapping cells with additional geometric costs. These geometric costs originate from the backward misalignment vector magnitude of each grid point in the overlapping cell magnitude subgraph. The overlapping cell with additional geometric costs includes at least the overlapping cell boundary, start and end points, the overlapping cell magnitude subgraph, the state of valid grid points, and the state of invalid grid points. Upon receiving this data, the S500 uses the overlapping cell boundary as the pathfinding space, the start and end points as path endpoints, and the backward misalignment vector magnitude as the cost of the candidate neighboring grids to perform path iteration.

[0091] S500: Based on the overlapping units with the additional geometric cost, the start and end intersections and the edge misalignment limit, perform path iteration according to the magnitude of the misalignment vector of the neighborhood grid to generate tessellation lines.

[0092] S500 is executed by the mosaicking line pathfinding unit. After receiving the overlapping cells with additional geometric costs from the S400 output, the mosaicking line pathfinding unit generates mosaicking lines for each overlapping cell. For each overlapping cell, the mosaicking line pathfinding unit reads the start and end points, the overlapping cell boundary, and the overlapping cell magnitude subgraph, and reads the edge offset misalignment limit. The edge offset misalignment limit is used to distinguish whether the magnitude of the offset misalignment vector of the neighboring grid is within the selectable range. This limit can be provided by DOM product edge requirements, image production task configuration, or manually input parameters.

[0093] At the start of path iteration, the mosaic line pathfinding unit designates the grid containing the starting intersection point or the nearest valid grid to the starting intersection point as the current pathfinding node, and the grid containing the ending intersection point or the nearest valid grid to the ending intersection point as the termination target. Subsequently, the mosaic line pathfinding unit determines the pathfinding trend direction based on the starting and ending intersection points and selects candidate neighboring grids near the current pathfinding node. The range of candidate neighboring grids does not extend beyond the overlapping cell boundary and does not include invalid grid points marked in S300 or S400. For candidate neighboring grids located on the boundary, the mosaic line pathfinding unit decides whether to retain them based on whether the grid's center point is located within the overlapping cell or on the boundary.

[0094] In the candidate neighborhood grids of the current pathfinding node, the tessellation line pathfinding unit first reads the magnitude of the backward misalignment vector of each candidate neighborhood grid. If there is a candidate neighborhood grid whose backward misalignment vector magnitude is not greater than the edge backward misalignment limit, then the grid with the smallest vertical distance to the line connecting the start and end points is selected as the next pathfinding node. Here, "closer to the line connecting the start and end points" means that the vertical distance from the center point of the candidate neighborhood grid to the line connecting the start and end points is small; when multiple candidate neighborhood grids have the same vertical distance, the tessellation line pathfinding unit selects the grid with the larger projected distance along the pathfinding trend direction as the next pathfinding node. After this processing, the path is constrained by the magnitude of the backward misalignment vector and also maintains the trend of advancing from the start point to the end point.

[0095] If the magnitude of the backward misalignment vector of all candidate neighboring grids is greater than the edge backward misalignment limit, the mosaic line pathfinding unit selects the grid with the smallest backward misalignment vector magnitude as the next pathfinding node. This branch is applicable to cases where the overall backward misalignment is large in a local area of ​​the overlapping unit; in this case, the path cannot completely avoid the grid exceeding the limit, but it still proceeds according to the local minimum magnitude. The mosaic line pathfinding unit writes the next pathfinding node into the path node sequence and updates the current pathfinding node based on the next pathfinding node, and continues to read new candidate neighboring grids. After the path node sequence reaches the termination point or a valid grid adjacent to the termination point, the mosaic line pathfinding unit converts the path node sequence into a linear object and generates the mosaic line corresponding to the overlapping unit.

[0096] During path iteration, if the path node sequence reaches the boundary of an overlapping unit and this boundary point is not a termination point, the mosaicking line pathfinding unit takes this boundary point as the new starting point and redetermines the pathfinding trend direction between this boundary point and the termination point. At this time, the formed path node sequence is retained, and the new path segment continues to iterate within the same overlapping unit. If all candidate neighborhood grids are invalid grid points, the mosaicking line pathfinding unit deletes the candidate neighborhood set and, in one implementation, expands the candidate neighborhood range; if there are still no valid grid points after expansion, the overlapping unit is marked as an overlapping unit that cannot automatically generate mosaic lines. If duplicate nodes appear in the path node sequence, the mosaicking line pathfinding unit deletes the candidate neighborhood grids corresponding to the duplicate nodes based on the visited node records to prevent the path from looping in local areas.

[0097] After the mosaic lines of each overlapping unit are generated, the mosaic line pathfinding unit merges multiple mosaic lines into a DOM mosaic line set. For cases where the endpoint distance at the boundary of adjacent overlapping units is less than one DOM resolution, the mosaic line pathfinding unit snaps the endpoints together to ensure the DOM mosaic line set remains continuous at the boundary. Subsequently, the DOM mosaic line set serves as the pixel selection boundary between the two DOM scenes and enters the mosaic result generation process. Specifically, the mosaic result generation unit divides the image capture range of the two DOM scenes within the overlapping area according to the DOM mosaic line set, selects pixels on both sides of the mosaic line, and in one embodiment, performs color balancing and feathering on the neighborhood of the mosaic line to generate the DOM mosaic result. This DOM mosaic result is saved together with the DOM mosaic line set. Subsequent image production tasks can directly read the mosaic line node sequence, overlapping unit number, and corresponding backward misalignment vector magnitude during result inspection or local recalculation.

Claims

1. A method for automatically mosaicking satellite remote sensing orthophotos that takes into account viewing angle and ground feature height, characterized in that, include: S100: Obtain the RPC model parameters, DEM, DSM, effective range of digital orthophoto DOM, and DOM resolution corresponding to the two satellite images; Based on the effective range of the DOM and the resolution of the DOM, construct a set of DOM grid coordinates for each scene, determine the overlapping grid index, and generate DOM grid mapping data; S200. Based on the DOM grid mapping data, perform DEM elevation interpolation, DSM elevation interpolation, and RPC coordinate calculation for each DOM grid point. Calculate the offset vector of the feature vertices relative to the ground bottom point, and generate DOM offset vector maps for each scene; S300. Based on the overlapping grid index and the DOM offset vector map of each scene, extract the offset vector of the corresponding grid point; perform vector difference calculation and magnitude calculation on the offset vectors of the two scenes to generate the reverse misalignment vector magnitude map. S400. Based on the effective range of the two DOMs, perform topological intersection and connected component partitioning to determine overlapping units and start and end intersection points; associate the backward misalignment vector magnitude graph with the overlapping units to generate overlapping units with additional geometric costs. S500: Based on the overlapping units with the additional geometric cost, the start and end intersections and the edge misalignment limit, perform path iteration according to the magnitude of the misalignment vector of the neighborhood grid to generate tessellation lines.

2. The method according to claim 1, characterized in that, The generation of the DOM grid mapping data includes: Based on the boundary coordinates of the effective range of each scene's DOM and the DOM resolution, generate a set of DOM grid coordinates for each scene; Transform the coordinate sets of each scene's DOM grid to the same plane coordinate system; Spatially overlay the effective ranges of the two DOM scenes, extract the grid point numbers located within the overlapping area of ​​the two DOM scenes, and generate an overlapping grid index.

3. The method according to claim 1, characterized in that, The DEM elevation interpolation, DSM elevation interpolation, and RPC coordinate calculation include: For each DOM grid point, the ground elevation is interpolated from the DEM, and the plane coordinates of the ground bottom point are calculated based on the ground elevation and RPC model parameters. For each scene's DOM grid points, interpolate the feature vertex elevations from the DSM, and calculate the feature vertex planar coordinates based on the feature vertex elevations and RPC model parameters; The difference between the planar coordinates of the vertex of the feature and the planar coordinates of the ground bottom point is used as the offset vector, and the offset vector is organized according to the DOM grid coordinate set to generate a DOM offset vector map.

4. The method according to claim 3, characterized in that, The generation of the DOM offset vector map also includes: The elevation difference of each scene's DOM grid points is marked based on the difference between the DSM elevation and the DEM elevation. Set the offset vector of the DOM grid points marked with height difference as low feature grid to zero vector; Mark DOM grid points with null DEM or DSM elevation values ​​as invalid grid points and remove them from the calculation range of the backward misalignment vector magnitude map.

5. The method according to claim 1, characterized in that, The generation of the backward misalignment vector magnitude map includes: Based on the overlapping grid index, the first offset vector is extracted from the first DOM offset vector map, and the second offset vector is extracted from the second DOM offset vector map; The vector difference between the first offset vector and the second offset vector is used as the reverse misalignment vector; Calculate the magnitude of the backward misalignment vector, and organize the magnitudes of the backward misalignment vector according to the overlapping grid positions to generate a backward misalignment vector magnitude map.

6. The method according to claim 1, characterized in that, The determination of the overlapping units and the start and end intersection points includes: Perform topological intersection on the effective ranges of the two DOMs to obtain the boundary of the overlapping region; The boundaries of the overlapping regions are divided into connected components to obtain multiple overlapping units; Extract the intersection points of each overlapping unit boundary and the effective range boundaries of the two DOMs. When there are two intersection points corresponding to the overlapping unit boundary, take one of the intersection points as the starting intersection point and the other intersection point as the ending intersection point. When there are more than two intersection points corresponding to the boundary of an overlapping unit, the intersection points are grouped into pairs according to the order in which they appear along the boundary of the effective range of the DOM, and each group of intersection points is used as the start and end intersection points of an overlapping unit.

7. The method according to claim 1, characterized in that, Associating the inverted misalignment vector magnitude map with the overlapping unit includes: The magnitude map of the inverted misalignment vector is cropped according to the boundary range of the overlapping unit to obtain the magnitude sub-map of the overlapping unit; The magnitude of the backward misalignment vector of each grid point in the overlapping unit magnitude subgraph is used as the geometric cost of the neighboring grid within the overlapping unit. The overlapping units, the start and end points, and the geometric cost are associated to generate overlapping units with additional geometric costs.

8. The method according to claim 1, characterized in that, The step of iterating the path according to the magnitude of the backward misalignment vector of the neighborhood grid to generate mosaic lines includes: The pathfinding trend direction is determined based on the start and end points, and candidate neighborhood grids are determined based on the current pathfinding node and the pathfinding trend direction. When there is a grid in the candidate neighborhood grid where the magnitude of the backward misalignment vector is not greater than the edge backward misalignment limit, the grid with the smallest vertical distance to the start and end intersection points is selected as the next pathfinding node. When the magnitude of the backward misalignment vector of the candidate neighborhood grid is greater than the edge backward misalignment limit, the grid with the smallest magnitude of the backward misalignment vector is selected as the next pathfinding node. Write the next pathfinding node into the path node sequence, and update the current pathfinding node based on the next pathfinding node, until the path node sequence reaches the termination point, generating a mosaic line.

9. The method according to claim 8, characterized in that, The path iteration also includes: When the path node sequence reaches the boundary of the overlapping unit and the boundary point is not the termination point, the boundary point is taken as the new starting point and the pathfinding trend direction is redefined. When a candidate neighborhood grid is an invalid grid point, the candidate neighborhood grid is removed from the selection range of the next pathfinding node; When duplicate nodes exist in the path node sequence, the candidate neighborhood grid corresponding to the duplicate node is deleted based on the visited node record.

10. The method according to claim 1, characterized in that, After generating the mosaic lines, the process further includes: Merge the mosaic lines corresponding to each overlapping unit to generate a DOM mosaic line set; The set of DOM mosaic lines is used as the pixel selection boundary between the two DOM scenes; Based on the pixel selection boundary, pixel selection, color balancing, and feathering are performed on the two DOM scenes to generate DOM mosaicking results.

Images