Spaceborne lidar point cloud and image feature joint extraction method, multi-level reference image construction method and electronic device

Abstract

The application provides a spaceborne laser radar point cloud and image feature joint extraction method, a multi-level reference image construction method and an electronic device, and relates to the technical field of remote sensing surveying and geographic information. The extraction method comprises the following steps: acquiring optical remote sensing images and radar data; a plurality of sliding windows are constructed in the along-track direction, and denoising processing is performed according to the plurality of sliding windows to obtain a plurality of candidate signal photons; the local terrain features of the candidate signal photons in each sliding window are calculated, and the candidate spot corresponding photons located in the smooth area are selected according to the local terrain features; it is judged whether the pixels of the plurality of candidate spot corresponding photons in the spot range satisfy a preset condition; the representative elevation of the plurality of signal photon elevation samples in the target spot range that satisfy the preset condition is selected, the pixel coordinates of the hard surface pixels in the target spot range in the optical remote sensing image are combined with the representative elevation, and the joint feature points are obtained.

Description

Technical Field

[0001] This application relates to the fields of remote sensing mapping and geographic information technology, and in particular to a method for joint extraction of point clouds and image features from spaceborne lidar, a method for constructing multi-level reference images, and an electronic device. Background Technology

[0002] Spaceborne photon-counting lidar utilizes a single-photon detection system to acquire high-temporal-frequency, high-spatial-density altimetry data. It features high precision, lightweight design, and wide coverage, and has been widely applied in fields such as ice sheet monitoring, forest structure retrieval, and sea surface and inland water elevation measurement. Meanwhile, high-resolution optical imagery offers significant advantages in ground feature identification, geometric texture representation, and spatial positioning. Synergistically utilizing these two technologies holds promise for constructing a joint control information system that combines planar texture representation capabilities with vertical elevation constraints.

[0003] In related technologies, multi-source data fusion methods typically spatially correlate laser altimetry data with image data. Based on the mapping relationship between the ground coordinates of the laser observation point and the image pixel coordinates, joint feature points containing elevation attributes and image location information are constructed to provide control information support for downstream tasks such as image geometric correction, multi-source data registration, and 3D positioning.

[0004] However, photon counting systems are extremely sensitive to weak echoes. The raw observation data contains atmospheric noise, surface scattering noise, and non-ground echo signals such as vegetation and building apex. Furthermore, the terrain undulations and surface cover types in the laser-covered areas vary. These factors are directly transmitted to the joint feature points during the spatial correlation of laser data and image data, making the constructed joint feature points susceptible to interference in terms of elevation accuracy and geometric stability. This limits the reliability of its application in high-precision remote sensing geometric correction and three-dimensional positioning. Summary of the Invention

[0005] In view of the above problems, this application provides a method for joint extraction of point cloud and image features from spaceborne lidar, a method for constructing multi-level reference images, and an electronic device.

[0006] The first aspect of this application provides a method for jointly extracting point clouds and image features from a spaceborne lidar, comprising: acquiring optical remote sensing images covering a target area and spaceborne photon-counting lidar data, wherein the spaceborne photon-counting lidar data includes multiple signal photons; constructing multiple sliding windows along the orbital direction, performing denoising processing on the multiple signal photons according to the multiple sliding windows to obtain multiple candidate signal photons within each sliding window; calculating the local terrain features of the surface area corresponding to the candidate signal photons within each sliding window, and selecting candidate light spot corresponding photons located in stable areas from the multiple candidate signal photons within each sliding window according to the local terrain features; determining whether the pixels of the multiple candidate light spot corresponding photons in the light spot range covered by the optical remote sensing image meet preset conditions, wherein the preset conditions are used to determine the surface type of the light spot range based on the spectral and texture features of each pixel; selecting the representative elevation of the elevation samples of multiple signal photons in the target light spot range that meet the preset conditions, and combining the pixel coordinates of the hard surface pixels in the optical remote sensing image of the target light spot range with the representative elevation to obtain joint feature points.

[0007] Furthermore, multiple sliding windows are constructed along the track direction, and multiple signal photons are denoised according to the multiple sliding windows to obtain multiple candidate signal photons in each sliding window. This includes: sorting each signal photon in chronological order, accumulating the distances between adjacent signal photons after sorting to obtain the cumulative distance along the track for each signal photon; dividing the cumulative distance along the track direction into multiple sliding windows, and counting the number of neighborhood points in the first preset neighborhood of each signal photon in each sliding window; determining a noise threshold based on the number of neighborhood points, and confirming signal photons in each sliding window whose number of neighborhood points is greater than or equal to the noise threshold as candidate signal photons.

[0008] Furthermore, local terrain features include slope and flatness; the local terrain features of the surface area corresponding to the candidate signal photons in each sliding window are calculated, and the candidate light spot corresponding photons located in the stable area are selected from the candidate signal photons in each sliding window based on the local terrain features, including: establishing a fitting curve of the elevation of the candidate signal photons in each sliding window as a function of distance along the track direction, and calculating the slope and flatness of the candidate signal photons along the track direction based on the fitting curve; the candidate signal photons whose absolute value of slope is not greater than a preset slope threshold and whose flatness is not greater than a preset flatness threshold are identified as the candidate light spot corresponding photons located in the stable area.

[0009] Furthermore, after selecting the candidate light spot corresponding photons located in the stable region from multiple candidate signal photons in each sliding window based on local terrain features, the method further includes: merging adjacent sliding windows that continuously meet the selection criteria along the track direction corresponding to each candidate light spot photon to form a track-smooth section window. The selection criteria are that the absolute value of the slope of the candidate signal photons in the sliding window is less than or equal to a preset slope threshold and the smoothness is less than or equal to a preset smoothness threshold; selecting the track-smooth section window according to preset continuity constraints, retaining the candidate light spot corresponding photons in the track-smooth section window that meet the preset continuity constraints, and updating the candidate light spot corresponding photons selected based on local terrain features. The preset continuity constraints include at least one of a continuous window number threshold, a continuous track length threshold, and an effective candidate photon density threshold.

[0010] Further, determining whether the pixels corresponding to multiple candidate light spots within the light spot range covered by the optical remote sensing image meet preset conditions includes: mapping the coordinates of the photons corresponding to multiple candidate light spots to the pixel coordinates of the optical remote sensing image to obtain the projection positions of the photons corresponding to multiple candidate light spots in the optical remote sensing image; acquiring the device parameters of the lidar, confirming the light spot radius of the lidar based on the device parameters, and determining the light spot range covered by the multiple candidate light spots within the optical remote sensing image based on each projection position and the light spot radius; acquiring multiple pixels within the light spot range, extracting the spectral and texture features of each pixel, and determining whether each pixel belongs to a hard surface pixel based on the spectral and texture features. The spectral features include at least grayscale, brightness, and image channel differences, and the texture features include at least contrast and the homogeneity of the pixel in the second preset neighborhood; calculating the ratio of the number of pixels belonging to hard surfaces to the total number of pixels within the light spot range. When the ratio is greater than a preset ratio threshold, it is determined that the pixels within the light spot range meet the preset conditions.

