Wave spectrum analysis method and system based on laser radar

Abstract

The invention belongs to the related technical field of ocean monitoring, and discloses a wave spectrum analysis method and system based on a laser radar. The method comprises the following steps: obtaining a plurality of groups of concentric circles in obtained wave field point cloud data; fitting discrete time-domain signals of the wave height and the derivative of the wave height corresponding to the concentric circles along with time change, and calculating a discrete value of the amplitude of the wave spectrum along with frequency change; calculating unknown parameters in the directional spread function; the discrete value of the amplitude of the wave spectrum along with the frequency change is multiplied by the spread function to obtain the discrete value of the direction spectrum in the frequency and direction, namely the analysis of the wave spectrum is realized. According to the invention, the problem of realizing wave field spectrum and direction spectrum analysis by using point cloud data in wave field analysis is solved.

Description

Technical Field

[0001] The present invention belongs to the technical field related to ocean monitoring, and more specifically, relates to a method and system for wave spectrum analysis based on lidar. Background Art

[0002] In many key fields such as ocean engineering, weather forecasting, and environmental science, obtaining the characteristic information of the wave field plays a crucial role. The core wave field characteristic elements include significant wave height, wave height spectrum, spectral peak period, direction moments (including the main wave propagation direction and direction divergence), and direction spectrum.

[0003] Commonly used spectrum estimation methods include the Walsh method, multi-taper method, and maximum entropy method. These different spectrum analysis methods are each adapted to specific application scenarios. Different spectrum analysis methods are adapted to different application scenarios. Kang et al. effectively extracted target features by using Walsh spectrum analysis combined with dynamic programming for feature selection when processing complex acoustic signals. Vahidreza et al. explored a new method for enhancing the directional spatial wave spectrum and obtained wave parameters such as significant wave height, peak period, and wavelength. Jeyaseelan and Balaji compared the spectral features obtained by the Fast Fourier Transform (FFT) with those obtained by the multi-taper method. The research shows that the wave spectrum obtained by the multi-taper method is superior to the FFT. Lygre and Krogstad introduced a general maximum entropy estimation method for estimating the directional distribution of ocean wave spectra.

[0004] The study of directional spectra is of great significance in ocean engineering and wave analysis, and can provide information on the distribution of wave energy in terms of frequency and direction. Commonly used methods for estimating directional spectra include the Maximum Entropy Method (MEM), the Extended Maximum Entropy Method (EMEP), the Bayesian Method (BDM), the Maximum Likelihood Method (MLM), and the Phase-Time Path Difference Method (PTPD), etc. Each method has its own advantages and limitations. The BDM method proposed by Hashimoto, N., & Kobune, K estimates the directional spectrum through the Akaike Bayesian Information Criterion (ABIC), assuming that the direction distribution function is a piecewise constant function. It can handle arbitrarily mixed instrument array measurements and is applicable to various wave measurement data. It performs well in processing a large amount of measurement data, can provide high-resolution directional spectrum estimates, and has high robustness to noise. Capon, J solves the directional spectrum through maximum likelihood estimation, assuming that the direction distribution function is a function of a specific form and is applicable to various wave measurement data. This method performs well in processing various wave measurement data, but is sensitive to noise. Draycott, S., Davey, T, etc. estimate the wave direction by measuring the phase difference of waves, which is applicable to the analysis of a single-direction wave field. The PTPD method performs well in processing a single-direction wave field, but has high requirements for the uniformity of the wave field. The EMEP proposed by Hashimoto, N., Nagai, T., & Asai, T is an extension of the MEM, considering the error of the cross-spectrum, and estimates the directional spectrum by minimizing the error. It can handle arbitrarily mixed instrument array measurements.

