High-precision DSM method for extracting stereoscopic image pair of island without measurement control point

Abstract

The invention provides a measurement control point-free island stereo image pair extraction high-precision DSM method, which comprises the following steps of: firstly, carrying out semi-global dense matching by adopting a semi-global dense matching algorithm to obtain a disparity map; secondly, performing coordinate three-dimensional positioning of the rational function model, and processing longitude, latitude and height coordinates to obtain a three-dimensional point cloud without deformation; and finally, realizing point cloud reconstruction DSM, and generating a DSM three-dimensional model. Experimental results show that the method for reconstructing the DSM by the high-resolution satellite stereo image pair is well applied to World-View3 satellite images, and the DSM three-dimensional model of the island reef can be constructed.

Description

Technical Field

[0001] The present invention relates to the field of geographic information systems and map making technology, and in particular to a method for extracting high-precision DSM from a stereo image pair of islands and reefs without measurement control points. Background Art

[0002] A digital surface model (DSM) is a digital model that describes the elevation of surface features, providing critical geographic information data for fields such as urban planning, resource management, and environmental monitoring. However, traditional aerial photogrammetry methods are currently not suitable for producing 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 achieves a full-color resolution of 0.3 meters. Its ultra-high resolution and multi-band data make it widely used in various fields. With the continuous development of my country's stereo mapping satellites, many constructive achievements have been made. In 2012, my country successfully launched the Ziyuan-3 satellite, its first independently developed civilian high-resolution stereo mapping satellite, and it has been widely used in the production of 5-meter resolution DSMs / DEMs. In addition to Ziyuan-3 01, 02, and 03 satellites, my country also has the Gaofen-7 satellite, launched in 2019, capable of producing 2-meter resolution DSMs / DEMs, and the Siwei Gaojing-1 and 2 satellites, scheduled for launch in 2022. This abundant satellite data source provides ample room for DSM data production, further strengthening the assurance of surveying and mapping products. Therefore, with the continuous development of high-resolution stereo mapping satellites, the use of high-resolution remote sensing imagery to produce surveying and 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 the three-dimensional point cloud, and how to use the three-dimensional point cloud to generate DSM are all hot and key issues that need to be solved urgently. Summary of the Invention

[0004] In view of the shortcomings of the above-mentioned prior art, the present invention proposes a method for extracting high-precision DSM from stereo image pairs of islands and reefs without measurement control points, which is used to generate DSM based on World-View3 satellite stereo image pairs.

[0005] To this end, the specific technical solutions adopted in the present invention are as follows:

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

[0007] Step 1: Obtain stereo images of remote islands and reefs from the World-View3 satellite;

[0008] Step 2: Based on the stereo image pairs of remote islands and reefs, semi-global dense matching is performed using an improved semi-global dense matching algorithm to obtain a disparity map;

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

[0010] Step 4: Pre-process the 3D point cloud obtained in step 3, and use the 3D point cloud reconstruction DSM algorithm to generate the remote island reef DSM, and then post-process the generated remote island reef DSM to construct a remote island reef DSM 3D model.

[0011] Furthermore, the specific steps of step 3 include:

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

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

[0014] Among them, r represents the right image and l represents the left image;

[0015] Step 3.2: Obtain the coordinates of the image points with the same name before cropping the left and right images according to the following formula:

[0016] Rn'=Rn+dR

[0017] Cn'=Cn+dC

[0018] Where dR and dC are the coordinate conversion 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 image points with the same name;

[0020] Step 3.4: Replace the longitude and latitude coordinates of the geodetic coordinate point with the corresponding pixel coordinates of 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 normalized translation parameters of the left and right images respectively, use it as the initial value, and convert these two initial values into normalized coordinates (X nl ,Y nl ,Z nl ) and (Xnr ,Y nr ,Z nr ), get the initial value of the geodetic coordinate (X 0 ,Y 0 ,Z 0 ), the specific calculation formula is:

[0024]

[0025] Among them, l and r are the left image and the right image respectively, X ol is the initial left image of the X-axis space rectangular coordinate, X or is the initial right image of the X-axis space rectangular coordinate, Y ol is the initial left image of the Y-axis spatial rectangular coordinate, Y or is the initial right image of the Y-axis space rectangular coordinate, Z ol is the initial left image of the Z-axis rectangular coordinate, Z or is the initial right image of the Z-axis rectangular coordinate, X sl is the left image after the iterative operation of the X-axis space rectangular coordinate, Y sl is the left image after the iterative operation of the Y-axis spatial rectangular coordinate, Z sl is the left image after the iterative operation of the rectangular coordinates of the Z axis space, X sr is the right image after the iterative operation of the X-axis space rectangular coordinate, Y sr is the right image after the iterative operation of the Y-axis space rectangular coordinate, Z sr It is the right image after the iterative calculation of the rectangular coordinates of the Z-axis space;

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

[0027]

[0028] Step 3.3.3: Based on the coordinates of the same-name image points in the left and right images (r l ,c l ) and (r r ,c r ), rewrite the above error equation as:

[0029]

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

[0031] Then the above formula is simplified into V = A△-l;

[0032] Step 3.3.4: Apply the least squares method to the above error equation and obtain the least squares solution of △ as:

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

