Field-free on-orbit calibration method for satellite-borne photon counting laser radar

Abstract

The invention provides a field-free on-orbit calibration method for a spaceborne photon counting laser radar, and relates to the technical field of spaceborne photon counting laser radars. According to the method, combined registration is carried out on a high-precision airborne laser point cloud and a high-resolution optical remote sensing image, and the plane positioning precision of the airborne laser point cloud is firstly improved by combining the plane positioning advantage of the high-resolution optical remote sensing image and the elevation precision advantage of the airborne laser point cloud; and the satellite-borne photon counting point cloud is matched with the corrected airborne laser point cloud, and the pointing error and the system ranging error of the satellite-borne photon counting laser radar are corrected, so that the field-free in-orbit geometric calibration of the satellite-borne photon counting laser radar is realized. According to the method, dependence on a ground calibration field laser detector is avoided, the method has the advantages of being high in adaptability, high in calibration precision and wide in application range, and a new technical means is provided for improving field-free on-orbit calibration of a spaceborne photon counting laser radar height measurement system.

Description

Technical Field

[0001] The present invention relates to the technical field of spaceborne photon counting laser radars, and in particular to a field-free on-orbit calibration method for spaceborne photon counting laser radars. Background Art

[0002] As an emerging Earth observation tool, spaceborne photon-counting lidar (LiDAR) boasts significant advantages such as micropulses, high repetition rates, multiple beams, and a small spot size. It has been widely used in fields such as global topography mapping, glacier monitoring, and ocean exploration. However, due to the complex and ever-changing attitude changes, orbital perturbations, and inherent precision limitations of the on-orbit environment, spaceborne photon-counting laser point cloud data often suffer from systematic ranging and pointing errors, which directly impact its spatial positioning accuracy and, in turn, restrict high-precision mapping applications based on spaceborne photon-counting laser data. Therefore, on-orbit calibration of the laser beam pointing error and systematic ranging error of spaceborne photon-counting LiDAR is of great significance.

[0003] Traditional methods for on-orbit calibration of high-precision spaceborne lidars rely primarily on ground-based laser detectors at a calibration site. However, the micropulse characteristics of photon-counting spaceborne lidars require extremely high detection sensitivity from ground-based laser detectors; their high repetition rate requires extremely high temporal synchronization accuracy; and their multi-beam, small spot size necessitates a high density of ground-based detectors. These characteristics make active calibration solutions based on calibration sites extremely labor-intensive, material-intensive, and financially expensive.

[0004] In response to the above shortcomings, a new technical approach is urgently needed to solve the problem of field-free on-orbit calibration of spaceborne photon counting lidar under the premise of limited reliance on ground information. Summary of the Invention

[0005] In view of the above problems, the present invention provides a field-free on-orbit calibration method for a space-borne photon counting lidar.

[0006] The present invention provides a field-free on-orbit calibration method for a spaceborne photon counting laser radar, comprising: step S1, obtaining a high-precision airborne laser point cloud, sequentially performing ground point cloud filtering and denoising on the airborne laser point cloud, and extracting a first structural feature from the denoised airborne laser point cloud; step S2, obtaining a high-resolution optical remote sensing image, preprocessing the optical remote sensing image, and extracting a second structural feature from the preprocessed optical remote sensing image; step S3, registering the first structural feature and the second structural feature after unifying the spatial reference, fusing the registered second structural feature with the first structural feature, and constructing a fused feature set; step S4, plane mapping the airborne laser point cloud based on the fused feature set. The method comprises the following steps: performing translation and plane rotation to correct two plane coordinates in the airborne laser point cloud except the elevation coordinates, and generating a corrected airborne laser point cloud; step S5, obtaining ranging data of the spaceborne photon counting laser radar, performing atmospheric delay correction and tidal correction on the ranging data in sequence, filtering the tidal corrected ranging data to form a spaceborne photon counting point cloud; unifying the spaceborne photon counting point cloud and the corrected airborne laser point cloud with a unified spatial reference, and performing registration to form a registered spaceborne photon counting point cloud; step S6, constructing a radar equation of the spaceborne photon counting laser radar, and determining the pointing error and system ranging error of the spaceborne photon counting laser radar based on the registered spaceborne photon counting point cloud and the radar equation.

[0007] Compared with existing technologies, the field-free on-orbit calibration method for spaceborne photon-counting lidar provided by this invention achieves high-precision estimation and correction of the geometric errors (including laser beam pointing error and ranging error) of the spaceborne photon-counting laser radar system (including laser beam pointing error and ranging error) by constructing a multi-source joint registration and error model framework between high-resolution optical remote sensing images, airborne laser point clouds, and spaceborne photon-counting laser data, without relying on a ground-based calibration field. This method has the following significant technical effects and application value:

[0008] 1. Achieve field-free on-orbit calibration, breaking through the limitations of ground dependence

[0009] This method completely eliminates reliance on costly and difficult-to-deploy ground-based laser receivers. Instead, it constructs a ground-based geometric reference using remote sensing imagery and airborne laser point clouds, enabling system calibration without the involvement of ground-based infrastructure. This approach significantly lowers the barrier to entry for on-orbit lidar geometric calibration, improving the deployment flexibility and operational independence of space application missions. It is particularly suitable for areas such as polar regions and mountainous terrain where traditional calibration facilities are difficult to deploy.

[0010] 2. Systematically solve the problem of joint estimation of ranging error and pointing error

[0011] This method constructs a complete laser error propagation model, which can establish a clear mathematical relationship between the actual spatial registration error and the system's internal parameter deviations (such as azimuth error Δθ, pitch angle error Δφ, and ranging error Δd). It also achieves quantitative estimation of system errors and parameter correction through the least squares method, effectively improving the geometric consistency and positioning accuracy of the laser altimeter system.

[0012] 3. Achieve consistent geometric calibration of multi-beam laser systems

[0013] Compared with the traditional method of calibrating a single beam or main beam, this method is applicable to multi-beam spaceborne photon counting lidar systems. It can jointly model, error compensate and align the data of multiple laser beams with different emission directions, thereby ensuring the spatial coordination and system stability of the entire altimeter system under multi-view conditions.

[0014] 4. Build a multi-source data collaborative registration framework to enhance error constraint capabilities