[0011] Furthermore, the determination of whether each pixel belongs to a hard surface pixel based on spectral and texture features includes: when each pixel simultaneously satisfies the following conditions: brightness is within a preset brightness range, image channel difference is within a preset difference range, contrast is less than or equal to a preset contrast threshold, and homogeneity is greater than or equal to a preset homogeneity threshold, the pixel is determined to belong to a hard surface.

[0012] Furthermore, representative elevations of multiple signal photon elevation samples within the target spot range that meet preset conditions are selected, including: using the median of each signal photon elevation sample as a reference elevation; calculating the deviation between each signal photon elevation sample and the reference elevation, and calculating the standard deviation of multiple deviations; correcting the standard deviations with a preset correction coefficient, removing signal photon elevation samples whose absolute deviation is greater than the corrected standard deviation, and calculating the mean of the remaining signal photon elevation samples to obtain the representative elevation.

[0013] The second aspect of this application provides a method for constructing a multi-level reference image, comprising: obtaining joint feature points, wherein the joint feature points are extracted according to the above-mentioned method for joint extraction of spaceborne lidar point cloud and image features; and constructing a multi-level reference image based on the joint feature points.

[0014] Furthermore, a multi-level reference image is constructed based on the joint feature points, including at least one of the following: taking the position of the joint feature points in the optical remote sensing image as the center, expanding the corresponding preset number of pixels along the row and column directions in the optical remote sensing image to obtain a point-level joint feature reference image; dividing the target area into multiple sub-regions, extracting multiple joint feature points in each sub-region to obtain a block-level joint feature reference image; and extracting multiple joint feature points in the target area to obtain a scene-level joint feature reference image.

[0015] A third aspect of this application provides an electronic device, comprising: one or more processors; and a memory for storing one or more programs, wherein when the one or more programs are executed by the one or more processors, the one or more processors cause the one or more processors to implement the above-described method.

[0016] Compared with existing technologies, the method for joint extraction of point clouds and image features from spaceborne lidar, the method for constructing multi-level reference images, and the electronic equipment provided in this application have at least the following advantages:

[0017] (1) By finely extracting signal photons and introducing dual constraints of track slope and flatness, it can be effectively ensured that the joint feature points mainly come from the local stable hard surface area, and the unstable areas such as undulating terrain, building facades and shading edges are automatically eliminated, which significantly improves the stability and reliability of elevation control.

[0018] (2) By explicitly considering the actual size of the spaceborne laser spot on the ground, the image spectrum and texture features are statistically analyzed within the spot scale, so that the laser altimetry information and the image texture / spectral response can establish a consistent correspondence within the same physical space, which can effectively avoid the spatial scale mismatch problem that is easy to occur in traditional "single point - single pixel" matching.

[0019] (3) By combining brightness, spectral differences and texture features to construct a hard surface discrimination strategy, and calculating the proportion of hard pixels within the light spot, only retaining targets suitable as long-term control and matching benchmarks such as roads, squares, and building roofs, the stability and reusability of joint features can be improved from the semantic level.

[0020] (4) By adopting a robust statistical strategy of “median-standard deviation elimination-mean” for the elevation samples of multiple signal photons in the light spot, the influence of isolated anomalous photons on the elevation control results can be effectively suppressed, making the elevation attributes of the joint feature points insensitive to noise and robust to local anomalies.

[0021] (5) By cropping fixed-size, multi-level reference images with joint feature points as the core, a correspondence of “pixel-latitude-longitude-elevation-attribute” can be established from point level, block level to scene level, forming a joint feature product system with unified structure, standard format, easy to port and batch process. Attached Figure Description

[0022] The above and other objects, features and advantages of this application will become clearer from the following description of embodiments with reference to the accompanying drawings, in which:

[0023] Figure 1 A flowchart illustrating the method for joint extraction of point cloud and image features from a spaceborne lidar according to an embodiment of this application is shown.

[0024] Figure 2 This illustration schematically shows a diagram of screening candidate light spots corresponding to photons based on slope and flatness according to an embodiment of this application;

[0025] Figure 3 The schematic diagram illustrates the principle of associating laser spot with image features according to an embodiment of this application;

[0026] Figure 4 A flowchart illustrating the determination of hard surfaces according to an embodiment of this application is shown schematically.

[0027] Figure 5 A flowchart illustrating a method for constructing multi-level reference images according to an embodiment of this application is shown schematically;

[0028] Figure 6 The diagram illustrates a block diagram of an electronic device suitable for implementing a method for jointly extracting point clouds and image features from spaceborne lidar or for constructing multi-level reference images, according to embodiments of this application. Detailed Implementation

[0029] The embodiments of this application will now 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 this application. In the following detailed description, numerous specific details are set forth to provide a thorough understanding of the embodiments of this application for ease of explanation. However, it will be apparent that one or more embodiments may be implemented without these specific details. Furthermore, descriptions of well-known structures and technologies are omitted in the following description to avoid unnecessarily obscuring the concepts of this application.

[0030] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the scope of this application. The terms “comprising,” “including,” etc., as 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.

[0031] 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 are to be interpreted in a manner consistent with the context of this specification, and not in an idealized or overly rigid way.

[0032] For the joint feature extraction problem of point clouds from spaceborne photon-counting lidar and high-resolution optical images, based on the principles and application effects of the multi-source data fusion methods mentioned in the background section, the shortcomings of existing technologies are summarized as follows:

[0033] (1) Lack of quantitative constraints on flat areas based on track slope and fitting residuals. When selecting control areas, existing technologies rely on external DEM (Digital Elevation Model) or experience-based judgment. They lack quantitative judgment on the flatness and smoothness of terrain based on the photon altimetry data itself. They do not use indicators such as track slope and fitting residuals to screen truly flat and stable hard surface areas, which easily introduces undulating terrain or local irregular areas into the control area, affecting the geometric stability of the control points.

[0034] (2) The lack of a "hard pixel ratio" constraint at the spot scale makes it impossible to finely screen stable features. Existing technologies mostly select control areas macroscopically within the image range, without using laser spots as units to statistically analyze multispectral grayscale, index, and texture features within the spot coverage area, and without defining constraints such as the "hard pixel ratio within the spot". This makes it difficult to remove unstable features such as vegetation, water bodies, and farmland in a timely manner, resulting in insufficient temporal stability and geometric reliability of the control targets.

[0035] (3) Robust statistical methods were not used to process photon elevations and assign them to image feature points. Existing technologies mainly rely on external DEMs to obtain elevation information. They do not perform robust statistical processing such as median-standard deviation elimination-mean at the spot level for multiple photons within the same pulse, nor do they uniformly assign the robustly estimated photon elevations to hard surface image feature points within the spot. Therefore, they are insufficient in terms of joint feature point elevation control accuracy and anti-anomaly capability.

[0036] (4) Lack of a multi-scale reference image system at the point, block, and scene levels. Existing technologies are relatively simple in their control film organization, usually generating only a single-scale control film or a small number of control points. They have not built a multi-scale product system around joint feature points, including "point-level reference images centered on feature points, block-level reference images with fixed spatial ranges, and whole-scene reference images," making it difficult to meet both the needs of fine geometric correction and the needs of large-scale engineering applications.

