A method for extracting high-precision DSM from stereo pairs of atolls without ground control points

Abstract

The application provides a high-precision DSM method for stereo image pair extraction of a no-control-point island reef, and first, a semi-global dense matching algorithm is used for semi-global dense matching to obtain a parallax map; second, a rational function model coordinate stereo positioning is performed to process longitude, latitude and height coordinates to obtain three-dimensional point clouds without deformation; finally, point cloud reconstruction DSM is realized to generate a DSM three-dimensional model. Experimental results show that the high-resolution satellite stereo image pair reconstruction DSM method is preferably applied to World-View3 satellite images and can generate a far-sea island reef DSM three-dimensional model.

Description

Technical Field

[0001] This invention relates to the fields of geographic information systems and cartography, and in particular to a high-precision DSM method for extracting stereo image pairs of islands and reefs without measurement control points. Background Technology

[0002] A digital surface model (DSM) is a digital model that describes the elevation features of the Earth's surface, providing crucial geographic information data for fields such as urban planning, resource management, and environmental monitoring. However, traditional aerial photogrammetry methods cannot currently be used to produce DSMs for remote islands and reefs.

[0003] Currently, various countries and regions are vigorously developing high-resolution stereo mapping satellites. Among them, the WorldView-3 satellite, with its panchromatic resolution of 0.3 meters, is widely used in various fields due to its ultra-high resolution and multi-band data advantages. With the continuous development of my country's stereo mapping satellites, many constructive achievements have also been made. In 2012, my country successfully launched the ZY-3 satellite, its first independently developed civilian high-resolution stereo mapping satellite, which has been widely used for 5-meter resolution DSM / DEM production. Currently, in addition to the ZY-3 01, 02, and 03 satellites, my country also has the Gaofen-7 satellite launched in 2019, which can be used for 2-meter resolution DSM / DEM production, and the Siwei Gaofen-1 and Gaofen-2 satellites launched in 2022. Abundant satellite data sources provide ample room for the production of DSM data, and the ability to guarantee mapping products will become increasingly stronger. Therefore, with the continuous development of high-resolution stereo mapping satellites, the use of high-resolution remote sensing imagery to produce mapping products is becoming increasingly common. However, in the process of generating DSM using satellite stereo image pairs, how to achieve high-precision dense matching, how to solve 3D point clouds, and how to use 3D point clouds to generate DSM are all hot and key issues that urgently need to be addressed. Summary of the Invention

[0004] To address the shortcomings of the prior art, this invention proposes a method for extracting high-precision DSMs from island and reef stereo image pairs without measurement control points, which is used to generate DSMs based on World-View3 satellite stereo image pairs.

[0005] Therefore, the specific technical solution adopted by the present invention is as follows:

[0006] A method for extracting high-precision DSM from stereo image pairs of islands and reefs without measurement control points includes the following steps:

[0007] Step 1: Acquire stereo image pairs of distant islands and reefs from World-View3 satellite;

[0008] Step 2: Based on the stereo image pairs of distant islands and reefs, use the improved semi-global dense matching algorithm to perform semi-global dense matching to obtain the disparity map;

[0009] Step 3: Based on the obtained disparity map, use a rational function model to perform stereo localization and obtain a 3D point cloud;

[0010] Step 4: Preprocess the 3D point cloud obtained in Step 3, and use the 3D point cloud reconstruction DSM algorithm to generate remote island and reef DSM. Then, postprocess the generated remote island and reef DSM to construct a 3D model of the remote island and reef DSM.

[0011] Furthermore, step 3 includes the following specific steps:

[0012] Step 3.1: Based on the disparity value d corresponding to each pixel in the disparity map, obtain the coordinates of the corresponding image points in the left and right images using the following formula:

[0013] r r =r l c r =c l -d

[0014] Where r represents the right image and l represents the left image;

[0015] Step 3.2: Obtain the coordinates of the corresponding image points before cropping the left and right images using the following formula:

[0016] Rn'=Rn+dR

[0017] Cn'=Cn+dC

[0018] Where dR and dC are coordinate transformation coefficients using the same pixels; (Rn,Cn) are the coordinates before image cropping, and (Rn',Cn') are the coordinates after image cropping;

[0019] Step 3.3: Based on the rational function model, calculate the corresponding geodetic coordinate points according to the coordinates of the corresponding image points;

[0020] Step 3.4: Replace the latitude and longitude coordinates of the geodetic points with the pixel coordinates corresponding to the left image;

[0021] Step 3.5: Output the final point cloud coordinates.

[0022] Furthermore, the specific steps of step 3.3 include:

[0023] Step 3.3.1: Calculate the average value of the standardized translation parameters for the left and right images respectively, use these as initial values, and convert these two initial values ​​into standardized coordinates (X). nl ,Y nl Z nl ) and (Xnr ,Y nr Z nr ), thus obtaining the initial value of the geodetic coordinates (X). 0 ,Y 0 Z 0 The specific calculation formula is as follows:

[0024]

[0025] Where l and r are the left and right images, respectively, and X ol The initial left image in X-axis spatial rectangular coordinates, X or The initial right image in Cartesian coordinates along the X-axis, Y... ol The initial left image in Y-axis spatial rectangular coordinates, Y or The initial right image in Y-axis spatial rectangular coordinates, Z ol The initial left image of the Z-axis spatial rectangular coordinate system, Z or The initial right image of the Z-axis spatial rectangular coordinate system, X sl The left image is obtained after iterative calculation using X-axis Cartesian coordinates, and Y... sl The left image is obtained after iterative calculation using Y-axis Cartesian coordinates, Z... sl The left image after iterative calculation using Z-axis Cartesian coordinates, X sr The right image is obtained after iterative calculation using X-axis Cartesian coordinates, and Y... sr The right image is obtained after iterative calculation using Y-axis spatial rectangular coordinates, Z... sr The right image is the result of iterative calculation using Z-axis spatial rectangular coordinates;

[0026] Step 3.3.2: Based on the standardized coordinates of the left and right images, construct the following error equation:

[0027]