[0015] This method constructs for the first time a spatial matching framework combining three sources: high-resolution imagery, airborne point cloud, and satellite-borne photon data. Through layer-by-layer registration of image edge features, point cloud structural features, and photon projection points, it forms a multi-layer feature fusion mechanism of surfaces, lines, and points, which enhances the stability of error residual modeling and observation redundancy, and greatly improves the reliability of error estimation and the accuracy of the results.

[0016] 5. High versatility, can be integrated into a variety of orbital remote sensing missions

[0017] This method is independent of specific equipment or platforms and is adaptable to various types of spaceborne photon-counting lidar systems. It can also be embedded in post-processing ground systems or data interpretation workflows, demonstrating its modularity and engineering scalability. Furthermore, this method can be used to perform post-processing geometric correction and accuracy improvement on historical data, imagery from areas lacking ground control points, or laser data. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] The above and other objects, features and advantages of the present invention will become more apparent through the following description of the embodiments of the present invention with reference to the accompanying drawings, in which:

[0019] Figure 1 Schematically shows one of the flow charts of the field-free on-orbit calibration method of a spaceborne photon counting lidar according to an embodiment of the present invention;

[0020] Figure 2 The second flowchart of the field-free on-orbit calibration method of a spaceborne photon counting lidar according to an embodiment of the present invention is schematically shown;

[0021] Figure 3A three-dimensional cloud diagram of observation data of a space-borne photon counting laser radar according to an embodiment of the present invention is schematically shown;

[0022] Figure 4 Schematically showing an elevation comparison diagram of a satellite-borne photon counting point cloud and an airborne laser point cloud before calibration according to an embodiment of the present invention;

[0023] Figure 5 The elevation comparison diagram of the spaceborne photon counting point cloud and the airborne laser point cloud after calibration according to an embodiment of the present invention is schematically shown. DETAILED DESCRIPTION

[0024] Hereinafter, embodiments of the present invention will be described with reference to the accompanying drawings. However, it should be understood that these descriptions are exemplary only and are not intended to limit the scope of the present invention. In the following detailed description, for ease of explanation, many specific details are set forth to provide a comprehensive understanding of embodiments of the present invention. However, it is apparent that one or more embodiments may also be implemented without these specific details. In addition, in the following description, descriptions of known structures and technologies are omitted to avoid unnecessary confusion of the concept of the present invention.

[0025] The terms used herein are only for describing specific embodiments and are not intended to limit the present invention. The terms "comprise," "include," etc. used herein indicate the presence of features, steps, operations, and / or components, but do not exclude the presence or addition of one or more other features, steps, operations, or components.

[0026] All terms used herein (including technical and scientific terms) have the meanings commonly understood by those skilled in the art unless otherwise defined. It should be noted that the terms used herein should be interpreted as having a meaning consistent with the context of this specification and should not be interpreted in an idealized or overly rigid manner.

[0027] An embodiment of the present invention provides a field-free on-orbit calibration method for a spaceborne photon-counting lidar, aiming to achieve high-precision geometric consistency registration between a multi-beam spaceborne photon-counting laser point cloud and a high-resolution optical remote sensing image, thereby supporting satellite image control point extraction and high-precision terrain mapping applications.

[0028] Figure 1 One of the flow charts of the field-free on-orbit calibration method of a spaceborne photon counting lidar according to an embodiment of the present invention is schematically shown. Figure 2 The second flowchart of the field-free on-orbit calibration method of a spaceborne photon counting lidar according to an embodiment of the present invention is schematically shown.

[0029] like Figure 1 and Figure 2As shown, the field-free on-orbit calibration method of the spaceborne photon counting lidar according to this embodiment may include steps S1 to S6.

[0030] Step S1: Acquire a high-precision airborne laser point cloud, perform ground point cloud filtering and denoising on the airborne laser point cloud in sequence, and extract a first structural feature from the denoised airborne laser point cloud.

[0031] It should be noted that airborne laser point clouds are three-dimensional laser point clouds, and their accuracy includes planar accuracy, elevation accuracy, and density. "High precision" in the present invention typically refers to a planar accuracy of 0.1-0.2 meters, an elevation accuracy of 3-5 cm, and a density of approximately 300 ± 50 points per square meter.

[0032] Furthermore, the above step S1 performs ground point cloud filtering and denoising processing on the airborne laser point cloud in sequence, including steps S11 to S12.

[0033] In step S11 , a cloth simulation filter (CSF) algorithm based on terrain fitting is used to remove ground point clouds from the airborne laser point cloud while retaining non-ground point clouds, wherein the non-ground point clouds represent regular artificial objects.

[0034] Regular artificial features such as buildings and road structures are often present in non-ground point clouds. By filtering the ground point cloud, the ground point cloud can be removed while retaining the non-ground point cloud. In other embodiments, for high-precision airborne laser point clouds, various other mature methods can be used to filter the ground point cloud.

[0035] Step S12: Eliminate local outlier noise points in the non-ground point cloud according to a preset number of neighborhood points.

[0036] After filtering the ground point cloud, some outlier noise points still exist in the non-ground point cloud. These noise points need to be removed to prevent interference with the subsequent extraction of regular ground feature points. A statistical analysis method is used to remove local outlier noise points from the non-ground point cloud. Since the point cloud contains multiple points, points in the non-ground point cloud that fall outside the standard deviation range are removed as noise points using a pre-set number of neighborhood points.

[0037] Furthermore, after the above-mentioned step S1 performs ground point cloud filtering and denoising on the airborne laser point cloud in sequence, steps S13 to S15 are also included.

[0038] In step S13 , the denoised airborne laser point cloud is divided into a plurality of clusters using a Euclidean clustering algorithm, wherein each cluster contains a plurality of point clouds.

[0039] After noise point removal in step S12, non-ground point clouds, such as building point clouds, appear as independent clusters. Using the Euclidean clustering algorithm, these independent clusters can be segmented into multiple clusters, thereby extracting multiple regular ground feature regions. Each cluster can be understood as an independent point set containing multiple point clouds.

[0040] Step S14: for each cluster, a random sample consistency algorithm is used to fit a plane equation to the local point cloud in the cluster to form a plane point set, each plane point set containing multiple point clouds.

[0041] After segmentation, a cluster may still contain point cloud data for objects such as buildings, vehicles, and vegetation. For each cluster, the Random Sample Consensus (RANSAC) algorithm is used to fit a plane equation to the local point cloud within the cluster for plane modeling. Each cluster is then transformed into a plane point set. This plane point set can, for example, be a point cloud of a building rooftop, facilitating the subsequent matching of airborne laser point cloud features with optical remote sensing imagery.