[0037] In view of this, embodiments of this application provide a method for joint extraction of point clouds and image features from spaceborne lidar, a method for constructing multi-level reference images, and an electronic device. The system can be designed from signal photon recognition, track flatness constraints, spot-scale image feature statistics, hard surface screening to joint feature product construction, so as to obtain a high-confidence and highly reusable three-dimensional joint feature set.

[0038] Figure 1 The flowchart illustrates a method for jointly extracting point clouds and image features from a spaceborne lidar according to an embodiment of this application.

[0039] like Figure 1 As shown, the method for joint extraction of point cloud and image features from spaceborne lidar in this embodiment includes operations S101 to S105.

[0040] In operation S101, optical remote sensing images and spaceborne photon counting lidar data covering the target area are acquired. The spaceborne photon counting lidar data includes multiple signal photons.

[0041] Optical remote sensing images can be panchromatic, color, or multispectral. Preferably, the spatial resolution of the optical remote sensing image is better than 0.5m to ensure sufficient spatial detail representation for subsequent image feature extraction and hard surface discrimination. The optical remote sensing image is accompanied by rational function model (RPC) parameters, which are used to map the ground 3D coordinates of the spaceborne lidar photon points to the image pixel coordinate system.

[0042] Spaceborne photon counting lidar data can be photon observation data acquired by single-photon altimetry systems such as ICESat-2, and it includes at least the following data: beam number, photon longitude. Photon latitude Photon observation time Photon Elevation .in, Indicates the first Photon. Optionally, the photon elevation is the ellipsoidal height under a reference ellipsoid consistent with the image's strictly geometric positioning model; when the input photon elevation is positive, it is converted into ellipsoidal height based on the geoid undulation value or elevation correction provided by geoid models such as EGM (Earth Gravitational Model) 2008, to ensure the consistency of subsequent laser data and image positioning model on the vertical reference.

[0043] To ensure joint analysis of the two types of data within a unified spatial framework, consistency verification and basic preprocessing can be performed. Consistency verification includes at least the following: horizontal coordinate baseline consistency verification, vertical baseline consistency verification, temporal coverage consistency check, spatial overlap area check, and preliminary checks for data integrity and outliers. Basic preprocessing includes: denoising, grayscale normalization, radiometric correction, or histogram stretching of the original imagery; and field filtering, invalid record removal, or beam-based sorting of the original photon data.

[0044] The temporal coverage consistency check refers to confirming whether the imaging time of the high-resolution optical image and the acquisition time of the spaceborne photon-counting lidar data meet the preset time interval requirements. Preferably, the time interval does not exceed 3 years; more preferably, the time interval does not exceed 1 year. For areas with rapid surface changes such as vegetation, water bodies, and farmland, the time interval can be further limited to no more than 6 months; for long-term stable hard surface areas such as roads, squares, and building roofs, the time interval can be appropriately relaxed according to the actual application accuracy requirements. The spatial overlap area check refers to confirming whether there is a valid spatial intersection between the coverage area of ​​the high-resolution optical image and the coverage area of ​​the spaceborne photon-counting lidar data. Only the local intersection area is used to construct joint feature products, and the actual spatial overlap area can be used as the scope for subsequent processing.

[0045] In operation S102, multiple sliding windows are constructed along the track direction, and multiple signal photons are denoised according to the multiple sliding windows to obtain multiple candidate signal photons in each sliding window.

[0046] In some embodiments, the above operation S102 may further include: sorting the signal photons in chronological order, and accumulating the distances between adjacent signal photons after sorting to obtain the cumulative distance along the track for each signal photon. The cumulative distance along the track is divided into multiple sliding windows along the track direction. The number of neighborhood points in the first preset neighborhood of the signal photon in each sliding window is counted. The noise threshold is determined based on the number of neighborhood points. Signal photons with a number of neighborhood points greater than or equal to the noise threshold in each sliding window are identified as candidate signal photons.

[0047] Specifically, for each photon sequence within a beam, according to the observation time Sort the data in ascending order, using the observation time of the first photon as the starting time. Calculate the cumulative distance along the orbit for each photon based on the ground velocity along the orbit from the nadir point or the ground distance between adjacent photons. A sliding window is constructed along the track direction. The photon elevation distribution and neighborhood statistical characteristics within each window are analyzed to identify candidate signal photons and initially eliminate noise photons. The length of the sliding window along the track can be 10m to 200m. If the window length is too small, the number of effective photons within the window may be insufficient, easily leading to unstable statistical characteristics and fitting results; if the window length is too large, it may cross different land features or different terrain change units, reducing the precision of local slope and flatness discrimination. The sliding step size can be set to 10% to 80% of the window length to maintain a certain overlap between adjacent windows. The number of photon neighborhood points or elevation histogram features can also be statistically analyzed within each sliding window. To ensure the reliability of statistical analysis and elevation fitting within the window, if the number of original photons within the window is less than 10, the observation sample within that window can be considered insufficient, and it will not participate in candidate signal photon identification or will be merged with adjacent windows; furthermore, if the number of original photons within the window is not less than 20, elevation histogram, multi-peak fitting, or neighborhood density statistical analysis will be performed.

[0048] Noisy photons typically exhibit a small number of neighborhood points, low local density, and discrete spatial distribution, while true signal photons typically exhibit a large number of neighborhood points, high local density, and are concentrated near major elevation peaks. Therefore, it is advisable to preferentially perform bi-Gaussian or multi-peak fitting on the neighborhood point count or local density statistics of each photon within the window. After fitting, the low-mean component can be considered as the noise photon component, and the high-mean component as the candidate signal photon component. The intersection point between the two components is preferentially selected as the candidate signal photon discrimination threshold. .

[0049] When bi-Gaussian fitting fails to converge, component separation is insufficient, or the photon distribution within the window does not satisfy the bimodal characteristic, candidate signal photons can be determined using methods such as elevation histogram main peak detection, neighborhood statistic quantile thresholding, robust statistical thresholding, or local density clustering. Alternatively, the window can be marked as invalid and skipped from subsequent processing. For the th A photon, if its neighborhood statistic (i.e., the first) The number of photons in a preset neighborhood (the number of photons in a given neighborhood) satisfies: If so, then it is marked as a candidate signal photon.

[0050] In other implementations, local density analysis, cluster analysis, histogram peak detection, machine learning classification models, or other identification methods suitable for distinguishing candidate signal photons from noise photons may also be used.

[0051] In operation S103, the local terrain features of the land area corresponding to the candidate signal photons in each sliding window are calculated, and the candidate light spot corresponding to the stable area is selected from multiple candidate signal photons in each sliding window based on the local terrain features.

[0052] In some embodiments, local terrain features include slope and smoothness. Specifically, a fitted curve is established to show the elevation of candidate signal photons within each sliding window as a function of distance along the track direction, and the slope and smoothness of the candidate signal photons along the track direction are calculated based on the fitted curve.

