Method for height correction of unmanned aerial vehicle oblique photogrammetry under complex terrain

Abstract

The application discloses a complex terrain under unmanned aerial vehicle oblique photogrammetry height correction method, belongs to the technical field of unmanned aerial vehicle photogrammetry. The method comprises the following steps: acquiring the multi-view oblique image of the unmanned aerial vehicle and the original height observation value; generating an initial three-dimensional point cloud and identifying the terrain mutation area and the non-mutation area according to the spatial distribution density of the point position; extracting the feature points on both sides of the boundary in the mutation area, calculating the parallax gradient change quantity, and inversely compensating the local error of the original height observation value according to the parallax gradient change quantity to generate a corrected height datum plane; and fitting and adjusting the original height observation value in the non-mutation area by using the height datum plane to complete the height correction of the whole area. The method can improve the problem that the height measurement precision is disturbed by the terrain mutation under the complex terrain condition, and realizes the differential accurate correction of the mutation area and the non-mutation area.

Description

Technical Field

[0001] This invention relates to the field of unmanned aerial vehicle (UAV) photogrammetry technology, specifically to a method for elevation correction in UAV oblique photogrammetry under complex terrain. Background Technology

[0002] When UAV oblique photogrammetry is applied in complex terrain areas, the acquired elevation data often contains significant errors due to factors such as severe surface undulations, local occlusion, and abrupt texture changes. In complex terrain, the parallax changes on both sides of abrupt boundaries such as steep slopes, gullies, and cliffs are extremely discontinuous. Conventional elevation difference calculation models struggle to recover the true surface elevation near these abrupt changes, leading to abrupt local elevation jumps. Existing elevation correction methods mostly rely on global filtering of the original point cloud or overall adjustment of elevation observations across the entire region based on a single adjustment model. These methods do not distinguish between abrupt and non-abrupt terrain areas, making it difficult to specifically address systematic distortions at abrupt changes. Global smoothing algorithms tend to smooth out the true abrupt features of terrain changes, while the correction effect in non-abrupt areas is affected by the drag from large errors at abrupt changes, resulting in a shift in the overall elevation datum. The cost of setting up field control points is high, and the coverage density is limited, making it impossible to apply dense constraints to all terrain abrupt change locations, leaving residual distortions difficult to eliminate. To address this issue, a method is needed that can automatically detect the location of terrain abrupt changes and apply targeted elevation compensation based on the differences in image features at the abrupt changes. At the same time, it is necessary to ensure that the elevation values ​​in non-abrupt areas smoothly transition under the guidance of the correction benchmark at the abrupt changes, thereby avoiding the propagation of global errors. Summary of the Invention

[0003] This invention provides a method for elevation correction in UAV oblique photogrammetry under complex terrain. The purpose is to solve the problems of insufficient suppression of distortion in areas with abrupt terrain changes and the influence of abrupt errors on the correction of non-abrupt terrain areas by abrupt errors in existing elevation correction methods. The method aims to achieve differentiated compensation for elevation errors based on terrain features, thereby improving the accuracy of elevation data and the authenticity of terrain representation across the entire region.

[0004] To achieve the above objectives, the present invention provides the following technical solution: The present invention provides a method for elevation correction in UAV oblique photogrammetry under complex terrain, the method comprising: The system acquires multi-view oblique image data collected by UAVs in complex terrain areas, and simultaneously acquires the original elevation observation values ​​corresponding to each image recorded by the airborne positioning system; it extracts feature points and performs stereo matching on the multi-view oblique image data to generate initial 3D point cloud data; based on the spatial distribution density of points in the initial 3D point cloud data, it identifies areas of abrupt terrain changes and areas of non-abrupt changes, accurately captures the location of abrupt terrain changes by utilizing the abrupt change characteristics of point cloud density, distinguishes different landform types, and provides a reliable basis for subsequent differential correction.

[0005] Within areas of abrupt topographic changes, image feature points located on both sides of the abrupt boundary are extracted, and their disparity gradient changes are calculated. This disparity gradient change quantitatively reflects the intensity of image disparity distortion caused by drastic elevation changes at the abrupt boundary. Based on the disparity gradient change, local error inversion compensation is performed on the original elevation observations to generate a corrected elevation datum, effectively eliminating local elevation distortion caused by abrupt topographic changes and making the elevation datum smooth and closely conforming to the actual terrain. In non-abrupt topographic areas, the corrected elevation datum is used to fit and adjust the original elevation observations. By transferring the correction datum constraint from the abrupt topographic area to the non-abrupt topographic area, the elevation correction of the entire region under complex terrain using UAV oblique photogrammetry is completed, ensuring the consistency and accuracy of elevation data across the entire region.

[0006] As a preferred technical solution of the present invention, the specific method for identifying abrupt terrain changes and non-abrupt terrain areas is as follows: The initial 3D point cloud data is divided into multiple grid cells of equal size according to horizontal projection coordinates; the number of points contained in each grid cell is counted, and the spatial distribution density value of the points in each grid cell is calculated; grid cells with spatial distribution density values ​​exceeding a preset density threshold are marked as candidate abrupt change units; four-connectivity clustering is performed on all marked candidate abrupt change units, and interconnected candidate abrupt change units are merged to form abrupt terrain changes, while the areas containing the remaining grid cells are marked as non-abrupt terrain areas. This method can automatically and robustly segment steep slopes, cliffs, and other abrupt terrain changes, as well as gentler areas, in complex terrain.

[0007] Preferably, the calculation process of the disparity gradient change is as follows: Multiple boundary control points are selected along the boundary line of the terrain abrupt change region. Each boundary control point extends a predetermined distance in a direction perpendicular to the boundary line towards the inside and outside of the abrupt change region, forming an inner sampling zone and an outer sampling zone. Corresponding image feature points are extracted within the inner and outer sampling zones. These corresponding image feature points exhibit the same texture features in tilted images from different viewpoints. The disparity value of each corresponding image feature point in the tilted images from different viewpoints is calculated. The average disparity value of all corresponding image feature points within the inner sampling zone is taken as the inner average disparity, and the average disparity value of all corresponding image feature points within the outer sampling zone is taken as the outer average disparity. The absolute value of the difference between the inner average disparity and the outer average disparity is taken as the disparity gradient change at that boundary control point. This disparity gradient change can fully reflect the degree of interference of elevation shift on the disparity of multi-view images caused by elevation changes on both sides of the abrupt change boundary.

[0008] Further, the step of generating the corrected elevation datum through local error inversion compensation includes: sorting the disparity gradient changes at each boundary control point in descending order, and selecting the top few boundary control points as strong distortion feature points; establishing a local inversion window centered on each strong distortion feature point, and collecting all original elevation observations within the local inversion window; calculating an inversion compensation coefficient for each original elevation observation within the local inversion window based on the disparity gradient change corresponding to the strong distortion feature point, wherein the inversion compensation coefficient is positively correlated with the disparity gradient change; multiplying each original elevation observation within the local inversion window by the corresponding inversion compensation coefficient to obtain the corrected local elevation value; and weighted fusion of the corrected local elevation values ​​from all local inversion windows to generate a corrected elevation datum covering areas with abrupt topographic changes.

