Field-free on-orbit calibration method for space-borne photon counting lidar

Abstract

The application provides a field-free on-orbit calibration method for a spaceborne photon counting laser radar, and relates to the technical field of the spaceborne photon counting laser radar.The method is combined with the plane positioning advantage of high-resolution optical remote sensing images and the height accuracy advantage of airborne laser point clouds, and the plane positioning accuracy of the airborne laser point clouds is improved first; then the pointing error and the system ranging error of the spaceborne photon counting laser radar are corrected by matching the spaceborne photon counting point clouds with the corrected airborne laser point clouds, so that the field-free on-orbit geometric calibration of the spaceborne photon counting laser radar is realized. The method avoids the dependence on the ground calibration field laser detector, has the advantages of strong adaptability, high calibration accuracy and wide application range, and provides a new technical means for improving the field-free on-orbit calibration of the spaceborne photon counting laser radar height measurement system.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of spaceborne photon counting lidar, and particularly relates to a field-free on-orbit calibration method of a spaceborne photon counting lidar. BACKGROUND

[0002] As a new emerging means of earth observation, spaceborne photon counting lidar has the advantages of micro-pulse, high repetition frequency, multi-beam and small spot, and has been widely used in global topographic mapping, glacier monitoring and ocean exploration. However, due to the complex and changeable attitude change, orbit disturbance and the precision limitation of the equipment itself in the on-orbit environment, the spaceborne photon counting lidar point cloud data usually has system ranging error and pointing error, which directly affects the spatial positioning accuracy, and further restricts the high-precision mapping application based on the spaceborne photon counting lidar data. Therefore, it is of great significance to calibrate the laser beam pointing error and system ranging error of the spaceborne photon counting lidar on-orbit.

[0003] The traditional high-precision spaceborne lidar on-orbit calibration method mainly depends on the calibration field ground laser detector. However, the micro-pulse characteristics of the spaceborne photon counting lidar require the ground laser detector to have extremely high detection sensitivity; the high repetition frequency characteristics require the ground detector to have extremely high time synchronization accuracy; the multi-beam and small spot characteristics require the ground detector to have a large layout density. These characteristics cause the active calibration scheme based on the calibration field to require a lot of manpower, material resources and financial resources.

[0004] In view of the above problems, a new technical path is needed to solve the field-free on-orbit calibration problem of the spaceborne photon counting lidar with limited dependence on ground information. SUMMARY

[0005] In view of the above problems, the present application provides a field-free on-orbit calibration method of a spaceborne 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] The method constructs a complete laser error propagation model, which can establish a clear mathematical relationship between the actual space registration error and the system internal parameter deviation (such as azimuth angle error Δθ, pitch angle error Δφ, and distance measurement error Δd), and realize quantitative estimation and parameter correction of system error through least square method, effectively improving the geometric consistency and positioning accuracy of the laser height measurement system.

[0012] 3. Realize the consistency geometric calibration of multi-beam laser system

[0013] Compared with the traditional method of calibrating single beam or main beam, the method is suitable for multi-beam spaceborne photon counting laser radar system, can jointly model, error compensate and align the data of multiple laser beams in different transmission directions, so as to ensure the spatial coordination and system stability of the whole height measurement system under multi-view conditions.

[0014] 4. Construct multi-source data collaborative registration framework to enhance error constraint ability

[0015] The method first constructs a three-source joint space matching framework of "high-resolution image - airborne point cloud - spaceborne photon data", through layer-by-layer registration of image edge features, point cloud structure features and photon projection points, forms a multi-layer feature fusion mechanism of surface, line and point, enhances the stability and observation redundancy of error residual modeling, and greatly improves the reliability and result accuracy of error estimation.

[0016] 5. Strong universality, can be integrated into various orbital remote sensing tasks

[0017] The method does not depend on specific equipment or platform, can be adapted to multiple types of spaceborne photon counting laser radar systems, and can be embedded in post-processing ground system or data interpretation process, has good modular implementation ability and engineering generalization. In addition, the method can also be used for post-geometric correction and precision improvement for historical data, area images or laser data lacking ground control points. BRIEF DESCRIPTION OF DRAWINGS

[0018] The above and other objects, features and advantages of the present application will become more apparent from the following description of embodiments of the present application, taken in conjunction with the accompanying drawings, in which:

[0019] Figure 1 One of the flowcharts of the on-orbit calibration method of the spaceborne photon counting laser radar according to the embodiment of the present application is schematically shown;

