An airborne laser radar depth sounding system intensity data extraction method and system

Abstract

The application provides a kind of airborne laser radar depth system intensity data extraction method and system, first to the intensity data obtained by airborne laser radar depth system;Through the feature extraction based on full waveform data and the setting of judging condition to identify the invalid waveform disturbed by noise, and then eliminate the corresponding invalid intensity data;Through the water-land segmentation based on satellite-borne multispectral image data auxiliary, the water-land range segmentation of airborne laser radar intensity data is realized;For water area range data, based on full waveform data decomposition, only the water surface reflection echo waveform of water surface echo is identified, and then the corresponding intensity data is eliminated, and finally the intensity data that can effectively represent the attributes of land and seabed objects is obtained without the interference of mixed intensity data, which is helpful for further object recognition, classification and segmentation.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of laser radar, in particular to a method and system for extracting intensity data of an airborne laser radar bathymetry system based on image data and full waveform data assistance. BACKGROUND

[0002] Airborne laser radar bathymetry technology is known as the dawn of island and coastal zone mapping, and is more and more applied to island and coastal zone geographic spatial information investigation due to its advantages of water-land integration, high precision and high resolution. The airborne laser radar bathymetry system not only provides geometric information (X, Y, Z), but also provides the return intensity (physical power) of the reflected laser pulse, which is regarded as a kind of spectral information and can be used to characterize the properties of seabed objects. The airborne laser radar bathymetry intensity data is of great significance to coastal object classification, seabed bottom survey, biological habitat mapping and coral reef monitoring.

[0003] However, due to the influence of environmental noise and sea surface reflection signals, the intensity data of the laser radar bathymetry is mixed with a large amount of invalid intensity data interference, making it difficult to express the object category attributes and further information mining. The intensity data recorded in the Las file generated by the airborne laser radar bathymetry system is not labeled or removed without related invalid intensity data. There is no effective solution to this problem in the existing research on airborne laser radar bathymetry intensity data correction and application. SUMMARY

[0004] In order to overcome the shortcomings of the prior art, the purpose of the present application is to provide a method and system for extracting intensity data of an airborne laser radar bathymetry system.

[0005] To achieve the above-mentioned purpose, the present application provides the following solutions:

[0006] In one aspect, the present application provides a method for extracting intensity data of an airborne laser radar bathymetry system, comprising:

[0007] parsing the airborne laser radar bathymetry Las file to obtain terrain data, intensity data and full waveform data;

[0008] performing feature extraction and threshold determination on the full waveform data to identify invalid waveforms and valid waveforms, and removing the intensity data corresponding to the invalid waveforms;

[0009] registering the terrain data and multispectral image data, and performing water-land segmentation on the multispectral image data based on a water body index to obtain a water-land segmentation result; the water-land segmentation result includes a water area range and a land area range;

[0010] According to the water-land segmentation result, the effective waveforms are classified to obtain water area waveforms and land area waveforms, and based on a constructed waveform decomposition model, waveform decomposition and Gaussian echo component number judgment are performed according to the water area waveforms to obtain water surface reflection waveforms and water body transmission waveforms, and intensity data corresponding to the water surface reflection echo waveforms are removed;

[0011] The intensity data corresponding to the remaining land range and the intensity data corresponding to the water body transmission waveforms are merged to obtain an extraction result.

[0012] Preferably, feature extraction and threshold judgment are performed on the full waveform data to identify invalid waveforms and effective waveforms, and intensity data corresponding to the invalid waveforms are removed, including:

[0013] The positions of the maximum and minimum intensity of the full waveform data are searched to obtain the minimum intensity and maximum intensity corresponding time and time deviation;

[0014] The minimum intensity and maximum intensity corresponding time and time deviation of the full waveform data obtained under noise interference are statistically analyzed to determine the best threshold;

[0015] According to the judgment condition set according to the best threshold, invalid waveforms and effective waveforms obtained under noise interference are identified according to the judgment condition, and intensity data corresponding to the invalid waveforms are removed.

[0016] Preferably, the terrain data and multispectral image data are registered, and water-land segmentation is performed on the multispectral image data based on a water body index to obtain a water-land segmentation result, including:

[0017] The multispectral image data and the terrain data are accurately registered based on a visualization method to obtain geometrically corrected multispectral image data;

[0018] Water body index calculation and image data water-land segmentation are performed on the geometrically corrected multispectral image data to obtain the water-land segmentation result of the terrain data.

[0019] Preferably, according to the water-land segmentation result, the effective waveforms are classified to obtain water area waveforms and land area waveforms, and based on a constructed waveform decomposition model, waveform decomposition and Gaussian echo component number judgment are performed according to the water area waveforms to obtain water surface reflection waveforms and water body transmission waveforms, and intensity data corresponding to the water surface reflection echo waveforms are removed, including:

[0020] A sea surface-water body-sea bottom-noise layered decomposition model is constructed;

[0021] Solving parameters of the sea surface-water body-sea bottom-noise layered decomposition model based on the LM algorithm, and sequentially fitting the full waveform data of the water area range to obtain a fitting result;

[0022] According to the fitting result parameters, classifying the water surface reflection waveform and the water body transmission waveform, and eliminating the intensity data corresponding to the water surface reflection echo waveform.

[0023] Preferably, the full waveform data is subjected to feature extraction and threshold determination to identify invalid waveforms and valid waveforms, and the intensity data corresponding to the invalid waveforms is eliminated, comprising:

[0024] A search function is written to sequentially calculate the minimum intensity and maximum intensity corresponding time and time deviation of each full waveform data, Wherein, f i (t) is the i-th waveform data, i is the waveform data identifier, i = 1, 2, …, N, N is the total number of full waveform data, t corresponds to the time, the waveform data recording start time is t n , is the intensity data recorded by the system at t is the maximum intensity corresponding to the acquisition time, is the minimum intensity corresponding to the acquisition time, and Δtime is the minimum intensity and maximum intensity acquisition time deviation;

[0025] By setting the judgment condition to screen the invalid waveform data that meets the condition, wherein if the maximum intensity of the waveform f i (t) is located in the front 4% of the entire recording time, that is The minimum intensity is located in the last 4% of the entire recording time, that is The minimum intensity and maximum intensity corresponding time deviation is not less than 90% of the entire recording time, that is Δtime > t n * 90%, f i (t) is considered to be an invalid waveform f i (t) 0 , otherwise it is judged to be a valid waveform f i (t) 1 ;

[0026] According to the determined category, the invalid waveform and the intensity data corresponding to the invalid waveform are eliminated, and the step of "writing a search function to sequentially calculate the minimum intensity and maximum intensity corresponding time and time deviation of each full waveform data" is returned until all the intensity data corresponding to the invalid waveform in the intensity data is eliminated.

[0027] Preferably, the topographic data and multispectral image data are registered according to the topographic data and multispectral image data, and water and land segmentation is performed on the multispectral image data based on a water body index to obtain a water and land segmentation result, comprising:

[0028] The topographic data corresponding to the effective waveform is registered with the multispectral image data collected in the same period through the registration processing, and a registration result is obtained;

[0029] The water body index Index is calculated NDWI , Index NDWI =(R Green -R NIR ) / (R Green +R NIR ); wherein R Green is the reflectivity recorded by the green band, R NIR is the reflectivity recorded by the near-infrared band, and if Index NDWI is greater than 0.15, it is a water body area, otherwise it is a land area;

[0030] Based on the registration result, the class label of the topographic data corresponding to the water body area of the multispectral image data is determined as water, and the corresponding full waveform data is and the class label of the topographic data corresponding to the land area is determined as land, and the corresponding full waveform data is

[0031] The full waveform data corresponding to the topographic data of the land area is retained The intensity data corresponding to the full waveform data is classified and removed.

[0032] Preferably, based on the constructed waveform decomposition model, waveform decomposition and Gaussian echo component number judgment are performed according to the water area waveform to obtain a water surface reflection waveform and a water body transmission waveform, and the intensity data corresponding to the water surface reflection echo waveform is removed, comprising:

[0033] A sea surface-water body-sea bottom-noise layered decomposition model is constructed, and the full waveform data corresponding to the water area is fitted based on the sea surface-water body-sea bottom-noise layered decomposition model; the expression of the sea surface-water body-sea bottom-noise layered decomposition model is:

[0034]

[0035] ​Wherein, y(t) is the superposition of sea surface echo, water scattering component, sea bottom echo and environmental noise, M is 1 for full waveform data directly reflected from the sea surface, only the sea surface echo component, refraction occurs at the sea surface, and then penetrates the water body to the sea bottom to reflect the full waveform data to the sensor, M is 2, when encountering phytoplankton, M is 3, wherein A j , μ j , ω j and γ j are the amplitude, peak position, half-wave width and waveform shape parameters of the jth Gaussian component, the water scattering component is defined as a triangle, represented by ∫(ax+b)+∫(cx+d), ax+b is the expression function of the rising waist of the triangle, and the start and end times are time1 and time2, respectively, cx+d is the expression function of the falling waist of the triangle, and the start and end times are time2 and time3, respectively, and the sum of the function integrals of the two in the duration is the water scattering component, is the noise component, and is the average noise level The integral in the entire recording time;

[0036] According to the preset target function, the parameters in the sea surface-water body-sea bottom-noise layered decomposition model are calculated;

[0037] According to the size of the fitted parameter M, the waveform category is classified, if M is equal to 1, the corresponding full waveform data is the water surface reflection echo waveform, otherwise the corresponding full waveform data is the water body transmission echo waveform;

[0038] The intensity data corresponding to the water surface reflection waveform identified is removed, and the remaining intensity data is retained.

[0039] Preferably, the expression of the target function is: min||f(t)-y(t)||; wherein f(t) is the original waveform data.

[0040] Preferably, the algorithm used to solve the target function is the LM algorithm.

[0041] Another aspect of the present application also provides an airborne laser radar depth sounding system intensity data extraction system, comprising:

[0042] A data analysis module is used to analyze the airborne laser radar depth sounding Las file to obtain terrain data, intensity data and full waveform data;

[0043] A waveform identification module is used to extract features and determine thresholds for the full waveform data to identify invalid waveforms and valid waveforms, and remove the intensity data corresponding to the invalid waveforms;

[0044] A water-land segmentation module is configured to register the terrain data and the multispectral image data, and perform water-land segmentation on the multispectral image data based on a water body index to obtain a water-land segmentation result, wherein the water-land segmentation result includes a water area range and a land area range.

[0045] A waveform classification module is configured to classify the effective waveforms based on the water-land segmentation result to obtain water area waveforms and land area waveforms, and perform waveform decomposition and Gaussian echo component quantity judgment on the water area waveforms based on a constructed waveform decomposition model to obtain water surface reflection waveforms and water body transmission waveforms, and eliminate intensity data corresponding to the water surface reflection waveforms.

[0046] A data extraction module is configured to combine intensity data corresponding to the remaining land area range and intensity data corresponding to the water body transmission waveforms to obtain an extraction result.

[0047] According to the embodiments of the present application, the following technical effects are achieved.

[0048] The present application provides a kind of airborne laser radar depth measurement system intensity data extraction method and system, first to the intensity data obtained by airborne laser radar depth measurement system;By the feature extraction and judgment condition setting based on full waveform data to identify the invalid waveform disturbed by noise, and then eliminate corresponding invalid intensity data;By the water-land segmentation of airborne laser radar intensity data assisted by multispectral image data based on satellite is realized;For water area range data, based on full waveform data decomposition identification only records the water surface echo of water surface reflection echo waveform, and then eliminate corresponding intensity data, finally obtain the intensity data that can effectively characterize the attribute of land and water bottom ground object without the interference of mixed intensity data, it is helpful to further carry out ground object identification, classification, segmentation. BRIEF DESCRIPTION OF DRAWINGS

[0049] In order to more clearly illustrate the technical solutions of the embodiments of the present application or the prior art, the drawings needed in the embodiments will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0050] Figure 1 The method schematic diagram provided for the embodiments of the present application;

[0051] Figure 2 The overall flow schematic diagram provided for the embodiments of the present application;

[0052] Figure 3 The invalid waveform distribution and morphology display diagram provided for the embodiments of the present application;

[0053] Figure 4 The image data and the water-land segmentation map of full waveform data provided for the embodiment of the present application;

