Image and laser height measurement combined adjustment method under water edge line constraint

Abstract

The present application belongs to the field of surveying and mapping, and more particularly relates to a method for combined adjustment of island and reef images and laser height measurement under water edge line constraint. The present application comprises the following steps: S1, input of multi-source satellite data; S2, extraction of water edge line points; S3, extraction of photon points; S4, preparation for combined adjustment, including S41, stereo image preprocessing, S42, generation of water edge line height control points, and S43, generation of laser point height control points; S5, combined adjustment and compensation, comprising S51, combined adjustment, and S52, affine transformation compensation on the adjustment result. The present application introduces elevation constraints such as water edge line and fuses spaceborne laser height measurement data to construct multi-source height control points, realizes high-precision three-dimensional mapping of sea reef areas under the condition of no ground control, and thus more efficiently and stably solves the mapping problems of difficult acquisition of field control points, insufficient precision of open source DEM, and uneven distribution of laser points.

Description

Technical Field

[0001] This invention belongs to the field of surveying and mapping technology, and more specifically, relates to a method for joint adjustment of island and reef images and laser altimetry under waterline constraints. Background Technology

[0002] Currently, with the rapid development of Earth observation technology, high-resolution optical satellite remote sensing imagery has become an indispensable data source for acquiring global geospatial information, including that of sea reefs, especially domestically produced high-resolution stereo mapping satellites such as ZY-3 and GF-7. Regional network adjustment combined with control points is a key step in improving positioning accuracy. However, control points for sea reefs are often difficult to obtain in the field. Therefore, how to improve positioning accuracy under controlless conditions has become a hot topic and a challenge in sea reef stereo mapping research.

[0003] Adding known elevation control constraints to stereo image pairs can effectively improve adjustment accuracy. Early methods extracted elevation information from open-source DEMs as elevation constraints for adjustment. However, the elevation accuracy of open-source DEMs is highly correlated with terrain type, achieving high accuracy only in flat terrain areas. Furthermore, data quality is unevenly distributed across different regions, making it difficult to apply to regional network adjustment in complex terrain areas. With the launch and operation of various satellites carrying spaceborne laser altimeters, the nominal elevation accuracy of the laser point data they acquire is generally high. For example, the accuracy of ICESat-2 laser altimeter points can reach 0.15m, and the absolute elevation accuracy of GF-7 laser altimeter points can reach 0.10~0.30m, providing a large number of elevation control points for regional network adjustment. Land terrain has rich textures and a large number of laser points, making the selection and establishment of laser elevation control points relatively simple. However, for small reefs far from the mainland, the number of photon points is usually sparse and unevenly distributed, making effective control of the reef difficult.

[0004] Chinese patent document CN118710941A describes a method for matching spaceborne laser altimetry data with stereo image pairs based on minimum elevation difference. This method finds the optimal matching point by generating candidate sub-points and calculating the overall elevation difference of the buffer zone, but it does not fully consider the constraint of consistent waterline elevation. In areas such as ocean reefs, where the number of laser points is sparse and their distribution is uneven, relying solely on laser altimetry data is insufficient to guarantee overall elevation accuracy.

[0005] Chinese patent document CN113124834A describes a method, system, and storage medium for regional network adjustment that combines multi-source data. It obtains water surface elevation data through lidar point cloud data and uses epipolar constraints to extract corresponding image points for regional network adjustment. However, it relies too heavily on lidar point cloud data for elevation control, which limits its application in areas such as oceans and reefs where it is difficult to obtain dense lidar point cloud data.

[0006] Chinese patent document CN104931022A proposes a satellite imagery stereo area network adjustment method based on spaceborne laser altimetry data. This method improves the compatibility of the adjustment algorithm with various types of control data by separating planar and elevation control data. However, it relies too heavily on the joint adjustment of data from the image control point library and the generalized elevation control point library. For areas such as oceans and reefs where it is difficult to obtain sufficient image control points, the application effect of this method is limited.

[0007] In summary, there is an urgent need to provide a joint adjustment method for island and reef imagery and laser altimetry under waterline constraints to solve the problems existing in the current technology. Summary of the Invention

[0008] The technical problem to be solved by this invention is to construct multi-source elevation control points by introducing elevation constraints such as the waterline and integrating satellite-borne laser altimetry data, thereby achieving high-precision three-dimensional mapping of sea reef areas under conditions without ground control, thus solving the mapping problems of difficulty in obtaining field control points, insufficient accuracy of open-source DEMs, and uneven distribution of laser points more efficiently and stably.

[0009] The technical terms related to this invention are explained as follows:

[0010] High-resolution optical satellites: Images of the Earth's surface taken by high-resolution optical satellites have high spatial resolution and can clearly display surface details. They are widely used in map making, urban planning, environmental monitoring, and other fields.

[0011] High-resolution stereo mapping satellite: Equipped with satellites capable of simultaneously acquiring multi-angle images of the same area, it generates a digital elevation model (DEM) through stereo image pairs, thereby realizing three-dimensional terrain mapping.

[0012] Control points: In photogrammetry and remote sensing, these are ground points with known precise geographic coordinates, used for geometric correction and positioning of images. Control points can be artificially established markers or natural feature points.

[0013] Spaceborne laser altimeter: A laser altimeter installed on a satellite that calculates the distance from the satellite to the ground by emitting laser pulses to the ground and measuring the echo time, thereby obtaining the ground elevation information.