[0028] Step 3.3.3: Based on the coordinates of corresponding image points (r) in the left and right images l ,c l ) and (r r ,c r The error equation described above can be rewritten as follows:

[0029]

[0030] Where r represents the row, c represents the column, and v rl v represents the row coordinates of the left image. cl v represents the column coordinates of the left image. rr v represents the row coordinates of the right image. cr Represents the column coordinates of the right image;

[0031] The above equation can then be simplified using symbols to represent V = A△ - l;

[0032] Step 3.3.4: Applying the least squares method to the above error equation, the least squares solution of Δ is obtained as follows:

[0033] △=[△X△Y△Z] T =(A T A) -1 A T l

[0034] Step 3.3.5: Determine whether the coordinate correction values ​​(△X, △Y, △Z) exceed the threshold of 1×10. －8 If so, then use the corrected geodetic coordinates (X). 1 ,Y 1 Z 1 Calculate the standardized coordinates of the left and right images and return to step 3.3.2 for iterative calculation; otherwise, complete the calculation, and the coordinates (X,Y,Z) at this time are the final geodetic coordinates.

[0035] Furthermore, step 4 includes the following specific steps:

[0036] Step 4.1: Perform point cloud preprocessing on the point cloud coordinates output in Step 3.3;

[0037] Step 4.2: Based on the preprocessed point cloud coordinates, perform DSM reconstruction using the point cloud reconstruction DSM algorithm;

[0038] Step 4.3: Use the CGAL library to post-process the DSM results generated in Step 4.2, and save the final DSM results in the graphical mesh model.

[0039] Furthermore, the point cloud preprocessing operation in step 4.1 includes:

[0040] 1) Deletion of mismatched points

[0041] Based on point cloud coordinates, the coordinate values ​​of mismatched points distributed outside a certain island are viewed after point cloud visualization using MATLAB. After obtaining the coordinate values, the coordinate points are deleted from the source file of the point cloud.

[0042] 2) Point cloud dilution

[0043] The point cloud after deleting mismatched points is diluted proportionally.

[0044] Furthermore, the specific steps of the point cloud reconstruction DSM algorithm in step 4.2 include:

[0045] Step 4.2.1: Read the point cloud coordinate data and create point set D1, edge set B1, triangle set T1 and baseline edge stack BaseStack. These sets contain three-dimensional coordinates and store information about edges and triangles.

[0046] Step 4.2.2: Select suitable coordinate points from point set D1 to create the convex closure of the triangulation, and store all edges on the convex closure into the baseline edge stack BaseStack; then construct point set D2 and store all points on the convex closure into it;

[0047] Step 4.2.3: Using the first edge in BaseStack as the baseline, find the expansion point that conforms to the Delaunay rule and perform LOP optimization, and construct the first triangle on the right; add the two newly generated edges to the edge set B1, add the newly generated triangle to the triangle set T1, and add these two new edges as new baseline edges to BaseStack;

[0048] Step 4.2.4: Using a point P in point set D2 i As a priority point, select the element containing P from BaseStack. i Using the edge as the baseline, construct a Delaunay triangulation on the right side and add the baseline edge to BaseStack; find the extension point that conforms to the Delaunay rule in the point set D1 and perform LOP optimization. After constructing the Delaunay triangle, add the two newly generated edges to the edge set B1 and the newly generated triangle to the triangle set T1. Add these two new edges as new baseline edges to BaseStack. If they are duplicates, they will not be added again.

[0049] Step 4.2.5: Check P i If a point is closed, remove it from point set D1 and remove the point containing P from BaseStack. i The edge;

[0050] Step 4.2.6: After creating the Delaunay triangle on the right, remove this baseline edge from the BaseStack and repeat steps 4.2.4 and 4.2.5;

[0051] Step 4.2.7: When BaseStack is empty, the DSM construction is complete.

[0052] Therefore, the present invention employs the above-mentioned method for extracting high-precision DSM from island and reef stereo images without measurement control points, which has the following beneficial effects:

[0053] First, this invention uses a rational function model for coordinate stereo positioning, obtains corresponding image points through disparity maps, constructs a set of error equations based on the initial coordinate values, uses the least squares method to obtain coordinate corrections, obtains the final geodetic coordinates through iteration, and processes the obtained latitude, longitude and altitude coordinates to obtain a three-dimensional point cloud without deformation.

[0054] Second, this invention uses TIN to represent DSM and employs an improved growth algorithm to generate TIN. Since the original growth algorithm is inefficient, this invention further optimizes the original growth algorithm through methods such as LOP optimization, redundant data removal, and convex closure construction, thereby improving the algorithm's efficiency. The DSM effect is further improved through point cloud preprocessing and DSM result post-processing.

[0055] Third, this invention utilizes an improved semi-global dense matching algorithm to obtain disparity maps, solving the problem of dense matching in cases of image blurring and noise interference. Finally, cost calculation optimization is performed, further improving the algorithm's running efficiency and matching accuracy.

[0056] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0057] Figure 1 This is a flowchart of the method proposed in this invention.

[0058] Figure 2 A detailed flowchart of DSM extraction based on satellite stereo image pairs.

[0059] Figure 3 The TIN triangulation network is constructed.

[0060] Figure 4 DSM algorithm for point cloud reconstruction.

[0061] Figures 5a-5b These are before-and-after comparison images of the DSM results after post-processing.

[0062] Figure 6 This is the visualization result of the point cloud.

[0063] Figure 7a and Figure 7b The DSM result diagrams generated by this invention and the DSM result diagrams generated by ENVI software are respectively;

[0064] Figure 8 Generate DEM results for Inpho software. Detailed Implementation

[0065] In the description of this invention, it should also be noted that, unless otherwise expressly specified and limited, these embodiments are for illustrative purposes only and are not intended to limit the scope of the invention. Furthermore, it should be understood that after reading the teachings of this invention, those skilled in the art can make various alterations or modifications to the invention, and these equivalent forms also fall within the scope defined by the appended claims.

[0066] This invention proposes a high-precision DSM extraction method for island and reef stereo image pairs without measurement control points, such as... Figure 1 As shown, it includes the following steps:

[0067] Step 1: Obtain stereo image pairs of distant islands and reefs from the World-View3 satellite;