[0005] In terms of the means of obtaining wave field information, currently, it mainly includes wave buoys and satellite remote sensing technology. Raghukumar et al. conducted a comprehensive evaluation of a new type of buoy, covering multiple aspects such as directional spectra, directional moments, and characteristic parameters (such as significant wave height, peak period, and mean wave direction). The buoy shows high reliability in the measurement of significant wave height and wave height spectra. However, its application has certain limitations, mainly reflected in the need for regular maintenance, high repair and replacement costs, and being fixed in a specific area, which greatly limits its application effect in large-scale sea area monitoring. Ozesmi and Bauer summarized the relevant literature on satellite remote sensing in wetland resource monitoring and management, emphasizing the significant advantages of satellite remote sensing in resource monitoring. Satellite remote sensing can cover a large area of the sea area and is suitable for global-scale wave information monitoring, especially performing well in the open ocean and deep sea areas. However, the spatial resolution of satellite measurements is relatively low, making it difficult to capture the fine wave characteristics of local sea areas, and it is easily affected by atmospheric conditions such as clouds and precipitation. Summary of the Invention

[0006] In view of the above deficiencies or improvement requirements of the prior art, the present invention provides a method and system for wave spectrum analysis based on lidar, which solves the problem of being unable to perform wave spectrum analysis using point cloud data.

[0007] To achieve the above object, according to one aspect of the present invention, a method for wave spectrum analysis based on lidar is provided. The method includes the following steps:

[0008] Select a plurality of fixed points from the acquired wave field point cloud data, and draw a set of concentric circles with each fixed point as the center to obtain multiple sets of concentric circles;

[0009] Fit the point cloud in each circle at each moment to obtain the wave height and wave guide at the center of the circle, so as to obtain the wave height and wave guide values at the center of the circle fitted by the point cloud in each circle at all moments, that is, obtain the discrete values of the wave height and wave guide at the center of the circle fitted by each circle in each set of concentric circles changing with time;

[0010] Optionally, compare the discrete values of the wave height and wave guide changing with time at the center of the circle fitted by each circle in a set of concentric circles with a preset standard value, and take the radius of the circle closest to the preset standard value as the optimal radius; use the discrete values of the wave height and wave guide changing with time at the center of the circle fitted with the optimal radius to calculate the discrete values of the wave spectrum amplitude changing with frequency;

[0011] Select at least three circles with the optimal radius from multiple sets of concentric circles, and use the discrete values of the wave height and wave derivative at the center of the circle with the optimal radius changing with time to calculate the unknown parameters in the direction diffusion function, so as to determine the diffusion function;

[0012] Multiply the discrete values of the wave spectrum amplitude changing with frequency by the diffusion function to obtain the discrete values of the direction spectrum in frequency and direction, that is, realize the analysis of the wave spectrum.

[0013] Further preferably, the discrete time domain signals of the wave height and the derivative of the wave height changing with time are obtained by fitting the point cloud data with a two-dimensional parabola.

[0014] Further preferably, the formula of the two-dimensional parabola is as follows:

[0015]

[0016] Where η xx (t), η yy (t) and η yx (t) are the time series of the second derivative of the wave height at the center of the circle with respect to x, the second derivative with respect to y, and the partial derivative with respect to y and x alternately, respectively. η(t), η y (t) and η x (t) are the time series of the wave height at the center of the circle and the first derivative of the wave height with respect to y and x, respectively. η i(t, x i , y i ) is a point cloud data set within a circular region, t is time, x i is the x-direction coordinate of the i-th point cloud, y i is the y-direction coordinate of the i-th point cloud.

[0017] Further preferably, the discrete values of the amplitude of the wave spectrum changing with frequency are obtained by performing wave spectrum analysis on the discrete values of the wave height and the derivative of the wave height changing with time using the multi-taper method.

[0018] Further preferably, the formula of the direction spreading function is as follows:

[0019]

[0020] where a n (f), b n (f), (n = 1,..., N) are unknown parameters, n is the order number, N is the total order, f is the wave frequency, and θ is the wave direction angle.

[0021] Further preferably, the unknown parameters in the direction spreading function are obtained by analysis using the extended maximum entropy method.

[0022] Further preferably, in the method, the main wave propagation direction (θ 1 , θ 2 ) and the direction divergence are calculated according to the following formula using the discrete values of the wave height and the wave guide changing with time

[0023]