[0054] Figure 5 The full waveform classification result map provided for the embodiment of the present application;

[0055] Figure 6 The original intensity data and the intensity image after removing invalid intensity interference provided for the embodiment of the present application. DETAILED DESCRIPTION

[0056] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative work fall within the protection scope of the present application.

[0057] The purpose of the present application is to provide an airborne laser radar sounding system intensity data extraction method and system, which can obtain airborne laser radar sounding system intensity data that can effectively express the properties of ground objects without invalid intensity interference.

[0058] In order to make the above-mentioned purposes, characteristics and advantages of the present application more obvious and easy to understand, the present application will be further described in detail below with reference to the drawings and specific embodiments.

[0059] Figure 1 The method schematic diagram provided for the embodiment of the present application is shown as Figure 1 The present application provides an airborne laser radar sounding system intensity data extraction method, which comprises:

[0060] Step 100: analyzing the airborne laser radar sounding Las file to obtain terrain data, intensity data and full waveform data;

[0061] Step 200: performing feature extraction and threshold determination on the full waveform data to identify invalid waveforms and valid waveforms, and removing the intensity data corresponding to the invalid waveforms;

[0062] Step 300: registering according to the terrain data and multispectral image data, and performing water-land segmentation on the multispectral image data based on a water body index to obtain a water-land segmentation result; the water-land segmentation result comprises a water area range and a land area range;

[0063] Step 400: classifying the effective waveform according to the water-land segmentation result to obtain water area waveform and land area waveform, and based on the constructed waveform decomposition model, performing waveform decomposition and Gaussian echo component number judgment according to the water area waveform to obtain water surface reflection waveform and water body transmission waveform, and removing the intensity data corresponding to the water surface reflection echo waveform;

[0064] Step 500: merging the intensity data corresponding to the remaining land range and the effective intensity data corresponding to the water body transmission waveform to obtain an extraction result.

[0065] Airborne laser radar sounding technology has become an ideal choice for island and coastal zone mapping due to its seamless water-land mapping advantage, and has received more and more attention. The intensity data recorded by the system can represent the attribute information of the ground object, which is of great significance for ground object classification and management. However, the intensity data recorded by the system is disturbed by environmental noise and sea surface return signals, resulting in a large amount of invalid intensity information, which causes great confusion in the use of intensity data. At present, neither commercial software nor scientific research has made relevant beneficial achievements. In order to solve this problem, the present application provides a method, and the specific flow chart is shown in Figure 2 which can be summarized as: Las file analysis to obtain original intensity data; invalid waveform identification and corresponding intensity removal based on waveform characteristics; water-land segmentation assisted by satellite image data; water surface reflection waveform identification and corresponding intensity data removal based on waveform decomposition; and finally obtaining effective airborne laser radar sounding intensity data. The specific details of each process are as follows:

[0066] S101, obtaining original intensity data by analyzing airborne laser radar sounding Las file

[0067] Specifically, the Las file analysis of airborne laser radar sounding is realized by writing corresponding code according to the Las file storage format compiled by the American Society for Photogrammetry and Remote Sensing (ASPRS), and the recorded coordinates XYZ, intensity data I and corresponding full waveform data FW are obtained, each full waveform data corresponding to a coordinate and an intensity.

[0068] S102, invalid echo identification and corresponding intensity data removal by full waveform feature extraction and threshold setting

[0069] Specifically, the process S102 includes:

[0070] Step one: by writing a maximum and minimum value search function to calculate the minimum intensity and maximum intensity corresponding time and time deviation of each full waveform data in turn:

[0071]

[0072]

[0073] f i (t) is the i-th waveform data, i = 1, 2, …, N, N is the total number of all waveform data, f i (t) ∈ FW, wherein t corresponds to time, the waveform data records the starting time t1 and the ending time t n , is the intensity data recorded by the system at t1, is the maximum intensity corresponding to the collection time, is the minimum intensity corresponding to the collection time, and Δtime is the time deviation of the minimum intensity and the maximum intensity.