[0068] Step 2: Perform semi-global dense matching using the improved semi-global dense matching algorithm to obtain the disparity map;

[0069] Step 3: Based on the obtained disparity map, use the rational function model for stereo localization to obtain the 3D point cloud;

[0070] Step 4: Preprocess the 3D point cloud obtained in Step 3, and use the 3D point cloud reconstruction DSM algorithm to generate the remote island and reef DSM. Then, postprocess the remote island and reef DSM to construct the 3D model of the DSM.

[0071] The key steps described above are explained below:

[0072] 1. Semi-global dense matching algorithm

[0073] Image dense matching algorithms are crucial for constructing a Digital Signal Module (DSM). Based on their optimization differences, dense matching algorithms can be categorized into locally optimal, globally optimal, and semi-global dense matching algorithms. Among these, semi-global dense matching algorithms currently offer higher matching efficiency and have proven effective in DSM reconstruction.

[0074] The semi-global dense matching algorithm was proposed primarily to address the excessive time and memory consumption of global matching. Existing technologies have improved upon the classic semi-global dense matching algorithm (SGM) by optimizing the calculation of matching costs, resulting in various methods based on Census transform. Additionally, there are object-side OSGM algorithms and tSGM algorithms employing a pyramid matching strategy, which have shown good performance in remote sensing image matching. In 2023, Ji Song et al. from the Information Engineering University proposed a semi-global constraint (Multi-View Vertical Line Locus, MVLL) matching method, which makes the matching results more reliable and improves image matching performance. To better apply the semi-global dense matching algorithm to DSM generation, Yang Xingbin et al. from Beijing University of Civil Engineering and Architecture proposed a DSM generation method based on improved semi-global matching in 2018, utilizing an image block strategy to constrain the initial disparity search range. The semi-global dense matching algorithm used in this invention is an improved version of the semi-global dense matching algorithm in the inventor's previous research result "A dense matching method for remote sensing imagesfused with CPS denoising". The original semi-global dense matching algorithm is based on mutual information (MI) and performs pixel-by-pixel matching. It adopts the idea of ​​optimizing the energy function, that is, finding the optimal disparity of each pixel to minimize the global energy function of the entire image. The entire matching process has four steps: initial cost calculation, cost aggregation, disparity calculation, and disparity optimization. From the DSM extraction process based on satellite stereo image pairs ( Figure 2 As can be seen from the diagram, the main steps of this improved semi-global dense matching algorithm are as follows: First, before matching, it is determined whether the peak signal-to-noise ratio (PSNR) is less than the threshold λ. If so, the CPS algorithm is used for image denoising. Second, the stereo image pair data is cropped and pixel coordinates are transformed. Third, based on the transformed image data, an epipolar image is constructed. Finally, the initial cost is calculated and the cost is aggregated based on the SGBM algorithm. Then, cost aggregation, disparity calculation, and disparity optimization are performed to finally output the disparity map.

[0075] Compared with the original semi-global dense matching algorithm, the improved algorithm includes the following improvements:

[0076] (1) Image cropping and pixel coordinate transformation

[0077] When performing semi-global dense matching on stereo image pairs, the original remote sensing images often contain excessive redundant regions. This significantly reduces the efficiency of image matching and consumes a large amount of computational space. Since the result of semi-global dense matching is a disparity map, which only contains information about the disparity between the left and right images, and the solution to disparity is related to the minimum solution of the global energy function, cropping the images greatly reduces the computational scale of the global energy function and improves the efficiency of disparity calculation. Therefore, to improve the efficiency of image matching, image cropping is performed to extract regions of interest.

[0078] After cropping an image, the pixel coordinates of the region of interest will change. If this is not addressed, it will lead to incorrect point cloud generation later. When cropping an image, it is necessary to ensure that the left and right images are the same size after cropping and that the same pixel coordinate transformation coefficients dR and dC are used. The relationship between the pixel coordinates (Rn', Cn') of the cropped image and the coordinates (Rn, Cn) before cropping is as follows:

[0079]

[0080] (2) Constructing epipolar images

[0081] Due to the parallel projection and orthogonal perspective projection along the direction of motion, linear motion causes the epipolar lines to become hyperbolic, while unavoidable nonlinear motion causes them to become general curves. Therefore, epipolar line correction is necessary before performing semi-global dense matching. When matching an image with corrected epipolar lines, only the corresponding pixels in each row need to be searched, thus reducing the search space from two dimensions to one dimension and significantly improving matching efficiency.

[0082] When constructing epipolar images, it is necessary to first establish an epipolar model. The traditional method is to divide the image into blocks, use affine transformation relationships to match image blocks, and then use the projection trajectory method to generate epipolar images. In addition, there is an epipolar image generation method that combines image-side and object-side methods. In order to minimize the vertical parallax of the epipolar images, this improved method adopts the epipolar model with the smallest vertical parallax to improve the accuracy and efficiency of semi-global dense matching.

[0083] The principle is as follows: for any pixel in the left image, let straight lines with different tilt angles pass through that point. Using the forward and inverse calculation model of a rational function model, calculate the corresponding image points of the points on the right image along these lines. When it is found that the distribution of corresponding image points of points on the right image along a straight line with a certain tilt angle is closest to a straight line, the left and right epipolar lines can be determined simultaneously. The specific steps for constructing the epipolar image are as follows:

[0084] Step 1: Let the angle of inclination of the line passing through pixel p on the left image be α. From this, we can obtain the equation of the left line as yy p =tanα(xx) p );

[0085] Step 2: Select n pixels at equal intervals on the left line, and determine their pixel coordinates (x, y, y). i ,y i ), i = 1, 2, ..., n and different geodetic heights Z i The corresponding geodetic coordinates (X) are calculated using the rational function inverse calculation model of the left image. i ,Y i );

[0086] Step 3: Calculate the geodetic coordinates (X, Y, F) of each point using the rational function forward model of the right image. i ,Y i Z i The coordinates (x, y) of the corresponding image point on the right image. i ',y i ');

[0087] Step 4: For (x) i ',y i By fitting a straight line and using the least squares method, the equation of the right-hand line is obtained as y'=kx'+b, with the coefficients of each term as follows:

[0088]

[0089] The right epipolar tilt angle α' = arctank', and the root mean square of the vertical disparity of each pixel is:

[0090]

[0091] Where Di is the root mean square of the vertical disparity of the i-th pixel, and

[0092]

[0093] Step 5: Record the root mean square of the vertical disparity σ obtained in the t-th test. If it is the current minimum value, then record σ. t =σ, α t =α;

[0094] Step 6: Within the range of epipolar dip angle values, repeat steps 1-5 until the test of each dip angle is completed;

[0095] Step 7: Calculate the left epipolar tilt angle α using the parabolic equation. min Then, by fitting the corresponding pixels of the left epipolar line on the right image, we obtain the right epipolar line f(α):

[0096] f(α) = u + v·α + w·α 2 (4)

[0097] Wherein, the parameters u, v, and w are represented by α t-1 α tα t+1 Solve the system of equations for the three corresponding root mean square equations of upper and lower parallax.

[0098] (3) Cost Calculation Optimization

[0099] When calculating costs based on mutual information, the calculation process is complex and requires a certain number of iterations, resulting in low efficiency. Currently, there are several improvement strategies for cost calculation. Among them, the SGBM (Semi-Global Block Matching) algorithm in the open-source computer vision library OpenCV is computationally simple and has high matching efficiency. Its BT (Block Truncation) strategy, starting from the information sampling perspective, can reduce errors caused by image sampling discretization. Furthermore, the SGBM algorithm adds a preprocessing step during cost calculation, using the Sobel operator for preprocessing, and then adding the BT cost of the preprocessed image to the BT cost of the original image as the initial cost. Therefore, this invention uses the SGBM algorithm for cost calculation, and its calculation formula is as follows:

[0100]

[0101] Where Sobel represents the horizontal Sobel operator, P represents the pixel value of each pixel in the image, and C(p,d) represents the cost of the current pixel p when the disparity is d. This represents the BT cost of the original image. This represents the BT cost of the preprocessed image.

[0102] (4) Image Denoising Based on Cauchy Proximal Split Algorithm

[0103] This improved algorithm employs the Cauchy Proximal Splitting (CPS) algorithm for image denoising, which boasts excellent convex optimization capabilities and good convergence performance, effectively addressing the L1 regularization problem. In image denoising, the CPS algorithm effectively constrains noise while preserving image details and texture information, exhibiting significant performance advantages. The main idea of ​​the CPS algorithm is to transform the image denoising problem into an equivalent convex optimization problem using the augmented Lagrangian method, and then iteratively calculates using proximal splitting algorithms (including Douglas Rachford splitting, forward-backward splitting, or alternating direction multiplier splitting). If the result converges, denoising is complete. This algorithm is applied to semi-global dense matching, primarily suitable for situations with low image quality and noise interference. It can improve image clarity and enhance object details, thereby improving matching accuracy. This improved algorithm adds a judgment condition before matching, using the Peak Signal-to-Noise Ratio (PSNR) to measure image quality; if the PSNR is less than a threshold, the CPS algorithm is required.

[0104] In summary, the improved algorithm adopted in this invention proposes an improved method that takes into account image quality based on the original SGM algorithm. It determines whether the CPS algorithm is needed based on image quality, and then further improves the algorithm by using image cropping and pixel coordinate transformation, constructing epipolar images, and cost calculation optimization. The improved algorithm has certain improvements in matching accuracy and algorithm efficiency, and performs well in image noise reduction. It can be well applied to dense matching of remote sensing images.

[0105] 2. Coordinate stereo positioning based on rational function model

[0106] (1) Rational function model

[0107] There are two main types of models for coordinate stereo positioning: strict geometric models and general geometric models. Traditional aerial photogrammetry uses strict geometric models for coordinate stereo positioning, which require the camera's interior and exterior orientation elements. However, commercial satellite sensor parameters are highly confidential, so satellite imagery typically uses general geometric models that are independent of sensor parameters. The rational function model is one such general geometric model, which calculates geodetic coordinates based on the RPC file (parameter file of the rational function model) provided by the satellite imagery. The rational function model is characterized by its independence from coordinate reference systems, independence from satellite sensors, and high confidentiality. It is an important imaging model for remote sensing satellite sensors, and it uses a set of polynomial ratios to represent the relationship between pixel coordinates and geodetic coordinates. The specific representation is as follows:

[0108]

[0109]

[0110]