[0034] Step 3.3.5: Determine whether the coordinate correction value (△X, △Y, △Z) exceeds the threshold value 1×10 －8 If so, use the corrected geodetic coordinates (X 1 ,Y 1 ,Z 1 ) calculates the standardized coordinates of the left and right images and returns to step 3.3.2 for iterative calculation; otherwise, the calculation is completed and the coordinates (X, Y, Z) at this time are the final geodetic coordinates.

[0035] Furthermore, the specific steps of step 4 include:

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

[0037] Step 4.2: Based on the pre-processed point cloud coordinates, DSM reconstruction is performed 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 a graphical mesh model.

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

[0040] 1) Deletion of incorrect matching points

[0041] Based on the point cloud coordinates, MATLAB is used to visualize the point cloud and check the coordinate values of the mismatched points distributed outside a certain island. 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 the 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 the point set D1, edge set B1, triangle set T1, and baseline edge stack BaseStack. These sets contain 3D coordinates and information about edges and triangles.

[0046] Step 4.2.2: Select appropriate coordinate points from point set D1 to create the convex closure of the triangulated network and store all edges on the convex closure in the baseline edge stack BaseStack. Then construct point set D2 and store all points on the convex closure in it.

[0047] Step 4.2.3: Using the first edge in BaseStack as the baseline, find an extension point that complies with the Delaunay theorem for LOP optimization and construct the first triangle on the right. Add the two newly generated edges to edge set B1 and the newly generated triangle to triangle set T1. Add these two new edges as new baseline edges to BaseStack.

[0048] Step 4.2.4: Take a point P in the point set D2 i As the priority point, select the one containing P from BaseStack i The edge of the triangle is used as the baseline, and the Delaunay triangulation is constructed on the right side. The baseline edge is added to the BaseStack. In the point set D1, an extension point that conforms to the Delaunay law is found for LOP optimization. After the Delaunay triangle is constructed, the two newly generated edges are added to the edge set B1 and the newly generated triangle is added to the triangle set T1. These two new edges are added as new baseline edges to the BaseStack. If there are duplicates, they are not added again.

[0049] Step 4.2.5: Check P i Is it a closed point? If so, delete it from the point set D1 and delete the point containing P from BaseStack. i edge;

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

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

[0052] Therefore, the present invention adopts 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, the present invention uses a rational function model for coordinate stereo positioning, obtains homonymous image points through a disparity map, constructs an error equation system based on the initial coordinate values, uses the least squares method to obtain the coordinate correction, and obtains the final geodetic coordinates through iteration. The obtained longitude and latitude coordinates are processed to obtain an undeformed three-dimensional point cloud;

[0054] Second, the present invention uses a TIN to represent the DSM and employs an improved growing algorithm to generate the TIN. Because the original growing algorithm is inefficient, this invention also improves its efficiency through methods such as LOP optimization, redundant data elimination, and convex closure construction. The DSM's effectiveness is further enhanced through point cloud preprocessing and DSM postprocessing.

[0055] Third, the present invention uses an improved semi-global dense matching algorithm to obtain the disparity map, which solves the dense matching problem in cases of image blur and noise interference. Finally, it performs cost calculation optimization to further improve the algorithm operation efficiency and matching accuracy.

[0056] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 1 The flowchart of the method proposed in the present invention is shown in FIG.

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

[0059] Figure 3 The constructed TIN triangulation network.

[0060] Figure 4 Reconstruct DSM algorithm for point clouds.

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

[0062] Figure 6 Visualize the results for the point cloud.

[0063] Figure 7a and Figure 7b Generate DSM result graphs for the present invention and ENVI software respectively;

[0064] Figure 8 Generate DEM result map for Inpho software. DETAILED DESCRIPTION

[0065] In the description of the present invention, it should be noted that, unless otherwise clearly specified and limited, these embodiments are only used to illustrate the present invention and are not intended to limit the scope of the present invention. In addition, it should be understood that after reading the content taught by the present invention, those skilled in the art will make various changes or modifications to the present invention, and these equivalent forms also fall within the scope defined by the appended claims of the application.

[0066] This paper proposes a method for extracting high-precision DSM from stereo image pairs of islands and reefs without measurement control points. Figure 1 As shown, it includes the following steps:

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

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

[0069] Step 3: Based on the obtained disparity map, use the rational function model for stereo positioning to obtain a three-dimensional point cloud;

[0070] Step 4: pre-process the three-dimensional point cloud obtained in step 3, and use the three-dimensional point cloud reconstruction DSM algorithm to generate the remote island reef DSM, and post-process the remote island reef DSM to construct a DSM three-dimensional model.

[0071] The following are the key steps described above:

[0072] 1. Semi-global dense matching algorithm

[0073] Image dense matching algorithms are crucial for constructing DSMs. Based on their optimization differences, these algorithms are categorized into locally optimal, globally optimal, and semi-global. The semi-global dense matching algorithm is currently the most efficient and has proven successful in DSM reconstruction.