[0053] For example, B-spline fitting, low-order polynomial fitting, local regression fitting, or other smoothing fitting methods can be used to fit the high-order sequence of candidate signal photons. When using B-spline fitting, the B-spline order can be 3; the B-spline nodes can be set at equal intervals along the track, with a node spacing of 5m to 50m. When using low-order polynomial fitting, the polynomial order can be 1 to 3, where 1st order is suitable for approximately flat or linearly changing regions, 2nd order is suitable for slowly undulating regions, and 3rd order is suitable for regions with slightly complex local changes but still continuous and smooth. It is preferable to avoid using excessively high-order polynomials to reduce the risk of overfitting.

[0054] Suppose that within a certain window, the first... The orbital distance of the candidate signal photons is Elevation is Then the fitted curve can be expressed as:

[0055]

[0056] in: Indicates the first The fitting curve obtained by fitting candidate signal photons. Indicates the first A candidate signal photon elevation fitting model is used. The slope is calculated based on this fitted curve. The average of the absolute values ​​of the first derivatives of the fitted curves within the window can be used as... This is used to characterize the overall slope level along the track of the window. It can be expressed as the ratio of the change in elevation to the change in distance along the track, and is a dimensionless ratio. It can also be converted into a percentage slope.

[0057] The residual sequence is obtained from the difference between the actual elevation and the fitted elevation:

[0058]

[0059] in: Indicates the first The fitting residuals of each candidate signal photon. Based on the residual sequence, the residual standard deviation, root mean square error, or other dispersion statistics can be further calculated, denoted as smoothness. Subsequently, candidate signal photons whose absolute slope value and flatness within each sliding window are not greater than a preset slope threshold and are not greater than a preset flatness threshold are identified as candidate spot photons located in the stable region. That is, assuming the slope threshold along the track is... The flatness threshold is When the following conditions are met simultaneously within a certain sliding window: and At that time, the candidate signal photons within the window are marked as the photons corresponding to the flat candidate light spot. Among them, It can be set according to the upper limit of the surface slope, with the preferred value range being 0.02 to 0.15. For flat and hard surfaces such as roads and squares, it can be 0.02 to 0.05. The upper limit of the residual standard deviation or root mean square error can be set, with a preferred value range of 0.05m to 1.00m. For the preset high-precision control area, a value of 0.05m to 0.20m can be used. Photons with smaller slopes and smaller fitting residuals are more suitable for use as photons corresponding to flat candidate spots in subsequent joint feature construction.

[0060] This screening mechanism can eliminate observation points located on steep slopes, building edges, obstructed edges, or areas with significant undulations from candidate signal photons, retaining candidate photons that are more likely to fall on geometrically stable surfaces.

[0061] Figure 2 This diagram schematically illustrates a method for selecting photons corresponding to candidate light spots based on slope and flatness according to an embodiment of this application. For example... Figure 2 As shown, the photons marked with circles represent the candidate light spots that are retained after screening and located in the flat and stable region. Discrete photons scattered at higher elevation positions are eliminated because they cannot meet the flatness constraints.

[0062] In some embodiments, after the above operation S103, the method may further include: merging adjacent sliding windows that continuously meet the screening conditions along the track direction corresponding to the photons of each candidate spot to form a track-smooth section window, wherein the screening condition is that the absolute value of the slope of the candidate signal photons in the sliding window is less than or equal to a preset slope threshold and the smoothness is less than or equal to a preset smoothness threshold; screening the track-smooth section window according to a preset continuity constraint condition, retaining the photons corresponding to the candidate spots in the track-smooth section window that meet the preset continuity constraint condition, and updating the photons corresponding to the candidate spots screened according to local terrain features, wherein the preset continuity constraint condition includes at least one of a continuous window number threshold, a continuous track length threshold, and an effective candidate photon density threshold.

[0063] In this context, "continuous" refers to adjacent sliding windows being physically adjacent or overlapping along the track direction, with the center distances of the windows distributed sequentially along the track. If there are no non-compliant windows between two adjacent windows that meet the conditions, they can be directly merged. If there is a non-compliant window in the middle, its length, the number of candidate photons, and the stability of adjacent windows can be used to determine whether it should be retained as a short discontinuity. To ensure the reliability of the flat section, in a preset high-precision control scenario, the non-compliant window in the middle can be designated as a break point.

[0064] The track-level smooth section window should have no fewer than three consecutive effective windows; the minimum track length is preferably 30m to 50m. Simultaneously, the track-level smooth section windows can be screened based on the effective candidate photon density. The effective candidate photon density can be defined as the ratio of the number of candidate signal photons within the track-level smooth section window to the track length of the window, preferably not less than 0.2 photons / m. Areas that do not meet the requirements for the number of consecutive windows, minimum track length, or effective candidate photon density can be eliminated to exclude undulating terrain, building facades, obstructed edges, and other unstable areas. Sliding windows that do not simultaneously meet the slope threshold and smoothness threshold, or areas within the sliding window where only some candidate photons or some sub-segments meet the smoothness conditions but the overall stability is insufficient, should be eliminated or designated as the break point for the track-level smooth section window to exclude undulating terrain, building facades, obstructed edges, and other unstable areas.

[0065] In operation S104, it is determined whether the pixels corresponding to the photons of multiple candidate light spots in the light spot range covered by the optical remote sensing image meet the preset conditions. The preset conditions are used to determine the land surface type where the light spot range is located based on the spectral and texture characteristics of each pixel.

[0066] Figure 3 The schematic diagram illustrates the principle of associating laser spot with image features according to an embodiment of this application.

[0067] like Figure 3 As shown, the laser spot has a real coverage area on the ground, which may include various ground features such as roads, buildings, and vegetation. To establish the spatial correspondence between spaceborne photon counting lidar observations and high-resolution optical images, the ground three-dimensional coordinates of the photons corresponding to the flat candidate laser spot need to be mapped to the image pixel coordinate system.

[0068] Figure 4 A flowchart illustrating the determination of a hard surface according to an embodiment of this application is shown schematically. Figure 4 As shown, in some embodiments, the above-described operation S104 may further include operations S401 to S404.

[0069] In operation S401, the coordinates of the photons corresponding to multiple candidate light spots are mapped to the pixel coordinates of the optical remote sensing image, thereby obtaining the projection positions of the photons corresponding to multiple candidate light spots in the optical remote sensing image.

[0070] Specifically, the transformation from ground latitude, longitude, and altitude to pixel coordinates can be performed using an image RPC model. Let the first... The ground coordinates of the photons corresponding to the candidate light spots are: First, normalize the latitude and longitude using the following formula:

[0071]

[0072]

[0073]

[0074] in: These represent the latitude normalized offset and the scale factor, respectively. These represent the longitude normalized offset and the scale factor, respectively. These represent the elevation normalization offset and the scale factor, respectively.

[0075] The aforementioned offset and scale factor are preferably the corresponding values ​​provided in the image RPC parameter file; if they need to be reset, they can be determined according to the latitude, longitude, and elevation range of the target image coverage or the target processing area, so that the normalized latitude, longitude, and elevation variables fall within the effective applicable range of the RPC model.