[0111] In the formula, P i (X n ,Y n Z n The form is as follows:

[0112]

[0113] Where i = 1, 2, 3, 4; (r n ,c n (X) represents the standardized pixel coordinates. n ,Y n Z n ) represents the standardized geodetic coordinates; r o c o X o Y o Z o All are standardized translation parameters, r s c s X s Y s Z s All are standardized proportional parameters, a i The coefficients are known polynomial coefficients.

[0114] In the rational function model, the first-order term is used to represent optical projection error, the second-order term is used to represent errors caused by factors such as the curvature of the earth, atmospheric refraction, and lens distortion, and the third-order term is used to represent errors with higher-order components such as camera shake.

[0115] (2) Coordinate stereo positioning algorithm

[0116] ① The steps for calculating geodetic coordinates are as follows: linearize the rational function model of the two images (i.e., step S1) to obtain four sets of error equations. Combine the four sets of error equations and use the least squares method to obtain the coordinate correction. Based on the initial solution, the final geodetic coordinates are obtained through iteration.

[0117] Coordinate stereo positioning, also known as spatial forward intersection in aerial photogrammetry, calculates the corresponding geodetic coordinates using corresponding image points of stereo image pairs. The final geodetic coordinates (X, Y, Z) can be obtained through the above steps. The specific steps include:

[0118] S1: Calculate the initial value of the geodetic coordinates (X) 0 ,Y 0 Z 0): Calculate the average value of the standardized translation parameters for the left and right images respectively as initial values, and convert these two initial values ​​into standardized coordinates (X). nl ,Y nl Z nl ) and (X nr ,Y nr Z nr The specific formula is as follows:

[0119]

[0120] Where l and r are the left and right images, respectively, and X ol The initial left image in X-axis spatial rectangular coordinates, X or The initial right image in Cartesian coordinates along the X-axis, Y... ol The initial left image in Y-axis spatial rectangular coordinates, Y or The initial right image in Y-axis spatial rectangular coordinates, Z ol The initial left image of the Z-axis spatial rectangular coordinate system, Z or The initial right image of the Z-axis spatial rectangular coordinate system, X sl The left image is obtained after iterative calculation using X-axis Cartesian coordinates, and Y... sl The left image is obtained after iterative calculation using Y-axis Cartesian coordinates, Z... sl The left image after iterative calculation using Z-axis Cartesian coordinates, X sr The right image is obtained after iterative calculation using X-axis Cartesian coordinates, and Y... sr The right image is obtained after iterative calculation using Y-axis spatial rectangular coordinates, Z... sr The right image is the result of iterative calculation using Z-axis spatial rectangular coordinates;

[0121] S2: Establish a system of error equations:

[0122] shilling Then r n and c n Substituting the standardized formula into it and simplifying, we get:

[0123]

[0124] Expanding the above formula at (X,Y,Z) into linear terms using Taylor's formula yields:

[0125]

[0126] Therefore, the established error equation is:

[0127]

[0128] according to The formulas for calculating the partial derivatives are as follows:

[0129]

[0130]

[0131] Then, the coordinates of the corresponding image points (r) of the left and right images l ,c l ) and (r r ,c r The following error equation can be obtained:

[0132]

[0133] In rational function models, r typically represents a row, c typically represents a column, and v rl v represents the row coordinates of the left image. cl v represents the column coordinates of the left image. rr v represents the row coordinates of the right image. cr Represents the column coordinates of the right image;

[0134] Equation (16) can be simplified using symbols as V = A△-l;

[0135] S3: Solving for coordinate corrections: Using the least squares method on the error equation, the least squares solution for △ is:

[0136] △=[△X△Y△Z] T =(A T A) -1 A T l (17)

[0137] S4: Perform iterative solution: If the coordinate correction values ​​(△X, △Y, △Z) exceed the threshold, then use the corrected geodetic coordinates (X... 1 ,Y 1 Z 1 Calculate the standardized coordinates of the left and right images, and return to step ② for iterative calculation; otherwise, complete the calculation, and the coordinates (X,Y,Z) at this time are the final geodetic coordinates.

[0138] The resulting coordinates (X, Y, Z) represent latitude, longitude, and elevation, respectively. Experiments show that the correction values ​​for X and Y are generally small, while the correction value for Z is large. Therefore, the threshold is set to 1 × 10⁻⁶. －8 When the time is right, the impact of subsequent iterations on the coordinate results can be ignored, and high-precision coordinate three-dimensional positioning can be achieved.

[0139] ② Obtain the coordinates of the point cloud from the disparity map;

[0140] Traditional methods for converting disparity maps to point clouds utilize disparity to calculate depth using triangulation, and then obtain the point cloud. This invention, based on a rational function model, can directly derive the point cloud coordinates from the disparity map. However, since the obtained coordinates are latitude, longitude, and altitude geodetic coordinates, the shape of the point cloud on the xoy plane will change during visualization. Furthermore, coordinate transformation of the point cloud is unreliable and distortion still occurs. Therefore, the latitude and longitude coordinates of the point cloud can be replaced with the pixel coordinates of the left image. This method avoids distortion, and the specific operation steps are as follows:

[0141] ① Based on the disparity value d corresponding to each pixel in the disparity map, the formula r r =r l c r =c l -d(r represents the right image, l represents the left image) calculates the coordinates of the corresponding image points of the left and right images;

[0142] ②The coordinates of the corresponding image points before image cropping are obtained according to formula (1);

[0143] ③ Based on the coordinates of the corresponding image points, the corresponding geodetic coordinate points are obtained using the above-mentioned geodetic coordinate calculation steps;

[0144] ④ Replace the latitude and longitude coordinates of the geodetic coordinate points with the pixel coordinates corresponding to the left image;

[0145] ⑤ Output the final point cloud coordinates;

[0146] In summary, the rational function model is a general geometric model with good confidentiality. It takes the form of a polynomial ratio, and linearizing it yields the error equation. By simultaneously solving the error equations of the left and right images and using the least squares method to solve for the coordinate correction, the geodetic coordinates can be obtained through iteration. Although this invention calculates the geodetic coordinates using the rational function model, the elevation error is still relatively large. Since the magnitude of the elevation error plays a decisive role in the accuracy of the point cloud, the positioning error of the rational function model is corrected, thereby reducing the coordinate positioning error.

[0147] 3. 3D Point Cloud Reconstruction DSM Algorithm

[0148] (1) Point cloud preprocessing

[0149] This invention performs point cloud preprocessing before reconstructing the 3D point cloud model (DSM). Point cloud preprocessing mainly involves two steps: removing mismatched points and point cloud dilution. For mismatched points in the point cloud, since they are mainly distributed outside a certain island and are relatively few in number, manual removal can be used. The principle of manual removal is quite simple: after visualizing the point cloud using MATLAB, the coordinate values ​​of these mismatched points can be viewed. After obtaining the coordinate values, the points can be directly deleted from the source file of the point cloud. Point cloud dilution is performed because the point cloud obtained from dense matching is too dense. Compared to real-world objects, most of the points are redundant, which will affect the program's performance and the final quality of the generated DSM. Therefore, the point cloud is diluted proportionally before generating the DSM, and the dilution ratio is adjusted according to specific circumstances.

[0150] (2) Principle of DSM Algorithm for Point Cloud Reconstruction

[0151] There are two main representation models for Digital Surface Models (DSMs): surface fitting models and Triangulated Irregular Networks (TINs). TINs offer two major advantages: ① less data redundancy and faster processing speed; ② higher computational efficiency, making them well-suited for high-precision modeling. This invention uses TINs to represent DSMs. While there are many methods for generating TINs, Delaunay triangulation is currently the most common.

[0152] There are three main algorithms for generating Delaunay triangulation networks: ① Divide and conquer algorithm: The basic idea is to decompose the problem, reducing the difficulty of generating the triangulation network. A large set of points is divided into smaller sets, and then the triangulation networks generated by these smaller sets are merged into the final triangulation network. Its advantage is high time efficiency, but the computation process requires a large number of recursive operations, consuming a significant amount of memory. ② Point-by-point insertion algorithm: The basic idea is to construct a convex hull containing all points, build an initial triangulation network within the hull, and then insert the remaining points one by one into the triangulation network. This algorithm is easy to implement but has low efficiency. ③ Triangulation growth algorithm: The basic idea of ​​the triangulation growth algorithm is to randomly select a point from the discrete elevation points, search for the nearest point to this point, and connect the two points as the initial baseline. Based on the empty circle property and the maximum and minimum angle characteristics of the Delaunay triangulation network, find the extension points that form a Delaunay triangle with the baseline; generate the Delaunay triangle, and then use the two newly generated sides of the triangle as the new baseline; repeat the search for extension points and the generation of Delaunay triangles until no extension points can be found.

[0153] The point cloud TIN generation algorithm of this invention employs a growth algorithm, and further optimizes the original triangular mesh growth algorithm through LOP optimization, redundant data elimination, and convex closure construction, making TIN generation more efficient. Its specific steps include:

[0154] ①LOP optimization

[0155] The Local Optimization Procedure (LOP) checks the LOP values ​​of a triangulated network. Four points form two triangles sharing a common edge. If the circumcircle of one triangle contains the fourth point, then the two diagonals are swapped. This algorithm can satisfy the Delaunay triangle assumption. The LOP optimization algorithm can use the maximum and minimum angle properties to determine whether swapping the diagonals BD and AC of a convex quadrilateral ABCD formed by two adjacent triangles sharing a common edge increases the angle values ​​of the six interior angles of these two triangles. This is done using common properties of inscribed angles of triangles. If:

[0156] If sin(∠A+∠C)=0, then point A lies on the circumcircle of triangle BCD, and diagonal swapping is not necessary.

[0157] If sin(∠A+∠C)>0, then point A is outside the circumcircle of triangle BCD, and diagonal swapping is not necessary.

[0158] If sin(∠A+∠C)<0, then point A is inside the circumcircle of triangle BCD, and the diagonals need to be swapped.

[0159] ② Remove redundant data

[0160] The Delaunay triangulation construction method centered on priority vertices is a further optimization of the method excluding closed vertices. The definition of a priority vertex: If, during the construction of the triangulation, vertex P... i If it is not a closed point, prioritize expanding the area containing P. i The edge until P i If P becomes a closed point, then i These are called priority points.

[0161] This invention employs a triangulation centered on a priority point and promptly removes closed points from the data point set. This reduces the number of traversals required to find triangle extension points, thus decreasing computation time. Whether a point is closed can be determined by checking if the sum of the included angles of all triangles formed by that point in the triangulation at that point equals 360° (except for points on the convex hull of the triangulation). Constructing the triangulation centered on the priority point involves prioritizing edges containing that priority point when selecting baselines. Once all edges containing that point in the edge set have participated in the triangulation construction, the priority point becomes a closed point and can then be removed from the point set, reducing the number of points that need to be traversed in subsequent triangulation construction. The removal of closed points does not affect subsequent triangulation construction and further reduces the number of data points in the data set, improving the efficiency of subsequent triangulation construction.

[0162] ③ Construct convex closure

[0163] The triangulation growth algorithm constructs a triangulation network by ordering the coordinates of the point set, requiring traversal of each point in the set to find expansion points. Using a convex closure, the original point set is divided into two point set regions, and the original baseline edge set is also divided into two edge set regions. Then, a priority point elimination method is used to exclude point sets and edge sets from different regions using different methods, which simplifies the algorithm.

[0164] When constructing a Delaunay triangulation, it's necessary to traverse the edge set storing the baselines to find the next baseline to construct. Each edge can serve as one side of at most two triangles, while an edge on a convex closure can only serve as a baseline once in the triangulation. Therefore, after determining that an edge has been used as a baseline twice, that edge can be deleted, thereby reducing the number of edges in the edge set and improving traversal efficiency. For example, in... Figure 3 In a triangulated network, an edge AB on the convex hull can only serve as a primary baseline of triangle ABC, while an edge BC inside the convex hull can serve as a baseline for both triangle ABC and triangle BCD. Therefore, edge AB can be removed from the edge set after it has been detected as a primary baseline, and edge BC can be removed from the edge set after it has been detected as a secondary baseline. This reduces the amount of data that needs to be processed subsequently.

[0165] (3) Point Cloud Reconstruction DSM Algorithm

[0166] Based on the above optimizations, the point cloud reconstruction DSM algorithm proposed in this invention, such as... Figure 4 As shown, the specific steps include the following:

[0167] Step 4.1: Read the point cloud coordinate data and create point set D1, edge set B1, triangle set T1 and baseline edge stack BaseStack. These sets contain three-dimensional coordinates and store information about edges and triangles.

[0168] Step 4.2: Select suitable coordinate points from point set D1 to create the convex closure of the triangulation, and store all edges on the convex closure into the baseline edge stack BaseStack; then construct point set D2 and store all points on the convex closure into it;

[0169] Specifically, when selecting suitable coordinate points, the origin set D1 is divided into two point set regions using convex closure, and the original baseline edge set B1 is divided into two edge set regions. Then, the method of excluding priority points is used to exclude each coordinate point in different point set regions and edge set regions to obtain the selected suitable coordinate points.

[0170] Step 4.3: Using the first edge in BaseStack as the baseline, find the expansion point that conforms to the Delaunay rule and perform LOP optimization, and construct the first triangle on the right side; add the two newly generated edges other than the first edge in BaseStack to the edge set B1, add the newly generated triangle to the triangle set T1, and add these two new edges as new baseline edges to BaseStack.

[0171] Step 4.4, take a point P in point set D2 i As a priority point, select the element containing P from BaseStack. i Using the edge as the baseline, construct a Delaunay triangulation on the right side and add the baseline edge to BaseStack; find the extension point that conforms to the Delaunay rule in the point set D1 and perform LOP optimization. After constructing the Delaunay triangle, add the two newly generated edges to the edge set B1 and the newly generated triangle to the triangle set T1. Add these two new edges as new baseline edges to BaseStack. If they are duplicates, they will not be added again.

[0172] Step 4.5, check P i If a point is closed, remove it from point set D1 and remove the point containing P from BaseStack. i The edge;

[0173] Step 4.6: After creating the Delaunay triangle on the right, remove this baseline edge from the BaseStack, and repeat steps 4.4 and 4.5.

[0174] Step 4.7: When BaseStack is empty, the DSM construction is complete;

[0175] (4) Post-processing of DSM results

[0176] After generating the initial TIN triangulation using the DSM reconstruction algorithm based on point clouds, the results in some areas are not ideal, containing many holes and incorrect triangular faces, thus requiring post-processing. First, the triangular faces in the TIN triangulation are processed by filtering out excessively large faces. Then, hole filling is performed. Hole identification is required during hole filling, and all holes except the outer shell are filled. The faces are then refined and smoothed to generate a better-shaped DSM. This processing flow uses the CGAL (Computational Geometry Algorithms Library), a third-party library widely used in computer graphics, geographic information systems, and other fields. It can save the final DSM result in a graphical mesh model (DSM 3D model) and generate a .ply file, facilitating the transfer and visualization of DSM results across platforms. A comparison of the post-processed DSM results is shown below. Figures 5a-5b As shown.

[0177] In summary, this invention uses TIN to represent DSM, and the algorithm for generating TIN from point clouds employs a growth algorithm. The original growth algorithm is inefficient; therefore, this invention optimizes the original growth algorithm through methods such as LOP optimization, redundant data removal, and convex closure construction, thereby improving the algorithm's efficiency. Before and after DSM reconstruction from point clouds, two additional steps—point cloud preprocessing and DSM result post-processing—are added, respectively, further improving the DSM's effectiveness.

[0178] Example

[0179] 1. Analysis of Point Cloud Generation Results

[0180] (1) Visualization and analysis of point cloud results

[0181] To verify the effectiveness of the point cloud generation method based on the rational number model in step 3, the point cloud results were visualized and analyzed. The final number of coordinate points generated from the disparity map was 16,250,442. MATLAB was used for visualization, and the final result is shown below. Figure 6 As shown, the point cloud visualization results reveal the approximate outline and elevation of the island. However, some mismatched points and sparse point clouds in certain areas are still visible, which will affect the final generated DSM result. Therefore, optimization measures are needed in subsequent processes. Furthermore, if the island is considered as a plane, the entire plane is tilted. While this does not affect the final visualization of the DSM 3D model, it will impact accurate positioning. The reason for this is that when using the least squares method to solve for coordinate corrections, inverting the matrix causes oscillations in the solution value.

[0182] (2) Point cloud result accuracy analysis

[0183] Accuracy analysis was performed on the generated point cloud, mainly analyzing its mean, standard deviation, and error percentage. The results are shown in Table 1. The table shows that the standard deviations for latitude and accuracy are relatively small, and their error percentages are also small, while the standard deviation for elevation is relatively large, and its error percentage is very high. Combined with the visualization results, it can be concluded that in the point cloud computing process, the magnitude of the elevation error plays a decisive role in the accuracy of the point cloud.

[0184] Table 1 Point Cloud Accuracy Analysis Table

[0185]

[0186] 2. Comparison and Accuracy Analysis of DSM Results

[0187] (1) Comparative analysis of DSM results

[0188] Generating a Directional Smart Point (DSM) using stereo image pairs is achievable in many commercial software programs with high accuracy. However, unlike the process described in this paper, commercial software requires aerial triangulation before DSM generation. Aerial triangulation involves first extracting tie points, then performing area network adjustment, and outputting relative orientation results. If ground control points are available, they can be used for adjustment to output absolute orientation results. With ground control points, both positioning and DSM accuracy are significantly improved. However, due to measurement limitations, the inability to obtain ground control points reduces the effectiveness of the DSM. A comparison of the DSM generated in this paper with the commercial software ENVI shows the following results: Figure 7a and Figure 7b As shown:

[0189] Depend on Figure 7a and Figure 7b It is evident that the DSM generated by ENVI describes the elevation undulations of ground features more accurately, but the elevation error of the sea surface is relatively large, with the elevation of most sea surface areas being higher than that of a certain island region, which is clearly incorrect. Although the DSM generated in this paper does not have the same visualization effect and description of ground feature elevation as ENVI, the overall DSM is smoother, with smaller elevation changes between ground features, which conforms to the actual topography of the island. Furthermore, the sea surface area is less affected by noise points, thus avoiding the influence of sea surface elevation errors on the DSM results.

[0190] (2) Accuracy analysis of DSM results

[0191] To analyze the accuracy of the DSM results, the DSM generated by the algorithm in this paper and the DSM generated by ENVI were compared with the DEM data generated by the professional software Inpho (e.g., ...). Figure 8The elevation errors are compared to those shown in the figure. The DEM data shows an average ground elevation of 26m, while the DSM generated by the algorithm in this paper has an average ground elevation of 37.9m, an error of +11.9m. The DSM generated by ENVI has an average ground elevation of 29.5m, an error of +3.5m. The elevation errors indicate that the DSM generated by ENVI is more accurate, but overall, the accuracy of both DSM results needs improvement.

[0192] The DSM error generated by the algorithm of this invention mainly comes from elevation calculation error, which is caused by four main reasons: ① the parameters in the RPC file are not accurate; ② the rational function model has certain systematic errors; ③ the geodetic coordinate solution method is not accurate; ④ the interference of noise points during the dense matching of images.

[0193] In summary, by comparing the DSM generated by this method with that generated by ENVI, it can be seen that this method provides a more realistic description of the surface of remote islands and reefs, and also reduces the impact of sea surface elevation errors on the DSM results.

[0194] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for extracting high-precision DSM from stereo image pairs of islands and reefs without measurement control points, characterized in that, Includes the following steps: Step 1: Acquire stereo image pairs of distant islands and reefs from World-View3 satellite; Step 2: Based on the stereo image pairs of distant islands and reefs, use the improved semi-global dense matching algorithm to perform semi-global dense matching to obtain the disparity map; Step 3: Based on the obtained disparity map, use a rational function model to perform stereo localization and obtain a 3D point cloud; Step 4: Preprocess the 3D point cloud obtained in Step 3, and use the 3D point cloud reconstruction DSM algorithm to generate remote island and reef DSM. Then, postprocess the generated remote island and reef DSM to construct a 3D model of the remote island and reef DSM. Step 3 includes the following specific steps: Step 3.1: Based on the disparity value d corresponding to each pixel in the disparity map, obtain the coordinates of the corresponding image points in the left and right images using the following formula: r r =r l ,c r =c l -d ; Among them, (r l ,c l (r) represents the coordinates of the corresponding pixel in the left image. r ,c r ) represents the coordinates of the corresponding pixel in the right image; Step 3.2: Obtain the coordinates of the corresponding image points before cropping the left and right images using the following formula: ； Where dR and dC are coordinate transformation coefficients using the same pixels; The coordinates of the image before cropping. The coordinates are those of the image after cropping; Step 3.3: Based on the rational function model, calculate the corresponding geodetic coordinate points according to the coordinates of the corresponding image points; Step 3.4: Replace the latitude and longitude coordinates of the geodetic points with the pixel coordinates corresponding to the left image; Step 3.5: Output the final point cloud coordinates; Furthermore, the specific steps of step 3.3 include: Step 3.3.1: Calculate the average value of the standardized translation parameters for the left and right images respectively, use them as initial values, and convert these two initial values ​​into standardized coordinates. and The initial values ​​of the geodetic coordinates are obtained. The specific calculation formula is as follows: ； in, The initial left image is the X-axis spatial rectangular coordinate. The initial right image is the rectangular coordinate system of the X-axis. The initial left image is the Y-axis spatial rectangular coordinate. The initial right image is the Y-axis spatial rectangular coordinate. The initial left image is the Z-axis spatial rectangular coordinate. The initial right image is the Z-axis spatial rectangular coordinate. The left image is the result of iterative calculation using X-axis spatial rectangular coordinates. The left image is the result of iterative calculation using Y-axis spatial rectangular coordinates. The left image is the result of iterative calculation using Z-axis Cartesian coordinates. The right image is the result of iterative calculation using X-axis spatial rectangular coordinates. The right image is the result of iterative calculation using Y-axis spatial rectangular coordinates. The right image is the result of iterative calculation using Z-axis spatial rectangular coordinates; Step 3.3.2: Based on the standardized coordinates of the left and right images, construct the following error equation: ； Step 3.3.3: Based on the coordinates of corresponding image points (r) in the left and right images l ,c l ) and (r r ,c r The error equation described above can be rewritten as follows: ； Among them, v rl v represents the row coordinates of the left image. cl v represents the column coordinates of the left image. rr v represents the row coordinates of the right image. cr Represents the column coordinates of the right image; The above formula can then be simplified using symbols as follows: ; Step 3.3.4: Apply the least squares method to the above error equation to obtain... The least squares solution is: ； Step 3.3.5: Determine the coordinate correction value Does it exceed the threshold of 1×10? －8 If so, then use the corrected geodetic coordinates. Calculate the standardized coordinates of the left and right images, and return to step 3.3.2 for iterative calculation; otherwise, complete the calculation, and the coordinates at this point are... This is the final geodetic coordinate.

2. The method for extracting high-precision DSM from island and reef stereo image pairs without measurement control points as described in claim 1, characterized in that, Step 4 includes the following specific steps: Step 4.1: Perform point cloud preprocessing on the point cloud coordinates output in Step 3.3; Step 4.2: Based on the preprocessed point cloud coordinates, perform DSM reconstruction using the point cloud reconstruction DSM algorithm; Step 4.3: Use the CGAL library to post-process the DSM results generated in Step 4.2, and save the final DSM results in the graphical mesh model.

3. The method for extracting high-precision DSM from island and reef stereo image pairs without measurement control points as described in claim 2, characterized in that, Step 4.1, the point cloud preprocessing operations, include: 1) Deletion of mismatched points Based on point cloud coordinates, point cloud visualization was performed using MATLAB to examine the distribution of errors outside a certain island. Match the coordinates of the point, obtain the coordinates, and then delete the point from the source file of the point cloud. 2) Point cloud dilution The point cloud after deleting mismatched points is diluted proportionally.