[0014] Water boundary line: A linear feature at the boundary between a body of water and land, serving as the dividing line between oceans, lakes, rivers, and the surrounding land. In remote sensing image processing, water boundary lines are often used to extract water body information or for geometric correction.

[0015] The technical problem to be solved by this invention is achieved by the following technical solution:

[0016] A method for joint adjustment of island and reef images and laser altimetry under waterline constraints includes the following steps:

[0017] S1. Input of multi-source satellite data:

[0018] Input high-resolution stereo mapping satellite stereo images as the optical stereo data source for sea reefs;

[0019] Input laser altimetry data from a satellite equipped with a spaceborne laser altimeter as the data source for measuring the elevation of sea reefs;

[0020] S2, Extraction of waterline points:

[0021] The instantaneous waterline data source of the sea reef is vectorized, rasterized and spatially intersected to obtain the spatial coordinates of the waterline points.

[0022] S3. Photon point extraction:

[0023] The data sources for sea reef elevation measurement are filtered by attributes and elevation differences, and the data is cleaned to obtain pure laser altimetry data.

[0024] S4. Preparation for joint adjustment:

[0025] This includes S41 stereo image preprocessing, S42 waterline elevation control point generation, and S43 laser point elevation control point generation.

[0026] S5. Joint Adjustment and Compensation:

[0027] Using the relatively oriented stereo image, the waterline elevation control point, and the laser point elevation control point as input, the S51 joint adjustment is performed, and the adjustment result is compensated by the S52 affine transformation to obtain the stereo image after affine transformation compensation, thus realizing high-precision three-dimensional reconstruction of the sea reef.

[0028] This technical solution combines the waterline elevation consistency constraint with the high-precision advantage of spaceborne laser altimetry data through multi-source data fusion and joint adjustment. First, it acquires the spatial locations of waterline points and photon points, then constructs elevation control points, and finally achieves high-precision reef mapping without ground control points through joint adjustment and affine transformation compensation. Specifically, the waterline of a reef naturally forms an approximately iso-elevation surface, which can serve as a virtual elevation benchmark; spaceborne laser altimetry data possesses centimeter-level elevation accuracy, providing accurate elevation values ​​for the waterline; the combination of these two technologies compensates for the scarcity of laser photon points and the insufficient accuracy of open-source DEMs in reef areas, constructing a high-precision elevation constraint system when no field control points are available. By establishing geometric relationships between images through relative orientation of stereo images, and then conducting joint adjustment with virtual waterline control points and laser points as dual constraints, image system errors and terrain undulation deviations can be eliminated, ultimately achieving high-precision three-dimensional mapping of reefs without ground control.

[0029] In addition, the method for joint adjustment of island and reef images and laser altimetry under waterline constraints proposed in this invention can also have the following additional technical features:

[0030] According to a preferred embodiment of the present invention, the input of multi-source satellite data in S1 is as follows:

[0031] The high-resolution stereo mapping satellite ZY-3 / GF-7 was selected, which provides stereo images that meet the requirements of consistent imaging time and complete coverage of the target sea reef waters and land areas, for obtaining sea reef topographic texture and performing planar position matching;

[0032] ATL03 data was acquired using a satellite equipped with a spaceborne laser altimeter. Laser photon point cloud data containing the target reef area was selected from this data as the data source for reef elevation measurement.

[0033] According to a preferred embodiment of the present invention, the extraction of the waterline points in S2 specifically includes:

[0034] S21. Using image grayscale segmentation and edge detection algorithms, vectorization extraction of the instantaneous waterline of the reef is completed, and the set of pixel coordinates of the waterline on a single scene image is obtained.

[0035] S22. Rasterize the vectorized waterline to unify the coordinate expression format and eliminate coordinate errors caused by differences in image resolution.

[0036] S23. Based on the RFM rational function model of stereo imagery, perform spatial intersection calculations of waterline points to solve for the three-dimensional object space coordinates of the waterline points, including the following calculations:

[0037] Assuming the reef image is a forward-looking image from a stereo image pair, for a point P in object space, its object space coordinates are: Point p is the image point corresponding to ground point P on the stereo image, and the coordinates of point p in the image coordinate system are: The RPC file included with the stereo image provides the 80 rational polynomial coefficients required for the rational function model, as well as the offsets needed for image-side and object-side coordinate normalization. , , , , ) and scaling factor ( , , , , );

[0038] Normalize the image and object coordinates:

[0039] (1)

[0040] Based on the normalized object coordinates ( , , Calculate the image-side coordinates. , ):

[0041] (2)

[0042] In the formula: Num and Den are the cubic polynomials of the numerator and denominator, respectively, containing a total of 80 rational polynomial coefficients; These are the polynomial coefficients, where n = 1, 2, 3, 4; polynomial coefficients Provided by the RPC file that comes with the stereoscopic image.

[0043] According to a preferred embodiment of the present invention, the extraction of photon points in S3 specifically includes:

[0044] S31. Filter the attribute fields that characterize the land surface type in the ATL03 laser altimetry data and remove ineffective surface photon points such as water bodies and noise.

[0045] S32. Combine ATL03 data to calculate the difference between the laser photon point elevation and the open-source DEM elevation, set an elevation difference threshold, and filter out abnormal photon points whose elevation deviation exceeds the threshold.