[0020] Figure 2 The second flowchart of the on-orbit calibration method of the spaceborne photon counting laser radar according to the embodiment of the present application is schematically shown;

[0021] Figure 3A three-dimensional cloud of observation data of a spaceborne photon counting lidar according to an embodiment of the present application is schematically shown;

[0022] Figure 4 A height comparison plot of a spaceborne photon counting point cloud before calibration and an airborne laser point cloud according to an embodiment of the present application is schematically shown.

[0023] Figure 5 A height comparison plot of a spaceborne photon counting point cloud after calibration and an airborne laser point cloud according to an embodiment of the present application is schematically shown. DETAILED DESCRIPTION

[0024] Hereinafter, embodiments of the present application will be described with reference to the accompanying drawings. It is to be understood, however, that these descriptions are merely exemplary and are intended to illustrate the scope of the present application, not to limit it. In the following detailed description of the embodiments of the present application, numerous specific details are set forth in order to provide a thorough understanding of the embodiments of the present application. However, it will be apparent to one skilled in the art that the embodiments of the present application can be practiced without these specific details. In other instances, well-known structures and

[0025] The terms used herein are merely used to describe specific embodiments, and are not intended to limit the present application. The terms "include" and "have" and the like used herein indicate the presence of the 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 same meanings as those generally understood by those skilled in the art, unless otherwise defined. It should be noted that the terms used herein should be interpreted as having meanings consistent with the context of the present specification, and should not be interpreted in an idealized or excessively formal manner.

[0027] Embodiments of the present application provide a field-free on-orbit calibration method for a spaceborne photon counting lidar, aiming to achieve high-precision geometric consistency registration between multi-beam spaceborne photon counting laser point clouds and high-resolution optical remote sensing images, thereby supporting satellite image control point extraction and high-precision topographic mapping applications.

[0028] Figure 1 A flowchart of a field-free on-orbit calibration method for a spaceborne photon counting lidar according to an embodiment of the present application is schematically shown. Figure 2 A flowchart of a field-free on-orbit calibration method for a spaceborne photon counting lidar according to an embodiment of the present application is schematically shown.

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

[0030] In step S1, a high-precision airborne laser point cloud is obtained, and ground point cloud filtering and denoising are sequentially performed on the airborne laser point cloud, and a first structural feature is extracted from the airborne laser point cloud after denoising.

[0031] It should be noted that the airborne laser point cloud belongs to a three-dimensional laser point cloud, and its precision includes plane precision, height precision, and density. The "high precision" of the embodiment of the present application generally refers to a plane precision of 0.1-0.2 meters, a height precision of 3-5 cm, and a density of about 300±50 points per square meter.

[0032] Further, the above step S1 sequentially performs ground point cloud filtering and denoising on the airborne laser point cloud, including steps S11-S12.

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

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

[0035] In step S12, according to a pre-set number of neighborhood points, local outlier noise points in the non-ground point clouds are removed.

[0036] There are still some outlier noise points in the non-ground point clouds after the ground point cloud filtering, and such noise points need to be removed to avoid interference with the extraction of regular feature points in the subsequent steps. A statistical analysis-based method is used to remove local outlier noise points from the non-ground point clouds. Since the point cloud contains multiple points, by pre-setting the number of neighborhood points, points in the non-ground point clouds that are outside the standard deviation range are deleted as noise points.

[0037] Further, after the above step S1 sequentially performs ground point cloud filtering and denoising on the airborne laser point cloud, steps S13-S15 are further included.

[0038] In step S13, a Euclidean clustering algorithm is used to segment the airborne laser point cloud after denoising into multiple clusters, wherein each cluster contains multiple point clouds.

[0039] After the step S12 of removing the noise points, the non-ground point cloud such as the building point cloud is in an independent cluster state, and the non-ground point cloud in the independent cluster state can be segmented into multiple clusters by using the Euclidean clustering algorithm, so as to extract multiple regular ground object regions. Each cluster can be understood as an independent point set containing multiple point clouds.

[0040] In step S14, for each cluster, a plane equation is fitted to the local point cloud in the cluster by using a random sample consistency algorithm, to form a plane point set, and each plane point set contains multiple point clouds.