[0074] The semi-global dense matching algorithm was proposed mainly to solve the problem of excessive time and memory consumption of global matching. The classic semi-global dense matching algorithm (Semi-Global Matching, SGM) proposed in the existing technology has been improved for the calculation of matching costs. A variety of methods based on Census transformation have emerged. In addition, there is the OSGM algorithm based on object space and the tSGM algorithm using a pyramid matching strategy, which are applied to remote sensing image matching with good results. In 2023, Ji Song et al. from the University of Information Engineering proposed a semi-global constrained (Multi-View Vertical Line Locus, MVLL) matching method to make the matching results more reliable and improve the image matching performance. In order to better apply the semi-global dense matching algorithm to generate DSM, Yang Xingbin et al. from Beijing University of Civil Engineering proposed a DSM generation method based on improved semi-global matching in 2018, using the image blocking strategy to constrain the initial disparity search range. The semi-global dense matching algorithm used in the present invention is an improved semi-global dense matching algorithm in the inventor's previous research results "A dense matching method for remote sensing images fused with CPS denoising". The original semi-global dense matching algorithm is based on mutual information (MI) pixel by pixel matching, and 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 ), the main steps of the improved semi-global dense matching algorithm are: first, before matching, it is determined whether the peak signal-to-noise ratio (PSNR) index is less than the threshold λ. If so, the CPS algorithm is used for image denoising; second, the stereo image data is cropped and pixel coordinates are converted; third, based on the converted image data, the epipolar image is constructed; finally, the initial cost calculation and cost aggregation are performed based on the SGBM algorithm, and then cost aggregation, disparity calculation and disparity optimization are performed, and finally the disparity map is output;

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

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

[0077] When performing semi-global dense matching on stereo image pair data, the original remote sensing images often contain excessive redundant regions, significantly reducing image matching efficiency and consuming significant computational space. Since the result of semi-global dense matching is a disparity map, which contains only the disparity between the left and right images, and disparity solution is dependent on minimizing the global energy function, cropping the images significantly reduces the computational scale of the global energy function and improves the efficiency of disparity solution. Therefore, to improve image matching efficiency, cropping the images is performed to extract the region of interest.

[0078] After cropping an image, the pixel coordinates of the region of interest will change. If not processed, it will result in an incorrect point cloud being generated later. When cropping an image, it is necessary to ensure that the left and right images are the same size after cropping, and use the same pixel coordinate conversion coefficients dR and dC. The relationship between the pixel coordinates of the cropped image (Rn', Cn') and the coordinates before cropping (Rn, Cn) is:

[0079]

[0080] (2) Constructing epipolar images

[0081] Due to parallel projection in the direction of motion and orthogonal perspective projection, linear motion causes the epipolar line to become a hyperbola, while inevitable nonlinear motion causes the epipolar line to become a normal curve. Therefore, before performing semi-global dense matching, epipolar correction is required on the image. When matching an image with epipolar correction, only the corresponding pixel points in each row need to be searched, reducing the search space from two dimensions to one, significantly improving matching efficiency.

[0082] When constructing an epipolar image, 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 the image blocks, and then use the projection trajectory method to generate the epipolar image. In addition, there is also an epipolar image generation method that combines image and object space. In order to minimize the upper and lower parallax of the epipolar image, this improved method uses an epipolar model with the smallest upper and lower parallax to improve the accuracy and efficiency of semi-global dense matching.

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

[0084] Step 1: Let the inclination angle of the line passing through the pixel point p on the left image be α, from which the equation of the left line can be obtained as yy p =tanα(xx p );

[0085] Step 2: Select n pixels at equal intervals on the left line and calculate the pixel coordinates (x i ,y i ), i=1,2,…,n and different geoid heights Z i , the corresponding geodetic coordinates (X i ,Y i );

[0086] Step 3: Calculate the geodetic coordinates (X i ,Y i ,Z i ) on the right image. i ',y i ');

[0087] Step 4: For (x i ',y i ') for straight line fitting, the least square method can be used to obtain the right line equation y'=kx'+b, and the coefficients are:

[0088]

[0089] The right core line inclination angle α'=arctank', the root mean square of the vertical disparity of each pixel is:

[0090]

[0091] Among them, Di is the root mean square of the upper and lower parallax of the i-th pixel, and

[0092]

[0093] Step 5: Record the root mean square of the upper and lower parallax obtained in the tth test. If it is the current minimum value, record σ t =σ,α t =α;

[0094] Step 6: Repeat steps 1 to 5 within the range of the epipolar line inclination until all inclination angles are tested.

[0095] Step 7: Use the parabola equation to calculate the left epipolar inclination angle α min , and then fit the left kernel line to the pixel points with the same name on the right image to obtain the right kernel line f(α):

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

[0097] Among them, the parameters u, v, and w are expressed by α t-1 , α t, α t+1 And the three corresponding upper and lower parallax root mean square equations are solved.

[0098] (3) Cost calculation optimization

[0099] When the cost is calculated based on mutual information, the calculation process of mutual information is relatively complicated, and the calculation process has requirements on the number of iterations, which is inefficient. There are currently a variety of improvement strategies for cost calculation, among which the SGBM (Semi-Global Block Matching) algorithm in the computer vision open source library OpenCV is simple to calculate and has high matching efficiency. The BT (Block Truncation) strategy it adopts can reduce the error caused by the discretization of image sampling from the perspective of information sampling. At the same time, the SGBM algorithm adds a preprocessing process when performing cost calculation, uses the Sobel operator for preprocessing, and then adds the BT cost of the preprocessed image and the BT cost of the original image as the initial cost. Therefore, the present invention adopts the SGBM algorithm for cost calculation, and its calculation formula is as follows:

[0100]

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

[0102] (4) Image denoising based on Cauchy proximal splitting algorithm

[0103] This improved algorithm uses the Cauchy Proximal Splitting (CPS) algorithm for image denoising. It boasts excellent convex optimization capabilities and good convergence performance, effectively solving the L1 regularization problem. The CPS algorithm effectively constrains noise while preserving image detail and texture information, offering significant performance advantages. The key idea behind the CPS algorithm is to transform the image denoising problem into an equivalent convex optimization problem using an augmented Lagrangian method. It then iteratively computes the problem using a proximal splitting algorithm (including Douglas Rachford splitting, forward-backward splitting, or alternating direction multiplier splitting). Noise reduction is considered complete if the results converge. This algorithm is primarily suitable for semi-global dense matching, particularly in situations with low image quality and the presence of noise. It can improve image clarity and enhance object detail, thereby improving matching accuracy. This improved algorithm adds a pre-matching condition: the Peak Signal-to-Noise Ratio (PSNR) metric measures image quality. If the PSNR is below a threshold, the CPS algorithm is used.

[0104] In summary, the improved algorithm adopted in the present invention proposes an improved method that takes image quality into consideration on the basis of the original SGM algorithm. The image quality is used to determine whether the CPS algorithm is needed. Then, the algorithm is further improved by using methods such as image cropping and pixel coordinate conversion, epipolar image construction, and cost calculation optimization. The improved algorithm has a certain degree of improvement in matching accuracy and algorithm efficiency, and has excellent performance in image denoising, and 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 categories of coordinate stereo positioning models: strict geometric models and general geometric models. Traditional aerial photogrammetry coordinate stereo positioning uses a strict geometric model, which requires the use of the camera's internal and external orientation elements. The sensor parameters of commercial satellites are highly confidential, so satellite images usually use a general geometric model that is independent of sensor parameters. The rational function model is a general geometric model that calculates geodetic coordinates based on the RPC file (i.e., the parameter file of the rational function model) provided by the satellite image. The rational function model has the characteristics of being independent of the coordinate reference system, independent of the satellite sensor, and having good confidentiality. It is an important imaging model for remote sensing satellite sensors. 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] Where, P i (X n ,Y n ,Z n ) is of the following form:

[0112]

[0113] Where i = 1, 2, 3, 4; (r n ,c n ) is the normalized pixel coordinate, (X n ,Y n ,Z n ) is the standardized geodetic coordinate; r o , c o , X o , Y o , Z o are all standardized translation parameters, r s , c s , X s , Y s , Z s are standardized proportion parameters, a i are known polynomial coefficients.

[0114] In the rational function model, the linear term is used to represent the optical projection error, the quadratic term is used to represent the errors caused by factors such as the earth's curvature, atmospheric refraction, and lens distortion, and the cubic term is used to represent errors with high-order components such as camera shake.

[0115] (2) Coordinate stereo positioning algorithm

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

[0117] Coordinate stereo positioning, also known as spatial intersection in aerial photogrammetry, uses the same-name image points of a stereo pair to calculate the corresponding geodetic coordinates. The above steps can be used to obtain the final geodetic coordinates (X, Y, Z). 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 of the left and right images as the initial value, 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] Among them, l and r are the left image and the right image respectively, X ol is the initial left image of the X-axis space rectangular coordinate, X or is the initial right image of the X-axis space rectangular coordinate, Y ol is the initial left image of the Y-axis spatial rectangular coordinate, Y or is the initial right image of the Y-axis space rectangular coordinate, Z ol is the initial left image of the Z-axis rectangular coordinate, Z or is the initial right image of the Z-axis space rectangular coordinate, X sl is the left image after the iterative operation of the X-axis space rectangular coordinate, Y sl is the left image after the iterative operation of the Y-axis spatial rectangular coordinate, Z sl is the left image after the iterative operation of the rectangular coordinates of the Z axis space, X sr is the right image after the iterative operation of the X-axis space rectangular coordinate, Y sr is the right image after the iterative operation of the Y-axis space rectangular coordinate, Z sr It is the right image after the iterative calculation of the rectangular coordinates of the Z-axis space;

[0121] S2: Establish error equations:

[0122] shilling Then r n and c n Substitute the standardized formula into it and we can get:

[0123]

[0124] Expand the above formula into a linear term at (X, Y, Z) using Taylor's formula to obtain:

[0125]

[0126] Therefore, the error equation established is:

[0127]

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

[0129]

[0130]

[0131] Then, the coordinates of the same-name image points of the left and right images (r l ,c l ) and (r r ,c r ) can be obtained as follows:

[0132]

[0133] In the rational function model, generally speaking, r usually represents row, c usually represents column, and v rl Indicates the row coordinates of the left image, v cl Indicates the column coordinates of the left image, v rr Indicates the row coordinates of the right image, v cr Represents the column coordinates of the right image;

[0134] After simplifying the equation (16) with symbols, it can be expressed as V = A△-l;