[0046] S33. Complete multiple rounds of data cleaning and quality verification to obtain clean, high-precision laser altimetry data for sea reefs.

[0047] According to a preferred embodiment of the present invention, the S41 stereoscopic image preprocessing specifically includes:

[0048] S411. Perform corresponding point matching on the ZY-3 / GF-7 stereo images, and use a robust feature matching algorithm to extract matching connection points between images to ensure that the connection points are evenly distributed and the matching accuracy meets the requirements.

[0049] S412. Based on the matching connection points, perform free network adjustment calculations, solve the relative orientation parameters between stereo images, complete the relative orientation processing of the images, and obtain the relative orientation stereo images.

[0050] According to a preferred embodiment of the present invention, the generation of the S42 waterline elevation control point specifically includes:

[0051] S421. Based on the image-space coordinate system, perform matching calculations between laser photon points and the shoreline, and select laser photon points that match the position of the shoreline.

[0052] S422: Replace the elevation of the three-dimensional object space position coordinates of the waterline point generated by the spatial intersection of S23 with the elevation of the laser photon point after matching in S421. Combining the elevation accuracy advantage of laser altimetry data with the planar position characteristics of the waterline, generate a waterline elevation control point that has both planar and elevation accuracy.

[0053] According to a preferred embodiment of the present invention, the generation of the S43 laser point elevation control point specifically includes:

[0054] S431. Select characteristic laser points from the clean laser altimetry data, and project the ground coordinates of the laser points onto the image space of the ZY-3 / GF-7 stereo image using the RPC rational function model to obtain the corresponding image point coordinates of the laser points on the image.

[0055] S432. Establish the relationship between the three-dimensional coordinates of the laser point on the ground and the coordinates of the image point, and complete the construction of the laser point elevation control point.

[0056] According to a preferred embodiment of the present invention, the S51 joint adjustment specifically includes:

[0057] S511. Using the matching connection points of the relative orientation stereo image, the waterline elevation control points, and the laser point elevation control points as the core observation values, construct a joint adjustment error equation model.

[0058] S512. Set reasonable weight allocation principles, carry out multi-constraint joint adjustment calculations, and solve the orientation parameters of stereo images and the high-precision coordinates of ground points.

[0059] According to a preferred embodiment of the present invention, the S52 affine transformation compensation specifically includes:

[0060] S521. Construct an image-side affine transformation model with 6 parameters, correct the image point coordinates calculated by RFM, and obtain the accurate image point coordinates of the adjusted observations.

[0061] S522. Using an affine transformation model to eliminate image system errors and coordinate deviations caused by topographic undulations, a affine transformation-compensated stereo image is obtained, achieving high-precision inversion of the planar position and elevation of the reef topography. This includes the following calculations:

[0062] The affine transformation model uses six parameters Image point coordinates are corrected using three types of parameters: translation, scaling, and coupling. The image point coordinate translation parameters are used to compensate for the overall image offset. The row and column aspect ratio parameters of the image points are used to compensate for differences in image scale. The row and column coupling parameters of the image points are used to compensate for image rotation and distortion; affine transformation formula:

[0063] (3)

[0064] In the formula: (r, c) are the image point coordinates obtained by RFM calculation; , ) is the image-side correction, obtained by compensating for the affine transformation parameters; , The coordinates of the image points are the precise image point coordinates after affine transformation compensation, which will be used as the observations for subsequent adjustment.

[0065] Regional network adjustment constructs error equations containing two types of unknowns: affine transformation parameters and ground point coordinates, and solves for the optimal solution using normal equations; Connectivity point error equations:

[0066] (4)

[0067] In the formula: the unknowns are the image-side affine transformation compensation parameters and the ground 3D coordinates of the connection points, and V1 is the image point coordinate residual vector of the connection points. , , , where is the measured image point coordinate of the connection point obtained through image matching. , These are the image point coordinates after affine transformation compensation and ground coordinate compensation.

[0068] This reflects the effect of changes in affine parameters on the coordinates of the image point. Let be the affine transformation parameter correction vector for the image point, representing the optimal solution corrections for the six affine parameters, which are unknowns to be solved; A is a 2×6 matrix, obtained by taking the partial derivatives of the RFM model with respect to the six affine parameters. The partial derivatives are calculated using numerical differentiation with a step size of 0.001. Then:

[0069]

[0070] In the formula: Reflects the influence of ground coordinates on image point coordinates; The ground three-dimensional coordinate correction vector for the connection point; Given a 2×3 matrix, obtained by taking the partial derivative of the RFM model with respect to ground coordinates and using numerical differentiation with a step size of 0.001, then:

[0071]

[0072] In the formula: = Let be the vector of constant terms in the error equation, where , The image point coordinates are calculated using the initial affine transformation parameters and the initial ground values. The weight matrix;

[0073] In summary, the matrix form of the connection point error equation is as follows:

[0074] (5)

[0075] Based on the connection point error equation (4) and the principle of least squares, the elevation control point error equation (6) is obtained:

[0076] (6)

[0077] Solving the normal equations yields the affine parameter correction t and the ground coordinate correction x1. Adding the corrections to the initial values ​​yields the optimal affine parameters and the ground coordinates of the connection points, enabling high-precision relative orientation of the stereo image pair.