[0024] where a 1 (f), a 2 (f), b 1 (f), and b 2 (f) are all Fourier coefficients related to the wave direction moment, f is the wave frequency, M 1 (f) and M 2 (f) are parameters related to the wave composed of Fourier coefficients.

[0025] Further preferably, when drawing a set of concentric circles with each fixed point as the center, for the moments containing incomplete data and outliers, the Kriging interpolation algorithm is used to fill the data missing area, and the interpolation weight function is:

[0026]

[0027] where d j is the distance from the known point to the point to be interpolated, d j≤10 m, where ρ is a power parameter, m is the distance in meters, and j is the number of the known point.

[0028] According to another aspect of the present invention, there is provided a lidar-based wave spectrum analysis system, which includes an actuator that executes a lidar-based wave spectrum analysis method as described above.

[0029] According to yet another aspect of the present invention, there is provided a computer-readable storage medium, on which a computer program is stored. When the computer program is executed by an actuator, it implements a lidar-based wave spectrum analysis method as described above.

[0030] Generally speaking, compared with the prior art by the above technical solutions conceived by the present invention, the following beneficial effects are achieved:

[0031] 1. The present invention realizes wave spectrum analysis by using the changing point cloud on the sea surface over a period of time, selecting fixed points in the point cloud data to draw concentric circles, and using the data in the concentric circles to fit the wave height and wave guide at the center of the circle, that is, at the fixed point, thereby realizing wave spectrum analysis using point cloud data.

[0032] 2. The present invention uses point cloud data for wave spectrum analysis. Compared with the existing method that only uses a limited number of wave buoys to collect points at the buoys for wave spectrum analysis, the data obtained by the present invention is more, the fitting accuracy is higher, and thus the accuracy of wave spectrum analysis is higher.

[0033] 3. The present invention uses time series of three or more wave heights and their derivatives to solve the unknown parameters in the diffusion function. Since the directional spectrum estimation needs to construct equations through the time series of wave heights at multiple fixed points, three or more fixed points can provide enough independent equation numbers M to solve the unknown parameters in the directional diffusion function.

[0034] 4. The present invention uses the Extended Maximum Entropy Method (EMEP) for directional spectrum estimation. By maximizing the entropy function and combining the cross-power spectrum constraint of the observed data, the directional spectrum of ocean waves is effectively estimated. Its implementation steps include initializing parameters, iteratively calculating the directional diffusion function, selecting the optimal model order, and numerical calculation.

[0035] 5. The present invention performs interpolation processing at the center of the circle for the measured lidar point cloud data in each frame to obtain discrete time-domain signals of the wave height and the derivative of the wave height changing with time. The interpolation processing uses a two-dimensional parabolic fitting method to ensure the integrity and continuity of the data, so as to obtain the time series of the wave height and the derivative at the center of the circle.

[0036] 6. Compared with technologies such as wave buoys, the present invention has lower deployment and maintenance costs, can achieve real-time monitoring of large areas of sea areas, is less affected by weather conditions such as clouds and precipitation, and has relatively high spatial resolution, which can meet the needs of fine wave feature monitoring, and is applicable to wave feature monitoring under various complex sea conditions, providing strong technical support for fields such as ocean engineering, weather forecasting, and environmental science. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 is a flowchart of a lidar-based wave spectrum analysis method and system constructed according to a preferred embodiment of the present invention;

[0038] Figure 2 is a schematic diagram of lidar collecting wave spectrum point cloud data constructed according to a preferred embodiment of the present invention;

[0039] Figure 3 is a 3D point cloud view of the collected lidar constructed according to a preferred embodiment of the present invention;

[0040] Figure 4 is a discrete time-domain signal of wave height and its derivative with respect to time constructed according to a preferred embodiment of the present invention, where (a) is the discrete time-domain signal of wave height with respect to time, and (b) is the discrete time-domain signal of the derivative of wave height in the x direction with respect to time;

[0041] Figure 5 is a flowchart of analyzing the directional spectrum by the extended maximum entropy method constructed according to a preferred embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0042] In order to make the objectives, technical solutions and advantages of the present invention clearer and more understandable, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0043] A lidar-based wave spectrum analysis method, the method comprising the following steps:

[0044] S1 A drone or other facility carries a lidar to hover at a set height H and dynamically scan the wave field of the sea surface, obtain high-precision point cloud data, and set the vertical angular resolution, horizontal angular scanning frequency, and sampling radius;

[0045] As Figure 2 shown, the lidar needs to hover at the set height H and dynamically scan the wave field to obtain high-precision point cloud data of the wave spectrum over a period of time. If a drone is used, the spatial position of the drone is fixed, and it is necessary to maintain an almost constant hovering position.

[0046] Since some of the point clouds will be absorbed by water when the actual lidar scans the water surface, to eliminate the influence of too little point cloud data on the subsequent process, for the moments with incomplete data and abnormal points, the Kriging interpolation algorithm is used to fill the missing data area, and the interpolation weight function is:

[0047]

[0048] Among them, d j represents the distance from the known point to the point to be interpolated, and ρ is a power parameter used to adjust the influence degree of the distance on the weight, usually taking a value of 2 or greater.

[0049] As Figure 3 shown, the point cloud data is divided into circular regions with multiple radii. Specifically, an arbitrary fixed spatial point in the point cloud data is selected as the center, and multiple concentric circles are drawn with a preset step size within a set radius range; different fixed spatial points are selected, and concentric circles are drawn with the same radius as above to obtain multiple groups of concentric circles;

[0050] S2 Fitting to obtain discrete time-domain signals of the wave height and its derivative with respect to time at different fixed points.

[0051] (1) Extract wave field parameters through fitting and spectral analysis. The specific steps are to fit the point cloud data within the circular sampling regions with different radii through two-dimensional parabolic fitting to obtain the time series of the wave height and the derivative of the wave height;

[0052] As Figure 5 shown, for the wave height time series fitting, for the point cloud data within each sampling region, the two-dimensional parabolic fitting method is used to fit the instantaneous wave surface elevation. η(t), η x (t) and η y (t) are the time series of the wave height at the center of the circle and its first-order derivatives with respect to x and y obtained by fitting, respectively.

[0053] Two-dimensional parabolic surface fitting:

[0054]

[0055] Among them, the parameters η xx (t), η yy (t) and η yx (t) are the time series of the second-order derivative of the wave height at the center of the circle with respect to x, the second-order derivative with respect to y, and the alternating partial derivative with respect to y and x, respectively. η(t), η y (t) and η x (t) are the time series of the wave height at the center of the circle and the first-order derivatives of the wave height with respect to y and x, respectively. η i (t,x i ,y i) is a point cloud data set within a circular area.

[0056] (2) Use the multi-taper method for spectral analysis, which are discrete values of the amplitude of the wave field spectrum varying with wave frequency; spectral analysis, perform spectral calculation on the wave height time series using the multi-taper method.

[0057] (3) Calculate the main wave propagation direction and direction divergence;

[0058] A. Direction moment calculation, based on the wave height spectrum and the spectra of the wave height derivative η x (t) and η y (t), calculate the first four Fourier coefficients:

[0059]

[0060] B. The main wave propagation direction can refer to two definitions, using the coefficients (a 1 , b 1 ) and (a 2 , b 2 ) respectively. The corresponding direction divergences are defined as follows:

[0061]

[0062] Where: M 1 and M 2 are shown as follows respectively,

[0063]

[0064] Where, is the main wave propagation direction obtained through energy-weighted Fourier coefficients (such as ), including the energy-weighted main wave propagation directions in the swell area and the sea wave area and Where, the energy-weighted main wave propagation direction in the swell area is shown as the following formula:

[0065]

[0066] The energy-weighted main wave propagation direction in the sea wave area is shown as the following formula:

[0067]

[0068] In the formula, the energy-weighted Fourier coefficient in the swell area is obtained from the following formula:

[0069]

[0070] The energy-weighted Fourier coefficient in the swell area is obtained from the following formula:

[0071]

[0072] The same applies to other energy-weighted Fourier coefficients.