[0074] Step two: filter the invalid waveform data that meets the condition by setting the judgment condition, wherein if the maximum intensity i of the waveform f (t) is recorded at the first 4% of the entire recording time, i.e. the minimum intensity is recorded at the last 4% of the entire recording time, i.e. the time deviation of the minimum intensity and the maximum intensity is not less than 90% of the entire recording time, i.e. Δtime > t n *90%, then f i (t) is considered as an invalid waveform f i (t) 0 , otherwise it is determined as an effective waveform f i (t) 1 ;

[0075] Step three: execute steps one to two on f i (t) in turn, f i (t) 0 with the determination category of invalid waveform constitutes FW 0 , f i (t) 1 with the determination category of effective waveform constitutes FW 1 , the corresponding intensity data I 0 and topographic data XYZ 0 of FW 1 are removed, and the corresponding intensity data I 1 and topographic data XYZ 1 of the effective waveform FW NDWI are temporarily retained, Figure 3 which shows the coordinate information corresponding to the invalid waveform identified by the airborne laser radar sounding data in the test area and the shape of 10 full waveforms.

[0076] S103, water-land segmentation of airborne LiDAR bathymetry data by registered spaceborne multispectral image data

[0077] The process S103 specifically includes:

[0078] Step one: precise registration of spaceborne multispectral image data and airborne LiDAR bathymetry data collected at the same time through visualization, so as to realize accurate correspondence of the two in coordinates;

[0079] Step two: calculation of water index Index based on green and near-infrared band information of spaceborne multispectral image data NDWI :

[0080] Index NDWI =(R Green -R NIR ) / (R Green +R NIR )

[0081] Wherein, R Green is the reflectivity recorded by the green band, R NIR is the reflectivity recorded by the near-infrared band, if Index NDWI is greater than 0.15, it is a water area, otherwise it is a land area;

[0082] Step three: based on the registration result of step one, the image data is divided into airborne LiDAR bathymetry data corresponding to the water area, and the category label is water, wherein the corresponding full waveform data is The terrain data is The intensity data is Similarly, the airborne LiDAR bathymetry data corresponding to the land area is assigned a category label of land The terrain data is The intensity data is The distribution of image data and airborne LiDAR data of the experimental area and the corresponding water-land segmentation result are shown in Figure 4 .

[0083] Step four: retaining the intensity data of airborne LiDAR bathymetry of land area The full waveform data and intensity data corresponding to the water area need to be further classified and removed;

[0084] S104, sea surface echo identification and corresponding intensity data removal by constructing waveform decomposition model to identify only one valid echo in water area

[0085] The process S104 specifically includes:

[0086] Step one: the effective water area full waveform data i = 1, 2, …, N', N' is the total number of effective water area full waveform data, through the construction of sea surface-water-seabed-noise layered decomposition model to the corresponding water area full waveform data Fitting, sea surface-water-seabed-noise layered decomposition model is defined as the superposition of sea surface echo, water scattering component, seabed echo and environmental noise:

[0087]

[0088] Wherein, y(t) is the sum of all components, sea surface and seabed echo component is defined as similar to the expression of Gaussian function, for the full waveform data from the sea surface directly reflected back to the sensor, M is 1, only sea surface echo component, refraction occurs at the sea surface, and then penetrate the water to the seabed reflected back to the sensor full waveform data M is generally 2, of course, when encountering phytoplankton transmission M is 3, A j , μ j , ω j and γ j are the amplitude, peak position, half wave width and waveform shape parameters of the jth Gaussian component, the water scattering component is defined as a triangle, which is represented by ∫(ax+b)+∫(cx+d), ax+b is the expression function of the triangle rising waist, the start and end time are timel and time2 respectively, cx+d is the expression function of the triangle falling waist, the start and end time are time2 and time3 respectively, the sum of the function integrals of the two in the duration is the water scattering component, is the noise component, which is the average noise level The integral in the whole recording time;

[0089] Step two: define the objective function min|f(t)-y(t)| 2 And calculate the parameters in the model, specifically:

[0090] First step: set the allowable error ε as the threshold condition for stopping iteration, set the maximum iteration number Max N As the iteration termination condition;