[0078] According to a preferred embodiment of the present invention, in step S432, based on the terrain undulation characteristics, corresponding connection points are searched within a certain distance range of the photon point. The image point coordinates of the corresponding connection points are spatially intersected to obtain ground coordinates. Neighboring connection points of the photon point are selected by setting a distance threshold. The elevation of the photon point is assigned to the neighboring connection points to construct laser elevation control points, including the following calculations:

[0079] Error equations for laser and waterline elevation control points:

[0080] (7)

[0081] In the formula: V2 is the image point coordinate residual vector; t is defined the same as the connection point and is the affine transformation parameter correction vector of the image point; The ground latitude and longitude correction vector is given. Since the elevation values ​​of the laser and the waterline elevation control points are true, no correction is needed. A and B2 are the corresponding coefficient matrices. L2 is the constant term obtained by substituting the initial values ​​into the observation equation. P2 is the weight matrix.

[0082] Error equation matrix of laser and waterline elevation control points:

[0083] (8)

[0084] The connection point error equation (4) and the laser and waterline elevation control point error equation (6) are combined and solved using the least squares method to solve for the unknown affine transformation correction t and ground point corrections x1 and x2. The initial value is added to the correction to obtain the optimal affine transformation parameters. In subsequent spatial intersection, the image point positions on the forward and backward images are compensated using the affine transformation parameters to obtain better image point coordinates, thus making the object elevation coordinates obtained from the intersection more accurate.

[0085] The one or more technical solutions provided by this invention have the following advantages compared with the prior art:

[0086] (1) Breaking the dependence on field control points: There is no need to set up control points on the ground. High-precision elevation constraints can be constructed based on the waterline and laser altimetry data, which is suitable for remote sea reef surveying scenarios.

[0087] (2) Improve mapping accuracy and reliability: By combining the elevation characteristics of the waterline and the high precision of laser altimetry, the defects of scarce photon points can be made up for, significantly improving the accuracy of elevation mapping in the reef area and providing stable support for marine geographic information;

[0088] (3) Adapt to complex terrain scenarios: The constraint logic is optimized for small-area sea reefs, which effectively solves the problem of uneven accuracy of open source DEM and improves the robustness and applicability of mapping under complex terrain. Attached Figure Description

[0089] Figure 1 This is a flowchart of the method of the present invention.

[0090] Figure 2 This is a schematic diagram for extracting the spatial coordinates of the waterline.

[0091] Figure 3 This is a schematic diagram of the fine processing of photon dots.

[0092] Figure 4 This is a schematic diagram of the adjustment results of the present invention. Detailed Implementation

[0093] Example 1

[0094] Please refer to Figures 1 to 4 This embodiment provides a method for joint adjustment of island and reef images and laser altimetry under waterline constraints, including the following steps:

[0095] S1. Input of multi-source satellite data: Input high-resolution stereo mapping satellite stereo imagery as the optical stereo data source for sea reefs; input laser altimetry data from satellites carrying onboard laser altimeters as the data source for sea reef elevation measurement.

[0096] S2. Extraction of waterline points: Please refer to... Figure 2In panchromatic imagery, water bodies and land exhibit distinctly different distribution characteristics in their grayscale histograms. This allows for the conversion of continuous grayscale display modes into discrete, categorized display modes. By suppressing grayscale gradations within similar land features, the contrast at the water-land interface is enhanced, facilitating the vectorization of water edges. However, the points obtained through vectorization represent the positions of water edge points within a single scene image, not the actual water edge positions. The Rational Function Model (RFM) is then used to project the vectorized reef water edge vertices back into the pixel coordinate system. RFM is a well-known satellite image positioning and projection model in photogrammetry, widely used for image-to-object coordinate transformation in optical satellite imagery. The steps for calculating pixel coordinates using RFM are shown below.

[0097] Suppose an image is the front view image in a stereo image pair. For a point P in object space, its object space coordinates are: Point p is the image point corresponding to ground point P in this image, and the coordinates of point p in the image coordinate system are: The image's built-in RPC file provides the 80 rational polynomial coefficients required for the RFM model, as well as the offsets needed for image-side and object-side coordinate normalization. , , , , ) and scaling factor ( , , , , ).

[0098] First, the image-side and object-side coordinates need to be normalized. The normalization formula is shown in (1).

[0099] (1)

[0100] Based on the normalized object coordinates ( , , Calculate the image-side coordinates. , The formula for ) is shown in (2).

[0101] (2)

[0102] In the formula: Num and Den are the numerator and denominator, respectively, cubic polynomials, containing a total of 80 rational polynomial coefficients. (n=1, 2, 3, 4) are polynomial coefficients, which are provided by the RPC file that comes with the image.

[0103] The closed waterline vector vertex is projected onto the image using RFM to obtain the vertex pixel coordinates. The sparse vertices are connected into a continuous pixel link using the Bresenham linear interpolation algorithm to obtain the row and column numbers of each pixel point passed through by the waterline, thus obtaining the waterline pixel coordinates on the stereo image pair. The above steps are performed on both the front-view and back-view images in the stereo image pair to obtain the pixel coordinates of the same closed waterline on the front-view and back-view images. Next, the spatial intersection method is used to obtain the object space coordinates of the true waterline point. The spatial intersection method uses multiple images covering the same area, and based on the observation values ​​of the same image points, constructs an error equation according to the generalized form of the collinearity condition equation. This error equation is a simplification of equation (4), that is, the only unknown is the object space coordinates. The three-dimensional spatial coordinates of the target point are solved by the least squares adjustment principle. Since the spatial intersection principle depends on the selection of the elevation position, an optimal elevation surface can be iteratively found so that the overlap between the front-view waterline and the back-view waterline is the highest after the front-view waterline is projected onto this elevation surface.

