A lidar-based wave spectrum analysis method and system

Abstract

The application belongs to the technical field of marine 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 concentric circles in the obtained wave field point cloud data; matching the discrete time domain signals of the wave height corresponding to the concentric circles and the derivative of the wave height with respect to time, and calculating the discrete values of the amplitude of the wave spectrum with respect to frequency; calculating unknown parameters in a directional spread function; and multiplying the discrete values of the amplitude of the wave spectrum with respect to frequency and the spread function to obtain the discrete values of the directional spectrum in frequency and direction, i.e., realizing the analysis of the wave spectrum. Through the application, the problem of realizing the analysis of the wave field spectrum and the directional spectrum by using point cloud data in the analysis of the wave field is solved.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field related to ocean monitoring, and more particularly, to a wave spectrum analysis method and system based on laser radar. BACKGROUND

[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, directional moment (including the main wave propagation direction and directional dispersion) and directional spectrum.

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

[0004] Directional spectrum research is of great significance in ocean engineering and wave analysis, as it provides information on the distribution of wave energy in frequency and direction. Common methods for estimating directional spectrum include Maximum Entropy Method (MEM), Extended Maximum Entropy Method (EMEP), Bayesian Method (BDM), Maximum Likelihood Method (MLM), and Phase Time Path Difference Method (PTPD). Each method has its advantages and limitations. Hashimoto, N., & Kobune, K proposed the BDM method, which estimates the directional spectrum using the Bayesian Information Criterion (ABIC) and assumes a piecewise constant function for the directional distribution function. This method can handle any mixed instrument array measurements and is suitable for various wave measurement data. It performs well in handling multiple measurement data, providing high-resolution directional spectrum estimates and is highly robust to noise. Capon, J used maximum likelihood estimation to solve the directional spectrum, assuming a specific form of function for the directional distribution function. This method performs well in handling multiple wave measurement data but is sensitive to noise. Draycott, S., Davey, T, et al estimated the direction of the wave by measuring the phase difference of the wave, which is suitable for analyzing single-direction wave fields. The PTPD method performs well in handling single-direction wave fields but requires high uniformity of the wave field. Hashimoto, N., Nagai, T., & Asai, T extended the MEM method to consider the error of the cross-spectrum and estimated the directional spectrum by minimizing the error, which can handle any mixed instrument array measurements.

[0005] Currently, the main methods for obtaining wave field information include wave buoys and satellite remote sensing technology. Raghukumar et al. conducted a comprehensive evaluation of a new type of buoy, involving directional spectrum, directional moment, and characteristic parameters such as significant wave height, peak period, and mean wave direction. The buoy showed high reliability in measuring significant wave height and wave height spectrum. However, its application has certain limitations, mainly reflected in the need for regular maintenance, high maintenance and replacement costs, and the need to be fixed in a specific area, which greatly limits its application effect in large-scale sea 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 monitoring resources. Satellite remote sensing can cover large areas of sea and is suitable for global wave information monitoring, especially in open ocean and deep sea areas. However, the spatial resolution of satellite measurements is relatively low, making it difficult to capture fine wave characteristics in local sea areas, and it is easily affected by atmospheric conditions such as clouds and precipitation. SUMMARY

[0006] In view of the above defects or improvement needs of the prior art, the present application provides a wave spectrum analysis method and system based on a laser radar, which solves the problem that point cloud data cannot be used for wave spectrum analysis.

[0007] To achieve the above object, according to one aspect of the present application, a wave spectrum analysis method based on a laser radar is provided, which comprises the following steps:

[0008] Selecting a plurality of fixed points in the acquired wave field point cloud data, and drawing a set of concentric circles with each fixed point as the center, thereby obtaining a plurality of sets of concentric circles;

[0009] Fitting the point cloud in each circle at each time to obtain the wave height and wave guide at the center of the circle, thereby obtaining the wave height and wave guide values at the center of the circle fitted by the point cloud in each circle at all times, i.e., 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, the discrete values of the wave height and wave guide at the center of the circle fitted by each circle in a set of concentric circles changing with time are compared with preset standard values, 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 wave guide at the center of the circle fitted by the optimal radius are used to calculate the discrete values of the wave spectrum amplitude changing with frequency;

[0011] At least three circles with the optimal radius are selected from the plurality of sets of concentric circles, and the discrete values of the wave height and wave guide at the center of the circle fitted by the optimal radius are used to calculate the unknown parameters in the directional spread function, thereby determining the spread function;

[0012] The discrete values of the wave spectrum amplitude changing with frequency are multiplied by the spread function to obtain the discrete values of the directional spectrum in the frequency and direction, i.e., the analysis of the wave spectrum is realized.

[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] wherein η xx (t), η yy (t) and η yx (t) are time series of the second-order derivative of the wave height at the center with respect to x, the second-order derivative with respect to y and the partial derivative with respect to y and x alternately, respectively, η y (t), η x (t) and η i(t, x i , y i ) is a point cloud dataset in a circular region, t is time, x i is the i-th point cloud x-direction coordinate, y i is the i-th point cloud y-direction coordinate.

[0017] Further preferably, the discrete values of the amplitude of the wave spectrum varying with frequency are obtained by using a multi-window method to analyze the discrete values of the wave height and the derivative of the wave height varying with time.

[0018] Further preferably, the formula of the directional diffusion function is as follows:

[0019]

[0020] wherein 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 directional diffusion function are obtained by using an extended maximum entropy method.

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

[0023]

[0024] wherein a1(f), a2(f), b1(f), and b2(f) are all Fourier coefficients related to the wave direction matrix, f is the wave frequency, and M1(f) and M2(f) are parameters related to the wave composed of Fourier coefficients.

[0025] Further preferably, when a set of concentric circles is drawn with each fixed point as the center, for the time containing incomplete data and abnormal points, the data missing area is filled by using a Kriging interpolation algorithm, and the interpolation weight function is as follows:

[0026]

[0027] wherein d j is the distance from the known point to the point to be interpolated, d j ≤ 10 m, ρ is a power parameter, m is the distance unit meter, and j is the number of known points.

[0028] According to another aspect of the present application, there is provided a laser radar-based wave spectrum analysis system, comprising an executor that executes the above-mentioned laser radar-based wave spectrum analysis method.

[0029] According to still another aspect of the present application, there is provided a computer-readable storage medium having stored thereon a computer program that, when executed by an executor, implements the above-mentioned laser radar-based wave spectrum analysis method.

[0030] Overall, the above technical solutions conceived by the present application have the following beneficial effects compared with the prior art:

[0031] 1. The present application realizes wave spectrum analysis by using the point cloud changing on the sea surface in a period of time and selecting fixed points in the point cloud data to draw concentric circles, and fitting the wave height and wave guide at the center of the concentric circles, i.e., the fixed points, to realize wave spectrum analysis using point cloud data.

[0032] 2. The present application uses point cloud data to perform wave spectrum analysis, which has more data and higher fitting accuracy compared with the prior art that only uses a limited number of wave buoys to collect points at the buoys to perform wave spectrum analysis, and thus the accuracy of wave spectrum analysis is higher.

[0033] 3. The present application uses more than three time series of wave height and its derivative to solve unknown parameters in the diffusion function, because direction spectrum estimation needs to construct an equation set through the time series of wave height of multiple fixed points, and more than three fixed points can provide sufficient number of independent equations M to solve unknown parameters in the direction diffusion function.

[0034] 4. The present application uses the extended maximum entropy method (EMEP) to estimate the direction spectrum, effectively estimates the direction spectrum of ocean waves by maximizing the entropy function and combining the cross-power spectrum constraint of observation data. The implementation steps include initializing parameters, iteratively calculating the direction diffusion function, selecting the optimal model order, and numerical calculation.

[0035] 5. The present application performs interpolation processing at the center of each frame of the measured radar point cloud data to obtain discrete time-domain signals of wave height and its derivative varying with time. The interpolation processing uses two-dimensional parabolic fitting, which ensures the integrity and continuity of the data to obtain the time series of wave height and its derivative at the center.

[0036] 6. Compared with the wave buoy technology, the deployment and maintenance cost is low, the real-time monitoring of large area sea area can be realized, the influence of weather conditions such as cloud layer and precipitation is small, the spatial resolution is relatively high, the demand of fine wave characteristic monitoring can be met, and the wave characteristic monitoring in various complex sea conditions is suitable, which provides strong technical support for the fields of ocean engineering, weather forecast, environmental science and the like. BRIEF DESCRIPTION OF DRAWINGS

[0037] Figure 1 is a flow chart of a laser radar-based wave spectrum analysis method and system constructed according to a preferred embodiment of the present application;

[0038] Figure 2 is a schematic diagram of laser radar collecting wave spectrum point cloud data according to a preferred embodiment of the present application;

[0039] Figure 3 is a 3D point cloud view of the collected laser radar according to a preferred embodiment of the present application;

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

[0041] Figure 5 is a flow chart of an extended maximum entropy method for analyzing directional spectrum according to a preferred embodiment of the present application. DETAILED DESCRIPTION

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

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

[0044] S1 The unmanned aerial vehicle or other facilities carries the laser radar to hover at a set height H and dynamically scan the wave field of the sea surface, to obtain high-precision point cloud data, and set the vertical angle resolution, horizontal angle scanning frequency and sampling radius;

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

[0046] Since the actual laser radar scanning water surface will have part of the point cloud absorbed by water, in order to eliminate the influence of too little point cloud data on the subsequent process, for the time containing incomplete data and abnormal points, the Kriging interpolation algorithm is used to fill the data missing area, and the interpolation weight function is:

[0047]

[0048] Where d j represents the distance from the known point to the point to be interpolated, and p is a power parameter for adjusting the influence degree of distance on weight, usually taking a value of 2 or more.

[0049] As shown in Figure 3 , the point cloud data is divided into multiple radius circular regions, specifically, any fixed spatial point in the point cloud data is selected as the center of the circle, and multiple concentric circles are drawn according to the preset step size; different fixed spatial points are selected, and concentric circles are drawn according to the same radius as above, to obtain multiple groups of concentric circles.

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

[0051] (1) Extract wave field parameters by fitting and spectral analysis, the specific steps are as follows: fit the point cloud data in the circular sampling region of different radii by two-dimensional parabolic fitting, to obtain the time series of wave height and the derivative of wave height;

[0052] As shown in Figure 5 , the wave height time series fitting, for the point cloud data in each sampling region, a 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 the first-order derivatives of x and y, respectively.

[0053] Two-dimensional parabolic fitting:

[0054]

[0055] Where 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 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-order derivatives of y and x, respectively. η i (t, x i ,y i) is a point cloud dataset in a circular region.

[0056] (2) Discrete values of the amplitude of the wave field spectrum as a function of wave frequency using a multi-taper spectral analysis; spectral analysis, the spectral computation of the wave height time series, using a multi-taper method.

[0057] (3) Calculation of the principal wave propagation direction and directional divergence;

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

[0059]

[0060] B. The principal wave propagation direction can be referenced by two definitions, using the coefficients (a1, b1) and (a2, b2) respectively. The corresponding directional divergence is as follows:

[0061]

[0062] where M1 and M2 are as follows,

[0063]

[0064] where, is the principal wave propagation direction obtained by energy-weighted Fourier coefficients (such as ), including the energy-weighted principal wave propagation direction in the swell region and the sea wave region and where the energy-weighted principal wave propagation direction in the swell region is as follows:

[0065]

[0066] The energy-weighted principal wave propagation direction in the sea wave region is as follows:

[0067]

[0068] where the energy-weighted Fourier coefficient in the swell region is obtained by the following formula:

[0069]

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

[0071]

[0072] Other energy-weighted Fourier coefficients can be obtained in the same way.

[0073] S3 calculates unknown parameters in the diffusion function

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

[0075] The corresponding wave height and wave height derivative of the point cloud in the circle with more than three optimal radii are input into the multiple window method for wave spectrum analysis to obtain unknown parameters a in the diffusion function n and b n .

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

[0077] As the order increases, iterative calculation. The directional diffusion function is calculated using the following formula:

[0078]

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

[0080] When observing data, the error contained in the cross power spectrum must be considered. Therefore, the existence of the error ε i can be considered to modify it:

[0081]

[0082] Where M is the number of independent equations left after removing meaningless co-spectrum, zero-orthogonal spectrum, etc. a n ,b n The independent equations in which a n ,b n are nonlinear, 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 residual between the solution and the approximate solution. Substitute ε i and rearrange to get the linearized equation:

[0085]

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

[0087]

[0088] The direction diffusion function G(0|f) is updated after iterative calculation, and the optimal model order N is selected by introducing 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 small after increasing the order, the order increasing is stopped. It is noted that in the iterative calculation process, the calculation is sometimes unstable in special cases, so a control parameter δ is introduced to stabilize the calculation, which may be unstable especially at high order. When the calculation is unstable, the value of δ is reduced, for example, δ = 0.5 k , and the calculation is restarted. The EMEP method effectively estimates the directional spectrum of ocean waves by maximizing the entropy function and combining the cross-power spectrum constraint of the observation data.

[0092] S4 uses the extended maximum entropy method (EMEP) to analyze the directional spectrum, and the directional spectrum function is S(f, θ) = S(f)G(θ|f)

[0093] The application will be further illustrated below with specific examples.

[0094] In this embodiment, a multi-line airborne laser radar is used for data sampling. The laser radar is hovered at a height H, the vertical angle resolution is 1°, the horizontal angle resolution is 0.5°, and the sampling frequency is 10 Hz, thereby obtaining N T frame radar point cloud images. The initial value of the radius range is 0.4≤R≤2.0m, and the radius interval is 0.2m. A plurality of fixed points are selected in the obtained wave field point cloud data, and a group of concentric circles is drawn with each fixed point as the center, thereby obtaining a plurality of groups of concentric circles.

[0095] The radar point cloud image used in this embodiment will inevitably have moments when the radar point cloud data is too small. The following formula is set to measure the proportion of the number of moments when the data is too small to the total number of moments:

[0096]

[0097] where N bad is the number of moments when the point cloud number is less than the cut-off point cloud number N c . For the moments containing incomplete data and abnormal points, the data of the two frames before and after the moment is obtained by Kriging interpolation, 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 p is a power parameter used to adjust the degree of influence of distance on weight, usually taking a value of 2 or more.

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

[0101] Methods commonly used in interpolation processing include plane fitting, two-dimensional parabolic fitting, and cubic fitting, etc. In this embodiment, the two-dimensional parabolic fitting method is used to process the data, and the least squares method is used to solve.

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

[0103]

[0104] where parameters η xx (t), η yy (t), and η yx (t) are time series of the second-order derivatives of the wave height at the center of the circle with respect to x, y, and the partial derivatives of the wave height with respect to y and x alternately, respectively, η(t), η y (t), and η x (t) are 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, and η i (t, x i , y i ) is a point cloud data set in a circular region.

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

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

[0107]

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

[0109]

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

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

[0112] The wave height spectrum and the spectrum of wave height derivative η x (t) and η y (t) obtained by the multiple window method are substituted into the formula to obtain the first four Fourier coefficients:

[0113]

[0114]

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

[0116]

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

[0118]

[0119] The direction moments (e.g. the principal wave propagation direction and the directional spread) of the actual response wave properties are solved: the principal wave propagation direction can be considered in two definitions, using the coefficients (a1, b1) and (a2, b2) respectively. The corresponding directional spread is as follows:

[0120]

[0121] where M1 and M2 are as follows:

[0122]

[0123] In the above formula, is the principal wave propagation direction obtained by energy-weighted Fourier coefficients (e.g. ) including the swell zone and the sea wave zone energy-weighted principal wave propagation directions and where the swell zone energy-weighted principal wave propagation direction is as follows:

[0124]

[0125] The sea wave zone energy-weighted principal wave propagation direction is as follows:

[0126]

[0127] The swell zone energy-weighted Fourier coefficients are obtained by the following formula:

[0128]

[0129] The sea wave zone energy-weighted Fourier coefficients are obtained by the following formula:

[0130]

[0131] Other energy-weighted Fourier coefficients can be obtained in the same way. According to the above formulas, the principal wave propagation direction and the directional spread can be plotted, and the analysis and rationality of the airborne laser radar in evaluating the wave directionality characteristics can be proved based thereon.

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

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

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

[0135] The directional spread function is calculated by the following formula:

[0136]

[0137] The directional spectrum of ocean waves is estimated by the extended maximum entropy method (EMEP) and is calculated from low order (N = 1) to high order. The initial approximate solution is set as and zero.

[0138] With the increase of order, the iterative calculation is performed. The directional spreading function usually takes a value greater than or equal to zero. However, in EMEP, the function is regarded as a function always taking a positive value, so the directional spreading function can be extended as

[0139]

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

[0141] When the observation data is considered, the error contained in the cross power spectrum must be considered. Therefore, the existence of the error ε i can be considered to modify it:

[0142]

[0143] where M is the number of independent equations left after removing meaningless co-spectrum, zero-orthogonal spectrum and other equations. a n ,b n The independent equations in which a n ,b n are non-linear, which is difficult to solve, which is why the local linearization and iteration technique 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 residual error between the solution and the approximate solution. Substituting ε i and rearranging can obtain the linearized equation:

[0147]

[0148] The solution is iteratively solved by assuming an appropriate approximate solution , in which:

[0149]

[0150]

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

[0152]

[0153] The directional spreading function G(θ|f) is updated after the iteration calculation, and the optimal model order N is selected by introducing 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 small after increasing the order, stop increasing the order. At the same time, it is noted that during the iterative calculation process, the calculation sometimes becomes unstable in special cases, so a control parameter δ is introduced to stabilize the calculation, especially the unstable phenomenon that may occur at high order.

[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 laser radar as an advanced measuring means has important practical significance for wave measurement and feature analysis. This method constructs a set of wave feature analysis process based on airborne laser radar synthetic data, through the completion of wave spectrum and wave direction moment, analyzes the fitting effect of two-dimensional parabolic fitting method in different radius circular sampling area, and the influence of various spectral analysis methods on wave features.

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

[0161] In the fitting of radar point cloud data within the radius, the wave spectrum, main wave propagation direction and direction dispersion obtained by the multiple window method have good comparison with the preset reference value, and the extended maximum entropy method effectively estimates the directional spectrum of ocean waves. The results show that the wave spectrum and direction moment feature extraction method based on airborne laser radar is feasible, and the airborne laser radar data is one of the effective tools for observing wave features.

[0162] It is to be understood that the above description is intended to be illustrative and not restrictive. Many other embodiments will be apparent to those of skill in the art upon reading and understanding the above description. The scope of the application should, therefore, be determined with reference to the appended claims, along with the full scope of equivalents to which such claims are entitled.

Claims

1. A laser radar-based wave spectrum analysis method, characterized by, The method comprises the following steps: Selecting a plurality of fixed points in the acquired wave field point cloud data, and drawing a set of concentric circles with each fixed point as the center to obtain a plurality of sets of concentric circles; Fitting the point cloud in each circle at each time to obtain the wave height and wave guide at the center of the circle, and thus obtaining the wave height and wave guide at the center of the circle fitted by each circle in each set of concentric circles, i.e. 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; Optionally, comparing 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 with preset standard values, taking the radius of the circle closest to the preset standard values as the optimal radius; and calculating the discrete values of the wave spectrum amplitude changing with frequency using the discrete values of the wave height and wave guide at the center of the circle fitted by the circle with the optimal radius changing with time; Selecting at least three circles with the optimal radius as the radius in the plurality of sets of concentric circles, and calculating unknown parameters in the directional spread function using the discrete values of the wave height and wave guide at the center of the circle fitted by the circle with the optimal radius changing with time, so as to determine the spread function; Multiplying the discrete values of the wave spectrum amplitude changing with frequency by the spread function to obtain the discrete values of the directional spectrum in frequency and direction, i.e. realizing the analysis of the wave spectrum.

2. A lidar-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 changing with time are fitted using a two-dimensional parabola to obtain the point cloud data.

3. A lidar-based wave spectrum analysis method as claimed in claim 2, characterized in that, The formula of the two-dimensional parabola is as follows: Where, η xx (t), η yy (t) and η yx η(t) represents the time series of the wave height at the center of the circle when the second derivative with respect to x, the second derivative with respect to y, and the partial derivatives with respect to y and x are taken alternately. y (t) and η x (t) represents the time series of wave height at the center of the circle and the first derivatives of the wave height with respect to y and x, respectively. η i (t,x i ,y i (x) represents a point cloud dataset within a circular region, where t is time and x is... i It is the x-coordinate of the i-th point cloud, and the y-coordinate is... i It is the y-coordinate of the i-th point cloud.

4. A lidar-based wave spectrum analysis method according to claim 1 or 3, characterized in that, The discrete values of the wave spectrum amplitude 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 a multiple window method.

5. A lidar-based wave spectrum analysis method according to claim 1 or 3, characterized in that, The formula of the directional spread function is as follows: 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.

6. A lidar-based wave spectrum analysis method according to claim 1 or 3, characterized in that, The unknown parameters in the directional spread function are obtained by analysis using an extended maximum entropy method.

7. A ladar-based wave spectrum analysis method according to claim 1 or 3, wherein, The method utilizes discrete values of wave height and waveguide variation with time to calculate the principal wave propagation direction (θ1, θ2) and directional divergence according to the following equations Wherein a1(f), a2(f), b1(f) and b2(f) are Fourier coefficients related to the wave direction matrix, f is the wave frequency, and M1(f) and M2(f) are parameters related to the wave composed of Fourier coefficients.

8. A ladar-based wave spectrum analysis method according to claim 1 or 3, wherein, When drawing a set of concentric circles with each fixed point as the center, for the time containing incomplete data and abnormal points, the data missing area is filled using a Kriging interpolation algorithm, and the interpolation weight function is as follows: where d j is the distance of the known point to the point to be interpolated, d j ≤ 10 m, p is a power parameter, m is the distance unit meter, and j is the number of the known point.

9. A lidar-based wave spectrum analysis system, characterized by, The system comprises an executor that executes the laser radar-based wave spectrum analysis method according to any one of claims 1-8.

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

Images