[0076] make This represents the basis function vector of the cubic rational polynomial of the RPC. Normalized row and column coordinates are calculated using the RPC coefficients, and then inverse normalized to obtain the nth... Each candidate light spot corresponds to the projection position of a photon in an optical remote sensing image. The pixel coordinates adopt the same indexing convention as the target image processing system, and can use a 0-based row and column coordinate system with the top left pixel as the origin.

[0077] In other implementations, a rigorous sensor model, a geometric correction model, or other equivalent positioning models can be used to complete the mapping from ground three-dimensional coordinates to image pixel coordinates.

[0078] In operation S402, the device parameters of the lidar are obtained, the spot radius of the lidar is confirmed according to the device parameters, and the spot range of multiple candidate spots corresponding to photons in the optical remote sensing image is determined according to each projection position and the spot radius.

[0079] Equipment parameters may include one or more of the following: orbital altitude, laser divergence angle, transmission and reception geometry, nominal spot size of the sensor system, etc. The approximate radius of the single-pulse spot on the ground can be calculated based on the orbital altitude and divergence angle; in other embodiments, publicly available sensor parameters, empirical nominal values, or calibrated spot radii can be directly used. In one embodiment of this application, when the spaceborne photon counting lidar data uses ICESat-2 / ATLAS data, the strong and weak beams can use their nominal ground spot diameter of approximately 17m as the spot diameter to construct the spot coverage pixel set. If calibration results or data products support this, it is also permissible to set parameters separately for each beam.

[0080] After obtaining the photons corresponding to the flat candidate light spots, one or more photons corresponding to the flat candidate light spots can be merged into candidate light spots based on the laser pulse, the continuity of photon observation time, spatial proximity, or projection position aggregation relationship. When a candidate light spot is represented by a single photon corresponding to a flat candidate light spot, the projection pixel position of that photon is taken as the projection center of the light spot; when a candidate light spot corresponds to multiple photons corresponding to flat candidate light spots, the median, mean, or weighted center of the coordinates of the multiple photon projection rows and columns is preferably taken as the projection center of the light spot.

[0081] The range of the light spot is constructed with the center of the light spot projection as the center and the radius of the light spot as the radius. :

[0082]

[0083] in: This represents the set of image pixels covered by a candidate light spot. Represents a set The Middle The row and column coordinates of each pixel This indicates the total number of pixels contained within the area of ​​the light spot. For pixels at the edge of the light spot, pixels whose centers fall within the area of ​​the light spot can be included in the set, or their inclusion can be determined based on the area overlap ratio between the edge pixels and the area covered by the light spot.

[0084] In operation S403, multiple pixels within the light spot range are acquired, and the spectral and texture features of each pixel are extracted. Based on the spectral and texture features, it is determined whether each pixel belongs to a hard surface pixel. The spectral features include at least grayscale, brightness, and image channel differences, and the texture features include at least contrast and the homogeneity of the pixel in the second preset neighborhood.

[0085] Specifically, for the first Each image channel, pixel The single-channel grayscale can be represented as When the input image contains When there are multiple channels, the brightness of a pixel can be defined as:

[0086]

[0087] in, Represents a pixel brightness, Indicates the number of channels. This represents the brightness of the c-th channel.

[0088] To enhance the spectral difference information between different channels, channel normalized difference features can be defined:

[0089]

[0090] in, Represents a pixel In the passage With channel Normalized difference characteristics between them For channel grayscale, For channel grayscale, This indicates an extremely small positive number to prevent the denominator from being zero. When the input image is a three-channel visible light image, the brightness can be represented by the average of the RGB three channels; when the input image is a multispectral image, the brightness can also be represented by a weighted combination of all channels or pre-selected channels.

[0091] For texture features, a second preset neighborhood is constructed centered on the pixel, and gray-level co-occurrence relationships are statistically analyzed within the second preset neighborhood. Let a certain pixel... The local neighborhood gray-level co-occurrence matrix is Then its contrast and homogeneity can be expressed as follows:

[0092]

[0093]

[0094] in, Represents a pixel Contrast Represents a pixel homogeneity, These represent the rows and columns of the gray-level co-occurrence matrix, which can be statistically analyzed and averaged in several preset directions, including 0°, 45°, 90°, and 135°. In other embodiments, entropy, energy, correlation, or other texture statistical features can also be further extracted.

[0095] Furthermore, the spectral and texture features of pixels within the light spot range can be statistically summarized to form light spot-level statistical features, which can be used to assist in analysis or threshold determination. For example, for the first... For each channel, the mean and standard deviation of the light spot level can be defined as follows:

[0096]

[0097] For brightness, the mean and standard deviation of spot-level brightness can also be defined:

[0098] in, These represent the first and second light spots within the light spot area, respectively. The mean and standard deviation of each channel. These represent the mean and standard deviation of the brightness within the light spot area, respectively.

[0099] In some embodiments, determining whether a pixel belongs to a hard surface pixel based on spectral and texture features may further include: determining that a pixel belongs to a hard surface when each pixel simultaneously satisfies the following conditions: brightness is within a preset brightness range, image channel difference is within a preset difference range, contrast is less than or equal to a preset contrast threshold, and homogeneity is greater than or equal to a preset homogeneity threshold. That is, the following conditions are simultaneously met:

[0100]

[0101]

[0102]

[0103] in, These represent the lower and upper limits of the preset brightness range, respectively. Represents a pixel brightness, These represent the lower and upper limits of the preset difference range, respectively. Indicates the preset contrast threshold. This represents a preset homogeneity threshold. This threshold can be determined based on image normalization results and statistical results of typical samples. When the brightness feature is normalized to the range of 0 to 1, The preferred value is 0.10 to 0.30. The preferred value is 0.70 to 0.95; the preferred discrimination interval for channel normalized difference features is -0.30 to 0.30; when texture contrast and homogeneity are normalized to the range of 0 to 1, The preferred value is 0.20 to 0.50. The preferred value is 0.50 to 0.80.

[0104] In practical applications, the above thresholds can be adaptively adjusted based on sensor type, image radiometric scale, land cover type, and sample statistical distribution. Brightness feature constraints are used to remove water bodies, strong shadows, or saturated bright areas; preset difference ranges are used to remove non-hard land cover such as vegetation and water bodies with significant differences in spectral response; texture features are used to retain stable hard land cover areas with relatively uniform texture and gentle changes.

[0105] In other implementations, hard surface discrimination can also be determined by the category labels or probability results output by supervised classification models, traditional machine learning models, or semantic segmentation models.

[0106] This allows for the generation of hard surface masks constrained by the light spot range. When pixel When it is located within the coverage area of ​​the candidate light spot and meets the hard surface discrimination condition. When pixel When the light spot is within the coverage area of ​​the candidate light spot, but does not meet the hard surface discrimination condition When pixel When not within the coverage area of ​​any candidate spot .in, Represents a pixel It was identified as a hard surface pixel; Represents a pixel Located within the light spot area, but classified as a non-hard surface pixel; Represents a pixel Located outside the light spot range, it does not participate in the hard surface discrimination.