[0104] Please refer to Figure 2 The left side shows a stereo image pair (two high-resolution satellite images, one fore-view and one for back-view) of the same reef, with the white instantaneous waterline outline extracted through edge detection. The middle section is the rasterization process: the vectorized waterline is converted into a raster (grid) form to unify the coordinate reference; the two images of the stereo image pair are rasterized separately to obtain an overlapping rasterized waterline layer; rasterization can eliminate image resolution differences and provide standardized pixel coordinates for subsequent spatial intersection. The right side shows the spatial intersection result: based on the RPC rational function model of the stereo image pair, forward intersection is performed on the rasterized waterline to calculate the three-dimensional object coordinates; finally, a high-precision waterline vector (thick black line) with geographic coordinates is generated and labeled with latitude and longitude (109°45′~109°46′E, 18°39′~18°40′N). This waterline can then be used as a virtual elevation control point for subsequent joint adjustment with laser altimetry data.

[0105] S3. Photon point extraction: Please refer to... Figure 3The data sources for reef elevation measurements were filtered by attributes and elevation differences, and the data was cleaned to obtain pure laser altimetry data. Initial screening of laser altimetry points was conducted using product quality evaluation indicators to ensure the reliability of photon elevation data. Based on two parameters, isurf (1-5 representing land, ocean, sea ice, land ice, and inland water, respectively) and conf (0-4 representing noise, background filling, low, medium, and high confidence signals, respectively), photon points with an isurf of 1 and conf of 3-4 were selected as high-confidence reef ground laser altimetry points. AW3D30 open-source DEM data was used to filter out photon points with large elevation deviations. AW3D30 DEM has the best elevation accuracy and can therefore be used as reference data for further screening, retaining laser altimetry points where the difference between the ground elevation of the photon point and the interpolated elevation of the DEM grid is less than 1m. Finally, statistical analysis of the photon point elevations was conducted to set elevation ranges, and photon points outside these ranges were excluded.

[0106] Please refer to Figure 3 The top left image shows the raw photon point cloud, displaying the original data distribution of the ICESat-2 ATL03 laser altimeter. Its elevation point cloud along the track distance is densely distributed, containing numerous discrete points, including noise points, water reflection points, and non-target feature points, making direct use inaccurate. The top right image shows the initial attribute screening, which filters the laser data's built-in attribute fields. The screening conditions are isurf=1 (surface type identified as land) and conf=4 (high confidence level), eliminating water reflections, noise points, and low-confidence invalid photons, retaining only points likely representing land features. The grayscale image in the top right corner is the open-source DEM for the corresponding area, used to assist in determining the terrain extent and aiding attribute screening. The bottom right image shows the precise elevation difference screening, which combines external reference data for elevation accuracy verification. Using an open-source DEM (such as the ALOS 30m DEM) as a benchmark, the difference between the laser point elevation and the open-source DEM elevation is calculated, filtering out outliers with excessively large deviations from the reference DEM elevation, further purifying the data. The elevation range control in the lower left corner is used for the final data boundary and range calibration. Based on the actual topographic features of the reef, a reasonable elevation range is set, and extreme values ​​that exceed the normal elevation range of the reef are eliminated. Finally, uniformly distributed, accurate, pure and high-quality laser altimetry data are obtained.

[0107] S4. Preparation for joint adjustment: including S41 stereo image preprocessing, S42 generation of waterline elevation control points, and S43 generation of laser point elevation control points.

[0108] S5. Joint Adjustment and Compensation: Please refer to... Figure 4 Using the relatively oriented stereo image, the waterline elevation control point, and the laser point elevation control point as input, the S51 joint adjustment is performed, and the adjustment result is compensated by the S52 affine transformation to obtain the stereo image after affine transformation compensation, thus realizing high-precision three-dimensional reconstruction of the sea reef. Figure 4 The method intuitively demonstrates the accuracy improvement effect of the joint adjustment method of island and reef imagery and laser altimetry under waterline constraints. The effectiveness of the method is verified by comparing the plane / elevation residual vectors of checkpoints before and after adjustment. Figure 4 Before the adjustment on the left, the black arrows are generally quite long, with some approaching or exceeding 10m, indicating large and unevenly distributed elevation errors; the gray arrows are also quite long, with obvious deviations in planar position; this indicates that the stereoscopic mapping results of the reef area have significant systematic and random errors, and the accuracy cannot meet the requirements of high-reliability mapping. Figure 4 After the adjustment on the right, the black arrow shortens significantly, almost to zero, indicating that the elevation error is significantly suppressed; the length of the gray arrow decreases significantly, indicating that the planar position deviation is effectively corrected; this shows that after the joint adjustment, the planar and elevation accuracy of the reef area is greatly improved, the error distribution is more uniform, and the convergence is better.

[0109] The principle is explained below:

[0110] To compensate for the systematic errors inherent in the RFM model and the relative distortions between images caused by factors such as attitude and terrain, a relative orientation step was performed on the fore- and back-view images participating in the joint adjustment. First, tie point matching was performed in the same region within the stereo image pair to obtain dense corresponding points, i.e., matching was performed within the sea, reef, and land area of ​​the test area; second, free-domain adjustment was performed using only corresponding tie points to obtain a high-precision relative-oriented stereo image pair.