[0009] As a technical solution of the present invention, in the local inversion compensation process, the radius of the local inversion window is dynamically and adaptively adjusted according to the boundary curvature of the terrain abrupt change region. The larger the boundary curvature, the smaller the radius of the local inversion window, so as to accurately track the details of terrain abrupt changes and avoid over-smoothing or under-compensation. Meanwhile, the specific method for calculating the inversion compensation coefficient for each original elevation observation value within the local inversion window is as follows: determine the coordinates of the center point of the local inversion window; calculate the planar distance between the coordinates of each original elevation observation value within the local inversion window and the coordinates of the center point; divide the planar distance corresponding to each coordinate by the radius of the local inversion window to obtain the normalized distance value of each coordinate; use the disparity gradient change of the strongly distorted feature point as the benchmark compensation intensity value; multiply the benchmark compensation intensity value by one minus the normalized distance value to obtain the inversion compensation coefficient at each coordinate, so that the closer the location is to the strongly distorted feature point, the stronger the error compensation.

[0010] The step of generating the corrected elevation datum using weighted fusion can be further designed as follows: For any target point located within the overlapping area of ​​multiple local inversion windows, collect multiple corrected local elevation values ​​corresponding to the target point within each local inversion window; obtain the spatial distance from the target point to the strongly distorted feature point corresponding to each local inversion window, and use the reciprocal of each spatial distance as the fusion weight value of that local inversion window for the target point; multiply each corrected local elevation value of the target point by its corresponding fusion weight value, sum the results, and then divide by the sum of all fusion weight values ​​to obtain the fused elevation value of the target point; traverse all target points within the terrain abrupt change area, and combine the fused elevation values ​​of all target points to form the corrected elevation datum. This fusion strategy based on distance reciprocal weighting ensures that the elevation transition in the overlapping area of ​​multiple inversion windows is natural and without abrupt changes.

[0011] In the fitting and adjustment process for non-abrupt terrain regions, a preferred approach is as follows: Extract a series of boundary elevation values ​​from the corrected elevation datum along the boundary line of the abrupt terrain region, while simultaneously selecting multiple original elevation observations within the non-abrupt region as control reference points; construct a radial basis function interpolation model using the boundary elevation values ​​as mandatory constraints and the control reference points as soft constraints; input each original elevation observation within the non-abrupt region into the radial basis function interpolation model, which outputs the elevation adjustment amount corresponding to each original elevation observation; add the corresponding elevation adjustment amount to each original elevation observation within the non-abrupt region to obtain the corrected elevation value for the non-abrupt region. This interpolation method, combining strong boundary constraints and internal soft constraints, can accurately transfer the corrected datum from abrupt terrain regions to non-abrupt regions without introducing unreasonable distortions in smooth terrain areas.

[0012] More preferably, the selection of control reference points adopts an adaptive point placement strategy. In non-abrupt regions, the closer the control reference points are to the boundary of abrupt terrain changes, the higher the density of control reference points. This better constrains elevation gradient changes near the boundary and improves the stability of the overall interpolation model. The specific steps for constructing the radial basis function interpolation model are as follows: The coordinates of the points corresponding to the boundary elevation values ​​are used as the central nodes of the radial basis function; at each central node, the boundary elevation value is set as a forced pass-through point of the interpolation model. The coordinates of the points corresponding to the control reference points are used as auxiliary nodes of the radial basis function; at each auxiliary node, the original elevation observation value of the control reference point is set as a soft constraint point of the interpolation model. A forced constraint weight value is assigned to each central node, and a soft constraint weight value is assigned to each auxiliary node, with the forced constraint weight value being greater than the soft constraint weight value. Based on the coordinates of all central and auxiliary nodes and their corresponding constraint values, the coefficient matrix of the interpolation model is constructed using the radial basis function. The model parameters of the radial basis function interpolation model are obtained by solving the coefficient matrix, thereby achieving accurate and continuous estimation of elevation adjustments in non-abrupt regions.

[0013] To obtain better elevation correction results, after completing the elevation correction of the entire area using UAV oblique photogrammetry under complex terrain, an iterative refinement step was set up, including: acquiring the corrected elevation data of the entire area; extracting elevation profile curves from vertical profiles in multiple different directions; performing second-order difference calculations on each elevation profile curve to obtain the elevation curvature change value at each point; identifying points where the elevation curvature change value exceeds a preset curvature threshold as residual error points, and collecting local point cloud data at the locations of the residual error points; re-matching feature points and recalculating elevations on the local point cloud data, and replacing the original elevation values ​​at the residual error points with the recalculated elevation values ​​to complete the iterative refinement of elevation correction. Through targeted reprocessing for residual errors, hidden minor errors are further eliminated, improving the accuracy of the results in complex terrains such as steep slopes and fault lines.

[0014] The specific process of re-matching feature points and recalculating elevations in the local point cloud data is as follows: Within a defined radius around the residual error points, all image feature points belonging to the local point cloud data are extracted; using a pixel-by-pixel displacement search method, the extracted image feature points are re-matched on oblique images from adjacent viewpoints to generate a densely matched disparity map for each image feature point; the disparity value of each image feature point is read from the densely matched disparity map, and converted into an elevation value according to the photogrammetric collinearity equation; the elevation values ​​of all image feature points are then subjected to median filtering to obtain the final corrected elevation value for each point in the local point cloud data. This iterative refinement process significantly enhances the adaptability of the elevation correction method to strongly abrupt terrain changes and minor terrain faults, ultimately obtaining high-precision elevation data that more faithfully reflects the true surface morphology.

[0015] The technical effects and advantages provided by the present invention in the above technical solution are as follows: By identifying abrupt and non-abrupt terrain regions based on the spatial distribution density of points in the initial 3D point cloud data, and calculating the disparity gradient change within abrupt regions using image feature points located on both sides of the abrupt boundary, the disparity changes on both sides of the abrupt boundary are drastic. This disparity gradient change can quantitatively reflect the intensity of image disparity distortion caused by terrain discontinuity. Based on this, point-by-point local error inversion compensation is performed on the original elevation observations within the abrupt region to generate a corrected elevation datum. This can specifically correct elevation anomalies at abrupt points without weakening the true abrupt terrain characteristics, avoiding the erroneous smoothing of terrain abrupt information by global filtering methods, and significantly reducing elevation jump errors at abrupt locations such as steep slopes and cliffs. In non-abrupt regions, the generated corrected elevation datum is used as a constraint boundary. By using the elevation values ​​on the boundary line of the abrupt terrain region as a mandatory constraint and the selected control reference points within the non-abrupt region as soft constraints, a radial basis function interpolation model is constructed to fit and adjust the original elevation observations within the non-abrupt region. This method enables a smooth transition of elevation values ​​in non-abrupt areas under the guidance of the calibration datum, eliminating the edge discontinuity effect caused by abrupt correction. At the same time, it preserves the topographic undulation details within non-abrupt areas, prevents residual errors at abrupt changes from spreading to flat areas, and ensures that the elevation data of the entire region are naturally connected and numerically reliable at abrupt and non-abrupt topographic changes. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0017] Figure 1This is a flowchart of the elevation correction method for UAV oblique photogrammetry in complex terrain; Figure 2 This is a flowchart for identifying areas of abrupt terrain change. Figure 3 This is a flowchart for calculating the parallax gradient change in areas of abrupt terrain change. Figure 4 This is a flowchart of the local error inversion compensation process for the elevation datum surface based on the parallax gradient change. Figure 5 This is a flowchart of the elevation fitting adjustment process for non-abrupt regions; Figure 6 This is a flowchart of the elevation correction iterative refinement process; Figure 7 This is a schematic diagram of the initial 3D point cloud data spatial distribution and the boundaries of terrain abrupt change areas; Figure 8 It is a curve comparing the original elevation observation value and the corrected elevation datum at the boundary control point of the area with abrupt topographic change; Figure 9 This is a comparison chart of the original elevation observations and corrected elevation fitting values ​​in the non-mutation region; Figure 10 This is a schematic diagram of the second-order difference elevation curvature variation and residual error points of the elevation profile curve. Detailed Implementation

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