[0041] In the multiple clusters after segmentation, a cluster can still include point cloud data of objects such as buildings, vehicles, and vegetation. For each cluster, a plane equation is fitted to the local point cloud in the cluster by using a random sample consistency algorithm (RANSAC), to perform plane modeling, and each cluster forms a plane point set. The plane point set can be, for example, a building roof plane point cloud, which is beneficial to subsequent matching of airborne laser point cloud features and optical remote sensing image features.

[0042] In step S15, at least one effective plane point set is selected from the multiple plane point sets of the multiple clusters, wherein the effective plane point set represents a regular shape roof plane point set.

[0043] The multiple plane point sets obtained through the above plane fitting contain a large number of plane segments. These plane point sets have the following characteristics: 1) some are large-scale planes corresponding to real building roofs; and 2) some are fragmented small planes caused by noise, vegetation tops, and residual shielding, and have no obvious geometric regularity. To ensure the accuracy of subsequent extraction of regular ground object boundary line segment features and corner point features, the effective plane point set can be further selected to remove meaningless small piece-shaped noise.

[0044] Further, the step S15 can further include steps S151 to S152.

[0045] In step S151, for each plane point set, the number of points, the area, and the rectangularity of the plane point set are calculated, wherein 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 minimum circumscribed rectangle area of the plane point set projected onto the fitting plane.

[0046] The point number, area and rectangularity of each of the plurality of plane point sets obtained by the traversal fitting are sequentially calculated. When calculating the area, the plurality of points included in the plane point set can be projected to the fitting plane to generate a convex hull in a two-dimensional space, and the area of the convex hull is calculated. The rectangularity is defined as the ratio of the area of the convex hull to the area of the minimum circumscribed rectangle, and can be expressed as :

[0047]

[0048] wherein, is the area of the convex hull, is the area of the minimum circumscribed rectangle.

[0049] In step S152, the point cloud number, area and rectangularity are compared with the set threshold values respectively, and at least one effective plane point set is selected from the plurality of plane point sets according to the comparison results.

[0050] Since the area of the plane is too small, the plane can be directly determined as a fragmented plane. The noise plane with a small area usually has a small number of points, and the rectangularity close to 1 indicates a regular shape, and a small rectangularity indicates a broken shape. Therefore, corresponding threshold values are set for the point cloud number, area and rectangularity, such as a point number threshold value, an area threshold value and a rectangularity threshold value. According to the set threshold values, the effective plane point set that meets all the rules is selected and retained, and the fragmented plane is removed. The effective plane point set retains the regular shape of the roof point cloud information.

[0051] For example, in the plurality of plane point sets, the fragmented plane point set with a point number lower than the point number threshold value, an area smaller than the area threshold value or a rectangularity smaller than the rectangularity threshold value can be removed, thereby selecting at least one effective plane point set.

[0052] After the effective plane point set selection in the above step S152, a set of roof point clouds with regular shapes is obtained, and these point cloud sets have obvious geometric structural features, including linear features and corner point features. In order to support the high-precision registration between the subsequent airborne laser point cloud and high-resolution optical remote sensing image, stable linear features and corner point features need to be extracted therefrom.

[0053] Further, in the above step S1, the first structural features include a plurality of first boundary line segments and a plurality of first corner points; and the first structural features extracted from the denoised airborne laser point cloud can include steps S161-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] In step S2, a high-resolution optical remote sensing image is acquired, and the optical remote sensing image is preprocessed, and second structural features are extracted from the preprocessed optical remote sensing image.

[0062] It should be noted that the "high resolution" of the embodiment of the present application generally refers to the resolution of the optical remote sensing image being 0-0.5 meters, for example, 0.1 meters or 0.5 meters.

[0063] The high-resolution optical remote sensing image has good geometric structure expression capability and can clearly display the boundary contour and corner point features of regular artificial objects (such as houses and roads). In order to realize accurate registration between the optical remote sensing image and the airborne laser point cloud, linear boundary features and corner point features corresponding to the point cloud structure are extracted from the high-resolution optical remote sensing image in this step.

[0064] Further, in the above step S2, the preprocessing of the optical remote sensing image can include steps S211-S212.

[0065] In step S211, the optical remote sensing image is denoised by using a median filtering algorithm and a bilateral filtering algorithm in sequence.

[0066] The optical remote sensing image is denoised by using the median filtering algorithm, which can smooth the small noise of the image, suppress local texture interference, and maintain clear edge contour. After the bilateral filtering, the texture region is further smoothed, and the edge structure is maintained.