[0111] Since the parameters of the RPC have no physical meaning, a face-space affine transformation is used for compensation. This affine transformation uses six parameters... Image point coordinates are corrected using three types of parameters: translation, scaling, and coupling. These are the image point coordinate translation parameters used to compensate for the overall image offset. The row and column aspect ratio parameters of the image points are used to compensate for differences in image scale. The row and column coupling parameters of the image points are used to compensate for image rotation and distortion. The affine transformation is shown in equation (3):

[0112] (3)

[0113] In the formula: (r, c) are the image point coordinates obtained by RFM calculation; , ) is the image-side correction, obtained by compensating for the affine transformation parameters; , The coordinates of the image points are the accurate image point coordinates after affine transformation compensation, which will be used as the observations for subsequent adjustment.

[0114] The core of regional network adjustment is to construct an error equation containing two types of unknowns: "image orientation parameters" (affine transformation parameters) and "ground point coordinates," and then solve for the optimal solution using the normal equation. The tie point error equation is shown in formula (4):

[0115] (4)

[0116] In the formula: the unknowns are the image-side affine transformation compensation parameters and the ground 3D coordinates of the connection points, and V1 is the image point coordinate residual vector of the connection points. , , , where is the measured image point coordinate of the connection point obtained through image matching. , These are the image point coordinates after affine transformation compensation and ground coordinate compensation.

[0117] It reflects the effect of changes in affine parameters on the coordinates of the image point. Let be the affine transformation parameter correction vector for the image point, representing the optimal solution corrections for the six affine parameters, which are unknowns to be solved; A is a 2×6 matrix, obtained by taking the partial derivatives of the RFM model with respect to the six affine parameters. The partial derivatives are calculated using numerical differentiation (step size 0.001).

[0118]

[0119] It reflects the influence of ground coordinates on image point coordinates. The ground three-dimensional coordinate correction vector for the connection point; Given a 2×3 matrix, obtained by taking the partial derivative of the RFM model with respect to ground coordinates (numerical differentiation method, step size 0.001), then...

[0120]

[0121] = Let be the vector of constant terms in the error equation, where , The image point coordinates are calculated using the initial affine transformation parameters and the initial ground values. It is a weight matrix.

[0122] In summary, the matrix form of the connection point error equation is as follows:

[0123] (5)

[0124] Based on the connection point error equation (4) and the principle of least squares, the elevation control point error equation (6) can be obtained:

[0125] (6)

[0126] Solving the normal equations yields the affine parameter correction t and the ground coordinate correction x1. Adding these corrections to the initial values ​​provides the optimal affine parameters and the ground coordinates of the connection points, achieving high-precision relative orientation of the stereo image pair. Since relative orientation only aligns corresponding image points and does not constrain the actual ground point positions, control points are needed for absolute orientation.

[0127] This invention combines the advantages of ZY-3 series data (possessing both orthographic imagery) and GF-7 data (high resolution), employing different strategies to establish a sufficient number of laser elevation control points. For ZY-3 imagery, photon points are projected onto the orthographic image using RFM to obtain photon pixel data. Then, the orthographic image is matched with corresponding points in the fore-sight and back-sight images to obtain their image coordinates, thus establishing the laser elevation control points. For GF-7 data, based on terrain undulation characteristics, connection points are searched within a certain distance of the photon points. The image coordinates of the matched connection points are spatially intersected to obtain the corresponding ground coordinates. An appropriate distance threshold is set to search for connection points near the photon points, and the photon point elevation is assigned to the nearest connection point to establish new elevation control points. The waterline elevation control points were established by projecting pre-selected photon points onto the forward and backward view images using RFM. The KD-Tree algorithm was then used to match the waterline pixels in the images. Considering the planar errors inherent in the images under the original uncontrolled adjustment conditions, a constraint was set for the search: the distance between the projected photon point and the waterline on both the forward and backward view images was less than 3 pixels. The elevations of the photon points meeting the constraint were extracted and averaged; this value represents the true elevation of the waterline points. Next, the true elevation values ​​were combined with the spatial coordinates of the waterline points obtained in the spatial coordinate extraction step to jointly construct the final virtual waterline elevation control points.

[0128] The error equation for the laser and the elevation control point of the waterline is shown in formula (7):

[0129] (7)

[0130] In the formula: V2 is the image point coordinate residual vector; t is defined the same as the connection point and is the affine transformation parameter correction vector of the image point; 1 is the ground latitude and longitude correction vector. Since the elevation values ​​of the laser and the waterline elevation control points are true values, no correction is needed. A and B2 are the corresponding coefficient matrices. L2 is the constant term obtained by substituting the initial values ​​into the observation equation. P2 is the weight matrix.

[0131] The error equation matrix form for the laser and the waterline elevation control points is as follows:

[0132] (8)

[0133] The connection point error equation (4) and the laser and waterline elevation control point error equation (6) are solved simultaneously using the least squares method to obtain the unknowns: the affine transformation correction t and the ground point corrections x1 and x2. Adding the corrections to the initial values ​​yields the optimal affine transformation parameters. During subsequent spatial intersection, the image point positions on the forward and backward views are compensated using the affine transformation parameters to obtain better image point coordinates, thus making the object-space elevation coordinates obtained from the intersection more accurate.