[0019] See Figure 1 This invention provides a method for elevation correction in UAV oblique photogrammetry under complex terrain. The overall implementation scheme is as follows: Acquire multi-view oblique image data collected by the UAV in a complex terrain area, and simultaneously acquire the original elevation observation values ​​corresponding to each image recorded by the airborne positioning system; extract feature points and perform stereo matching on the multi-view oblique image data to generate initial three-dimensional point cloud data; identify terrain abrupt change areas and non-abrupt change areas based on the spatial distribution density of points in the initial three-dimensional point cloud data; in terrain abrupt change areas, extract image feature points located on both sides of the abrupt change boundary and calculate their disparity gradient changes; perform local error inversion compensation on the original elevation observation values ​​based on the disparity gradient changes to generate a corrected elevation datum; in non-abrupt change areas, use the corrected elevation datum to fit and adjust the original elevation observation values, completing the UAV oblique photogrammetry elevation correction for the entire area under complex terrain.

[0020] Example 1:

[0021] In specific implementation, please refer to Figure 2 The process of identifying abrupt and non-abrupt terrain regions based on the spatial distribution density of points in the initial 3D point cloud data is as follows.

[0022] The initial 3D point cloud data is divided into multiple equally sized grid cells according to the horizontal projection coordinates. The horizontal projection coordinates refer to the X and Y coordinates of each point in the 3D point cloud on the horizontal plane, using the local horizontal coordinate system adopted by UAV oblique photogrammetry. When dividing the grid cells, the minimum and maximum values ​​Xmin and Ymax in the X direction, and the minimum and maximum values ​​Ymin and Ymax in the Y direction of the initial 3D point cloud data are first obtained, forming the horizontally bounding rectangular region of the point cloud. The side length L of the grid cells is set, determined based on the overall average point spacing of the initial 3D point cloud data. The overall average point spacing is obtained by averaging the nearest neighbor distances of all points, and the side length L is 2 to 5 times the overall average point spacing. The horizontally bounding rectangular region is divided into multiple columns along the X direction starting from Xmin with a step size L, and into multiple rows along the Y direction starting from Ymin with a step size L, resulting in several grid cells of size L×L. Each grid cell has a coverage range of [Xmin+i×L,Xmin+(i+1)×L] in the X direction and [Ymin+j×L,Ymin+(j+1)×L] in the Y direction, where i and j are integer indices.

[0023] Count the number of points contained in each grid cell. Iterate through all points in the initial 3D point cloud data. For each point, determine the grid cell it falls into based on its horizontal projection coordinates, and increment the point count within that grid cell by 1. After completing the iteration, obtain the number of points in each grid cell. Calculate the spatial distribution density of points in each grid cell. The spatial distribution density is equal to the number of points in that grid cell divided by the area of ​​the grid cell. The area of ​​each grid cell is fixed at L×L.

[0024] Grid cells with a spatial distribution density value exceeding a preset density threshold are marked as candidate mutation cells. The preset density threshold is determined as follows: calculate the average value μ and standard deviation σ of the spatial distribution density values ​​of all grid cells, and set the preset density threshold to μ+2σ. For each grid cell, if the spatial distribution density value of the grid cell is greater than μ+2σ, then the grid cell is marked as a candidate mutation cell.

[0025] Four-connected clustering is performed on all labeled candidate mutation units. A four-connected component is a grid cell that is connected to its four horizontally adjacent grid cells above, below, to the left, and to the right. The clustering process is as follows: A label matrix corresponding to the grid division is established, and the matrix elements corresponding to all candidate mutation units are initialized to an unvisited state. All candidate mutation units are traversed. When an unvisited candidate mutation unit is encountered, a new connected component label is created, and four-connected region growth is performed starting from this candidate mutation unit. During region growth, the current candidate mutation unit is marked as visited and assigned a connected component label. The adjacent grid cells above, below, to the left, and to the right are checked. If an adjacent grid cell is also a candidate mutation unit and has not been visited, the adjacent grid cell is added to the growth queue, and processing continues until the queue is empty. After visiting all candidate mutation units, candidate mutation units with the same connected component label form a connected region. The interconnected candidate mutation units are merged, and each connected region forms a topographic mutation region. Of all grid cells, those that are classified as terrain abrupt change regions are marked as non-abrupt change regions.

[0026] See Figure 7 The figure shows the distribution of the initial 3D point cloud data on the horizontal projection plane. The horizontal axis represents the horizontal projection X-coordinate, ranging from 0 meters to 100 meters, and the vertical axis represents the horizontal projection Y-coordinate, also ranging from 0 meters to 100 meters. The point cloud data is presented as black scattered dots, with a relatively sparse overall distribution. Three areas of dense point cloud clusters are clearly visible in the figure. These three dense point cloud areas are located at approximate coordinates (20,25), (65,30), and (45,70), respectively. Their boundaries are marked with black dashed rings in the figure, representing the boundaries of terrain abrupt change areas identified in Example 1 based on spatial point distribution density statistics and threshold determination.

[0027] Specifically, the point cloud density inside the three circular dashed boundaries is significantly higher than that in the outer areas, indicating local point cloud aggregation, which aligns with the process in Example 1 of identifying abnormally dense grid cells by dividing the area into horizontal grid cells and counting the number of points. The shape and size of the dashed boundaries reflect local variations in the spatial distribution of the point cloud, identifying abrupt terrain changes in complex terrain. The points outside the dashed boundaries are more dispersed, representing non-abrupt areas.

[0028] This figure effectively illustrates how grid cells are divided based on the horizontal projection coordinates of the point cloud, the density of points is statistically analyzed, the mean and standard deviation of the density are calculated, and a threshold (mean plus twice the standard deviation) is set to screen candidate abrupt change units. Subsequently, a four-connected domain clustering algorithm is used to merge adjacent high-density grid cells, ultimately forming a specific illustration of the terrain abrupt change region. This provides a basis for spatial region division for elevation correction based on terrain abrupt change regions in subsequent embodiments.

[0029] Example 2:

[0030] In specific implementation, please refer to Figure 3 The process of extracting image feature points located on both sides of the abrupt change boundary within the terrain change area and calculating the disparity gradient change is as follows.

[0031] Multiple boundary control points were selected along the boundary line of the terrain abrupt change region. The boundary line of the terrain abrupt change region consists of closed polygons formed by the merging of four connected domains in Example 1, and the boundary line is composed of a series of connected grid cell edge segments. The boundary control points were selected as follows: points were collected sequentially along the boundary line from the starting point at a set arc length interval. The arc length interval was set based on the overall average point spacing of the initial 3D point cloud data. The overall average point spacing was obtained by calculating the nearest neighbor distance of all points and taking the average value. The arc length interval was twice the overall average point spacing. The collected points were used as boundary control points.

[0032] For each boundary control point, determine the direction perpendicular to the boundary line. This direction is determined as follows: Obtain two adjacent boundary points on either side of the boundary control point on the boundary line. Connect these two points to form a chord. Calculate the direction vector of this chord on the horizontal plane. Rotate this direction vector counterclockwise by 90 degrees to obtain a reference direction vector perpendicular to the boundary line and pointing outwards from the abrupt change area. Rotate this direction vector clockwise by 90 degrees to obtain a reference direction vector perpendicular to the boundary line and pointing inwards from the abrupt change area. Extend a predetermined distance from the boundary control point in the direction perpendicular to the boundary line and pointing inwards from the abrupt change area to form an inner sampling zone. Extend a predetermined distance from the boundary control point in the direction perpendicular to the boundary line and pointing outwards from the abrupt change area to form an outer sampling zone. The predetermined distance is determined based on the scale of the abrupt change area, taking one-tenth of the equivalent radius of the abrupt change area, which is the radius of a circle with the same area as the abrupt change area. The inner sampling zone is a rectangular strip that starts from the boundary control point, has the same width as the set distance, and extends inward along a direction perpendicular to the boundary line. The outer sampling zone is a rectangular strip that starts from the boundary control point, has the same width as the set distance, and extends outward along a direction perpendicular to the boundary line.

[0033] Image feature points with the same name are extracted within both the inner and outer sampling bands. The extraction process is as follows: Within the spatial range covered by the inner sampling band, all 3D points located within this spatial range are selected from the initial 3D point cloud data, and the corresponding image feature points are obtained. These image feature points are pixels with feature descriptors determined during feature point extraction and stereo matching in the multi-view oblique image data. Similarly, within the spatial range covered by the outer sampling band, all 3D points located within this spatial range are selected from the initial 3D point cloud data, and the corresponding image feature points are obtained. Image feature points with the same name refer to pixels corresponding to the same 3D point on oblique images from different viewpoints. These pixels have the same texture features in the oblique images from different viewpoints, and these texture features are characterized by feature descriptors. For each 3D point within the inner sampling band, record the pixel coordinates of the 3D point on at least two different angled tilted images to form a set of image feature points with the same name for the inner sampling band; for each 3D point within the outer sampling band, record the pixel coordinates of the 3D point on at least two different angled tilted images to form a set of image feature points with the same name for the outer sampling band.

[0034] Calculate the disparity value of each corresponding image feature point on tilted images at different viewing angles. The disparity value is calculated as follows: For a set of corresponding image feature points, select a reference viewing angle tilted image, calculate the Euclidean distance between the pixel coordinates of every two different viewing angles for each corresponding image feature point in the set, and take the average of all Euclidean distances as the disparity value of that corresponding image feature point. Calculate the arithmetic mean of the disparity values ​​of all corresponding image feature points within the inner sampling band as the inner average disparity. Calculate the arithmetic mean of the disparity values ​​of all corresponding image feature points within the outer sampling band as the outer average disparity.

[0035] The absolute value of the difference between the inner average disparity and the outer average disparity is taken as the disparity gradient change at the boundary control point, expressed by the formula:

[0036] in, This represents the disparity gradient change at the b-th boundary control point, where b is the index of the boundary control point. The value of b is a positive integer, ranging from 1 to the total number of boundary control points. The inner average disparity represents the average disparity of all image feature points with the same name within the inner sampling band. The inner average disparity is obtained by summing the disparity values ​​of all image feature points with the same name within the inner sampling band and then dividing by the number of image feature points with the same name within the inner sampling band. The average disparity value represents the sum of the disparity values ​​of all image feature points with the same name within the outer sampling band, i.e., the outer average disparity. The outer average disparity is obtained by summing the disparity values ​​of all image feature points with the same name within the outer sampling band and then dividing by the number of image feature points with the same name within the outer sampling band.

[0037] Example 3

[0038] In specific implementation, please refer to Figure 4 The process of performing local error inversion compensation on the original elevation observations based on the parallax gradient change and generating the corrected elevation datum is as follows.

[0039] The disparity gradient changes at each boundary control point are sorted in descending order. The disparity gradient changes are calculated in Example 2. Let 'b' be a positive integer, ranging from 1 to the total number of boundary control points. The boundary control points are sorted in descending order, with the boundary control points exhibiting the largest disparity gradient change at the beginning and the smallest at the end, resulting in a sorted sequence. The top few boundary control points are selected as strong distortion feature points. The number of selected points is determined by a preset proportion of the total number of boundary control points, set at 20%. This is based on the following reasoning: in complex terrain scenes, elevation distortion is mainly caused by the most drastic disparity changes in abrupt boundary changes. Selecting the top 20% of boundary control points can cover the main distortion source areas while avoiding the introduction of too many non-dominant distortion points, which could lead to excessive overlap of the inversion window and computational redundancy. When the total number of boundary control points is N, the number of strong distortion feature points is ceil(0.2×N), where ceil represents rounding up.

[0040] A local inversion window is established centered on each strongly distorted feature point. The local inversion window is a circular region on the horizontal projection plane centered on the strongly distorted feature point and with an inversion radius as its radius. The radius of the local inversion window is dynamically and adaptively adjusted according to the boundary curvature of the terrain abrupt change region. The boundary curvature is calculated as follows: on the boundary line of the terrain abrupt change region, using the boundary control point corresponding to the strongly distorted feature point as the reference point, an arc length interval is extended to both sides along the boundary line. This arc length interval is consistent with the arc length interval used when selecting the boundary control points in Example 2, i.e., twice the overall average point spacing, resulting in two adjacent boundary points. A circular arc is fitted using the reference point and the two adjacent boundary points, and the curvature of the arc is calculated as the boundary curvature at the strongly distorted feature point. When the boundary curvature is greater, the radius of the local inversion window is smaller. The radius is inversely proportional to the boundary curvature. The radius of the local inversion window is the product of the reciprocal of the boundary curvature at the strong distortion feature point and a reference radius coefficient. The reference radius coefficient is one-tenth of the equivalent radius of the terrain abrupt change region. The equivalent radius of the terrain abrupt change region is the radius of a circle with the same area as the terrain abrupt change region.

[0041] All raw elevation observations are collected within the local inversion window. These raw elevation observations are the elevation values ​​corresponding to each image recorded by the airborne positioning system, and each raw elevation observation corresponds to a 3D point coordinate. The raw elevation observations corresponding to all points whose horizontal projection coordinates fall within the circular area of ​​the local inversion window are collected into the dataset of that local inversion window.