[0067] In step S212, the denoised optical remote sensing image is histogram equalized to enhance the brightness and contrast of the optical remote sensing image.

[0068] The histogram equalization processing can improve the brightness and contrast of the optical remote sensing image, enhance the gray difference between the boundary of regular objects (such as buildings) and the background, and ensure that the boundary of regular objects is strengthened in edge detection.

[0069] Further, in the above step S2, the second structural features include a plurality of second boundary line segments and a plurality of second corner points; and the extraction of the second structural features from the preprocessed optical remote sensing image can include steps S22-S23.

[0070] In step S22, a two-dimensional Hough transform algorithm is used to extract the second boundary line segment from the preprocessed optical remote sensing image.

[0071] The two-dimensional Hough transform algorithm is used to extract the second boundary line segment from the preprocessed optical remote sensing image.

[0072] Optionally, for the extracted second boundary line segments, high-low double thresholds of the slope and the intercept can also be set adaptively based on the image statistics, and the features of the partial boundary line segments that meet both high-low double thresholds are screened out from the second boundary line segments, so as to ensure the edge continuity and accuracy of the extracted second boundary line segments.

[0073] The start and end pixel coordinates, the direction angle and the line segment length of each second boundary line segment are recorded. Optionally, the too short line segments and the non-structural boundary line segments with a length lower than a set threshold can also be discarded.

[0074] In step S23, the included angle between two second boundary line segments at the intersection point of the two second boundary line segments is detected, and if the included angle is a preset regular angle, the intersection point is taken as a second corner point.

[0075] The included angle being a preset regular angle can be that the included angle is close to a right angle (80°-100°) or a regular angle (such as 45° or 135°), and in this case, the intersection point can be extracted as a second corner point.

[0076] It should be noted that the above step S1 and the above step S2 do not have a strict sequence, and the two steps can be performed synchronously or sequentially, and the sequence is not limited.

[0077] In step S3, the first structural feature and the second structural feature are unified with a spatial reference, and then registered, and the registered second structural feature is fused with the first structural feature to construct a fused feature set.

[0078] On the basis that the airborne laser point cloud and the high-resolution optical remote sensing image have the same geographic coordinate system, in order to realize high-precision registration of the two, the first structural feature extracted in the above step S1 and the second structural feature extracted in the above step S2 need to be spatially matched. In the embodiment of the present application, the boundary line segment feature and the corner point feature are used as the matching primitives.

[0079] Further, the above step S3 can include steps S31-S35.

[0080] In step S31, the first structural feature is projected to the image coordinate system in which the second structural feature is located.

[0081] For the first structural feature extracted from the airborne laser point cloud, including a plurality of first boundary line segments (such as roof edges) and a plurality of first corner points (such as corner points), using a preliminary geographic coordinate projection relationship, the first structural feature is projected to the image coordinate system in which the optical remote sensing image is located, and the image coordinate system is a plane coordinate system.

[0082] For example, an initial coarse alignment model (such as unit conversion or ground reference plane simplification) known in the art can be used to project the three-dimensional points of the first structural feature to the two-dimensional image coordinate system. ) mapping to the plane pixel coordinates in the image coordinate system ).

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

[0084] In step S32, the projected first boundary line segments and the second boundary line segments are matched according to the angle difference value and the line segment length ratio to form a line segment matching set.

[0085] For each projected first boundary line segment, a set of second boundary line segments close in position is first queried from the plurality of second boundary line segments of the optical remote sensing image, the direction included angle between the projected first boundary line segment and any second boundary line segment in the set is calculated, and second boundary line segments similar in direction are screened out from the set according to a set angle threshold. For example, the difference in the direction included angle between the two can be set to be less than an angle threshold value , and the threshold value may be 10°.

[0086] Then, the line segment length ratio is compared, and candidate line segments with too large length difference are screened out. The part of the second boundary line segments that satisfy both the angle difference value and the line segment length ratio and the projected first boundary line segment are recorded as a line segment matching pair to form a line segment matching set , which can be expressed as:

[0087]

[0088] wherein, represents the i-th projected first boundary line segment; represents the j-th second boundary line segment.

[0089] In step S33, the projected first corner points and the second corner points are matched according to the spatial distance and the direction between the adjacent two boundary line segments to form a corner point matching set.