[0134] This invention utilizes the characteristics of short intervals between the forward and backward views of stereo mapping satellite images, consistent elevation of all points on the instantaneous waterline, and the presence of ICESat-2 stripes to provide a large number of evenly distributed elevation control points for joint adjustment. This can greatly eliminate the elevation system errors in the images and provide highly reliable basic geographic information support for various types of reefs in a vast marine area without the need for ground control.

[0135] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the technical solutions of the present invention, and are not intended to limit the specific implementation of the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the claims of the present invention should be included within the protection scope of the claims of the present invention.

Claims

1. A method for joint adjustment of island and reef images and laser altimetry under waterline constraints, characterized in that, Includes the following steps: S1. Input of multi-source satellite data: Input high-resolution stereo mapping satellite stereo images as the optical stereo data source for sea reefs; Input laser altimetry data from a satellite equipped with a spaceborne laser altimeter as the data source for measuring the elevation of sea reefs; S2, Extraction of waterline points: The instantaneous waterline data source of the sea reef is vectorized, rasterized and spatially intersected to obtain the spatial coordinates of the waterline points. S3. Photon point extraction: The data sources for sea reef elevation measurement are filtered by attributes and elevation differences, and the data is cleaned to obtain pure laser altimetry data. S4. Preparation for joint adjustment: This includes S41 stereo image preprocessing, S42 waterline elevation control point generation, and S43 laser point elevation control point generation. S5. Joint Adjustment and Compensation: Using the relatively oriented stereo image, the waterline elevation control point, and the laser point elevation control point as input, the S51 joint adjustment is performed, and the adjustment result is compensated by the S52 affine transformation to obtain the stereo image after affine transformation compensation, thus realizing high-precision three-dimensional reconstruction of the sea reef. The S52 affine transformation compensation specifically includes: S521. Construct an image-side affine transformation model with 6 parameters, correct the image point coordinates calculated by RFM, and obtain the accurate image point coordinates of the adjusted observations. S522. Using an affine transformation model to eliminate image system errors and coordinate deviations caused by topographic undulations, a affine transformation-compensated stereo image is obtained, achieving high-precision inversion of the planar position and elevation of the reef topography. This includes the following calculations: The affine transformation model uses six parameters Image point coordinates are corrected using three types of parameters: translation, scaling, and coupling. The image point coordinate translation parameters are used to compensate for the overall image offset. The row and column aspect ratio parameters of the image points are used to compensate for differences in image scale. The row and column coupling parameters of the image points are used to compensate for image rotation and distortion; affine transformation formula: （3） In the formula: (r, c) are the image point coordinates calculated by RFM; , ) is the image-side correction, obtained by compensating for the affine transformation parameters; , The coordinates of the image points are the precise image point coordinates after affine transformation compensation, which will be used as the observations for subsequent adjustment. Regional network adjustment constructs error equations containing two types of unknowns: affine transformation parameters and ground point coordinates, and solves for the optimal solution using normal equations; Connectivity point error equations: （4） In the formula: the unknowns are the image-side affine transformation compensation parameters and the ground 3D coordinates of the connection points, and V1 is the image point coordinate residual vector of the connection points. , , , where is the measured image point coordinate of the connection point obtained through image matching. , These are the image point coordinates after affine transformation compensation and ground coordinate compensation. This reflects the effect of changes in affine parameters on the coordinates of the image point. Let be the affine transformation parameter correction vector for the image point, representing the optimal solution corrections for the six affine parameters, which are unknowns to be solved; A is a 2×6 matrix, obtained by taking the partial derivatives of the RFM model with respect to the six affine parameters. The partial derivatives are calculated using numerical differentiation with a step size of 0.

001. Then: In the formula: Reflects the influence of ground coordinates on image point coordinates; The ground three-dimensional coordinate correction vector for the connection point; Given a 2×3 matrix, obtained by taking the partial derivative of the RFM model with respect to ground coordinates and using numerical differentiation with a step size of 0.001, then: In the formula: = Let be the vector of constant terms in the error equation, where , The image point coordinates are calculated using the initial affine transformation parameters and the initial ground values. The weight matrix; In summary, the matrix form of the connection point error equation is as follows: （5） Based on the connection point error equation (4) and the principle of least squares, the elevation control point error equation (6) is obtained: （6） Solving the normal equations yields the affine parameter correction t and the ground coordinate correction x1. Adding the corrections to the initial values ​​yields the optimal affine parameters and the ground coordinates of the connection points, enabling high-precision relative orientation of the stereo image pair.

2. The method for joint adjustment of island and reef images and laser altimetry under waterline constraints as described in claim 1, characterized in that: The input of the multi-source satellite data in S1: The high-resolution stereo mapping satellite ZY-3 / GF-7 was selected, which provides stereo images that meet the requirements of consistent imaging time and complete coverage of the target sea reef waters and land areas, for obtaining sea reef topographic texture and performing planar position matching; ATL03 data was acquired using a satellite equipped with a spaceborne laser altimeter. Laser photon point cloud data containing the target reef area was selected from this data as the data source for reef elevation measurement.