4. The method for extracting high-precision DSM from stereo image pairs of islands and reefs without measurement control points as described in claim 3, characterized in that, The specific steps of the DSM algorithm for point cloud reconstruction in step 4.2 include: Step 4.2.1: Read the point cloud coordinate data and create point set D1, edge set B1, triangle set T1 and baseline edge stack BaseStack. These sets contain three-dimensional coordinates and store information about edges and triangles. Step 4.2.2: Select suitable coordinate points from point set D1 to create the convex closure of the triangulation, and store all edges on the convex closure into the baseline edge stack BaseStack; then construct point set D2 and store all points on the convex closure into it; Step 4.2.3: Using the first edge in BaseStack as the baseline, find the expansion point that conforms to the Delaunay rule and perform LOP optimization, and construct the first triangle on the right; add the two newly generated edges to the edge set B1, add the newly generated triangle to the triangle set T1, and add these two new edges as new baseline edges to BaseStack; Step 4.2.4: Using a point P in point set D2 i As a priority point, select the element containing P from BaseStack. i Using the edge as the baseline, construct a Delaunay triangulation on the right side and add the baseline edge to BaseStack; find the extension point that conforms to the Delaunay rule in the point set D1 and perform LOP optimization. After constructing the Delaunay triangle, add the two newly generated edges to the edge set B1 and the newly generated triangle to the triangle set T1. Add these two new edges as new baseline edges to BaseStack. If they are duplicates, they will not be added again. Step 4.2.5: Check P i If a point is closed, remove it from point set D1 and remove the point containing P from BaseStack. i The edge; Step 4.2.6: After creating the Delaunay triangle on the right, remove this baseline edge from the BaseStack and repeat steps 4.2.4 and 4.2.5; Step 4.2.7: When BaseStack is empty, the DSM construction is complete.

Images