[0107] In operation S404, the ratio of the number of pixels belonging to hard surfaces to the total number of pixels within the light spot area is calculated. When the ratio is greater than a preset ratio threshold, the pixels within the light spot area are determined to meet a preset condition. Specifically, the proportion of pixels within the light spot area that are determined to be hard surfaces is statistically analyzed:

[0108]

[0109] in, This represents the proportion of pixels on a hard surface. This represents the total number of pixels within the area of ​​the light spot. This represents the total number of hard surface pixels within the light spot area.

[0110] Set hard pixel ratio threshold When satisfied At that time, it is determined that the pixels within the light spot range meet the preset conditions. Among them, It can be set according to the actual application scenario, and the preferred value range is 0.6 to 0.8.

[0111] Next, for the light spot on the hard surface, within its corresponding spot range Inside, based on a hard surface mask Image feature extraction operators are used to extract image feature points located within hard surface regions, forming a set of hard surface image feature points. These operators can be corner operators, edge feature operators, locally salient feature operators, or other operators suitable for extracting stable image feature points from high-resolution images. Only those satisfying the specified conditions are retained. Image feature points extracted at the pixel location or within its local neighborhood.

[0112] In operation S105, representative elevations of multiple signal photon elevation samples within the target spot range that meet preset conditions are selected. The pixel coordinates of the hard surface pixels in the optical remote sensing image within the target spot range are combined with the representative elevations to obtain joint feature points.

[0113] In some embodiments, selecting representative elevations from multiple signal photon elevation samples within a target spot range that meet preset conditions may further include:

[0114] The median of the elevation samples of each signal photon is used as the reference elevation. ; Calculate the deviation between each signal photon elevation sample and the reference elevation: And calculate the standard deviation of multiple deviations. ; with preset correction coefficient The standard deviation is corrected, and signal photon elevation samples whose absolute deviation is greater than the corrected standard deviation are removed. That is, signal photon elevation samples that meet the following conditions will be identified as outliers and removed: ,in, This represents the preset correction coefficient, with a preferred value range of 1 to 3.

[0115] Calculate the mean of the remaining signal photon elevation samples to obtain the representative elevation. In other implementations, the representative elevation may also be obtained using the median, weighted average, or other robust location estimation methods.

[0116] Finally, the pixel coordinates of the hard surface pixels in the optical remote sensing image are combined with the representative elevation to make it appear in the image coordinates. Based on this, elevation attributes provided by spaceborne lidar are added to form joint feature points. For multiple image feature points belonging to the same spot range, they can share the same control elevation, thus forming a set of joint feature points that simultaneously contain image two-dimensional position and laser elevation information.

[0117] Figure 5 A flowchart illustrating a method for constructing multi-level reference images according to an embodiment of this application is shown schematically.

[0118] like Figure 5 As shown, this application also provides a method for constructing multi-level reference images, including operations S501 to S502.

[0119] In operation S501, joint feature points are obtained, which are extracted according to the above-mentioned method for joint extraction of spaceborne lidar point cloud and image features.

[0120] This operation S501 is applied to the joint extraction method of spaceborne lidar point cloud and image features in any of the above embodiments, and the technical effect achieved is also similar. For specific details, please refer to the embodiment section of the joint extraction method of spaceborne lidar point cloud and image features, which will not be repeated here.

[0121] In operation S502, a multi-level reference image is constructed based on the joint feature points.

[0122] In some embodiments, the above operation S502 may include at least one of the following.

[0123] The first method involves taking the location of the joint feature points in the optical remote sensing image as the center, and expanding the corresponding preset number of pixels along the row and column directions in the optical remote sensing image to obtain a point-level joint feature reference image.

[0124] Specifically, using joint feature points Image coordinates Centered on the original image, a predetermined number of pixels are expanded to both sides along the row and column directions to construct a fixed-size point-level cropping window, resulting in a point-level joint feature reference image. The predetermined number of pixels is preferably 32–256 pixels per side; in high-resolution optical image scenes with a resolution better than 1 meter, 128 pixels per side is preferable. In practical applications, this can be adjusted according to image resolution, texture richness, and subsequent matching requirements. The point-level joint feature reference image includes at least one or more of the following: a point-level image block centered on the joint feature point; the image row and column coordinates of the joint feature point; the latitude and longitude coordinates of the joint feature point; the control elevation of the joint feature point; source image identifier, sensor information, and imaging time; quality assessment information or land cover type information.

[0125] The second approach involves dividing the target area into multiple sub-regions, extracting multiple joint feature points within each sub-region, and obtaining a block-level joint feature reference image.

[0126] Specifically, the ground space can be divided into several block-level units according to a regular grid. The ground coverage area of ​​each block-level unit can be a fixed size, such as 1km × 1km, but is not limited to this. For each block-level unit, the corresponding block-level image is obtained by cropping, and all joint feature points falling within the range of that block-level unit are retrieved to form a block-level joint feature reference image. To ensure that the block-level joint feature product has effective geometric control capabilities, each valid block-level unit should contain no less than 3 joint feature points. The block-level joint feature reference image should include at least one or more of the following: the block-level image corresponding to the range of the block-level unit; a list of all joint feature points within the block-level range; the image location, geographical location, and control elevation of each joint feature point; block-level unit boundary information; source image identification and block-level index information; and quality statistics within the block-level range.

[0127] The third method involves extracting multiple joint feature points within the target area to obtain a scene-level joint feature reference image.

[0128] Specifically, based on the original high-resolution image of the entire scene, all joint feature points within the image range are integrated to construct a scene-level joint feature reference image. In one embodiment, the scene-level joint feature reference image does not perform local cropping on the original image, but retains the entire scene image range, while simultaneously establishing an index of all joint feature points within the coverage area of ​​the scene image. The scene-level joint feature reference image includes at least one or more of the following: the original image of the entire scene or reference information for the entire scene image; an index of all joint feature points within the entire scene range; two-dimensional image positions and three-dimensional control information of each joint feature point; scene metadata, including row and column dimensions, spatial resolution, coordinate system, sensor information, and imaging time; and a scene-level description file or a global index file.

[0129] Furthermore, a corresponding description file can be generated for each point-level joint feature reference image, block-level joint feature reference image, or scene-level joint feature reference image. The description file can adopt XML, JSON, text tables, or other structured metadata formats to record the spatial extent of the reference image, image index information, and joint feature point attribute information.

[0130] For point-level, block-level, and scene-level joint feature products, a unified naming convention, directory organization method, and metadata indexing system can be established to achieve standardized management of joint feature products. The file naming convention can include one or more of the following information: source image number; product level identifier; joint feature point number; block-level grid number; scene-level image number; time identifier or track identifier. Simultaneously, a directory organization method can be established according to "scene-level-type" or other hierarchical structures to maintain clear grouping of point-level, block-level, and scene-level joint feature products in physical storage.

[0131] Based on the aforementioned structured organization and metadata indexing information, joint feature products can be retrieved, located, and invoked. In one embodiment, joint feature products can be retrieved based on one or more of the following conditions: specifying a spatial range; specifying the location of joint feature points; specifying a block-level grid number; specifying a scene-level image number; specifying a product level; specifying a quality threshold or attribute condition. After retrieving the target joint feature product, the corresponding image data, joint feature point list, description file, and its attribute information can be output for geometric correction, image registration, control information analysis, 3D positioning, or other downstream applications. In other embodiments, an interface-based invocation method can also be provided, allowing access and transmission of joint feature products through file interfaces, database interfaces, network service interfaces, or other data exchange methods.