[0090] For each projected first corner point, a set of second corner points closest in spatial distance is first queried from the plurality of second corner points of the optical remote sensing image, the direction between the adjacent two first boundary line segments in which the projected first corner point is located and the direction between the adjacent two second boundary line segments in which the second corner point is located are calculated, and it is determined whether the two directions form the same angle structure such as “L type” and “cross type”. If yes, the projected first corner point and the second corner point are recorded as a corner point matching pair to form a corner point matching set , which can be expressed as:

[0091] ​​

[0092] wherein, represents the first projected first corner point; represents the first projected second corner point. represents the first projected second corner point. represents the first projected second corner point. is a three-dimensional corner point of the airborne laser point cloud, is a two-dimensional corner point of the optical remote sensing image.

[0093] Step S34, each second corner point in the corner point matching set is assigned to the corresponding first corner point, and the assigned first corner point is formed.

[0094] For each pair of matching pairs in the corner point matching set The mapping operation is performed. The coordinates of the two-dimensional corner point of the optical remote sensing image are assigned as additional attribute values to the corresponding three-dimensional corner point of the airborne laser point cloud , the feature of the three-dimensional corner point of the airborne laser point cloud is extended, and the following five-dimensional vector is formed:

[0095]

[0096] wherein, represents the five-dimensional vector of the th three-dimensional corner point of the airborne laser point cloud.

[0097] Thus, 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] Step S35, the plurality of assigned first corner points are summarized to form a fusion feature set.

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

[0100]

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

[0102] Up to now, after the airborne laser point cloud is registered with the high-resolution optical remote sensing image, the fusion feature set is obtained.

[0103] Step S4, based on the fusion feature set, the airborne laser point cloud is subjected to planar translation and planar rotation to correct two planar coordinates in the airborne laser point cloud except the elevation coordinate, and a corrected airborne laser point cloud is generated.

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

[0105] Further, the above step S4 can include steps S41-S44.

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

[0107] The three-dimensional coordinates of the airborne laser point cloud in the fusion feature set are represented by , and the two-dimensional plane coordinates of the optical remote sensing image corresponding to the airborne laser point cloud are represented by .

[0108] It should be noted that is the coordinates of the two-dimensional corner point of the optical remote sensing image. is obtained by using its corresponding rational polynomial coefficient (RPC) conversion. The correction target of the correction model is to correct the airborne laser point cloud in the fusion feature set to the two-dimensional coordinates of the image coordinate system of the optical remote sensing image, so that and the two-dimensional plane coordinates of the optical remote sensing image corresponding to the airborne laser point cloud in the fusion feature set are as close as possible.

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

[0110]

[0111] wherein, is a small-angle rotation matrix, is a two-dimensional plane translation vector, and the translation amount in the elevation direction is set to 0.

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

[0113]

[0114] Therefore, the three-dimensional coordinates of the corrected airborne laser point cloud are:

[0115]

[0116] wherein, is the micro-rotation angle to be solved, is unchanged.

[0117] Step S42, a first objective function of the correction model is constructed by using the plane re-projection error.

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

[0119]

[0120] wherein, N is the number of matched pairs in the fusion feature set, and is also the number of matched pairs in the corner point matching set in the foregoing step S35.

[0121] After substituting the two-dimensional coordinates of the above-mentioned corrected airborne laser point cloud into , the following is obtained:

[0122]

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

[0124] In order to minimize the plane re-projection error , the partial derivatives of are respectively taken and set to zero, and the following normal equation set is obtained:

[0125]

[0126] Taking as an example, the partial derivative of is taken:

[0127]

[0128] Setting it to zero, the following is obtained:

[0129]

[0130] After arrangement, the following is obtained:

[0131]

[0132] By analogy, the partial derivative of Take the partial derivative and make it zero, and arrange to get another two formulas. The whole can be represented as a small linear equation group, which is solved by using the least square method to obtain the optimal solution of .

[0133] Step S44, the plane translation vector and the small angle rotation matrix are used to transform the airborne laser point cloud, and the corrected airborne laser point cloud is obtained.

[0134] The optimal solution of is applied to all the airborne laser point clouds in the above step S1, the plane precision of the whole point cloud is improved, and the height information is not affected:

[0135]

[0136] Wherein, is the three-dimensional coordinate of the corrected airborne laser point cloud.

[0137] So far, the corrected airborne laser point cloud generated can be simply represented as .