[0042] Step S15 , screening out at least one valid plane point set from the multiple plane point sets of the multiple clusters, wherein the valid plane point set represents a roof plane point set with a regular shape.

[0043] The multiple plane point sets obtained through the above plane fitting process contain a large number of plane fragments. These plane point sets have the following characteristics: 1) Some are large planes corresponding to real building roofs; 2) Some are small, fragmented planes caused by noise, vegetation tops, and occlusion residues, lacking clear geometric regularity. To ensure the accuracy of subsequent extraction of regular feature boundary segment features and corner point features, valid plane point sets can be further filtered to remove meaningless small pieces of noise.

[0044] Furthermore, the above step S15 may also include steps S151 and S152.

[0045] Step S151: For each plane point set, calculate the number of points, area and rectangularity of the plane point set, where the number of points is the number of point clouds contained in the plane point set. , the area is the convex hull area of ​​the plane point set projected onto the fitting plane, and the rectangularity is the ratio of the convex hull area to the area of ​​the minimum circumscribed rectangle of the plane point set projected onto the fitting plane.

[0046] Traverse each plane point set in the multiple plane point sets obtained by fitting, and calculate the three properties of the plane point set: number of points, area, and rectangularity. When calculating the area, the multiple point clouds contained in the plane point set can be projected onto the fitting plane, and the convex hull (Convex Hull) can be generated in two-dimensional space to calculate the area of ​​the convex hull. The rectangularity is defined as the ratio of the convex hull area to the area of ​​the minimum circumscribed rectangle, which can be expressed as :

[0047]

[0048] in, is the convex hull area, is the minimum enclosing rectangle area.

[0049] In step S152 , the number, area, and rectangularity of the point cloud are compared with set thresholds respectively, and at least one valid plane point set is selected from the multiple plane point sets based on the comparison results.

[0050] Since planes with very small areas can be directly identified as fragmented planes, small noisy planes typically have very few points. A rectangularity close to 1 indicates a regular shape, while a smaller value indicates a fragmented shape. Therefore, corresponding thresholds are set for the number of points, area, and rectangularity, such as the point count threshold, area threshold, and rectangularity threshold. Attributes are compared one by one according to the set thresholds, and a valid plane point set that meets all the rules is retained, while fragmented planes are eliminated. This valid plane point set retains the regularly shaped roof point cloud information.

[0051] For example, among multiple plane point sets, fragmented plane point sets with a point count lower than a point count threshold, an area smaller than an area threshold, or a rectangularity smaller than a rectangularity threshold can be eliminated, thereby screening out at least one valid plane point set.

[0052] After filtering the valid planar point sets in step S152, a collection of regularly shaped rooftop point clouds is obtained. These point clouds possess distinct geometric structural features, including linear and corner features. To support the subsequent high-precision registration of airborne laser point clouds with high-resolution optical remote sensing imagery, stable linear and corner features must be extracted from these points.

[0053] Furthermore, in the above step S1, the first structural feature includes multiple first boundary line segments and multiple first corner points; extracting the first structural feature from the denoised airborne laser point cloud may include steps S161 and S162.

[0054] Step S161: For each valid plane point set in at least one valid plane point set, project the local coordinates of the valid plane point set onto the best fit principal plane according to the plane equation of the valid plane point set to form a two-dimensional projection point set; extract multiple boundary points of each two-dimensional projection point set, and determine a first boundary line segment from the multiple boundary points.

[0055] For each valid plane point set, the local coordinates of the valid plane point set are projected onto the best-fit principal plane according to the plane equation fitted in step S14 above, forming a two-dimensional projected point set to facilitate linear feature detection. The best-fit principal plane is different from the fitted plane equation and is an approximate horizontal plane formed after fine-tuning the plane equation.

[0056] use The algorithm extracts multiple boundary points from each two-dimensional projection point set, ensuring that each boundary point closely encloses the building edge while avoiding excessive concavity. Preferably, the multiple boundary points are arranged in spatial order, and the local direction between two adjacent boundary points is calculated. If the local direction changes within a certain tolerance, the two boundary points are demarcated as belonging to the same line segment; if the local direction changes suddenly beyond a threshold, the two boundary points are determined to be segment separators. In this way, the first boundary segment can be determined from the multiple boundary points.

[0057] For each extracted first boundary segment, record the following characteristic parameters: 1) start and end coordinates; 2) segment length; 3) segment orientation. Preferably, valid segments with lengths greater than a set threshold are retained, while short segments are discarded.

[0058] Step S162: At the intersection of any two first boundary line segments, detect the angle between the two boundary line segments. If the angle is at a preset regular angle, remap the plane coordinates of the intersection to three-dimensional space, restore the elevation coordinates of the intersection according to the plane equation, and use the restored intersection as the first corner point.

[0059] The angle is a preset regular angle, which can be close to a right angle (80°-100°) or a regular angle (such as 45° or 135°). For the intersection of the two first boundary segments, the precise plane coordinates of the intersection are first determined. The intersection is then remapped to three-dimensional space. The corresponding elevation coordinates are restored based on the fitted plane equation. The intersection of the plane coordinates and the elevation coordinates is then used as the first corner point. This first corner point can be, for example, a corner of a house.

[0060] In this way, the first boundary line segment and the first corner point are extracted from each valid plane point set for subsequent registration of the high-precision airborne laser point cloud with the high-resolution optical remote sensing image.

[0061] Step S2: acquiring a high-resolution optical remote sensing image, preprocessing the optical remote sensing image, and extracting a second structural feature from the preprocessed optical remote sensing image.

[0062] It should be noted that the “high resolution” in the embodiments of the present invention generally refers to a resolution of the optical remote sensing image of 0-0.5 meters, such as 0.1 meters or 0.5 meters.

[0063] High-resolution optical remote sensing images have excellent geometric structure representation capabilities, clearly displaying the boundary contours and corner features of regular artificial features (such as houses and roads). To achieve precise registration between optical remote sensing images and airborne laser point clouds, this step extracts linear boundary features and corner features from the high-resolution optical remote sensing images that correspond to the point cloud structure.