[0132] In summary, this application provides a method for jointly extracting features from spaceborne lidar point clouds and imagery. By acquiring optical remote sensing images covering the target area and spaceborne photon-counting lidar data, multiple sliding windows are constructed along the orbital direction to denoise the signal photons. Candidate photons corresponding to light spots located in stable areas are selected by calculating local terrain features. The surface type is determined by judging whether the pixels within the light spot range meet preset conditions. Finally, robust statistical processing is used to obtain representative elevations, which are then combined with the pixel coordinates of hard surface pixels to obtain joint feature points. Based on these joint feature points, multi-level reference images at the point, block, and scene levels are further constructed, forming a unified, standardized, easily portable, and batch-processable joint feature product system.

[0133] As can be seen from the above description, the embodiments of this application achieve at least the following technical effects:

[0134] (1) Regarding the flatness and geometric stability of the control area, this application performs B-spline fitting on the candidate signal photon high-order sequence, calculates the slope along the track and statistically analyzes the standard deviation or root mean square error of the fitting residual, and adopts the joint constraint of slope threshold and flatness threshold to retain only the flat section along the track with small slope and low fitting residual, and automatically eliminates unstable areas such as undulating terrain, building facades and shading edges, thereby ensuring the geometric flatness and stability of the control area and providing high-quality regional support for subsequent geometric constraints.

[0135] (2) Regarding the consistency of the spatial scale of laser and image, this application explicitly introduces the actual size of the spaceborne laser spot on the ground. Using the spot as the basic spatial unit, a set of pixels is constructed within the coverage area of ​​the spot, and the multi-channel grayscale, normalized difference index and texture features are statistically analyzed. This establishes a correspondence between the laser altimetry information and the image spectral / texture features on the "spot scale", effectively avoiding the spatial scale mismatch problem that is easy to occur in the traditional "single photon - single pixel" matching, and improving the consistency and physical interpretability of the joint features of laser and image.

[0136] (3) Regarding the semantic stability and temporal reusability of ground features, this application constructs a hard surface discrimination strategy by combining brightness, spectral differences and texture features, and generates a hard surface mask at the light spot scale within the entire image range; further, it calculates the proportion of hard pixels within the light spot, retains only hard surfaces with stable geometric shapes and relatively slow temporal changes, such as roads and squares, and automatically removes areas with unstable geometric and temporal characteristics, such as vegetation, water bodies, and cultivated land, thereby improving the long-term stability and cross-temporal reusability of joint feature points from a semantic perspective.

[0137] (4) Regarding the robustness of elevation control, this application calculates the median of the elevation samples of multiple signal photons within the spot as the initial elevation for each joint feature candidate point. Then, based on the deviation from the median and its standard deviation, abnormal elevation points are eliminated. Finally, the remaining samples are averaged to obtain the final control elevation value, and this control elevation is uniformly assigned to all hard surface image feature points within the spot. This effectively suppresses the influence of isolated abnormal photons, making the elevation attributes of the joint feature points insensitive to noise and robust to local anomalies.

[0138] (5) Regarding the form of joint feature products and data standardization, this application constructs three types of reference images—point-level, block-level, and scene-level—on the original image based on joint feature points, and generates a matching description file for each type of reference image. It establishes a correspondence between “pixel-latitude-longitude-elevation-attribute” from point-level, block-level to scene-level, forming a joint feature product system with unified structure, standard format, easy portability, and batch processing.

[0139] Figure 6 The diagram illustrates a block diagram of an electronic device suitable for implementing a method for jointly extracting point clouds and image features from spaceborne lidar or for constructing multi-level reference images, according to embodiments of this application. Figure 6 The electronic device shown is merely an example and should not impose any limitation on the functionality and scope of use of the embodiments of this application.

[0140] like Figure 6 As shown, an electronic device 600 according to an embodiment of this application includes a processor 601, which can perform various appropriate actions and processes according to a program stored in a read-only memory (ROM) 602 or a program loaded from a storage portion 608 into a random access memory (RAM) 603. The processor 601 may include, for example, a general-purpose microprocessor (e.g., a CPU), an instruction set processor and / or an associated chipset and / or a special-purpose microprocessor (e.g., an application-specific integrated circuit (ASIC)), etc. The processor 601 may also include onboard memory configured for caching purposes. The processor 601 may include a single processing unit or multiple processing units configured to perform different actions of the method flow according to an embodiment of this application.

[0141] RAM 603 stores various programs and data required for the operation of electronic device 600. Processor 601, ROM 602, and RAM 603 are interconnected via bus 604. Processor 601 executes various operations of the method flow according to embodiments of this application by executing programs in ROM 602 and / or RAM 603. It should be noted that the programs may also be stored in one or more memories other than ROM 602 and RAM 603. Processor 601 may also execute various operations of the method flow according to embodiments of this application by executing programs stored in said one or more memories.

[0142] According to embodiments of this application, the electronic device 600 may further include an input / output (I / O) interface 605, which is also connected to a bus 604. The electronic device 600 may also include one or more of the following components connected to the input / output (I / O) interface 605: an input section 606 including a keyboard, mouse, etc.; an output section 607 including a cathode ray tube (CRT), liquid crystal display (LCD), etc., and a speaker, etc.; a storage section 608 including a hard disk, etc.; and a communication section 609 including a network interface card such as a LAN card, modem, etc. The communication section 609 performs communication processing via a network such as the Internet. A drive 610 is also connected to the input / output (I / O) interface 605 as needed. A removable medium 611, such as a disk, optical disk, magneto-optical disk, semiconductor memory, etc., is installed on the drive 610 as needed so that computer programs read from it can be installed into the storage section 608 as needed.

[0143] The flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of this application. In this regard, each block in a flowchart or block diagram may represent a module, segment, or portion of code containing one or more executable instructions configured to perform a specified logical function. It should also be noted that in some alternative implementations, the functions indicated in the blocks may occur in a different order than those indicated in the drawings. For example, two consecutively indicated blocks may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. It should also be noted that each block in a block diagram or flowchart, and combinations of blocks in a block diagram or flowchart, may be implemented using a dedicated hardware-based system that performs the specified function or operation, or using a combination of dedicated hardware and computer instructions.

[0144] Those skilled in the art will understand that the features described in the various embodiments of this application can be combined and / or combined in various ways, even if such combinations or combinations are not explicitly described in this application. In particular, the features described in the various embodiments of this application can be combined and / or combined in various ways without departing from the spirit and teachings of this application. All such combinations and / or combinations fall within the scope of this application.

[0145] The embodiments of this application have been described above. However, these embodiments are merely illustrative and not intended to limit the scope of this application. Although various embodiments have been described above, this does not mean that the measures in the various embodiments cannot be used advantageously in combination. Without departing from the scope of this application, those skilled in the art can make various substitutions and modifications, all of which should fall within the scope of this application.