[0042] Based on the disparity gradient change corresponding to the strongly distorted feature point, an inversion compensation coefficient is calculated for each original elevation observation within the local inversion window. The calculation process for the inversion compensation coefficient is as follows: Determine the coordinates of the center point of the local inversion window, where the center point coordinates are the horizontal projection coordinates of the boundary control point corresponding to the strongly distorted feature point; calculate the planar distance between the coordinates of each original elevation observation within the local inversion window and the coordinates of the center point; divide the planar distance corresponding to each coordinate by the radius of the local inversion window to obtain the normalized distance value for each coordinate; using the disparity gradient change of the strongly distorted feature point as the baseline compensation intensity value, multiply the baseline compensation intensity value by one and subtract the normalized distance value to obtain the inversion compensation coefficient at each coordinate, expressed by the formula:

[0043] in, This represents the inversion compensation coefficient of the k-th original elevation observation within the local inversion window corresponding to the q-th strong distortion feature point. The value of q is a positive integer, ranging from 1 to the total number of strong distortion feature points, and the value of k is a positive integer, ranging from 1 to the total number of original elevation observations within the q-th local inversion window. The disparity gradient change corresponding to the qth strongly distorted feature point is used as the baseline compensation intensity value. The disparity gradient change is obtained by the disparity gradient change calculation method at the boundary control point in Example 2. This represents the planar distance between the coordinates of the point corresponding to the k-th original elevation observation within the q-th local inversion window and the coordinates of the center point of the local inversion window. It is obtained by calculating the Euclidean distance between the two points on the horizontal projection plane. Represents the radius of the q-th local inversion window. Inversion compensation coefficients. With disparity gradient change There is a positive correlation; the greater the change in the disparity gradient, the greater the magnitude of the inversion compensation coefficient. At the same time, the inversion compensation coefficient decreases linearly with the increase of the distance from the point to the center point within the local inversion window.

[0044] Each original elevation observation within the local inversion window is multiplied by its corresponding inversion compensation coefficient to obtain the corrected local elevation value. The corrected local elevation value represents the compensated elevation value of that point within the local inversion window.

[0045] Weighted fusion of the corrected local elevation values ​​from all local inversion windows is performed to generate a corrected elevation datum covering the terrain abrupt change area. For any target point within the terrain abrupt change area that lies within the overlapping area of ​​multiple local inversion windows, multiple corrected local elevation values ​​corresponding to that target point in each local inversion window are collected. The target point refers to the grid node or original point cloud point within the terrain abrupt change area that requires elevation correction. The spatial distance from the target point to the strongly distorted feature point corresponding to each local inversion window is obtained. The spatial distance is the Euclidean distance between the 3D coordinates of the target point and the 3D coordinates of the boundary control point corresponding to the strongly distorted feature point. The reciprocal of each spatial distance is used as the fusion weight value for the target point in that local inversion window. The corrected local elevation values ​​of the target point are multiplied by their corresponding fusion weight values, summed, and then divided by the sum of all fusion weight values ​​to obtain the fused elevation value of the target point. Traverse all target points within the terrain abrupt change area, calculate the fused elevation value for each target point, and combine the fused elevation values ​​of all target points to form a corrected elevation datum covering the terrain abrupt change area.

[0046] See Figure 8 In the figure, the horizontal axis represents the spatial horizontal position in meters (m), ranging from 0 to 200 meters; the vertical axis represents the elevation value of the corresponding position in meters (m), ranging from approximately 148 meters to 183 meters. The legend shows that the gray curve represents the "original elevation observation value," and the black curve represents the "corrected elevation datum." The curves, connected by continuous points, reflect the changing trend of terrain elevation within that horizontal range.

[0047] As can be seen from the figure, the original elevation observation curve exhibits a distinct elevation peak at approximately 85 meters, with a peak value approaching 182 meters, indicating a certain degree of elevation distortion. A secondary peak also occurs at approximately 135 meters, with a peak value of about 166 meters. In contrast, the peak values ​​of the corrected elevation datum curve within the same location range are significantly suppressed. The peak value at 85 meters decreases to approximately 174 meters, while the peak value at 135 meters increases to approximately 172 meters. The overall curve is smoother, eliminating the local anomalies in the original observations. Furthermore, the corrected curve maintains a high degree of consistency with the original elevation observations in the 0–50 meter and 150–200 meter ranges, indicating that the correction method reasonably fitted and adjusted the elevation in non-abrupt regions, preserving the original topographic trend.

[0048] The data in this figure reflects the effect of generating a corrected elevation datum after local error inversion compensation based on the disparity gradient change in Example 3. By identifying strong distortion feature points, establishing a local inversion window, calculating the inversion compensation coefficients, and weighted fusion, effective correction of the original elevation observations in abrupt change areas is achieved, significantly reducing local elevation distortion and improving the accuracy and continuity of elevation data.

[0049] Example 4

[0050] In specific implementation, please refer to Figure 5 The process of fitting and adjusting the original elevation observations using the corrected elevation datum in non-abrupt regions is as follows.

[0051] A series of boundary elevation values ​​are extracted from the corrected elevation datum on the boundary line of the terrain abrupt change area. The corrected elevation datum is obtained from Example 3. The boundary line of the terrain abrupt change area is consistent with the boundary line defined in Example 2, and is a closed polygon composed of a series of connected grid cell edge segments. Sampling points are collected along the boundary line at uniform arc length intervals. The arc length interval is the same as the arc length interval when selecting boundary control points in Example 2, and is twice the overall average point spacing. For each sampling point, the corresponding elevation value is obtained in the corrected elevation datum through bilinear interpolation based on the horizontal projection coordinates of the sampling point. This elevation value is used as the boundary elevation value. At the same time, the three-dimensional coordinates of the sampling point are recorded. The three-dimensional coordinates consist of the horizontal projection coordinates and the boundary elevation value.

[0052] Multiple original elevation observations were selected as control reference points within the non-abrupt terrain region. An adaptive point placement strategy was employed, with higher density of control reference points closer to the boundary of the abrupt terrain change zone within the non-abrupt region. The specific implementation of the adaptive point placement strategy involved calculating the shortest planar distance from each original elevation observation point within the non-abrupt region to the boundary line of the abrupt terrain change zone, denoted as 'd'. A minimum point spacing was then set. and maximum spacing between dots , Take 1 times the overall average point spacing. Take 5 times the overall average point spacing. For each original elevation observation point, generate a random number r, which is uniformly distributed between 0 and 1. Calculate the acceptance probability p of this original elevation observation point being selected as a control reference point, where p takes the value of:

[0053] in, This represents the shortest planar distance from the original elevation observation point to the boundary line of the terrain abrupt change area; Indicates the minimum spacing between dots; This represents the maximum spacing between dots. When d is less than... When d is greater than 1, the acceptance probability p is set to 1, meaning that the point will definitely be selected as a control reference point; when d is greater than 1, the acceptance probability p is set to 1. dmaxWhen the acceptance probability p is set to 0, the point is not selected as a control reference point. For each original elevation observation point, if the random number r is less than the acceptance probability p, then the original elevation observation point is selected as a control reference point. Through the above probability screening, in areas closer to the boundary, d is smaller and p is larger, resulting in a higher probability of the control reference point being selected, thus forming a distribution with a high density of control reference points near the boundary and a low density far from the boundary.