[0091] Second step: consider y(t) = F(t, P) as the function of undetermined coefficients Wherein A j , μ j , ω j , γ j are the parameters of the jth echo component, m is the total number of parameters to be estimated, initialize the parameter vector P 0 = {p1 0 , p20 , p3 0 , …, P m 0}, and F (t, P) is expanded in Taylor series, and for the sake of calculation simplicity, the second and higher order terms are omitted, to obtain:

[0092]

[0093] The third step is to estimate the parameter vector to be estimated by driving the LM algorithm, specifically including:

[0094] (1) The damping coefficient ζ is initialized, and the Jacobian matrix J(t) is calculated, specifically:

[0095] Assume the function ∏:

[0096]

[0097] ζ≥0 is the damping coefficient, in order to make the first-order partial derivative of the function ∏ with respect to P be 0, that is, to minimize ∏, the Jacobian matrix J(t) is introduced:

[0098]

[0099] (2) According to the Jacobian matrix, calculate the Hessian matrix f = J(t) T J(t) ;

[0100] (3) Calculate the residual ξ of the original waveform and the fitting waveform corresponding to the Nth iteration N :

[0101] ξ N = |f(t)-F(t, P) N | 2

[0102] (4) Calculate the step size ΔP:

[0103] ΔP = -[Γ+ζD]J(t)ξ N

[0104] Where D is the unit matrix.

[0105] The fourth step is to update the parameters according to ΔP if ξ N > ε, and the iteration number N < Max N , re-execute step two, otherwise stop iteration and output the final parameter vector P;

[0106] Step three: the size of the number of Gaussian echo components M obtained after executing step two is used to classify the waveform, if M is equal to 1, the corresponding full waveform data is the water surface reflection echo waveform The corresponding intensity data is Otherwise, the corresponding full waveform data is water transmission echo waveform The corresponding intensity data is Figure 5 The full waveform classification result image of the experimental area is shown, including invalid waveform, land waveform, water surface reflection waveform and water transmission waveform;

[0107] Step four: remove the intensity data corresponding to the water surface reflection echo waveform identified, and the remaining intensity data is retained;

[0108] S105, the remaining intensity data is used as the intensity data of the airborne laser radar depth sounding system which can effectively express the attribute of the object

[0109] The process S105 specifically includes:

[0110] The intensity data corresponding to the land waveform The water transmission echo waveform The corresponding intensity data is Merging constitutes the intensity data set of the airborne laser radar depth sounding system which can effectively express the attribute of the object, and the corresponding Figure 6 The original intensity data image recorded by the airborne laser radar depth sounding system in the experimental area and the effective intensity data image after removing the intensity corresponding to the invalid waveform and the water surface reflection waveform are shown.

[0111] Corresponding to the above method, the embodiment also provides an airborne laser radar depth sounding system intensity data extraction system, comprising:

[0112] The data analysis module is used for analyzing the airborne laser radar depth sounding Las file to obtain terrain data, intensity data and full waveform data;

[0113] The waveform identification module is used for feature extraction and threshold determination on the full waveform data to identify invalid waveform and valid waveform, and remove the intensity data corresponding to the invalid waveform;

[0114] The water-land segmentation module is used for registration according to the terrain data and multispectral image data, and water-land segmentation based on water body index to obtain water-land segmentation result; The water-land segmentation result includes water area range and land area range;

[0115] The waveform classification module is used for classifying the valid waveform according to the water-land segmentation result to obtain water area waveform and land area waveform, and based on the constructed waveform decomposition model, according to the water area waveform, waveform decomposition and Gaussian echo component number judgment are carried out to obtain water surface reflection waveform and water transmission waveform, and the intensity data corresponding to the water surface reflection echo waveform is removed;

[0116] The data extraction module is configured to combine the intensity data corresponding to the remaining land range and the effective intensity data corresponding to the water body transmission waveform to obtain an extraction result.

[0117] The various embodiments are described in a progressive manner in the specification, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other. For the system disclosed in the embodiments, the description is relatively simple because it corresponds to the method disclosed in the embodiments. The relevant parts can be referred to the description of the method.

[0118] The principles and implementation manners of the present application are described by using specific examples in the specification. The above description of the embodiments is only used to help understand the method of the present application and its core idea. Meanwhile, for those skilled in the art, the specific implementation manners and application ranges can be changed according to the idea of the present application. In conclusion, the content of the specification should not be understood as a limitation of the present application.