Claims

1. A method for jointly extracting point clouds and image features from a spaceborne lidar system, characterized in that, include: Acquire optical remote sensing images and spaceborne photon counting lidar data covering the target area, wherein the spaceborne photon counting lidar data includes multiple signal photons; Multiple sliding windows are constructed along the track direction, and the multiple signal photons are denoised according to the multiple sliding windows to obtain multiple candidate signal photons in each sliding window. Calculate the local terrain features of the land area corresponding to the candidate signal photons in each sliding window, and filter out the candidate light spot corresponding to the stable region from multiple candidate signal photons in each sliding window based on the local terrain features; Determine whether the pixels corresponding to the photons of the multiple candidate light spots in the light spot range covered by the optical remote sensing image meet the preset conditions. The preset conditions are used to determine the land surface type where the light spot range is located based on the spectral and texture features of each pixel. Representative elevations of multiple signal photon elevation samples within the target spot range that meet the preset conditions are selected. The pixel coordinates of the hard surface pixels in the optical remote sensing image within the target spot range are combined with the representative elevations to obtain joint feature points. The local terrain features include slope and flatness. Calculating the local terrain features of the surface area corresponding to the candidate signal photon within each sliding window includes: A fitting curve is constructed based on the orbital distance of the candidate signal photons within each sliding window, and the average value of the absolute value of the first derivative of the fitting curve is calculated to obtain the slope of the land surface area corresponding to the candidate signal photons within each sliding window. Let the first sliding window be... The orbital distance of the candidate signal photons is Elevation is Then the fitted curve can be expressed as: in: Indicates the first The fitting curve obtained by fitting candidate signal photons. Indicates the first A candidate signal photon along-track elevation fitting model; The residual sequence is calculated based on the fitted curve and the actual elevation of the candidate signal photons, and the residual standard deviation is calculated based on the residual sequence to obtain the smoothness. The residual sequence can be represented as: in: Indicates the first The fitting residuals of the candidate signal photons Indicates the first The actual elevation of each candidate signal photon.

2. The method according to claim 1, characterized in that, The process involves constructing multiple sliding windows along the track direction, and performing denoising processing on the multiple signal photons according to the multiple sliding windows to obtain multiple candidate signal photons within each sliding window, including: The signal photons are sorted in chronological order, and the distances between adjacent signal photons after sorting are accumulated to obtain the cumulative distance along the track for each signal photon. The cumulative distance along the track is divided into multiple sliding windows along the track direction, and the number of neighborhood points in the first preset neighborhood of the signal photon in each sliding window is counted. A noise threshold is determined based on the number of neighboring points in each of the aforementioned windows, and signal photons whose number of neighboring points in each of the aforementioned sliding windows is greater than or equal to the noise threshold are identified as candidate signal photons.

3. The method according to claim 2, characterized in that, The step of filtering candidate light spot photons located in stable regions from candidate signal photons within each sliding window based on the local terrain features includes: Candidate signal photons whose absolute slope value within each sliding window is not greater than a preset slope threshold and whose flatness is not greater than a preset flatness threshold are identified as candidate spot photons located in the stable region.

4. The method according to claim 3, characterized in that, After filtering out the candidate light spot corresponding to the stable region from multiple candidate signal photons within each sliding window based on the local terrain features, the method further includes: The adjacent sliding windows that continuously meet the screening conditions along the track direction corresponding to the photons of each candidate spot are merged to form a smooth section window along the track. The screening conditions are that the absolute value of the slope of the candidate signal photons in the sliding window is less than or equal to a preset slope threshold and the smoothness is less than or equal to a preset smoothness threshold. The track-side smooth section window is filtered according to a preset continuity constraint condition. The photons corresponding to candidate light spots in the track-side smooth section window that satisfy the preset continuity constraint condition are retained. The photons corresponding to candidate light spots filtered according to the local terrain features are updated. The preset continuity constraint condition includes at least one of the following: a continuous window number threshold, a continuous track length threshold, and an effective candidate photon density threshold.

5. The method according to claim 1, characterized in that, The step of determining whether the pixels corresponding to the photons of the plurality of candidate light spots in the light spot range covered by the optical remote sensing image meet the preset conditions includes: The coordinates of the photons corresponding to the multiple candidate light spots are mapped to the pixel coordinates of the optical remote sensing image to obtain the projection position of the photons corresponding to the multiple candidate light spots in the optical remote sensing image. The device parameters of the lidar are obtained, the spot radius of the lidar is determined according to the device parameters, and the spot range of the photons corresponding to the multiple candidate spots in the optical remote sensing image is determined according to each projection position and the spot radius. Multiple pixels within the light spot range are acquired, and the spectral and texture features of each pixel are extracted. Based on the spectral and texture features, it is determined whether each pixel belongs to a hard surface pixel. The spectral features include at least grayscale, brightness, and image channel differences, and the texture features include at least contrast and pixel homogeneity in a second preset neighborhood. Calculate the ratio of the number of pixels belonging to the hard surface to the total number of pixels within the light spot range. When the ratio is greater than a preset ratio threshold, determine that the pixels within the light spot range meet the preset condition.

6. The method according to claim 5, characterized in that, The step of determining whether each pixel belongs to a hard surface pixel based on the spectral features and the texture features includes: When the spectral and texture features of each pixel simultaneously satisfy the following conditions: the brightness is within a preset brightness range, the image channel difference is within a preset difference range, the contrast is less than or equal to a preset contrast threshold, and the homogeneity is greater than or equal to a preset homogeneity threshold, the pixel is determined to belong to a hard surface.

7. The method according to claim 1, characterized in that, The representative elevations of multiple signal photon elevation samples within the target spot range that meet the preset conditions include: The median of each of the aforementioned signal photon elevation samples is used as the reference elevation; Calculate the deviation between each of the signal photon elevation samples and the reference elevation, and calculate the standard deviation of multiple deviations; The standard deviation is corrected using a preset correction factor. Signal photon elevation samples whose absolute value of the deviation is greater than the corrected standard deviation are removed. The mean of the remaining signal photon elevation samples is calculated to obtain the representative elevation.

8. A method for constructing multi-level reference images, characterized in that, include: The joint feature points are obtained, wherein the joint feature points are extracted by the method for joint extraction of star-borne lidar point cloud and image features according to any one of claims 1-7; A multi-level reference image is constructed based on the joint feature points.

9. The method according to claim 8, characterized in that, A multi-level reference image is constructed based on the joint feature points, including at least one of the following: Using the position of the joint feature point in the optical remote sensing image as the center, the corresponding preset number of pixels are extended along the row and column directions in the optical remote sensing image to obtain a point-level joint feature reference image. The target region is divided into multiple sub-regions, and multiple joint feature points are extracted from each sub-region to obtain a block-level joint feature reference image. Multiple joint feature points are extracted within the target area to obtain a scene-level joint feature reference image.

10. An electronic device, comprising: One or more processors; Memory, used to store one or more computer programs. The characteristic is that the one or more processors execute the one or more computer programs to implement the steps of the method according to any one of claims 1 to 7 or the method according to any one of claims 8 to 9.

Images