[0135] S3: Solve the coordinate correction value: Use the least squares method for the error equation, and the least squares solution of △ can be obtained as:

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

[0137] S4: Perform iterative solution: If the coordinate correction value (△X,△Y,△Z) exceeds the threshold, 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, the calculation is completed and the coordinates (X, Y, Z) at this time are the final geodetic coordinates;

[0138] The coordinates (X, Y, Z) obtained at this time represent latitude, longitude, and altitude respectively. After experiments, the correction values of X and Y are generally small, while the correction value of Z is large. The threshold is set to 1×10 －8 When , the influence of subsequent iterations on the coordinate results can be ignored, and high-precision coordinate stereo positioning can be achieved.

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

[0140] The traditional method of converting a disparity map into a point cloud is to use disparity to obtain depth through triangulation, and then obtain a point cloud. The present invention can directly obtain the coordinates of the point cloud from the disparity map based on the rational function model. However, since the longitude and latitude coordinates are obtained, the shape on the xoy plane will change when the point cloud is visualized. Moreover, the result of the coordinate conversion of the point cloud is not reliable and will still be deformed. Therefore, the longitude and latitude coordinates of the point cloud can be replaced with the pixel coordinates of the left image. This method will not cause deformation. The specific operation steps are as follows:

[0141] ① According to 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) to obtain the coordinates of the image points with the same name in the left and right images;

[0142] ② Obtain the coordinates of the image points with the same name before image cropping according to formula (1);

[0143] ③According to the coordinates of the image points with the same name, use the above-mentioned calculation steps of the geodetic coordinates to obtain the corresponding geodetic coordinate points;

[0144] ④Replace the longitude and latitude coordinates of the geodetic coordinate point 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 universal geometric model with good confidentiality. It is in the form of a ratio of polynomials. After linearization, the error equation can be obtained. The error equations of the left and right images are combined, and the coordinate correction amount is solved using the least squares method. The geodetic coordinates can be solved by iteration. After the present invention calculates the geodetic coordinates through the rational function model, the elevation error is still too large, and the elevation error plays a decisive role in the accuracy of the point cloud. Therefore, 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] The present invention performs point cloud preprocessing before reconstructing the DSM from the three-dimensional point cloud. There are two main steps in point cloud preprocessing, namely, eliminating mismatched points and diluting the point cloud. For the mismatched points in the point cloud, since they are mainly distributed outside a certain island and are relatively small in number, they can be manually eliminated. The principle of manual elimination is relatively simple. After visualizing the point cloud using MATLAB, the coordinate values of these mismatched points can be viewed. After obtaining the coordinate values, the coordinate points can be directly deleted from the source file of the point cloud. The reason for diluting the point cloud is that the point cloud obtained by dense matching is too dense. Compared with real objects, most of the points are redundant points, which will affect the performance of the program and the effect of the final generation of the DSM. Therefore, the point cloud is diluted proportionally before generating the DSM, and the dilution ratio is adjusted according to the specific situation.

[0150] (2) Point cloud reconstruction DSM algorithm principle

[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: 1) less data redundancy and faster processing; 2) high computational efficiency, making them well-suited for high-precision modeling. This paper uses TINs to represent DSMs. While there are many methods for generating TINs, the Delaunay triangulation is currently the most common.

[0152] There are three main algorithms for generating Delaunay triangulations: 1. The divide-and-conquer algorithm, whose basic idea is to decompose the problem and reduce the difficulty of generating a triangulation. A large point set is divided into smaller sets, and then the triangulations generated by each set are merged into the final triangulation. This algorithm has the advantage of being time-efficient, but the computational process requires a large number of recursive operations and consumes a lot of memory. 2. The point-by-point insertion algorithm, whose basic idea is to construct a convex hull containing all points, establish an initial triangulation within the convex hull, and then insert the remaining points one by one into the triangulation. Although this algorithm is easy to implement, it is relatively inefficient. 3. The triangulation growth algorithm. The basic idea of the triangulation growth algorithm is to randomly select a point from discrete elevation points, search for the closest 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, the algorithm finds an extension point that forms a Delaunay triangle with the baseline. A Delaunay triangle is generated, and the two newly generated edges of the triangle are used as the new baseline. The search for extension points and the generation of Delaunay triangles are repeated until no extension points are found.

[0153] The algorithm for generating TIN from point cloud in the present invention adopts a growth algorithm and optimizes the original triangulation growth algorithm by LOP optimization, eliminating redundant data, and constructing a convex closure, making TIN generation more efficient. The specific steps include:

[0154] ①LOP optimization

[0155] The comprehensive local optimization procedure (LOP) performs an LOP check on the constructed triangulated network. Four points form two colateral triangles. If the circumcircle of a triangle contains the fourth point, the two diagonals are swapped. This algorithm can meet the requirements of Delaunay triangles. The LOP optimization algorithm can use the maximum and minimum angle characteristics to make judgments. That is, in the convex quadrilateral ABCD composed of two adjacent colateral triangles, the two diagonals BD and AC are swapped to determine whether the angle values of the six internal angles of the two triangles will increase. The properties of the commonly used triangle circumference angle are used for judgment. If:

[0156] If sin(∠A+∠C)=0, then point A is on the circumcircle of triangle BCD, and there is no need to swap diagonals.