[0054] A radial basis function (RBF) interpolation model is constructed using boundary elevation values ​​as mandatory constraints and control reference points as soft constraints. The specific steps for constructing the RBF interpolation model are as follows: The coordinates of the points corresponding to the boundary elevation values ​​are used as the center nodes of the RBF, and the boundary elevation value is set at each center node as a mandatory pass-through point for the interpolation model. The coordinates of the points corresponding to the control reference points are used as auxiliary nodes of the RBF, and the original elevation observation value of the control reference point is set at each auxiliary node as a soft constraint point for the interpolation model. A mandatory constraint weight value is assigned to each center node, and a soft constraint weight value is assigned to each auxiliary node, with the mandatory constraint weight value being greater than the soft constraint weight value. The mandatory constraint weight value is uniformly set to 1.0, and the soft constraint weight value is uniformly set to 0.1. The rationale is as follows: the mandatory constraint requires the interpolation surface to pass precisely through the center node, thus assigning it the largest weight; the soft constraint allows for a certain deviation of the interpolation surface near the auxiliary nodes to avoid interference from local point cloud noise on the overall fitting of non-abrupt regions. The weight value of 0.1 was determined experimentally, which can smooth out random errors while preserving the terrain trend.

[0055] Based on the coordinates of all central and auxiliary nodes and their corresponding constraint values, a coefficient matrix for the interpolation model is constructed using radial basis functions. Solving the coefficient matrix yields the model parameters of the radial basis function interpolation model. The radial basis functions are selected as quadratic functions, in the form of... ,in, The shape parameter is determined based on the average spacing between the center node and auxiliary nodes, and is taken as 0.5 times the average spacing. Let there be... Each central node and There are auxiliary nodes, and the total number of nodes is N= + The interpolation model is expressed as follows: for any spatial point coordinates Its elevation adjustment amount Calculated using the following formula:

[0056] in, This represents the horizontal projection coordinates corresponding to the i-th center node, where i ranges from 1 to... , The total number of central nodes; This represents the horizontal projected coordinates of the j-th auxiliary node, where j ranges from 1 to... , This represents the total number of auxiliary nodes. Indicated by point To the central node The Euclidean distance is the value of a quadratic function of the independent variable. Indicated by point To auxiliary node The Euclidean distance is the value of a quadratic function of the independent variable; This represents the model coefficients to be determined for the i-th center node. Let represent the model coefficients to be determined for the j-th auxiliary node. Substitute the coordinates of all central and auxiliary nodes, and enforce the boundary elevation adjustment (the difference between the boundary elevation value and the original elevation observation value) for the central nodes and the soft constraint value for the auxiliary nodes (the soft constraint value is set to 0, indicating that the original elevation observation value should be as close as possible to the original value after elevation adjustment). Construct a system of linear equations containing N equations and N unknowns. After adding weighting factors, solve for all the unknowns. and The weighting factors are added by multiplying the left side of the equation corresponding to the central node in the equation system by the mandatory constraint weight value of 1.0, and the left side of the equation corresponding to the auxiliary node by the soft constraint weight value of 0.1. The model parameters are obtained by solving the linear equation system.

[0057] Each original elevation observation within the non-mutation region is input into the radial basis function interpolation model. For any original elevation observation within the non-mutation region, the corresponding point coordinates... ,Will Substituting into the above interpolation model formula, the elevation adjustment amount is calculated. The original elevation observations within the non-mutation region. Add the corresponding elevation adjustment The elevation values ​​after correction for non-mutation regions were obtained. .

[0058] See Figure 9 In the figure, the horizontal axis represents the horizontal position in meters, ranging from 0 to 150 meters, and the vertical axis represents the elevation in meters, ranging from approximately 98 to 136 meters. The gray solid line in the figure represents the original elevation observation value, and the black dashed line and triangular marker line represent the corrected elevation value after fitting and adjusting through the radial basis function interpolation model.

[0059] As can be seen from the curve trend in the figure, the original elevation observations show a relatively stable upward trend overall, starting at approximately 100 meters and gradually rising to around 135 meters as the horizontal position increases. There are some fluctuations in the original elevation observations, especially in the 100-140 meter range, where the fluctuation amplitude is relatively significant, indicating that the elevation data contains a certain amount of random error and noise.

[0060] The corrected elevation curves maintain the overall trend of the original elevation observations, but with significantly reduced fluctuations and smoother curves. In particular, in the region between 40 and 80 meters, the corrected curves effectively suppress short fluctuations in the original data, demonstrating the effective fitting and smoothing capabilities of the radial basis function interpolation model for elevation observations in non-abrupt regions.

[0061] This figure illustrates the effectiveness of using a corrected elevation datum along the boundary line of a terrain abrupt change area as a mandatory constraint within a non-abrupt region, combined with control reference points selected using an adaptive point placement strategy, to construct a radial basis function interpolation model for fitting and adjusting the original elevation observations. This method effectively filters out random noise within the non-abrupt region while ensuring boundary elevation constraints, achieving smooth and accurate correction of elevation data. This aligns with the technical solution described in Example 4, which uses a radial basis function interpolation model to fit and adjust elevation values ​​in non-abrupt regions.

[0062] Example 5

[0063] In specific implementation, please refer to Figure 6 After completing the elevation correction of the entire area under complex terrain using UAV oblique photogrammetry, the process of iterative refinement of the elevation correction is as follows.

[0064] Obtain the corrected elevation data for the entire region. The full-region elevation data is obtained by merging the corrected elevation values ​​of the non-abrupt terrain areas with the corrected elevation datum surface generated in Example 3 within the terrain abrupt change areas, covering the entire complex terrain area measured by UAV oblique photogrammetry. The corrected full-region elevation data is stored in the form of a regular grid or point cloud, with each elevation data point containing horizontal projected coordinates and the corresponding elevation value.

[0065] Elevation profile curves are extracted from vertical profiles in multiple directions. The vertical profiles are set in four directions: due north, due east, northeast, and northwest, covering terrain features with different orientations. The vertical profiles are laid out as follows: parallel profile lines are evenly distributed at preset intervals on a horizontal plane covering the entire area's elevation data. The interval is 10 times the overall average point spacing, which is obtained by averaging the nearest neighbor distances of all points. For each vertical profile line, all elevation data points traversed by the line are arranged in order of distance along the profile line, and the elevation value of each data point is extracted to form an elevation profile curve. The elevation profile curve is represented as a function of the distance variable t along the profile line. t is the cumulative planar distance along the profile line starting from the starting point.