[0064] Furthermore, in the above step S2, pre-processing the optical remote sensing image may include steps S211 to S212.

[0065] Step S211 , using a median filter algorithm and a bilateral filter algorithm in sequence to reduce noise on the optical remote sensing image.

[0066] Using median filtering to reduce noise in optical remote sensing images can smooth out small noise in the image, suppress local texture interference, and maintain clear edge contours. Bilateral filtering can further smooth texture areas and maintain edge structure.

[0067] Step S212 : performing histogram equalization processing on the optical remote sensing image after noise reduction to enhance the brightness and contrast of the optical remote sensing image.

[0068] Histogram equalization processing can improve the brightness and contrast of optical remote sensing images, enhance the grayscale difference between the boundaries of regular objects (such as buildings) and the background, and ensure that the boundaries of regular objects are strengthened in edge detection.

[0069] Furthermore, in the above step S2, the second structural feature includes a plurality of second boundary line segments and a plurality of second corner points; extracting the second structural feature from the pre-processed optical remote sensing image may include steps S22 to S23.

[0070] Step S22 : extracting a second boundary line segment from the pre-processed optical remote sensing image using a two-dimensional Hough transform algorithm.

[0071] The second boundary line segment is extracted from the preprocessed optical remote sensing image using a two-dimensional Hough transform algorithm.

[0072] Optionally, for the multiple extracted second boundary line segments, the high and low dual thresholds of the slope and the high and low dual thresholds of the intercept can be adaptively set based on image statistics, and the features of some boundary line segments that meet both high and low dual thresholds can be screened out from the multiple second boundary line segments to ensure that the edges of the extracted second boundary line segments are coherent and accurate.

[0073] The starting and ending pixel coordinates, direction angle, and segment length of each second boundary line segment are recorded. Optionally, short line segments with lengths below a set threshold and non-structural boundary line segments can be discarded.

[0074] Step S23: At the intersection of any two second boundary line segments, detect the angle between the two second boundary line segments. If the angle is within a preset regular angle, use the intersection as a second corner point.

[0075] The included angle is at a preset regular angle, which can be close to a right angle (80°~100°) or a regular angle (such as 45°, 135°). In this case, the intersection point can be extracted as the second corner point.

[0076] It should be noted that there is no strict order between the above step S1 and the above step S2. The two steps can be performed simultaneously or sequentially, and the order is not specifically limited.

[0077] Step S3: registering the first structural feature and the second structural feature after unifying the spatial reference, fusing the registered second structural feature with the first structural feature to construct a fusion feature set.

[0078] To achieve high-precision registration of the airborne laser point cloud and the high-resolution optical remote sensing image, given that they share the same geographic coordinate system, the first structural features extracted in step S1 and the second structural features extracted in step S2 must be spatially matched. This embodiment of the present invention uses boundary line segment features and corner point features as matching primitives.

[0079] Furthermore, the above step S3 may include steps S31 to S35.

[0080] Step S31 : Projecting the first structural feature to the image coordinate system where the second structural feature is located.

[0081] The first structural features extracted from the airborne laser point cloud, including multiple first boundary segments (such as roof edges) and multiple first corner points (such as corner points), are projected into the image coordinate system of the optical remote sensing image using a preliminary geographic coordinate projection relationship. The image coordinate system is a plane coordinate system.

[0082] For example, a known initial coarse alignment model (such as unit conversion or ground reference plane simplification) can be used to convert the 3D points of the first structural feature ( ) is mapped to the plane pixel coordinates in the image coordinate system ( ).

[0083] Optionally, a spatial index structure may be established for the second structural features extracted from the optical remote sensing image, including multiple second boundary line segments and multiple second corner points, to support fast query of a candidate feature set within a specific spatial range.

[0084] Step S32 : matching the projected first boundary line segment and the second boundary line segment according to the angle difference and the line segment length ratio to form a line segment matching set.

[0085] For each first boundary line segment after projection, first search for a set of second boundary line segments with similar positions from multiple second boundary line segments in the optical remote sensing image, calculate the direction angle between the first boundary line segment after projection and any second boundary line segment in the set, and filter out second boundary line segments with similar directions from the set according to the set angle threshold. For example, the difference between the direction angles of the two can be set to be less than the angle threshold. , threshold The value can be 10°.

[0086] Then, compare the line segment length ratios and filter out candidate line segments with large length differences. The second boundary line segments that meet both the angle difference and line segment length ratio and the projected first boundary line segments are recorded as line segment matching pairs to form a line segment matching set. , which can be expressed as:

[0087]

[0088] in, Indicates the The first boundary segment after projection; Indicates the A second boundary segment.

[0089] Step S33 : matching the projected first corner point and the second corner point according to the spatial distance and the direction between two adjacent boundary line segments to form a corner point matching set.

[0090] For each first corner point after projection, first find the set of second corner points with the closest spatial distance from multiple second corner points of the optical remote sensing image, calculate the direction between the two adjacent first boundary segments where the first corner point is located after projection, and the direction between the two adjacent second boundary segments where the second corner point is located, and determine whether the two directions form the same angular structure such as "L-shaped" or "cross-shaped". If so, the first corner point after projection and the second corner point are recorded as a corner point matching pair to form a corner point matching set. , expressed as:

[0091]

[0092] in, Indicates the The first corner point after projection; Indicates the The second corner point. is the 3D corner point of the airborne laser point cloud, is the two-dimensional corner point of the optical remote sensing image.

[0093] Step S34: assign each second corner point in the corner point matching set to the corresponding first corner point to form an assigned first corner point.

[0094] For each pair of matching pairs in the corner matching set Perform mapping operations. The coordinates of the two-dimensional corner points of the optical remote sensing image are As additional attribute values ​​assigned to the corresponding 3D corner points of the airborne laser point cloud , the three-dimensional corner point features of the airborne laser point cloud are expanded to form a five-dimensional vector as shown below:

[0095]

[0096] in, Represents the first A five-dimensional vector of three-dimensional corner points.

[0097] In this way, the three-dimensional spatial position of the airborne laser point cloud and the planar position of the optical remote sensing image can be integrated.

[0098] In step S35 , the multiple assigned first corner points are aggregated to form a fusion feature set.

[0099] For example, construct the following fusion feature set:

[0100]

[0101] in, Indicates the number of matching pairs in the corner point matching set.