[0157] If sin(∠A+∠C)>0, then point A is outside the circumcircle of triangle BCD, and there is no need to swap diagonals.

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

[0159] ② Eliminate redundant data

[0160] The Delaunay triangulation method based on priority points is a further optimization of the method of excluding closed points. Definition of priority points: If, in the process of constructing the triangulation, the vertex P i Not a closed point, priority is given to expanding to include P i The edge until P i Become a closed point, then P i It is called the priority point.

[0161] The present invention adopts the method of constructing a triangulated network with a priority point as the center, and timely excludes closed points from the data point set, which can reduce the number of traversals in the process of finding triangle extension points and reduce the calculation time. Whether a point is a closed point can be determined by judging whether the sum of the angles of all triangles formed by the point in the triangulated network at the point is 360° to determine whether it is a closed point (except for points on the convex closure of the triangulated network). Constructing a triangulated network with a priority point as the center is to give priority to selecting the edge containing the priority point when selecting the baseline. When all the edges containing the point in the edge set participate in the construction of the triangulated network, the priority point will become a closed point. After that, the point can be eliminated from the point set, which can reduce the number of points that need to be traversed in the subsequent construction of the triangulated network. Deleting the closed point will not affect the subsequent construction of the triangulated network, and can also reduce the number of data points in the set, thereby improving the efficiency of the subsequent construction of the triangulated network.

[0162] ③Construct convex closure

[0163] The triangulation growth algorithm constructs a triangulation in order of point coordinate size, requiring a 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 divided into two edge set regions. The algorithm is then simplified by using a method called priority point exclusion to eliminate point and edge sets in different regions using different methods.

[0164] When constructing a Delaunay triangulation, it is necessary to traverse the edge set that stores the baselines to find the next baseline to be constructed. Each edge can be used as a side of at most two triangles, while an edge on a convex closure can only be used as a baseline once to construct the triangulation. Therefore, after an edge is determined to be a baseline twice, it 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 such as AB on the convex closure can only serve as the primary baseline of △ABC, while an edge BC inside the convex closure can serve as the baseline of both △ABC and △BCD. Therefore, we can remove the edge AB from the edge set after detecting it as the primary baseline, and remove the edge BC from the edge set after detecting it as the secondary baseline. This can reduce the data that needs to be processed later.

[0165] (3) Point cloud reconstruction DSM algorithm

[0166] Based on the above optimization, the point cloud reconstruction DSM algorithm proposed in this invention is as follows: Figure 4 As shown, the following specific steps are included:

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

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

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

[0170] Step 4.3: Using the first edge in BaseStack as the baseline, find an extension point that complies with the Delaunay theorem for LOP optimization and construct the first triangle on the right. Add the two newly generated edges, excluding the first edge in BaseStack, to edge set B1 and the newly generated triangle to triangle set T1. These two new edges are then added to BaseStack as new baseline edges.

[0171] Step 4.4, take a point P in the point set D2 i As the priority point, select the one containing P from BaseStack i The edge of the triangle is used as the baseline, and the Delaunay triangulation is constructed on the right side. The baseline edge is added to the BaseStack. In the point set D1, an extension point that conforms to the Delaunay law is found for LOP optimization. After the Delaunay triangle is constructed, the two newly generated edges are added to the edge set B1 and the newly generated triangle is added to the triangle set T1. These two new edges are added as new baseline edges to the BaseStack. If there are duplicates, they are not added again.

[0172] Step 4.5, check P i Is it a closed point? If so, delete it from the point set D1 and delete the point containing P from BaseStack. i edge;

[0173] Step 4.6, after creating the Delaunay triangle on the right, delete this baseline edge from BaseStack and repeat steps 4.4 and 4.5;

[0174] In step 4.7, when BaseStack is empty, the construction of DSM is completed;

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

[0176] After the initial TIN triangulation is generated according to the point cloud reconstruction DSM algorithm, the effect in some areas is not ideal, and there are still many holes and erroneous triangular faces, so post-processing is required. First, the triangular faces in the TIN triangulation are processed, and the faces that are too large are filtered out, and then the holes are filled. When filling the holes, hole identification is required, and all holes except the outer shell are filled, and the faces are refined and smoothed to generate a DSM with a better shape. The above processing flow is implemented using the CGAL library (Computational Geometry Algorithms Library). This third-party library is widely used in computer graphics, geographic information systems and other fields. The final generated DSM results can be saved in a graphic mesh model (DSM three-dimensional model) and generated into a .ply format file, which facilitates the transmission and visualization of DSM results between platforms. The comparison effect of post-processing of DSM results is shown below. Figure 5a-5b shown.

[0177] In summary, the present invention uses a TIN to represent the DSM, and the algorithm for generating the TIN from a point cloud employs a growing algorithm. The original growing algorithm is inefficient, so the present invention optimizes it through methods such as LOP optimization, redundant data elimination, and convex closure construction, improving its efficiency. Adding point cloud preprocessing and DSM result postprocessing before and after point cloud DSM reconstruction improves the DSM performance.

[0178] Example

[0179] 1. Analysis of point cloud generation results

[0180] (1) Visual analysis of point cloud results