Claims

1. An airborne lidar bathymetry system intensity data extraction method, characterized by, The method comprises the following steps: parsing a laser radar sounding (Las) file to obtain terrain data, intensity data, and full waveform data; performing feature extraction and threshold determination on the full waveform data to identify invalid waveforms and valid waveforms, and eliminating the intensity data corresponding to the invalid waveforms; registering the terrain data and multispectral image data, and performing water-land segmentation on the multispectral image data based on a water body index to obtain a water-land segmentation result; the water-land segmentation result comprises a water area range and a land area range; classifying the valid waveforms according to the water-land segmentation result to obtain water area waveforms and land area waveforms, performing waveform decomposition and Gaussian echo component quantity judgment on the water area waveforms based on a constructed waveform decomposition model, obtaining water surface reflection waveforms and water body transmission waveforms, and eliminating the intensity data corresponding to the water surface reflection waveforms; merging the intensity data corresponding to the remaining land area range and the intensity data corresponding to the water body transmission waveforms to obtain an extraction result.

2. The airborne lidar depth system intensity data extraction method of claim 1, wherein, The method comprises the following steps: searching for the positions of the maximum and minimum values of the intensity of the full waveform data to obtain the minimum intensity and maximum intensity corresponding time and time deviation; statistically analyzing the minimum intensity and maximum intensity corresponding time and time deviation of the full waveform data disturbed by noise to determine an optimal threshold value; setting a judgment condition according to the optimal threshold value, identifying invalid waveforms and valid waveforms disturbed by noise according to the judgment condition, and eliminating the intensity data corresponding to the invalid waveforms.

3. The airborne lidar depth system intensity data extraction method of claim 1, wherein, The method comprises the following steps: accurately registering the multispectral image data and the terrain data based on a visualization method to obtain geometrically corrected multispectral image data; performing water body index calculation and image data water-land segmentation on the geometrically corrected multispectral image data to obtain the water-land segmentation result of the terrain data.

4. The airborne lidar depth system intensity data extraction method of claim 1, wherein, The method comprises the following steps: constructing a sea surface-water body-sea bottom-noise layered decomposition model; solving the parameters of the sea surface-water body-sea bottom-noise layered decomposition model based on an LM algorithm, and sequentially fitting the full waveform data of the water area range to obtain a fitting result; classifying the waveforms into the water surface reflection waveforms and the water body transmission waveforms according to the fitting result, and eliminating the intensity data corresponding to the water surface reflection waveforms.

5. The airborne lidar depth system intensity data extraction method of claim 1, wherein, The method comprises the following steps: The search function is programmed to calculate the minimum intensity and maximum intensity corresponding time and time deviation of each full waveform data in turn, wherein f i (t) is the i-th waveform data, i is the waveform data identifier, i = 1, 2, …, N, N is the total number of all waveform data, t corresponds to time, the waveform data record start time is t1, and the end time is t n , is the intensity data recorded by the system at t1, is the maximum intensity corresponding acquisition time, is the minimum intensity corresponding acquisition time, and Δtime is the minimum intensity and maximum intensity acquisition time deviation; The invalid waveform data meeting the condition is filtered by setting a judgment condition, wherein if the waveform f i (t) is located at the first 4% of the whole recording time, i.e. The minimum intensity is located at the last 4% of the whole recording time, i.e. The time deviation of the minimum intensity from the maximum intensity is not less than 90% of the whole recording time, i.e. Δtime > t n *90%, the waveform f i (t) is considered as an invalid waveform f i (t) 0 , otherwise, it is determined as a valid waveform f i (t) 1 ; According to the determined category, the invalid waveform and the intensity data corresponding to the invalid waveform are removed, and the step of "writing a search function to calculate the minimum intensity and the maximum intensity corresponding time and time deviation of each full waveform data in turn" is returned until all the intensity data corresponding to the invalid waveform in the intensity data is removed.