[0102] At this point, after the airborne laser point cloud and the high-resolution optical remote sensing image are registered, the fused feature set is obtained.

[0103] Step S4: performing plane translation and plane rotation on the airborne laser point cloud based on the fusion feature set to correct two plane coordinates in the airborne laser point cloud except the elevation coordinates, thereby generating a corrected airborne laser point cloud.

[0104] This step performs planar translation and planar rotation on the airborne laser point cloud based on the fusion feature set to correct the X and Y coordinate components in the airborne laser point cloud and keep the Z (elevation) coordinate component unchanged.

[0105] Furthermore, the above step S4 may include steps S41 to S44.

[0106] Step S41: Based on the fusion feature set, the predefined plane translation vector and the small-angle rotation matrix, a correction model of the airborne laser point cloud before and after correction is established, and the elevation coordinates of the airborne laser point cloud before and after correction are kept unchanged.

[0107] The 3D coordinates of the airborne laser point cloud in the fusion feature set are used Indicates that the two-dimensional plane coordinates of the optical remote sensing image corresponding to the airborne laser point cloud are expressed as express.

[0108] It should be noted that The coordinates of the two-dimensional corner points of the aforementioned optical remote sensing image are The correction target of the correction model is to transform the airborne laser point cloud in the fusion feature set into 2D coordinates corrected to the image coordinate system , making The two-dimensional plane coordinates of the optical remote sensing image corresponding to the airborne laser point cloud in the fusion feature set As close as possible.

[0109] Definition: The correction model for the rigid transformation from the airborne laser point cloud to the image coordinate system of the optical remote sensing image is:

[0110]

[0111] in, is the small-angle rotation matrix, is a two-dimensional plane translation vector, For the overall translation compensation to be solved, the translation amount in the elevation direction is set to 0.

[0112] Small angle rotation is defined as: a small angle rotation around the Z axis , this small angle satisfies Based on this, the small angle rotation matrix It can be approximately expanded as:

[0113]

[0114] Therefore, the 3D coordinates of the corrected airborne laser point cloud are for:

[0115]

[0116] in, is the small rotation angle compensation to be solved, Remain unchanged.

[0117] Step S42: constructing a first objective function of the correction model using the plane reprojection error.

[0118] It is The position of the corrected airborne laser point cloud, It is The position of the optical remote sensing image corresponding to the airborne laser point cloud in the fusion feature set. Definition: The first objective function of the correction model is the plane reprojection error , which is expressed as follows:

[0119]

[0120] Wherein, N is the number of matching pairs in the fused feature set, which is also the number of matching pairs in the corner point matching set in the aforementioned step S35.

[0121] The two-dimensional coordinates of the corrected airborne laser point cloud use After substituting, we get:

[0122]

[0123] Step S43: solving the plane translation vector and the small-angle rotation matrix by minimizing the first objective function.

[0124] To minimize the planar reprojection error , respectively Taking the partial derivatives and setting them to zero, we obtain the following normal equations:

[0125]

[0126] by For example, Find the partial derivative:

[0127]

[0128] Setting it to zero gives:

[0129]

[0130] After finishing, we can get:

[0131]

[0132] By analogy, we can also Find the partial derivative and set it to zero, and then sort out the other two equations. The whole can be expressed as a small linear equation system, which can be solved using the least squares method to obtain The optimal solution of .

[0133] Step S44 , transforming the airborne laser point cloud using the plane translation vector and the small-angle rotation matrix to obtain a corrected airborne laser point cloud.

[0134] Use the following formula to The optimal solution is applied to all airborne laser point clouds in step S1 above to improve the overall plane accuracy of the point cloud while keeping the elevation information unaffected:

[0135]

[0136] in, It is The 3D coordinates of the calibrated airborne laser point cloud.

[0137] At this point, the generated corrected airborne laser point cloud can be simplified as .

[0138] Step S5: Acquire ranging data from a spaceborne photon counting lidar, perform atmospheric delay correction and tidal correction on the ranging data in sequence, filter the tidal-corrected ranging data to form a spaceborne photon counting point cloud, and align the spaceborne photon counting point cloud and the corrected airborne laser point cloud with a unified spatial reference to form a registered spaceborne photon counting point cloud.

[0139] This step aims to spatially register the spaceborne photon counting point cloud with the corrected airborne laser point cloud generated in step S4 above, thereby establishing a spatial correspondence between the two types of laser sensors, and providing basic support for the subsequent correction of the pointing error and system ranging error of the spaceborne photon counting lidar.

[0140] It should be noted that the "spaceborne photon counting lidar" in the embodiment of the present invention can be single-beam or multi-beam, and the number of beams can reach six or more.

[0141] Figure 3 A three-dimensional cloud map of observation data of a space-borne photon counting lidar according to an embodiment of the present invention is schematically shown.

[0142] The embodiment of the present invention constructs the observation data of multi-beam spaceborne photon counting lidar based on the characteristics of polar orbit satellites and the original airborne lidar point cloud data, such as Figure 3 As shown, there are six beams in total. The simulation introduces a known laser pointing error ( ) and system ranging error ( meters), and taking into account atmospheric delays and tidal effects.

[0143] Figure 4 The elevation comparison diagram of a satellite-borne photon counting point cloud and an airborne laser point cloud before calibration according to an embodiment of the present invention is schematically shown.

[0144] Comparative analysis of the elevation of the six-beam satellite-borne photon counting point cloud and the airborne laser point cloud before calibration Figure 4 As shown. Figure 4 It can be seen that there is a certain difference in the elevation of the two data due to the laser pointing error and the system ranging error.

[0145] Furthermore, in the above step S5, atmospheric delay correction and tidal correction are performed on the ranging data of the spaceborne photon counting lidar to eliminate the ranging deviation caused by atmospheric delay and earth tidal effect during the laser ranging process.

[0146] For example, the Saastamoinen model can be used to calculate atmospheric interference delay and wet term delay The input data used include atmospheric pressure, temperature and water vapor pressure obtained from the ERA5 dataset. Atmospheric delay correction as follows:

[0147]

[0148] Raw height measurement for each photon point , perform atmospheric delay correction for laser ranging:

[0149]

[0150] in, It is the ranging data after atmospheric delay correction.

[0151] For example, GOT4.8 and IERS global tidal models can be used for tidal corrections, which include: solid tide , tide , Extreme Tide Corrected, the calculation formula is as follows:

[0152]

[0153] in, is the tidal correction.

[0154] For each photon point, according to its latitude and longitude (lat, lon) and timestamp , interpolate and calculate the above tidal correction. To make tidal corrections:

[0155]

[0156] in, The distance measurement data after tide correction.

[0157] Furthermore, in the above step S5, the tidal-corrected ranging data is filtered to form a satellite-borne photon counting point cloud.

[0158] For example, the tidal-corrected ranging data is first calculated using the LiDAR equation to form the initial spaceborne photon counting point cloud data. Then, the initial spaceborne photon counting point cloud data is filtered based on the linear feature photon data filtering algorithm (LFPSE) to remove noise and low-confidence photon points and retain signal photons, thereby forming a spaceborne photon counting point cloud distributed along the track. .

[0159] Furthermore, in the above step S5, the corrected airborne laser point cloud and the spaceborne photon counting point cloud are aligned with a unified spatial reference to form a registered spaceborne photon counting point cloud, which may include steps S53 to S55.

[0160] Step S53 , unifying the spatial coordinate references of the calibrated airborne laser point cloud and the satellite-borne photon counting point cloud into the geocentric coordinate system.

[0161] For example, the spaceborne photon counting point cloud Compared with the rectified airborne laser point cloud The spatial coordinate reference is unified to the geocentric coordinate system where the satellite-borne photon counting point cloud is located, such as the WGS84 coordinate system.

[0162] Furthermore, after step S53, the airborne laser point cloud and the satellite-borne photon counting point cloud, after the spatial reference is unified, are post-processed. For example, the satellite-borne photon counting point cloud can be divided into a plurality of equally spaced segments according to the trajectory direction in the geocentric coordinate system, and the corrected airborne laser point cloud can be cropped according to the trajectory range of the satellite-borne photon counting point cloud in each equally spaced segment.

[0163] Specifically, the satellite-borne photon counting point cloud can be divided into multiple equally spaced segments, each 500 meters long, for subsequent segmented registration. This approach is applicable to areas with large-scale terrain variations and improves local fitting quality. Within each equally spaced segment, the airborne laser point cloud is clipped based on the trajectory range of the satellite-borne photon counting point cloud, retaining only the overlapping areas with the satellite-borne photon counting point cloud, improving computational efficiency and matching accuracy.

[0164] In step S54, the centroid coordinates of the corrected airborne laser point cloud and the satellite-borne photon counting point cloud are calculated respectively, and the corrected airborne laser point cloud and the satellite-borne photon counting point cloud are roughly registered in the elevation profile according to the difference between the two centroid coordinates.

[0165] Computing spaceborne photon counting point clouds and rectified airborne laser point cloud The center of gravity coordinates of , get the difference between the two barycentric coordinates . Use the difference as the initial translation vector and and Perform initial translation and extract cross-section elevation distribution profile to and Complete the coarse registration.

[0166] Step S55: For each satellite point in the satellite-borne photon counting point cloud, search for the nearest neighbor point corresponding to the satellite point in the calibrated airborne laser point cloud according to the distance threshold; establish a second objective function regarding the translation vector and the rotation matrix based on each satellite point and the corresponding nearest neighbor point; solve the translation vector and the rotation matrix by minimizing the second objective function; and transform the satellite-borne photon counting point cloud using the translation vector and the rotation matrix to obtain the registered satellite-borne photon counting point cloud.

[0167] Specifically, for the coarse registration Each satellite point in , according to the distance threshold in the airborne laser point cloud after coarse registration Search for the satellite point The corresponding nearest neighbor , filtered by a distance threshold to avoid false matches.

[0168] Satellite and the matched nearest neighbor Based on this, we establish the following rotation matrix and translation vectors The second objective function:

[0169]

[0170] By minimizing the second objective function, the rotation matrix can be solved and translation vectors The optimal solution of .

[0171] Then, the optimal solution is used to calculate the spaceborne photon counting point cloud Each satellite point Apply the following transformation model:

[0172]

[0173] in, Indicates the satellite point Position after registration.

[0174] So far, for each satellite point , you can refer to the airborne laser point cloud Nearest neighbor Align it and obtain the satellite-borne point after alignment The registered satellite photon counting point cloud can be simplified as .

[0175] To further improve the registration accuracy, after the above step S55, the following step may be performed: the point set after the registration between the onboard laser point cloud and the satellite-borne photon counting point cloud is calculated. The Euclidean residual of all matching points in :

[0176]

[0177] in, Represents the i-th pair of matching points in the point set after registration.

[0178] According to the Euclidean residual and the preset accuracy index, the abnormal satellite points in the satellite-borne photon counting point cloud are eliminated to obtain the registered satellite-borne photon counting point cloud.

[0179] For example, by using statistical accuracy indicators such as the root mean square error (RMSE), maximum deviation, and standard deviation, satellite-borne points in the satellite-borne photon counting point cloud with errors exceeding three times the standard deviation are eliminated to avoid extreme points affecting subsequent error modeling, and weighted nearest neighbor iterative optimization is used to suppress the impact of low-quality areas.

[0180] Step S6: constructing a radar equation for the spaceborne photon counting lidar, and determining the pointing error and system ranging error of the spaceborne photon counting lidar based on the aligned spaceborne photon counting point cloud and the radar equation.

[0181] Furthermore, the above step S6 may include steps S61 to S63.

[0182] Step S61: construct radar equations for multiple satellite-borne points of a satellite-borne photon counting laser radar.

[0183] For multiple photon points of spaceborne photon counting lidar (also called satellite point), use Indicates the The one-way flight distance of a photon. The unit vector of the emission direction is According to the azimuth of the launch direction and pitch angle Further defined as:

[0184]

[0185] The measured position of a spaceborne photon counting laser can be expressed by the following radar equation:

[0186]

[0187] in, is the satellite center position at the moment of photon emission, is the satellite placement matrix at the time of photon emission, is the measured distance of the i-th photon.

[0188] Step S62, defining the pointing error and system ranging error of the spaceborne photon counting lidar; and correcting the radar equation according to the pointing error and the system ranging error.

[0189] In an embodiment of the present invention, the pointing error includes an azimuth error and pitch angle error , the system ranging error is also the systematic deviation of ranging The corrected unit direction vector after the error is introduced for:

[0190]

[0191] Among them, small disturbance Expands to:

[0192]