[0073] S3 calculates the unknown parameters in the diffusion function

[0074] Select the optimal radius. Select a set of discrete values of wave height and its derivative with respect to time obtained at the center of a concentric circle of any fixed point, and compare them with the preset discrete values of wave height and its derivative with respect to time. The radius with the closest data is the optimal radius;

[0075] Input the relationships of wave height and wave height derivative with respect to time corresponding to the point clouds within more than 3 optimal radius circles into the multi-window method for spectral analysis to obtain the unknown parameters a n and b n .

[0076] Set the initial approximate solution and to be zero.

[0077] As the order increases, perform iterative calculations. Use the following formula to calculate the directional diffusion function:

[0078]

[0079] where a n (f), b n (f), (n = 1,..., N) are unknown parameters.

[0080] When observing data, the errors contained in the cross-power spectrum must be considered. Therefore, the existence of the error ε i can be considered and modified as follows:

[0081]

[0082] where M is the number of independent equations remaining after removing meaningless co-spectra, zero orthogonal spectra, etc. The independent equations where a n , b n are located are non-linear and difficult to solve, which is why local linearization and iterative techniques are introduced to solve this problem. If the approximate solution is known, the solution can be written as:

[0083]

[0084] a′ n , b′ n can be regarded as the residuals between the solution and the approximate solution. Substituting ε i and rearranging gives the linearized equation:

[0085]

[0086] The relationship between the cross-power spectrum and the directional spectrum satisfies:

[0087]

[0088] After iterative calculation, the directional spreading function G(θ|f) is updated, and at the same time, the optimal model order N is selected by introducing the Akaike Information Criterion (AIC):

[0089]

[0090] where M is the number of independent equations, is the variance estimate of ε i

[0091] If the AIC increases or the difference is very small after increasing the order, stop increasing the order. At the same time, it should be noted that during the iterative calculation process, in special cases, the calculation sometimes becomes unstable. Therefore, a control parameter δ is introduced to stabilize the calculation, especially for the instability phenomenon that may occur at high orders. When the calculation is unstable, reduce the value of δ, for example, δ = 0.5 k , and recalculate. The EMEP method effectively estimates the directional spectrum of ocean waves by maximizing the entropy function and combining the cross-power spectrum constraints of the observed data.

[0092] S4 uses the Extended Maximum Entropy Method (EMEP) to analyze and obtain the directional spectrum, and the directional spectrum function: S(f,θ) = S(f)G(θ|f)

[0093] The present invention will be further described below in conjunction with specific embodiments.

[0094] In this embodiment, a multi-line airborne lidar is used for data sampling. The lidar hovers at a height H, the vertical angular resolution of the radar is 1°, the horizontal angular resolution is 0.5°, and the sampling frequency is 10 Hz. Thus, N T frames of radar point cloud images are obtained. The initial value of the radius range is 0.4 ≤ R ≤ 2.0 m, and the radius interval is 0.2 m. Multiple fixed points are selected from the obtained wave field point cloud data, and a set of concentric circles is drawn with each fixed point as the center to obtain multiple sets of concentric circles.

[0095] For the radar point cloud images used in this embodiment, there will inevitably be moments when the radar point cloud data is too small. The following formula is set to measure the proportion of the number of poor data moments in the total number of moments:

[0096]

[0097] where N bad is the number of point clouds less than the cumulative number of point clouds N c ​The number of moments. For moments containing incomplete data and outliers, the data at that moment is obtained by Kriging interpolation using the data from the two adjacent frames before and after. The algorithm fills the data missing area, and the interpolation weight function is:

[0098]

[0099] where d j represents the distance from the known point to the point to be interpolated, and ρ is a power parameter used to adjust the influence degree of the distance on the weight, usually taking a value of 2 or greater.

[0100] In this embodiment, the N T frame radar point cloud data is interpolated at the center of the circle for each frame to obtain the time series of η, η x and η y changing with time t. η(t), η x (t) and η y (t) are the time series of the wave height at the center of the circle and its first-order derivatives with respect to x and y respectively, which are obtained by fitting.