[0138] Step S5, the ranging data of the spaceborne photon counting laser radar is obtained, the atmospheric delay correction and the tidal correction are sequentially performed on the ranging data, and the ranging data after the tidal correction is filtered to form a spaceborne photon counting point cloud; the spaceborne photon counting point cloud and the corrected airborne laser point cloud are unified in space reference and then registered 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 the above step S4, so as to establish the spatial correspondence between the two types of laser sensors, thereby providing a basic support for the subsequent correction of the pointing error and the system ranging error of the spaceborne photon counting laser radar.

[0140] It should be noted that the "spaceborne photon counting laser radar" of the embodiment of the application can be a single beam or a multi-beam, and the number of beams can reach six or more.

[0141] Figure 3 A three-dimensional cloud chart of the observation data of the spaceborne photon counting laser radar according to the embodiment of the application is schematically shown.

[0142] The embodiment of the application is based on the characteristics of the polar orbit satellite and the original airborne laser radar point cloud data, and constructs the observation data of the multi-beam spaceborne photon counting laser radar, as shown in Figure 3 , which includes six beams in total. The known laser pointing error ( ) and the system ranging error ( m) and taking into account atmospheric delay and tidal effects.

[0143] Figure 4 The elevation comparison of the spaceborne photon counting point cloud and the airborne laser point cloud before calibration is schematically shown.

[0144] The elevation comparison of the spaceborne photon counting point cloud and the airborne laser point cloud before calibration is schematically shown. Figure 4 As can be seen from Figure 4 , the elevations of the two data sets have some differences due to laser pointing errors and system ranging errors.

[0145] Further, in the above step S5, the ranging data of the spaceborne photon counting lidar is corrected for atmospheric delay and tides to eliminate the ranging deviation caused by atmospheric delay and earth tide effects in the laser ranging process.

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

[0147]

[0148] The original height of each photon point is corrected for atmospheric delay in laser ranging:

[0149]

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

[0151] For example, GOT4.8 and IERS global tide models can be used for tidal correction, which includes solid tide , ocean tide , polar tide correction, and the calculation formula is as follows:

[0152]

[0153] wherein, is the tidal correction amount.

[0154] For each photon point, the above tidal correction amount is calculated by interpolation according to its latitude (lat), longitude (lon) and timestamp . The ranging data after the above atmospheric delay correction is corrected for tides:

[0155]

[0156] wherein, is the ranging data after tide correction.

[0157] Further, in the step S5, the ranging data after tide correction is filtered to form the spaceborne photon counting point cloud.

[0158] For example, the ranging data after tide correction is firstly calculated to form the initial spaceborne photon counting point cloud data using the laser radar equation. Then, the initial spaceborne photon counting point cloud data is filtered based on the linear feature photon data filtering algorithm (LFPSE) to remove the noise and low-confidence photon points and retain the signal photons, thereby forming the spaceborne photon counting point cloud distributed along the orbit. .

[0159] Further, after the step S5, the corrected airborne laser point cloud and the spaceborne photon counting point cloud are unified in the spatial reference to perform registration, and the registered spaceborne photon counting point cloud can include the steps S53-S55.

[0160] In the step S53, the spatial coordinate references of the corrected airborne laser point cloud and the spaceborne photon counting point cloud are unified to the geocentric coordinate system.

[0161] For example, the spatial coordinate reference of the spaceborne photon counting point cloud is unified to the geocentric coordinate system, such as the WGS84 coordinate system, of the corrected airborne laser point cloud .

[0162] Further, after the step S53, the corrected airborne laser point cloud and the spaceborne photon counting point cloud unified in the spatial reference are further processed. For example, in the geocentric coordinate system, the spaceborne photon counting point cloud is divided into multiple equal-interval segments according to the orbit direction, and the corrected airborne laser point cloud is cropped according to the orbit range of the spaceborne photon counting point cloud in each equal-interval segment.

[0163] Specifically, the spaceborne photon counting point cloud can be divided into multiple equal-interval segments according to every 500 meters, which is used for subsequent segmented registration to adapt to large-scale terrain change areas and improve the local fitting quality. In each equal-interval segment, the airborne laser point cloud is cropped according to the orbit range of the spaceborne photon counting point cloud, and only the overlapping area with the spaceborne photon counting point cloud is retained to improve the calculation efficiency and matching accuracy.

[0164] In the step S54, the barycentric coordinates of the corrected airborne laser point cloud and the spaceborne photon counting point cloud are respectively calculated, and the corrected airborne laser point cloud and the spaceborne photon counting point cloud are coarsely registered in the elevation profile according to the difference between the two barycentric coordinates.