[0193] The partial derivatives are:

[0194]

[0195] Corrected measured distance for:

[0196]

[0197] Therefore, the measurement position after introducing error correction into the above radar equation is for:

[0198]

[0199] Should The formula is also the modified radar equation.

[0200] Step S63 , solving the pointing error and the system ranging error according to the corrected radar equation and the registered satellite-borne photon counting point cloud.

[0201] Define the registered spaceborne photon counting point cloud The measured position after error correction The residual between for:

[0202]

[0203] Substitute the pointing error and system ranging error into the equation and expand the residual approximation:

[0204]

[0205] Construct the residual matrix for all satellite points (N photon points) and use the Jacobian matrix express:

[0206]

[0207] The residual can be further expressed as:

[0208]

[0209] in, is a 3N×1 residual vector, is a 3N×3 Jacobian matrix, is the error parameter vector to be estimated.

[0210] Solve using the least squares method:

[0211]

[0212] At this point, the pointing error and system ranging error can be obtained .

[0213] Then, you can use Measuring the position of a spaceborne photon-counting laser Make corrections, update the direction angle (azimuth, pitch angle), update the ranging parameters, recalculate the three-dimensional coordinates of the satellite-borne photon counting point cloud, and obtain the updated measurement position .

[0214] The field-free on-orbit calibration method for spaceborne photon counting lidar proposed in the embodiment of the present invention is used to perform error modeling and calculate the direction angle (azimuth, elevation) and ranging parameters. The results show that this method achieves extremely high-precision estimation of the three error parameters. The estimated value is -4.003 arc seconds, and the pitch angle Estimated to be -2.002 arc seconds, systematic ranging error The estimated value is -0.400 meters, which is almost completely consistent with the actual set values ​​of 4 arc seconds, 2 arc seconds, and 0.4 meters, respectively. The errors are only within the order of 0.003 arc seconds and 1 mm, verifying the high robustness and error sensitivity of this method.

[0215] Figure 5 The elevation comparison diagram of the spaceborne photon counting point cloud and the airborne laser point cloud after calibration according to an embodiment of the present invention is schematically shown.

[0216] like Figure 5 As shown in the figure, the three-dimensional coordinates of the satellite-borne photon counting point cloud recalculated after updating the direction angle (azimuth, pitch angle) and ranging parameters are highly consistent with the point cloud data of the airborne lidar.

[0217] This method has been fully validated for its applicability in multi-beam spaceborne photon-counting lidar. Relying solely on the spaceborne platform's own data and auxiliary imagery, it can achieve high-precision estimation and correction of the multi-beam photon laser system's ranging and directional errors without relying on ground-based calibration sites or additional equipment. This method demonstrates strong engineering feasibility and widespread application, particularly in remote sensing missions where deploying ground-based calibrators is difficult.

[0218] In summary, the embodiments of the present invention provide a method for field-free on-orbit calibration of satellite-borne photon counting lidar, which combines the high-precision airborne laser point cloud with the high-resolution optical remote sensing image for joint registration, and combines the planar positioning advantage of the high-resolution optical remote sensing image with the elevation accuracy advantage of the airborne laser point cloud to first improve the plane positioning accuracy of the airborne laser point cloud; then, by matching the satellite-borne photon counting point cloud with the corrected airborne laser point cloud, the pointing error and system ranging error of the satellite-borne photon counting lidar are corrected, thereby realizing field-free on-orbit geometric calibration of the satellite-borne photon counting lidar. This method avoids dependence on ground calibration field laser detectors, has the advantages of strong adaptability, high calibration accuracy, and a wide range of applications, and provides a new technical means for improving the field-free on-orbit calibration of satellite-borne photon counting lidar altimetry systems.

[0219] The flowcharts and block diagrams in the accompanying drawings illustrate the possible implementation architecture, functions and operations of the systems, methods and computer program products according to various embodiments of the present invention. In this regard, each box in the flowchart or block diagram can represent a module, program segment, or a part of code, and the above-mentioned module, program segment, or a part of code contains one or more executable instructions for implementing the specified logical function. It should also be noted that in some alternative implementations, the functions marked in the box can also occur in an order different from that marked in the accompanying drawings. For example, two boxes represented in succession can actually be executed substantially in parallel, and they can sometimes be executed in the opposite order, depending on the functions involved. It should also be noted that each box in the block diagram or flowchart, and the combination of boxes in the block diagram or flowchart, can be implemented with a dedicated hardware-based system that performs the specified function or operation, or can be implemented with a combination of dedicated hardware and computer instructions.

[0220] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly specifying the number of the technical features indicated. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include one or more of the features. In the description of the present invention, "plurality" means at least two, such as two, three, etc., unless otherwise specifically defined. Furthermore, the word "a" or "an" preceding an element does not exclude the presence of a plurality of such elements.

[0221] It will be understood by those skilled in the art that the features described in the various embodiments of the present invention may be combined and / or coupled in various ways, even if such combinations or couplings are not explicitly described in the present invention. In particular, the features described in the various embodiments of the present invention may be combined and / or coupled in various ways without departing from the spirit and teachings of the present invention. All such combinations and / or couplings fall within the scope of the present invention.

[0222] The above describes embodiments of the present invention. However, these embodiments are for illustrative purposes only and are not intended to limit the scope of the present invention. Although each embodiment has been described separately above, this does not mean that the measures in each embodiment cannot be advantageously used in combination. Without departing from the scope of the present invention, those skilled in the art may make various substitutions and modifications, which should all fall within the scope of the present invention.

Claims

1. A field-free on-orbit calibration method for a spaceborne photon counting lidar, characterized in that: include: Step S1, obtaining a high-precision airborne laser point cloud, performing ground point cloud filtering and denoising on the airborne laser point cloud in sequence, and extracting a first structural feature from the denoised airborne laser point cloud; Step S2, acquiring a high-resolution optical remote sensing image, preprocessing the optical remote sensing image, and extracting a second structural feature from the preprocessed optical remote sensing image; Step S3, registering the first structural feature and the second structural feature after unifying the spatial reference, fusing the registered second structural feature with the first structural feature to construct a fused feature set; Step S4, performing a plane translation and a plane rotation on the airborne laser point cloud based on the fusion feature set to correct two plane coordinates in the airborne laser point cloud except the elevation coordinates, thereby generating a corrected airborne laser point cloud; Step S5: acquiring ranging data from a spaceborne photon counting lidar, sequentially performing atmospheric delay correction and tidal correction on the ranging data, filtering the tidal-corrected ranging data to form a spaceborne photon counting point cloud; and registering the spaceborne photon counting point cloud and the corrected airborne laser point cloud with a unified spatial reference to form a registered spaceborne photon counting point cloud. Step S6: construct a radar equation for the spaceborne photon counting lidar, and determine the pointing error and system ranging error of the spaceborne photon counting lidar based on the aligned spaceborne photon counting point cloud and the radar equation.