[0181] In order to verify the effect of the method of generating point cloud based on rational number model in step 3, the point cloud results are visualized and analyzed. The number of coordinate points generated based on the disparity map is 16250442. The final result is visualized using MATLAB. Figure 6 As shown in the figure, the point cloud visualization shows the general outline and elevation of the island, but some mismatched points and sparse point clouds in some areas are still visible. These will affect the final DSM results, so optimization measures will be required in the subsequent process. Furthermore, if the island is viewed as a plane, the entire plane is tilted. While this does not affect the visualization of the final DSM 3D model, it does affect the precise positioning. Analysis shows that when using the least squares method to solve the coordinate correction, the matrix inversion causes the solution value to oscillate.

[0182] (2) Accuracy analysis of point cloud results

[0183] The generated point cloud was analyzed for accuracy, focusing on its mean, mean square error, and error percentage. The results are shown in Table 1. As can be seen from the table, the mean square error of latitude and accuracy is small, and the error percentage is also small, while the mean square error of elevation is large, and the error percentage is very high. Combined with the visualization results, it can be seen that the error size of elevation plays a decisive role in the accuracy of point clouds during point cloud computing.

[0184] Table 1 Point cloud accuracy analysis table

[0185]

[0186] 2. DSM results comparison and accuracy analysis

[0187] (1) Comparative analysis of DSM results

[0188] Using stereo image pairs to generate DSM can be achieved in many commercial software with high accuracy. Unlike the processing flow in this article, commercial software needs to perform aerial triangulation before generating DSM. The aerial triangulation process first extracts the connection points, and then performs regional network adjustment to output relative orientation results. If there are ground control points, the control points can also be used for adjustment processing to output absolute orientation results. With ground control points, the positioning accuracy and DSM accuracy will be greatly improved. However, due to measurement limitations, ground control points cannot be obtained, so the effect of DSM will be reduced. The DSM generated in this article is compared with the commercial software ENVI. The results are as follows: Figure 7a and Figure 7b As shown:

[0189] Depend on Figure 7a and Figure 7b As can be seen, the DSM generated by ENVI more accurately describes the elevation fluctuations of the land features, but the sea surface elevation error is large. The elevation of most sea surface areas is higher than that of the island area, which is obviously incorrect. Although the visualization effect and description of the land feature elevation of the DSM generated in this paper are not as good as ENVI, the entire DSM is smoother, the elevation changes between land features are smaller, which is consistent with the actual topography of the island. The sea surface area is less affected by noise points, which avoids the influence of sea surface elevation errors on the DSM results.

[0190] (2) DSM result accuracy analysis