[0165] Computing a spaceborne photon-counting point cloud and a corrected airborne laser point cloud of the barycenter coordinates , obtaining a difference between the two barycenter coordinates . Taking the difference as an initial translation vector, respectively performing initial translation on and , extracting cross-section elevation profile curves, and performing coarse registration on and based on the elevation profile.

[0166] Step S55, for each spaceborne point in the spaceborne photon-counting point cloud, searching for the nearest neighbor point corresponding to the spaceborne point in the corrected airborne laser point cloud according to a distance threshold, and establishing a second objective function about the translation vector and the rotation matrix according to each spaceborne point and the corresponding nearest neighbor point; solving the translation vector and the rotation matrix by minimizing the second objective function; and transforming the spaceborne photon-counting point cloud using the translation vector and the rotation matrix to obtain a registered spaceborne photon-counting point cloud.

[0167] Specifically, for each spaceborne point in the coarse-registered , searching for the nearest neighbor point corresponding to the spaceborne point in the coarse-registered airborne laser point cloud according to a distance threshold, and filtering by the distance threshold to avoid false matching.

[0168] Taking the spaceborne point and the matched nearest neighbor point as the basis, a second objective function about the rotation matrix and the translation vector is established as follows:

[0169]

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

[0171] Then, the optimal solution is used to apply the following transformation model to each spaceborne point of the spaceborne photon-counting point cloud :

[0172]

[0173] wherein represents the registered position of the spaceborne point .

[0174] So far, for each satellite-borne point , it can be registered with the nearest neighbor point in the airborne laser point cloud to obtain the registered satellite-borne point . The registered satellite-borne photon counting point cloud can be simply represented as .

[0175] To further improve the registration accuracy, after the above step S55, it can further include: calculating the Euclidean residual of all matching points in the point set registered by the airborne laser point cloud and the satellite-borne photon counting point cloud:

[0176]

[0177] wherein, represents the ith pair of matching points in the point set after registration.

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

[0179] For example, the accuracy indicators such as root mean square error (RMSE), maximum deviation, and standard deviation are counted, and the satellite-borne points in the satellite-borne photon counting point cloud whose error exceeds 3 times the standard deviation are removed to avoid the influence of extreme points on subsequent error modeling, and a weighted nearest neighbor iterative optimization is used to suppress the influence of low-quality areas.

[0180] Step S6, constructing a radar equation of the satellite-borne photon counting lidar, determining the pointing error and system ranging error of the satellite-borne photon counting lidar according to the registered satellite-borne photon counting point cloud and the radar equation.

[0181] Further, the above step S6 can include steps S61-S63.

[0182] Step S61, constructing a radar equation about multiple satellite-borne points of the satellite-borne photon counting lidar.

[0183] For multiple photon points (also referred to as satellite-borne points) of the satellite-borne photon counting lidar, the one-way flight distance of the ith photon is denoted as . The unit vector of the emission direction can be further defined as: and the pitch angle .

[0184] ​

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

[0186]

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

[0188] Step S62, define the pointing error and the system ranging error of the spaceborne photon counting laser radar; correct the radar equation according to the pointing error and the system ranging error.

[0189] In the embodiment of the present application, the pointing error includes an azimuth angle error and a pitch angle error , and the system ranging error is a systematic deviation of ranging . The modified unit direction vector after introducing the errors is:

[0190]

[0191] wherein, the small perturbation is expanded as:

[0192]

[0193] wherein, the partial derivatives are respectively:

[0194]

[0195] The measured distance after correction is:

[0196]

[0197] Therefore, the measured position after introducing the error correction to the above radar equation is:

[0198]

[0199] The formula of the is the corrected radar equation.

[0200] Step S63, solve the pointing error and the system ranging error according to the corrected radar equation and the registered spaceborne photon counting point cloud.

[0201] Define the registered spaceborne photon counting point cloud and the measured position after error correction Residual error between is:

[0202]

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

[0204]

[0205] Construct a residual error matrix for all the spaceborne points (N photon points), and use the Jacobian matrix to represent:

[0206]

[0207] The further residual error can be represented as:

[0208]

[0209] wherein, is a residual error vector of 3N x 1, is a Jacobian matrix of 3N x 3, is an error parameter vector to be estimated.

[0210] Solve by using the least square method:

[0211]

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