[0101] Commonly used methods for interpolation processing include plane fitting, two-dimensional parabola fitting, and cubic fitting, etc. In this embodiment, the two-dimensional parabola fitting method is used to process the data, and the least squares method is used for solution.

[0102] The two-dimensional parabola fitting formula is as follows:

[0103]

[0104] where the parameters η xx (t), η yy (t) and η yx (t) are the time series of the second-order derivative of the wave height at the center of the circle with respect to x, the second-order derivative with respect to y, and the alternating partial derivative with respect to y and x respectively. η(t), η y (t) and η x (t) are the time series of the wave height at the center of the circle and the first-order derivatives of the wave height with respect to y and x. η i (t,x i ,y i ) is the point cloud data set within the circular region.

[0105] Figure 4 Taking the 40-second time series as an example, the time series images of the wave height η and the wave height derivative η x at R = 0.8m and R = 2.0m obtained by the two-dimensional parabola fitting method are plotted. It can be seen that the wave height time series images at the two radii are basically coincident with the reference value images, and only small burrs appear at the wave peaks and wave troughs.

[0106] To quantitatively measure the effectiveness of the fitting method at different radii, the significant wave height H s (R) and the mean square wave steepness are statistically analyzed and compared with the preset reference values. The significant wave height H s (R) can be defined from η(t) as follows:

[0107]

[0108] The mean square wave steepness is the average of the sum of the squares of η x (t) and η y (t), used to measure the steepness of the wave surface, and the formula is as follows:

[0109]

[0110] Wave spectrum feature extraction: Based on the point cloud data within the radii R = 0.8m and R = 2.0m, the time series η(t) of the wave height at the center (0,0) is obtained through two-dimensional parabolic fitting, and the wave height spectrum S η (f) is obtained through the multi-taper method.

[0111] For the convenience of subsequent analysis of the wave spectrum and direction moment in different frequency ranges, in this method, the frequency range 0.04 ≤ f ≤ 0.1Hz is called the swell band, then the frequency range 0.1 ≤ f ≤ 0.4Hz is called the sea band, and the frequency range 0.4 ≤ f ≤ 1Hz is called the chop band.

[0112] Substitute the wave height spectrum and the spectra of the wave height derivatives η x (t) and η y (t) obtained by the multi-taper method into the formula to obtain the first four Fourier coefficients:

[0113]

[0114]

[0115] In the formula, E(f,θ) is the frequency-direction spectrum. is the imaginary part of the cross-spectrum between the time series η(t) of the signal wave height and the time series η x (t) of the partial derivative of the signal wave height with respect to x. The cross-spectrum is represented by the real part and the imaginary part, as shown in the following formula:

[0116]

[0117] Similarly, and are respectively the imaginary part and the real part of the cross-spectrum of the corresponding two signals,

[0118]

[0119] Solving the direction moments (such as the main wave propagation direction and direction divergence) that actually reflect the wave nature: Two definitions can be considered for the main wave propagation direction, using the coefficients (a 1 , b 1 ) and (a 2 , b 2 ) respectively. The corresponding direction divergences are defined as follows:

[0120]

[0121] where M 1 and M 2 are shown as follows respectively:

[0122]

[0123] In the above formula, is the main wave propagation direction obtained through energy-weighted Fourier coefficients (such as ), including the energy-weighted main wave propagation directions in the swell area and the sea wave area and where the energy-weighted main wave propagation direction in the swell area is shown as follows:

[0124]

[0125] The energy-weighted main wave propagation direction in the sea wave area is shown as follows:

[0126]

[0127] In the formula, the energy-weighted Fourier coefficient in the swell area is obtained from the following formula:

[0128]

[0129] The energy-weighted Fourier coefficient in the swell area is obtained from the following formula:

[0130]

[0131] Other energy-weighted Fourier coefficients can be obtained in the same way. Based on the above formulas, the main wave propagation direction and direction divergence can be plotted, and based on this, the analyzability and rationality of airborne lidar in evaluating the wave directional characteristics can be proved.