2. The method according to claim 1, characterized in that In step S1, the airborne laser point cloud is sequentially subjected to ground point cloud filtering and denoising, including: Step S11, using a cloth filtering algorithm based on terrain fitting to remove ground point clouds from the airborne laser point cloud and retain non-ground point clouds, wherein the non-ground point clouds represent regular artificial objects; Step S12: removing local outlier noise points in the non-ground point cloud according to a preset number of neighborhood points.

3. The method according to claim 1, characterized in that In step S1, after the airborne laser point cloud is subjected to ground point cloud filtering and denoising processing, the following steps are further included: Step S13, using a Euclidean clustering algorithm to segment the denoised airborne laser point cloud into multiple clusters, where each cluster contains multiple point clouds; Step S14: for each of the clusters, using a random sample consistency algorithm to fit a plane equation to the local point cloud in the cluster to form a plane point set, each of the plane point sets containing multiple point clouds; Step S15 , selecting at least one valid plane point set from the plurality of plane point sets of the plurality of clusters, wherein the valid plane point set represents a roof plane point set with a regular shape.

4. The method according to claim 3, characterized in that The step S15 includes: Step S151: For each of the planar point sets, calculate the number of points, area, and rectangularity of the planar point set, where the number of points is the number of point clouds contained in the planar point set, the area is the area of ​​the convex hull of the planar point set projected onto the fitting plane, and the rectangularity is the ratio of the area of ​​the convex hull to the area of ​​the minimum circumscribed rectangle of the planar point set projected onto the fitting plane; Step S152 : Compare the number, area, and rectangularity of the point cloud with set thresholds respectively, and select at least one valid plane point set from the plurality of plane point sets based on the comparison results.

5. The method according to claim 3, characterized in that In step S1, the first structural feature includes a plurality of first boundary line segments and a plurality of first corner points; and extracting the first structural feature from the denoised airborne laser point cloud includes: Step S161: for each valid plane point set in the at least one valid plane point set, project the local coordinates of the valid plane point set onto a best-fit principal plane according to the plane equation of the valid plane point set to form a two-dimensional projected point set; extract multiple boundary points from each of the two-dimensional projected point sets, and determine the first boundary line segment from the multiple boundary points; Step S162: At the intersection of any two of the first boundary line segments, detect the angle between the two first boundary line segments. If the angle is at a preset regular angle, remap the plane coordinates of the intersection to three-dimensional space, restore the elevation coordinates of the intersection according to the plane equation, and use the restored intersection as the first corner point.

6. The method according to claim 1, characterized in that In step S2, the second structural feature includes a plurality of second boundary line segments and a plurality of second corner points; and extracting the second structural feature from the pre-processed optical remote sensing image includes: Step S22, extracting the second boundary line segment from the pre-processed optical remote sensing image using a two-dimensional Hough transform algorithm; Step S23: At the intersection of any two second boundary line segments, detect the angle between the two second boundary line segments; if the angle is within a preset regular angle, use the intersection as the second corner point.

7. The method according to claim 6, characterized in that The step S3 comprises: Step S31, projecting the first structural feature to the image coordinate system where the second structural feature is located; Step S32, matching the projected first boundary line segment and the second boundary line segment according to the angle difference and the line segment length ratio to form a line segment matching set; Step S33, matching the projected first corner point and the second corner point according to the spatial distance and the direction between two adjacent boundary segments to form a corner point matching set; Step S34, assigning each second corner point in the corner point matching set to the corresponding first corner point to form an assigned first corner point; Step S35 , aggregating the multiple assigned first corner points to form the fused feature set.

8. The method according to claim 1, characterized in that The step S4 comprises: Step S41, establishing a correction model of the airborne laser point cloud before and after correction based on the fusion feature set, the predefined plane translation vector and the small-angle rotation matrix, and keeping the elevation coordinates of the airborne laser point cloud before and after correction unchanged; Step S42, constructing a first objective function of the correction model using a plane reprojection error; Step S43, solving the plane translation vector and small-angle rotation matrix by minimizing the first objective function; Step S44 : transforming the airborne laser point cloud using the plane translation vector and the small-angle rotation matrix to obtain a corrected airborne laser point cloud.

9. The method according to claim 1, characterized in that In step S5, the corrected airborne laser point cloud and the satellite-borne photon counting point cloud are aligned with each other in a unified spatial reference to form an aligned satellite-borne photon counting point cloud, including: Step S53, unifying the spatial coordinate references of the corrected airborne laser point cloud and the satellite-borne photon counting point cloud into a geocentric coordinate system; Step S54, calculating the centroid coordinates of the corrected airborne laser point cloud and the satellite-borne photon counting point cloud in the geocentric coordinate system, and performing a coarse registration of the corrected airborne laser point cloud and the satellite-borne photon counting point cloud in the elevation profile based on the difference between the two centroid coordinates; Step S55: for each satellite point in the satellite-borne photon counting point cloud, search for the nearest neighbor point corresponding to the satellite point in the corrected airborne laser point cloud according to a distance threshold; establish a second objective function regarding a translation vector and a rotation matrix based on each satellite point and the corresponding nearest neighbor point; solve the translation vector and the rotation matrix by minimizing the second objective function; and transform the satellite-borne photon counting point cloud using the translation vector and the rotation matrix to obtain a registered satellite-borne photon counting point cloud.

10. The method according to claim 1, characterized in that The step S6 comprises: Step S61, constructing radar equations for multiple satellite-borne points of the satellite-borne photon counting laser radar; Step S62, defining the pointing error and system ranging error of the spaceborne photon counting laser radar; and correcting the radar equation according to the pointing error and system ranging error; Step S63 , solving the pointing error and the system ranging error according to the corrected radar equation and the registered satellite-borne photon counting point cloud.

Images