[0213] Then, the measurement position of the spaceborne photon counting laser can be corrected, the direction angle (azimuth angle, elevation angle) is updated, the ranging parameter is updated, the three-dimensional coordinates of the spaceborne photon counting point cloud are recalculated, and the updated measurement position is obtained.

[0214] By using the spaceborne photon counting laser radar field-free on-orbit calibration method proposed in the embodiments of the present application, error modeling is performed, and the direction angle (azimuth angle, elevation angle) and the ranging parameter are calculated. The results show that the present method realizes extremely high precision estimation for the three error parameters. Among them, the estimated value of the azimuth angle is-4.003 angular seconds, the estimated value of the elevation angle is-2.002 angular seconds, and the estimated value of the system ranging error is-0.400 meters, which are almost completely consistent with the true set values of 4 angular seconds, 2 angular seconds and 0.4 meters respectively, and the errors are only within the order of magnitude of 0.003 angular seconds and 1 millimeter, which verifies the high robustness and error sensitivity of the present method.

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

[0216] As shown in Figure 5 The three-dimensional coordinates of the spaceborne photon counting point cloud recalculated by updating the direction angle (azimuth angle, pitch angle) and the ranging parameter have high consistency with the point cloud data of the airborne laser radar.

[0217] So far, the applicability of the method in the multi-beam spaceborne photon counting laser radar is comprehensively verified. The method does not need to rely on the ground calibration field or additional equipment, but only relies on the data of the spaceborne platform itself and auxiliary image resources, so that the ranging error and the direction error of the multi-beam photon laser system can be estimated and corrected with high precision, and the method has strong engineering feasibility and popularization value, and is especially suitable for remote sensing task scenes which are difficult to arrange ground calibrators.

[0218] In summary, the embodiment of the application provides a field-free on-orbit calibration method for a spaceborne photon counting laser radar. The airborne laser point cloud with high precision and the high-resolution optical remote sensing image are jointly registered, the advantages of the planar positioning of the high-resolution optical remote sensing image and the height accuracy of the airborne laser point cloud are combined, and the planar positioning accuracy of the airborne laser point cloud is improved. Then, the pointing error and the system ranging error of the spaceborne photon counting laser radar are corrected by matching the spaceborne photon counting point cloud with the corrected airborne laser point cloud, so that the field-free on-orbit geometric calibration of the spaceborne photon counting laser radar is realized. The method avoids the dependence on the ground calibration field laser detector, has the advantages of strong adaptability, high calibration accuracy and wide application range, and provides a new technical means for improving the field-free on-orbit calibration of the spaceborne photon counting laser radar height measurement system.

[0219] The flowcharts and block diagrams in the drawings illustrate the possible implementation architecture, function and operation of the system, method and computer program product according to various embodiments of the application. In this regard, each block in the flowchart or block diagram can represent a module, a program segment, or a part of code containing one or more executable instructions for implementing a specified logical function. It should also be noted that in some alternative implementations, the functions noted in the blocks can occur in different order than that shown in the figure. For example, two blocks that are shown in succession can actually be executed substantially in parallel, and sometimes they can be executed in reverse order, depending on the function involved. It should also be noted that each block in the block diagram or flowchart, and the combination of blocks in the block diagram or flowchart, can be implemented by a dedicated hardware-based system that performs the specified function or operation, or can be implemented by a combination of special-purpose hardware and computer instructions.

[0220] Furthermore, the terms "first", "second", etc. are used herein for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Therefore, the features defined with "first", "second" can explicitly or implicitly include one or more of the features. In the description of the present application, the meaning of "a plurality of" is at least two, for example two, three, etc., unless otherwise explicitly specified. In addition, the word "one" or "an" before an element does not exclude the presence of a plurality of such elements.

[0221] It will be appreciated by those skilled in the art that features recited in the various embodiments of the present application can be combined and / or integrated in a variety of ways, even if such combinations or integrations have not been expressly specified in the present application. In particular, features recited in the various embodiments of the present application can be combined and / or integrated in a variety of ways without departing from the spirit and teachings of the present application. All such combinations and / or integrations are within the scope of the present application.

[0222] The above describes embodiments of the present application. However, these embodiments are merely for illustrative purposes and are not intended to limit the scope of the present application. Although each embodiment is described above separately, this does not mean that the measures in each embodiment cannot be used advantageously in combination. Those skilled in the art can make various substitutions and modifications without departing from the scope of the present application, and these substitutions and modifications should all fall within the scope of the present application.

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 fusion 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