[0191] In order 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 (such as Figure 8The DEM data shows an average ground elevation of 26 meters on the island, while the DSM generated by our algorithm shows an average ground elevation of 37.9 meters, with an error of +11.9 meters. The DSM generated by ENVI shows an average ground elevation of 29.5 meters, with an error of +3.5 meters. The ground elevation errors indicate that the ENVI-generated DSM is more accurate, but overall, the accuracy of both DSMs needs to be improved.

[0192] The error source of the DSM generated by the algorithm of the present invention is mainly the elevation calculation error. There are four main reasons for this error: ① The parameter accuracy in the RPC file is not high; ② There are certain systematic errors in the rational function model; ③ The accuracy of the geodetic coordinate solution method is not high; ④ The interference of noise points in the image dense matching process.

[0193] In summary, by comparing the DSM generated by this method with that generated by ENVI, it can be seen that this method describes the surface of remote islands and reefs more realistically 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 rather than to limit the same. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solutions of the present invention with equivalents, and these modifications or equivalent replacements 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: The following steps are involved: Step 1: Obtain stereo images of remote islands and reefs from the World-View3 satellite; Step 2: Based on the stereo image pairs of remote islands and reefs, semi-global dense matching is performed using an improved semi-global dense matching algorithm to obtain a disparity map; Step 3: Based on the obtained disparity map, use the rational function model to perform stereo positioning and obtain a 3D point cloud; Step 4: Pre-process the 3D point cloud obtained in step 3, and use the 3D point cloud reconstruction DSM algorithm to generate the remote island reef DSM, and then post-process the generated remote island reef DSM to construct a remote island reef DSM 3D model.

2. The method for extracting high-precision DSM from stereo image pairs of islands and reefs without measurement control points according to claim 1, characterized in that: The specific steps of step 3 include: Step 3.1: Based on the disparity value d corresponding to each pixel in the disparity map, the coordinates of the same-name image points in the left and right images are obtained using the following formula: r r =r l ,c r =c l -d Among them, r represents the right image and l represents the left image; Step 3.2: Obtain the coordinates of the image points with the same name before cropping the left and right images according to the following formula: Rn'=Rn+dR Cn'=Cn+dC Where dR and dC are the coordinate conversion coefficients using the same pixels; (Rn, Cn) are the coordinates before image cropping, and (Rn', Cn') are the coordinates after image cropping; Step 3.3: Based on the rational function model, calculate the corresponding geodetic coordinate points according to the coordinates of the image points with the same name; Step 3.4: Replace the longitude and latitude coordinates of the geodetic coordinate point with the corresponding pixel coordinates of the left image; Step 3.5: Output the final point cloud coordinates.

3. The method for extracting high-precision DSM from stereo image pairs of islands and reefs without measurement control points as claimed in claim 2, characterized in that: The specific steps of step 3.3 include: Step 3.3.1: Calculate the average value of the normalized translation parameters of the left and right images respectively, use it as the initial value, and convert these two initial values into normalized coordinates (X nl ,Y nl ,Z nl ) and (X nr ,Y nr ,Z nr ), get the initial value of the geodetic coordinate (X 0 ,Y 0 ,Z 0 ), the specific calculation formula is: Among them, l and r are the left image and the right image respectively, X ol is the initial left image of the X-axis space rectangular coordinate, X or is the initial right image of the X-axis space rectangular coordinate, Y ol is the initial left image of the Y-axis spatial rectangular coordinate, Y or is the initial right image of the Y-axis space rectangular coordinate, Z ol is the initial left image of the Z-axis rectangular coordinate, Z or is the initial right image of the Z-axis space rectangular coordinate, X sl is the left image after the iterative operation of the X-axis space rectangular coordinate, Y sl is the left image after the iterative operation of the Y-axis spatial rectangular coordinate, Z sl is the left image after the iterative operation of the rectangular coordinates of the Z axis space, X sr is the right image after the iterative operation of the X-axis space rectangular coordinate, Y sr is the right image after the iterative operation of the Y-axis space rectangular coordinate, Z sr It is the right image after the iterative calculation of the rectangular coordinates of the Z-axis space; Step 3.3.2: Based on the normalized coordinates of the left and right images, construct the following error equation: Step 3.3.3: Based on the coordinates of the same-name image points in the left and right images (r l ,c l ) and (r r ,c r ), rewrite the above error equation as: Among them, r represents the row, c represents the column, and v rl Indicates the row coordinates of the left image, v cl Indicates the column coordinates of the left image, v rr Indicates the row coordinates of the right image, v cr Represents the column coordinates of the right image; Then the above formula is simplified into V = A△-l; Step 3.3.4: Apply the least squares method to the above error equation and obtain the least squares solution of △ as: △=[△X△Y△Z] T =(A T A) -1 A T l Step 3.3.5: Determine whether the coordinate correction value (△X, △Y, △Z) exceeds the threshold value 1×10 －8 If so, use the corrected geodetic coordinates (X 1 ,Y 1 ,Z 1 ) calculates the standardized coordinates of the left and right images and returns to step 3.3.2 for iterative calculation; otherwise, the calculation is completed and the coordinates (X, Y, Z) at this time are the final geodetic coordinates.

4. The method for extracting high-precision DSM from stereo image pairs of islands and reefs without measurement control points as claimed in claim 3, characterized in that: The specific steps of step 4 include: Step 4.1: Perform point cloud preprocessing on the point cloud coordinates output from step 3.3; Step 4.2: Based on the pre-processed point cloud coordinates, DSM reconstruction is performed 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 a graphical mesh model.

5. The method for extracting high-precision DSM from stereo image pairs of islands and reefs without measurement control points as claimed in claim 4, characterized in that: The point cloud preprocessing operations in step 4.1 include: 1) Deletion of incorrect matching points Based on the point cloud coordinates, MATLAB is used to visualize the point cloud and check the coordinate values of the mismatched points distributed outside a certain island. After obtaining the coordinate values, the coordinate points are deleted from the source file of the point cloud. 2) Point cloud dilution The point cloud after deleting the mismatched points is diluted proportionally.

6. The method for extracting high-precision DSM from stereo image pairs of islands and reefs without measurement control points according to claim 5, characterized in that: The specific steps of the point cloud reconstruction DSM algorithm in step 4.2 include: Step 4.2.1: Read the point cloud coordinate data and create the point set D1, edge set B1, triangle set T1, and baseline edge stack BaseStack. These sets contain 3D coordinates and information about edges and triangles. Step 4.2.2: Select appropriate coordinate points from point set D1 to create the convex closure of the triangulated network and store all edges on the convex closure in the baseline edge stack BaseStack. Then construct point set D2 and store all points on the convex closure in it. Step 4.2.3: Using the first edge in BaseStack as the baseline, find an extension point that complies with the Delaunay theorem for LOP optimization and construct the first triangle on the right. Add the two newly generated edges to edge set B1 and the newly generated triangle to triangle set T1. Add these two new edges as new baseline edges to BaseStack. Step 4.2.4: Take a point P in the point set D2 i As the priority point, select the one containing P from BaseStack i The edge of the triangle is used as the baseline, and the Delaunay triangulation is constructed on the right side. The baseline edge is added to the BaseStack. In the point set D1, an extension point that conforms to the Delaunay law is found for LOP optimization. After the Delaunay triangle is constructed, the two newly generated edges are added to the edge set B1 and the newly generated triangle is added to the triangle set T1. These two new edges are added as new baseline edges to the BaseStack. If there are duplicates, they are not added again. Step 4.2.5: Check P i Is it a closed point? If so, delete it from the point set D1 and delete the point containing P from BaseStack. i edge; Step 4.2.6: After creating the Delaunay triangle on the right side, delete this baseline edge from BaseStack and repeat steps 4.2.4 and 4.2.5; Step 4.2.7: When BaseStack is empty, DSM construction is complete.

Images