[0132] The flow chart of obtaining the direction spectrum by the extended maximum entropy method is as Figure 5 shown, and its specific principle is as follows:

[0133] The directional spectrum S(f,θ) describes the distribution of wave energy over frequency f and propagation direction θ. The directional spectrum can be expressed as the product of the frequency spectrum S(f) and the directional spreading function G(θ|f):

[0134] S(f,θ) = S(f)G(θ|f)

[0135] The directional spreading function is calculated using the following formula:

[0136]

[0137] The directional spectrum of ocean waves is estimated by solving through the Extended Maximum Entropy Method (EMEP), calculating from low order (N = 1) to high order. Set the initial approximate solution and to zero.

[0138] As the order increases, iterative calculations are performed. The directional spreading function usually takes values greater than or equal to zero. However, in EMEP, the function is considered to always take positive values, so the directional spreading function can be extended to:

[0139]

[0140] where a n (f), b n (f), (n = 1,..., N) are unknown parameters.

[0141] When observing data, the errors contained in the cross-power spectrum must be considered. Therefore, the existence of the error ε i can be considered and it is modified as follows:

[0142]

[0143] where M is the number of independent equations remaining after removing meaningless co-spectra, zero orthogonal spectra, etc. The independent equations where a n , b n are located are non-linear and difficult to solve, which is why local linearization and iterative techniques are introduced to solve this problem.

[0144] If the approximate solution is known, the solution can be written as:

[0145]

[0146] a′ n , b′ n can be regarded as the residuals between the solution and the approximate solution. Substituting ε i and rearranging gives the linearized equation:

[0147]

[0148] By assuming an appropriate approximate solution to solve iteratively, where in the above formula:

[0149]

[0150]

[0151] The constraint conditions to be satisfied include the relationship between the cross-power spectrum and the directional spectrum:

[0152]

[0153] After iterative calculation, the directional spreading function G(θ|f) is updated, and at the same time, the optimal model order N is selected by introducing the Akaike Information Criterion (AIC):

[0154]

[0155] where M is the number of independent equations, is the variance estimate of ε i .

[0156] If the AIC increases or the difference is very small after increasing the order, stop increasing the order. At the same time, it should be noted that during the iterative calculation process, in special cases, the calculation sometimes becomes unstable. Therefore, a control parameter δ is introduced to stabilize the calculation, especially the unstable phenomenon that may occur at high orders.

[0157]

[0158] When the calculation is unstable, reduce the value of δ, for example, δ = 0.5 k , and recalculate. The EMEP method effectively estimates the directional spectrum S(f,θ) of ocean waves by maximizing the entropy function and combining the cross-power spectrum constraints of the observed data.

[0159] In the field of ocean monitoring technology, airborne lidar, as an advanced measurement means, has important practical significance for wave measurement and feature analysis. This method constructs a set of wave feature analysis processes based on airborne lidar synthetic data, obtains the wave spectrum and wave direction moments, analyzes the fitting effect of the two-dimensional parabola fitting method on the point cloud data in different radius circular sampling areas, and the influence of various spectrum analysis methods on wave features.

[0160] The analysis of the results provides an important theoretical basis for the application of airborne lidar in wave measurement, further proves its high efficiency and accuracy in ocean monitoring, and provides strong technical support for ocean engineering and ocean environment research.

[0161] In the radar point cloud data within the fitting radius, the wave spectrum, main wave propagation direction, direction divergence and other sea wave characteristics obtained by the multiple window method have a good comparison with the preset reference values. The extended maximum entropy method is used to effectively estimate the directional spectrum of ocean waves. The results show the feasibility of the method for extracting wave spectrum and direction moment features based on airborne lidar. Using airborne lidar data is one of the effective tools for observing wave characteristics.

[0162] It is easy for those skilled in the art to understand that the above are only preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention should be included within the protection scope of the present invention.

Claims