[0066] For each elevation profile curve, a second-order difference calculation is performed to obtain the elevation curvature change value at each point. The second-order difference calculation method is as follows: for three consecutive points arranged in sequence on the elevation profile curve, with point numbers m-1, m, and m+1, the corresponding elevation values ​​are respectively... , , The corresponding distances along the section line are respectively , , Calculate the forward difference interval. and backward differential spacing Location The first-order difference value is calculated using the central difference formula, and the second-order difference value is calculated using the difference of the first-order difference. The change in elevation curvature is the absolute value of the second-order difference value, expressed by the formula:

[0067] in, This indicates that the distance along the profile line is... The change in elevation curvature at the point; Indicates the distance along the profile line as The elevation value at the point; Indicates the distance along the profile line as The elevation value at the point; Indicates the distance along the profile line as The elevation value at the point; Indicates from arrive Forward differential spacing, The value is the actual planar distance between two adjacent points along the profile line; Indicates from arrive backward differential spacing, The value of is the actual planar distance between two adjacent points along the profile line; m is the index of the point on the elevation profile curve, and the value of m is a positive integer greater than or equal to 2. For the first and last endpoints of the elevation profile curve, since a three-point continuous sequence cannot be formed, second-order difference calculation is not performed.

[0068] Points whose elevation curvature changes exceed a preset curvature threshold are identified as residual error points. The preset curvature threshold is determined by statistically analyzing the distribution of elevation curvature changes across the entire region's elevation profile curves, using the 95th percentile of all elevation curvature changes in the entire region as the preset threshold. For each point where the elevation curvature change is calculated, if the change exceeds the preset curvature threshold, the three-dimensional coordinates of that point are recorded as a residual error point.

[0069] Collect local point cloud data at the locations of residual error points. The local point cloud data is collected as follows: Using the horizontal projection coordinates of the residual error points as the center, define a radius around which all points belonging to the initial 3D point cloud data are collected. The radius is taken as 5 times the overall average point spacing.

[0070] Feature point matching and elevation calculation are performed again on the local point cloud data. Within a radius defined around the residual error points, all image feature points belonging to the local point cloud data are extracted. Image feature points are pixels with feature descriptors determined during feature point extraction and stereo matching in multi-view oblique image data. Each image feature point corresponds to a three-dimensional point in the local point cloud data, and the pixel coordinates of each image feature point are recorded on at least two different view oblique images.

[0071] A pixel-wise displacement search method is employed to re-match extracted image feature points on oblique images from adjacent viewpoints. The specific implementation of the pixel-wise displacement search method is as follows: For each image feature point in the local point cloud data, a reference image is determined. A template window of size W×W is extracted from the reference image, centered on the pixel coordinates of the image feature point. The side length W of the template window is set to 21 pixels. This setting is based on the fact that, at typical UAV oblique photogrammetry resolutions, a 21-pixel window can contain sufficient texture information to distinguish corresponding points while avoiding excessive background interference. On oblique images from adjacent viewpoints that overlap with the reference image, a one-dimensional search is performed along the epipolar direction with a search step size of 1 pixel. The search range is from negative S pixels to positive S pixels along the epipolar direction, and the search range S is set to 30 pixels. This setting is based on the fact that the parallax range between adjacent viewpoints in UAV oblique photogrammetry typically does not exceed 30 pixels. At each search location, a candidate window of the same size as the template window is extracted. The normalized cross-correlation coefficient between the template window and the candidate window is calculated, with the value ranging from -1 to 1. The search location with the largest normalized cross-correlation coefficient is taken as the matching location, and the pixel coordinates of the matching location are used as the re-matching coordinates of the image feature point on adjacent viewpoint images. For each image feature point in the local point cloud data, the above pixel-by-pixel displacement search and matching is performed on the tilted images of each adjacent viewpoint, generating a dense matching disparity map for each image feature point. The dense matching disparity map records the pixel coordinates of each image feature point after re-matching on image pairs of different viewpoints.

[0072] The disparity value of each image feature point is read from the densely matched disparity map. The disparity value is calculated in the same way as in Example 2. For a set of image feature points with the same name, a tilted image of a reference viewpoint is selected, and the Euclidean distance between the pixel coordinates of each two different viewpoint images in the set of image feature points with the same name is calculated. The average value of all Euclidean distances is taken as the disparity value of the image feature point.

[0073] The parallax values ​​are converted into elevation values ​​using the photogrammetric collinearity equation. The photogrammetric collinearity equation is a fundamental equation in photogrammetry, establishing the mathematical relationship between image point coordinates, object-space coordinates, and camera interior / exterior orientation elements. For each image feature point, given the rematched parallax value, the object-space 3D coordinates corresponding to the feature point are solved using spatial forward intersection, combined with the interior orientation elements obtained from camera calibration and the exterior orientation elements recorded by the airborne positioning system. The elevation value is then extracted from these object-space 3D coordinates. This transformation is performed on all image feature points to obtain the elevation value of each image feature point in the local point cloud data.

[0074] Median filtering is applied to the elevation values ​​of all image feature points. The specific implementation of median filtering is as follows: For each 3D point in the local point cloud data, a spherical neighborhood is determined with the point as its center and a radius of three times the overall average point spacing. Elevation values ​​of all image feature points within this spherical neighborhood are collected and sorted numerically. The median of these sorted values ​​is taken as the median-filtered elevation value for that 3D point. If the number of image feature points within the spherical neighborhood is odd, the median is the elevation value located in the middle after sorting; if the number of image feature points is even, the median is the arithmetic mean of the two middle elevation values ​​after sorting. After median filtering, the final corrected elevation value for each point in the local point cloud data is obtained. The original elevation values ​​at residual error points are replaced with the recalculated final corrected elevation values, completing the iterative refinement of elevation correction.

[0075] See Figure 10 In the figure, the horizontal axis represents the distance along the set profile line, in meters, ranging from 0 to 200 meters. The vertical axis represents the change in elevation curvature, in units of 1 / m, with a value range of approximately 0 to 0.018. The curve shows the trend of elevation curvature change with the distance from the profile line. The solid black line shows the continuous change in elevation curvature. The scattered "×" marks indicate the locations of identified residual error points, and the elevation curvature change values ​​corresponding to these points are all large outliers.

[0076] As can be seen from the figure, the elevation curvature changes are mostly concentrated in the lower range below 0.0025, and the curve generally exhibits small fluctuations, reflecting that the elevation changes in most areas are gentle and the curvature changes are small. However, at approximately 50 meters, 105 meters, 155 meters, and 200 meters, the curve shows obvious peaks, corresponding to multiple residual error points. The highest peak reaches 0.018, which is significantly higher than the curvature changes in other areas, indicating that there are abrupt changes in the terrain or large residual errors in elevation correction at these locations.

[0077] The distribution of residual error points is consistent with the peak height of elevation curvature change, indicating that the elevation curvature change value calculated by second-order difference in this embodiment can effectively locate points where there are still large errors after elevation correction. The elevation curvature change value and residual error points shown in this figure are the key judgment criteria in the elevation correction iterative refinement step in Example 5. Based on this result, local point cloud data around the residual error points are subsequently selected for re-feature point matching and elevation calculation, thereby further improving the elevation correction accuracy in complex terrain areas.