6. The airborne lidar depth system intensity data extraction method of claim 1, wherein, According to the terrain data and the multispectral image data, registration is performed, and water and land segmentation is performed on the multispectral image data based on a water body index to obtain a water and land segmentation result, including: Through the registration processing, the terrain data corresponding to the effective waveform is registered with the multispectral image data collected in the same period to obtain a registration result; Calculate water index NDWI , Index NDWI =(R Green -R NIR ) / (R Green +R NIR ); where R Green is the reflectance recorded in the green band, R NIR is the reflectivity recorded in the near-infrared band. If Index NDWI If it is greater than 0.15, it is a water area, otherwise it is a land area; Based on the registration result, a category label of terrain data corresponding to a water body region of multispectral image data is determined as water, and corresponding full waveform data is and a category label of terrain data corresponding to a land region is determined as land, and corresponding full waveform data is full waveform data corresponding to topographic data of a terrestrial region corresponding intensity data, corresponding intensity data corresponding full waveform data is classified and rejected.

7. The airborne lidar bathymetry system intensity data extraction method of claim 1, wherein, Based on the constructed waveform decomposition model, waveform decomposition and Gaussian echo component quantity judgment are performed according to the water area waveform to obtain a water surface reflection waveform and a water body transmission waveform, and intensity data corresponding to the water surface reflection waveform is removed, including: A sea surface-water body-seabed-noise layered decomposition model is constructed, and based on the sea surface-water body-seabed-noise layered decomposition model, full waveform data corresponding to a water area is decomposed fitting is performed; and an expression of the sea surface-water body-seabed-noise layered decomposition model is: Wherein, y(t) is the superposition of sea surface echo, water scattering component, sea bottom echo and environmental noise, M is 1 for the full waveform data directly reflected from the sea surface to the sensor, only the sea surface echo component, refraction occurs at the sea surface, and then penetrates the water to the sea bottom to reflect the full waveform data to the sensor, M is 2, when encountering phytoplankton, M is 3, wherein A j , μ j , ω j and γ j are the amplitude, peak position, half-wave width and waveform shape parameters of the jth Gaussian component respectively, the water scattering component is defined as a triangle, represented by ∫(ax+b)+∫(cx+d), ax+b is the expression function of the rising waist of the triangle, and the start and end times are time1 and time2 respectively, cx+d is the expression function of the falling waist of the triangle, and the start and end times are time2 and time3 respectively, and the sum of the function integrals of the two in the duration is the water scattering component, is the noise component, and is the average noise level Integral in the entire recording time; According to a preset target function, parameters in the sea surface-water body-sea bottom-noise layered decomposition model are calculated; According to the size of the fitted parameters M, the waveform category is classified, if M is equal to 1, the corresponding full waveform data is a water surface reflection waveform, otherwise the corresponding full waveform data is a water body transmission echo waveform; Intensity data corresponding to the water surface reflection waveform identified is removed, and the remaining intensity data is retained.

8. The airborne lidar depth system intensity data extraction method of claim 7, wherein, The expression of the target function is: min||f(t)-y(t)||; wherein f(t) is the original waveform data.

9. The airborne lidar depth system intensity data extraction method of claim 7, wherein, The algorithm used to solve the target function is the LM algorithm.

10. An airborne lidar bathymetry system intensity data extraction system characterized by, Including: The data analysis module is configured to analyze the airborne laser radar sounding Las file to obtain terrain data, intensity data, and full waveform data; The waveform recognition module is configured to perform feature extraction and threshold determination on the full waveform data to identify invalid waveforms and effective waveforms, and remove intensity data corresponding to the invalid waveforms; The water and land segmentation module is configured to register the terrain data and the multispectral image data, and perform water and land segmentation on the multispectral image data based on a water body index to obtain a water and land segmentation result; the water and land segmentation result includes a water area range and a land area range; The waveform classification module is configured to classify the effective waveforms according to the water and land segmentation result to obtain a water area waveform and a land area waveform, and based on a constructed waveform decomposition model, perform waveform decomposition and Gaussian echo component quantity judgment on the water area waveform to obtain a water surface reflection waveform and a water body transmission waveform, and remove intensity data corresponding to the water surface reflection waveform; The data extraction module is configured to merge the intensity data corresponding to the remaining land area range and the effective intensity data corresponding to the water body transmission waveform to obtain an extraction result.

Images