1. A wave spectrum analysis method based on laser radar, characterized in that: The method comprises the following steps: Select multiple fixed points in the acquired wave field point cloud data, draw a group of concentric circles with each fixed point as the center, so as to obtain multiple groups of concentric circles; The wave height and waveguide at the center of the circle are obtained by fitting the point cloud in each circle at each time, so as to obtain the wave height and waveguide values ​​at the center of the circle fitted by the point cloud in each circle at all times, that is, to obtain the discrete values ​​of the wave height and waveguide at the center of the circle fitted by each circle in each group of concentric circles over time; The discrete values ​​of the wave height and waveguide at the center of each circle of a group of concentric circles that are fitted over time are compared with the preset standard value, and the radius of the circle closest to the preset standard value is taken as the optimal radius; the discrete values ​​of the wave height and waveguide at the center of the circle of the optimal radius that are fitted over time are used to calculate the discrete values ​​of the wave spectrum amplitude that changes with frequency; At least three circles with the optimal radius as the radius are selected from the multiple groups of concentric circles, and the unknown parameters in the directional diffusion function are calculated using the discrete values ​​of the wave height and wave derivative at the center of the circle with the optimal radius that vary with time, so as to determine the diffusion function; The discrete values ​​of the directional spectrum in frequency and direction are obtained by multiplying the discrete values ​​of the amplitude of the wave spectrum with the frequency with the diffusion function, thereby realizing the analysis of the wave spectrum.

2. A laser radar based wave spectrum analysis method as claimed in claim 1, characterized in that: The discrete time domain signals of the wave height and the derivative of the wave height varying with time are obtained by fitting the point cloud data with a two-dimensional parabola.

3. A laser radar based wave spectrum analysis method as claimed in claim 2, characterized in that: The formula for the two-dimensional parabola is as follows: Among them, η xx (t), η yy (t) and η yx (t) are the time series of the wave height at the center of the circle for the second-order derivative of x, the second-order derivative of y, and the partial derivatives of y and x alternately, η(t), η y (t) and η x (t) are the time series of the wave height at the center of the circle and the first-order derivatives of the wave height with respect to y and x, respectively, η i (t,x i ,y i ) is the point cloud dataset in the circular area, t is the time, x i is the x-coordinate of the ith point cloud, y i is the y-coordinate of the i-th point cloud.

4. A laser radar-based wave spectrum analysis method as claimed in claim 1 or 3, characterized in that: The discrete values ​​of the amplitude of the wave spectrum that vary with frequency are obtained by performing spectral analysis on the discrete values ​​of the wave height and the derivative of the wave height that vary with time using a multiple window method.

5. A laser radar based wave spectrum analysis method as claimed in claim 1 or 3, characterized in that: The formula of the directional diffusion function is as follows: Among them, a n (f), b n (f), (n=1,...,N) are unknown parameters, n is the order number, N is the total order, f is the wave frequency, and θ is the wave direction angle.

6. A laser radar-based wave spectrum analysis method as claimed in claim 1 or 3, characterized in that: The unknown parameters in the directional diffusion function are obtained by using the extended maximum entropy method.

7. A laser radar-based wave spectrum analysis method as claimed in claim 1 or 3, characterized in that: The method uses the discrete values ​​of wave height and waveguide varying with time to calculate the main wave propagation direction (θ1, θ2) and directional divergence according to the following formula: Where a1(f), a2(f), b1(f) and b2(f) are the Fourier coefficients related to the wave directional moment, f is the wave frequency, and M1(f) and M2(f) are wave-related parameters composed of Fourier coefficients.

8. A laser radar based wave spectrum analysis method as claimed in claim 1 or 3, characterized in that: When a set of concentric circles are drawn with each fixed point as the center, for moments with incomplete data and abnormal points, the Kriging interpolation algorithm is used to fill in the data missing area, and the interpolation weight function is: Among them, d j is the distance from the known point to the point to be interpolated, d j ≤10m, ρ is a power parameter, m is the distance unit meter, and j is the number of the known point.

9. A wave spectrum analysis system based on laser radar, characterized in that: The system includes an actuator, which executes a laser radar-based wave spectrum analysis method as described in any one of claims 1-8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by an executor, the laser radar-based wave spectrum analysis method described in any one of claims 1 to 8 is implemented.

Images