[0078] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for elevation correction in UAV oblique photogrammetry under complex terrain, characterized in that, The steps include the following: Acquire multi-view oblique image data collected by UAVs in complex terrain areas, and simultaneously acquire the original elevation observation values ​​corresponding to each image recorded by the airborne positioning system; Feature point extraction and stereo matching are performed on multi-view oblique image data to generate initial 3D point cloud data; Based on the spatial distribution density of points in the initial 3D point cloud data, regions with abrupt terrain changes and regions without abrupt changes are identified. Within the terrain abrupt change area, image feature points located on both sides of the abrupt change boundary are extracted and their disparity gradient changes are calculated; Local error inversion compensation is performed on the original elevation observations based on the parallax gradient change to generate a corrected elevation datum. In non-abrupt areas, the corrected elevation datum is used to fit and adjust the original elevation observations to complete the elevation correction of the entire area under complex terrain by UAV oblique photogrammetry. This includes: extracting a series of boundary elevation values ​​of the corrected elevation datum on the boundary line of the terrain abrupt change area, and selecting multiple original elevation observations as control reference points in non-abrupt areas. A radial basis function interpolation model is constructed using boundary elevation values ​​as mandatory constraints and control reference points as soft constraints. Each original elevation observation value within the non-mutation region is input into the radial basis function interpolation model, and the radial basis function interpolation model outputs the elevation adjustment amount corresponding to each original elevation observation value; Add the corresponding elevation adjustment amount to each original elevation observation value in the non-mutation region to obtain the corrected elevation value for the non-mutation region.

2. The method for elevation correction in UAV oblique photogrammetry under complex terrain according to claim 1, characterized in that, The specific steps for identifying abrupt and non-abrupt terrain regions based on the spatial distribution density of points in the initial 3D point cloud data are as follows: The initial 3D point cloud data is divided into multiple grid cells of equal size according to the horizontal projection coordinates; Count the number of points contained in each grid cell and calculate the spatial distribution density value of points in each grid cell; Grid cells whose spatial distribution density values ​​exceed a preset density threshold are marked as candidate mutation cells; Four-connected domain clustering is performed on all labeled candidate mutation units. Interconnected candidate mutation units are merged to form terrain mutation regions, and the regions where the remaining grid units are located are labeled as non-mutation regions.

3. The method for elevation correction in UAV oblique photogrammetry under complex terrain according to claim 2, characterized in that, The specific steps for extracting image feature points located on both sides of the abrupt change boundary and calculating their disparity gradient change within the terrain abrupt change area are as follows: Multiple boundary control points are selected along the boundary line of the terrain change zone. Each boundary control point is extended a set distance in a direction perpendicular to the boundary line to the inside and outside of the change zone, respectively, forming an inner sampling zone and an outer sampling zone. Image feature points with the same name are extracted in the inner and outer sampling bands respectively. These image feature points have the same texture features in tilted images from different viewpoints. Calculate the disparity value of each corresponding image feature point on tilted images at different viewpoints. Take the average disparity value of all corresponding image feature points in the inner sampling band as the inner average disparity, and take the average disparity value of all corresponding image feature points in the outer sampling band as the outer average disparity. The absolute value of the difference between the inner average disparity and the outer average disparity is taken as the disparity gradient change at the boundary control point.

4. The method for elevation correction in UAV oblique photogrammetry under complex terrain according to claim 3, characterized in that, The specific steps for performing local error inversion compensation on the original elevation observations based on the disparity gradient change to generate the corrected elevation datum are as follows: The disparity gradient changes at each boundary control point are sorted in descending order, and the top few boundary control points are selected as strong distortion feature points. A local inversion window is established with each strongly distorted feature point as the center, and all original elevation observations are collected within the local inversion window; Based on the disparity gradient change corresponding to the strong distortion feature point, an inversion compensation coefficient is calculated for each original elevation observation within the local inversion window. The inversion compensation coefficient is positively correlated with the disparity gradient change. Each original elevation observation within a local inversion window is multiplied by its corresponding inversion compensation coefficient to obtain a corrected local elevation value. The corrected local elevation values ​​from all local inversion windows are then weighted and fused to generate a corrected elevation datum covering areas with abrupt topographic changes.

5. The method for elevation correction in UAV oblique photogrammetry under complex terrain according to claim 4, characterized in that, The radius of the local inversion window is dynamically and adaptively adjusted according to the boundary curvature of the terrain abrupt change region; the greater the boundary curvature, the smaller the radius of the local inversion window.

6. The method for elevation correction in UAV oblique photogrammetry under complex terrain according to claim 4, characterized in that, The specific steps for calculating an inversion compensation coefficient for each original elevation observation within a local inversion window are as follows: Determine the coordinates of the center point of the local inversion window, and calculate the planar distance between the coordinates of the point corresponding to each original elevation observation value within the local inversion window and the coordinates of the center point. Divide the planar distance corresponding to each point coordinate by the radius of the local inversion window to obtain the normalized distance value of each point coordinate; Using the disparity gradient change of strongly distorted feature points as the baseline compensation intensity value, the inversion compensation coefficient at each point coordinate is obtained by multiplying the baseline compensation intensity value by one minus the normalized distance value.

7. The method for elevation correction in UAV oblique photogrammetry under complex terrain according to claim 4, characterized in that, The specific steps for weighted fusion of the corrected local elevation values ​​from all local inversion windows to generate a corrected elevation datum covering areas of abrupt topographic changes are as follows: For any target point located within the overlapping area of ​​multiple local inversion windows, collect multiple corrected local elevation values ​​corresponding to the target point in each local inversion window. Obtain the spatial distance from the target point to the strong distortion feature point corresponding to each local inversion window, and use the reciprocal of each spatial distance as the fusion weight value of the local inversion window for the target point; The fused elevation value of the target point is obtained by multiplying each corrected local elevation value of the target point by its corresponding fusion weight value, summing the results, and then dividing by the sum of all fusion weight values. Traverse all target points within the area of ​​abrupt terrain change, and combine the merged elevation values ​​of all target points to form a corrected elevation datum.

8. The method for elevation correction in UAV oblique photogrammetry under complex terrain according to claim 1, characterized in that, The selection of control reference points adopts an adaptive deployment strategy, with a higher density of control reference points located closer to the boundary of terrain change areas within non-change areas.

9. The method for elevation correction in UAV oblique photogrammetry under complex terrain according to claim 8, characterized in that, The specific steps for constructing the radial basis function interpolation model are as follows: The coordinates of the points corresponding to the boundary elevation values ​​are used as the center nodes of the radial basis function, and the boundary elevation value is set at each center node as the forced pass point of the interpolation model; The coordinates of the control reference point are used as auxiliary nodes of the radial basis function, and the original elevation observation value of the control reference point is set at each auxiliary node as a soft constraint point of the interpolation model. Assign a mandatory constraint weight value to each central node and a soft constraint weight value to each auxiliary node, with the mandatory constraint weight value being greater than the soft constraint weight value. Based on the coordinates of all central and auxiliary nodes and their corresponding constraint values, the coefficient matrix of the interpolation model is constructed using radial basis functions, and the model parameters of the radial basis function interpolation model are obtained by solving the coefficient matrix.

Images