3. The method for joint adjustment of island and reef images and laser altimetry under waterline constraints as described in claim 1, characterized in that: The extraction of the waterline points in S2 specifically includes: S21. Using image grayscale segmentation and edge detection algorithms, vectorization extraction of the instantaneous waterline of the reef is completed, and the set of pixel coordinates of the waterline on a single scene image is obtained. S22. Rasterize the vectorized waterline to unify the coordinate expression format and eliminate coordinate errors caused by differences in image resolution. S23. Based on the RFM rational function model of stereo imagery, perform spatial intersection calculations of waterline points to solve for the three-dimensional object space coordinates of the waterline points, including the following calculations: Assuming the reef image is a forward-looking image from a stereo image pair, for a point P in object space, its object space coordinates are: Point p is the image point corresponding to ground point P on the stereo image, and the coordinates of point p in the image coordinate system are: The RPC file included with the stereo image provides the 80 rational polynomial coefficients required for the rational function model, as well as the offsets needed for image-side and object-side coordinate normalization. , , , , ) and scaling factor ( , , , , ); Normalize the image and object coordinates: （1） Based on the normalized object coordinates ( , , Calculate the image-side coordinates. , ): （2） In the formula: Num and Den are the cubic polynomials of the numerator and denominator, respectively, containing a total of 80 rational polynomial coefficients; The coefficients are polynomial coefficients, where n = 1, 2, 3, 4; Provided by the RPC file that comes with the stereoscopic image.

4. The method for joint adjustment of island and reef images and laser altimetry under waterline constraints as described in claim 2, characterized in that: The extraction of the photon point in S3 specifically includes: S31. Filter the attribute fields that characterize the land surface type in the ATL03 laser altimetry data, and remove water bodies, noise and ineffective photon points; S32. Combine ATL03 data to calculate the difference between the laser photon point elevation and the open-source DEM elevation, set an elevation difference threshold, and filter out abnormal photon points whose elevation deviation exceeds the threshold. S33. Complete multiple rounds of data cleaning and quality verification to obtain clean, high-precision laser altimetry data for sea reefs.

5. The method for joint adjustment of island and reef images and laser altimetry under waterline constraints as described in claim 2, characterized in that: The S41 stereoscopic image preprocessing specifically includes: S411. Perform corresponding point matching on the ZY-3 / GF-7 stereo images, and use a robust feature matching algorithm to extract matching connection points between images to ensure that the connection points are evenly distributed and the matching accuracy meets the requirements. S412. Based on the matching connection points, perform free network adjustment calculations, solve the relative orientation parameters between stereo images, complete the relative orientation processing of the images, and obtain the relative orientation stereo images.

6. The method for joint adjustment of island and reef images and laser altimetry under waterline constraints as described in claim 5, characterized in that: The generation of the S42 waterline elevation control point specifically includes: S421. Based on the image-space coordinate system, perform matching calculations between laser photon points and the shoreline, and select laser photon points that match the position of the shoreline. S422: Replace the elevation of the three-dimensional object space position coordinates of the waterline point generated by the spatial intersection of S23 with the elevation of the laser photon point after matching in S421. Combining the elevation accuracy advantage of laser altimetry data with the planar position characteristics of the waterline, generate a waterline elevation control point that has both planar and elevation accuracy.

7. The method for joint adjustment of island and reef images and laser altimetry under waterline constraints as described in claim 6, characterized in that: The generation of the S43 laser point elevation control point specifically includes: S431. Select characteristic laser points from the clean laser altimetry data, and project the ground coordinates of the laser points onto the image space of the ZY-3 / GF-7 stereo image using the RPC rational function model to obtain the corresponding image point coordinates of the laser points on the image. S432. Establish the relationship between the three-dimensional coordinates of the laser point on the ground and the coordinates of the image point, and complete the construction of the laser point elevation control point.

8. The method for joint adjustment of island and reef images and laser altimetry under waterline constraints as described in claim 7, characterized in that: The S51 joint adjustment specifically includes: S511. Using the matching connection points of the relative orientation stereo image, the waterline elevation control points, and the laser point elevation control points as the core observation values, construct a joint adjustment error equation model. S512. Set reasonable weight allocation principles, carry out multi-constraint joint adjustment calculations, and solve the orientation parameters of stereo images and the high-precision coordinates of ground points.

9. The method for joint adjustment of island and reef images and laser altimetry under waterline constraints as described in claim 7, characterized in that: In step S432, based on the terrain undulation characteristics, corresponding connection points are searched within a certain distance range of the photon point. The image point coordinates of the corresponding connection points are spatially intersected to obtain the ground coordinates. By setting a distance threshold, neighboring connection points of the photon point are selected, and the elevation of the photon point is assigned to the neighboring connection points to construct laser elevation control points. This includes the following calculations: Error equations for laser and waterline elevation control points: （7） In the formula: V2 is the image point coordinate residual vector; t is defined the same as the connection point and is the affine transformation parameter correction vector of the image point; The ground latitude and longitude correction vector is given. Since the elevation values ​​of the laser and the waterline elevation control points are true, no correction is needed. A and B2 are the corresponding coefficient matrices. L2 is the constant term obtained by substituting the initial values ​​into the observation equation. P2 is the weight matrix. Error equation matrix of laser and waterline elevation control points: （8） The connection point error equation (4) and the laser and waterline elevation control point error equation (6) are combined and solved using the least squares method to solve for the unknown affine transformation correction t and ground point corrections x1 and x2. The initial value is added to the correction to obtain the optimal affine transformation parameters. In subsequent spatial intersection, the image point positions on the forward and backward images are compensated using the affine transformation parameters to obtain better image point coordinates, thus making the object elevation coordinates obtained from the intersection more accurate.

Images