1. A wave spectrum calculation method based on inertial measurement is characterized by comprising the following steps:

acquiring N groups of horizontal system accelerations within preset time, wherein N is more than or equal to 2;

carrying out vector decomposition on the acceleration of the X axis and the acceleration of the Y axis of the horizontal system to obtain acceleration components of the geographical coordinate system in the east direction and the north direction, which are respectively defined as an east acceleration A _ E and a north acceleration A _ N;

performing frequency domain numerical integration on the Z-axis acceleration of the horizontal system for 2 times to obtain a horizontal system vertical displacement sequence D _ Z, and calculating an inbred spectrum of the horizontal system vertical displacement sequence D _ Z to obtain a sequence C₁₁；

Calculating the orthogonal spectrum of the horizontal system vertical displacement sequence D _ Z and the eastern acceleration A _ E to obtain a sequence Q₁₂；

Calculating the orthogonal spectrum of the horizontal system vertical displacement sequence D _ Z and the north direction acceleration A _ N to obtain a sequence Q₁₃；

Calculating the self-cross spectrum of the acceleration A _ E in the east direction to obtain a sequence C₂₂；

Calculating the self-cross spectrum of the acceleration A _ N in the north direction to obtain a sequence C₃₃；

Calculating a covariance spectrum of the east acceleration A _ E and the north acceleration A _ N to obtain a sequence C23;

according to sequence C₁₁、Q₁₂、Q₁₃、C₂₂、C₃₃And C₂₃Calculating a direction spectrum S (f, theta)The first M coefficients of Fourier decomposition (M) are more than or equal to 2, the M coefficients are synthesized into the direction spectrum S (f, theta), and the direction spectrum S (f, theta) is smoothed.

2. A wave spectrum calculation method based on inertial measurement according to claim 1, characterized in that fast Fourier transform is performed on N sets of horizontal system vertical displacement sequences D _ Z to obtain N sets of real parts and imaginary parts, then the square sum of the real parts and the imaginary parts is obtained and divided by the number N of sampling points to obtain a rough power spectrum of the horizontal system vertical displacement sequences D _ Z, and the rough power spectrum is smoothed to obtain a smoothed power spectrum S (f).

3. An inertial measurement based wave spectrum calculation method according to claim 2 wherein the spectral peak frequency f of the power spectrum s (f) of the horizontal series of vertical displacements D _ Z is calculated_pTo obtain the wave peak value period T_pThe bandwidth of the adaptive adjustment filter is 1/(T)_p+5)<f_p<1/(T_p-5) only the wave components within K seconds before and after the dominant wave period are retained, K is more than or equal to 5 and less than or equal to 15.

4. An inertial measurement based wave spectrum calculation method as claimed in claim 3, wherein the formula is a discrete integralPerforming reserved signal frequency band integration on the power spectrum S (f) to obtain wave zero order moment m₀Calculating the spectral characteristic effective wave height of the power spectrum S (f)

5. An inertial measurement based wave spectrum calculation method as claimed in claim 4, wherein the formula is a discrete integralPerforming reserved frequency band integration on the power spectrum S (f) to obtain a first moment m of the wave₁Calculating the average wave period T of the spectral characteristics of the power spectrum S (f)_a＝m₀/m₁Wherein, S (f)_n) The power spectrum s (f) corresponding to the nth frequency is shown as a function value.

6. An inertial measurement based wave spectrum calculation method as claimed in claim 5, wherein θ is given by the formula_1n＝arctan(B₁/A₁) Obtaining the time sequence of the average direction of the wave, and calculating the spectrum peak value frequency f of the time sequence of the average direction of the wave_pCorresponding angle D_mean＝θ_1n(f_p) And D is_meanAs the average direction of the wave corresponding to the peak frequency of the wave, wherein B₁And A₁Is the first order fourier coefficient of the directional spectrum.

7. An inertial measurement based wave spectrum calculation method as claimed in claim 6, wherein θ is given by the formula_2n＝0.5*arctan(B₂/A₂) Obtaining a wave dominant wave direction time sequence, and calculating the spectrum peak frequency f of the wave dominant wave direction time sequence_pCorresponding angle D_m＝θ_2n(f_p) And D is_mAs the principal direction of the wave corresponding to the peak frequency of the wave, where B₂And A₂Is the second order fourier coefficient of the directional spectrum.

Background

The sea wave is researched, the wave is generally regarded as a superposition of a plurality of groups of random sine waves, the frequency of each group of sine waves is different, but the average value of the heights of the groups of sine waves is zero, and the average value of the heights of the wave surfaces obtained by superposition is zero, namely the waves do not change along with time. Generally, the wave changes slowly in the height variance within a short observation time interval (10-30 minutes), and the change of the sea wave during this period can be regarded as a quasi-stationary random process. Therefore, describing sea waves in a random process using a spectral method is a very important research approach. The wave spectrum includes a power spectrum S (f) and a direction spectrum S (f, θ). The ocean wave power spectrum s (f) is a density division of the reaction energy, giving the distribution of energy provided by the different frequency component waves with respect to frequency. Generally, in the wave component, energy is mainly provided by the component waves concentrated in a narrow frequency band, and the energy provided by the component waves with small or large frequencies at two ends of the power spectrum is small. The wave spectrum characteristics are reflected to the macroscopic objective expression of the ocean: the sea wave height in a certain sea area is different, but the wave height of small frequency (long period gravity wave) and high frequency (short period capillary wave, broken clutter and the like) is very small, and some large waves with obvious periodicity occupy the main position. The direction spectrum S (f, theta) is a spectrum reflecting the internal directional structure of the sea wave, and can provide the distribution of the energy of each component wave in different directions relative to the frequency. For a given wave frequency, the directional spectrum may give the energy in different directional intervals, i.e. the distribution of the constituent wave energy of a given frequency with respect to the wave direction.

Wave measurement systems (inertial wave sensors, pressure wave sensors and the like) are mainly installed in buoys at home and abroad to obtain wave characteristic parameters, and a direction spectrum, an average wave direction and a main wave direction (called as a PRZ method) are obtained through cross spectrum calculation among 3 time sequences of roll angle (roll), pitch angle (pitch) and vertical direction displacement (D _ Z). However, when the size of the float is large, the float is easy to shake, and the inclination angle is small or even no inclination occurs, so that errors occur in calculation of the direction spectrum.

Disclosure of Invention

Aiming at the problems, the invention provides a wave spectrum calculation method based on inertial measurement, which mainly solves the problem that the direction spectrum calculation has errors due to small shaking inclination angle or no inclination of a large-size buoy.

In order to solve the technical problems, the technical scheme of the invention is as follows:

a wave spectrum calculation method based on inertial measurement comprises the following steps:

acquiring N groups of horizontal system accelerations within preset time, wherein N is more than or equal to 2; carrying out vector decomposition on the acceleration of the X axis and the acceleration of the Y axis of the horizontal system to obtain acceleration components of the geographical coordinate system in the east direction and the north direction, which are respectively defined as an east acceleration A _ E and a north acceleration A _ N; performing frequency domain numerical integration on the Z-axis acceleration of the horizontal system for 2 times to obtain a horizontal system vertical displacement sequence D _ Z, and calculating an inbred spectrum of the horizontal system vertical displacement sequence D _ Z to obtain a sequence C₁₁(ii) a Calculating the orthogonal spectrum of the horizontal system vertical displacement sequence D _ Z and the eastern acceleration A _ E to obtain a sequence Q₁₂(ii) a Calculating the orthogonal spectrum of the horizontal system vertical displacement sequence D _ Z and the north direction acceleration A _ N to obtain a sequence Q₁₃(ii) a Calculating the self-cross spectrum of the acceleration A _ E in the east direction to obtain a sequence C₂₂(ii) a Calculating the self-cross spectrum of the acceleration A _ N in the north direction to obtain a sequence C₃₃(ii) a Calculating the covariance spectrum of the east acceleration A _ E and the north acceleration A _ N to obtain a sequence C₂₃(ii) a According to sequence C₁₁、Q₁₂、Q₁₃、C₂₂、C₃₃And C₂₃Calculating the first M coefficients of Fourier decomposition of the direction spectrum S (f, theta), wherein M is more than or equal to 2, synthesizing the M coefficients into the direction spectrum S (f, theta), and smoothing the direction spectrum S (f, theta).

In some embodiments, the N sets of horizontal system vertical displacement sequences D _ Z are subjected to fast fourier transform to obtain N sets of real and imaginary parts, and then the square sum of the real and imaginary parts is obtained and divided by the number N of sampling points to obtain a rough power spectrum of the horizontal system vertical displacement sequences D _ Z, and the rough power spectrum is smoothed to obtain a smoothed power spectrum s (f).

In some embodiments, calculating a spectral peak frequency f of the power spectrum s (f) of the horizontal series of vertical shift sequences D _ Z_pTo obtain the wave peak value period T_pThe bandwidth of the adaptive adjustment filter is 1/(T)_p+5)<f_p<1/(T_p-5) only the wave components within K seconds before and after the dominant wave period are retained, K is more than or equal to 5 and less than or equal to 15.

In some embodiments, formulated as discrete integralsPerforming reserved signal frequency band integration on the power spectrum S (f) to obtain wave zero order moment m₀Calculating the spectral characteristic effective wave height of the power spectrum S (f)

In some embodiments, formulated as discrete integralsPerforming reserved frequency band integration on the power spectrum S (f) to obtain a first moment m of the wave₁Calculating the average wave period T of the spectral characteristics of the power spectrum S (f)_a＝m₀/m₁Wherein, s (fn) represents the power spectrum s (f) corresponding to the nth frequency.

In some embodiments, according to the formula θ_1n＝arctan(B₁/A₁) Obtaining the time sequence of the average direction of the wave, and calculating the spectrum peak value frequency f of the time sequence of the average direction of the wave_pCorresponding angle D_mean＝θ_1n(f_p) And D is_meanAs the average direction of the wave corresponding to the peak frequency of the wave.

In some embodiments, according to the formula θ_2n＝0.5*arctan(B₂/A₂) Obtaining a wave dominant wave direction time sequence, and calculating the spectrum peak frequency f of the wave dominant wave direction time sequence_pCorresponding angle D_m＝θ_2n(f_p) And D is_mAs the main wave direction corresponding to the peak wave frequency.

The invention has the beneficial effects that: the method adopts the cross spectrum among 3 time sequences of the eastern acceleration, the northern acceleration and the horizontal series vertical displacement sequence to obtain the direction spectrum, the average wave direction and the main wave direction, and compared with the traditional method, the method can solve the problem of direction spectrum calculation error caused by the condition that the buoy shakes and inclines less or even has no inclination due to the larger size shape of the buoy.

Drawings

Fig. 1 is a schematic flow chart of a method for calculating a wave spectrum based on inertial measurement according to an embodiment of the present invention.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention clearer and clearer, the following detailed description of the present invention is provided with reference to the accompanying drawings and detailed description. It is to be understood that the specific embodiments described herein are merely illustrative of the invention and are not limiting of the invention. It should be further noted that, for the convenience of description, only some but not all of the relevant aspects of the present invention are shown in the drawings.

As shown in fig. 1, the present embodiment provides a wave spectrum calculation method based on inertial measurement, which can be used for on-board calculation to improve the performance of a wave measurement system, and mainly includes the following steps:

101, acquiring N groups of horizontal system accelerations within preset time, wherein N is more than or equal to 2.

In this embodiment, the acquisition of the carrier system 3-axis acceleration, the 3-axis angular acceleration and the 3-axis magnetometer is performed according to a set sampling frequency, and the real-time attitude calculation is performed by synchronously performing a complementary fusion filtering algorithm based on 9-axis parameters to acquire 3-axis attitude angles (pitch, roll, yaw), thereby acquiring the horizontal system 3-axis acceleration (a _ X, A _ Y, A _ Z). The sampling frequency of 2/4/8Hz and the sampling time length of 4/8/16 minutes which are set in advance are selectable, and when the sampling time is up, the wave calculation program is entered.

102, performing vector decomposition on the acceleration of the X axis (A _ X) and the acceleration of the Y axis (A _ Y) of the horizontal system by combining the heading angle yaw to obtain acceleration components of the geographical coordinate system in the east direction and the north direction, which are respectively defined as an east acceleration A _ E and a north acceleration A _ N;

in the embodiment, FFT conversion is carried out on the east acceleration A _ E, and the real part/imaginary part is respectively stored in a preset storage space; performing FFT (fast Fourier transform) on the north direction acceleration A _ N, and respectively storing a real part/an imaginary part in a preset storage space;

calculating a covariance spectrum C between variables_xySelf-cross-bred spectrum Q_xyObtaining corresponding sequences, the calculation formulas are respectively

C_xy＝(Re[x]*Re[y]-Im[x]*Im[y])/N，Q_xy＝(Im[x]*Re[y]+Re[x]*Im[y]) The specific steps are as follows in step 103-108, without the following sequence:

103, performing frequency domain numerical integration on the Z-axis acceleration of the horizontal system for 2 times to obtain a horizontal system vertical displacement sequence D _ Z, and calculating an inbred spectrum of the horizontal system vertical displacement sequence D _ Z to obtain a sequence C₁₁；

104, calculating the orthogonal spectrum of the horizontal system vertical displacement sequence D _ Z and the east acceleration A _ E to obtain a sequence Q₁₂；

105, calculating the orthogonal spectrum of the horizontal system vertical displacement sequence D _ Z and the north direction acceleration A _ N to obtain a sequence Q₁₃；

106, calculating the self-cross spectrum of the east acceleration A _ E to obtain a sequence C₂₂；

107, calculating the self-cross spectrum of the acceleration A _ N in the north direction to obtain a sequence C₃₃；

108, calculating a covariance spectrum of the acceleration A _ E in the east direction and the acceleration A _ N in the north direction to obtain a sequence C₂₃；

109 according to sequence C₁₁、Q₁₂、Q₁₃、C₂₂、C₃₃And C₂₃Calculating the first M coefficients of Fourier decomposition of the direction spectrum S (f, theta), wherein M is more than or equal to 2, synthesizing the M coefficients into the direction spectrum S (f, theta), and smoothing the direction spectrum S (f, theta).

In this embodiment, M is 5, so the first 5 coefficients are a0, a1, a2, B1, and B2, and the calculation formulas are: a. the₀＝C₁₁；A₁＝Q₁₂；B₁＝Q₁₃；A₂＝C₂₂+C₃₃；B₂＝2C₂₃And the coefficient k is the number of waves obtained by each measurement, and can be obtained by performing wave time domain zero crossing statistics in the early stage of a program. According to a Fourier series decomposition formula: s (f, theta) ═ A₀+A₁ cosθ+B₁sinθ+A₂ cos2θ+B₂sin2 theta are respectively substituted into 5 Fourier coefficients, and the wave direction spectrum S (f, theta) is obtained through calculation. Because only the first 5 terms are used to approximate the fourier series instead of the directional spectral function, a certain error is generated on the real result. Therefore, the direction spectrum needs to be smoothed, and after a plurality of experiments, it is found that the best effect is obtained when the weight coefficients are 2/5, 3/4 and 1/5, and the smoothing formula is as follows:

the on-board computation of the direction spectrum S (f, θ) is a two-dimensional variable function. If the calculation is substituted according to 4096 sets of frequency components (4Hz, with 4/4096 frequency intervals) and 360 sets of direction components (the 360-degree range is divided according to 1-degree intervals), 5760K of memory is required, which is far beyond 256K of the STM32F 427. Moreover, the calculation takes a long time, and the real-time calculation requirement cannot be met. By combining the physical characteristics of sea waves, the resolution of a computational grid is properly reduced in the onboard computation process, frequency components f are dispersed into 128 groups (4Hz, the frequency interval is 4/128), direction components theta are dispersed into 72 groups (the 360-degree range is divided according to 5-degree intervals), 9216 groups of S (f, theta) values are finally obtained, and only 36K memory is needed.

The method adopts the cross spectrum among 3 time sequences of the eastern acceleration, the northern acceleration and the horizontal series vertical displacement sequence to obtain the direction spectrum, the average wave direction and the main wave direction, and compared with the traditional method, the method can solve the problem of direction spectrum calculation error caused by the condition that the buoy shakes and inclines less or even has no inclination due to the larger size shape of the buoy. Specifically, the principle is as follows: the buoy which carries out direction spectrum calculation through two inclination angles (rolling and pitching) and vertical direction displacement has a good effect on the small-sized wave carrier buoy, and because the wave following performance of the small-sized buoy is obvious, the following wave can obviously shake at two angles, so that the direction can be well identified. However, in some large buoys, the stability is good, and the obvious shaking on two axes does not necessarily occur along with the waves, so the method for providing the shaking angle on the large buoy has poor effect on calculating the direction spectrum (the large waves can cause a small inclination angle of the large buoy, and the large buoy is not sensitive to the change of the inclination angle). However, no matter on the large/small buoy, the sensitivity of two accelerations in the horizontal direction is higher than the angle, and when a small wave impacts, obvious acceleration changes can be measured in the east and north directions of the large buoy, so that the wave spectrum is calculated through the spectrums of the accelerations in the east/north directions and the displacement in the vertical direction, and the platform adaptability is better.

Further, the wave characteristic parameter T is inverted through a power spectrum S (f) and a direction spectrum S (f, theta)_p,T_a，H_m0,D_mean,D_mThe following are:

201, performing fast fourier transform on N sets of horizontal system vertical displacement sequences D _ Z to obtain N sets of real parts and imaginary parts, then solving the square sum of the real parts and the imaginary parts and dividing by the number N of sampling points to obtain a rough power spectrum of the horizontal system vertical displacement sequences D _ Z, and smoothing the rough power spectrum to obtain a smoothed power spectrum s (f).

202, calculating the spectrum peak frequency f of the power spectrum S (f) of the horizontal vertical shift sequence D _ Z_pTo obtain the wave peak value period T_pThe power spectrum S (f) adaptively adjusts the bandwidth of the filter to be 1/(T)_p+5)<f_p<1/(T_p-5) only the wave components within K seconds before and after the dominant wave period are retained, K is more than or equal to 5 and less than or equal to 15.

In this example, K is 10.

In this embodiment, the rough power spectrum obtained in step 201 is subjected to 3-point moving average, and after a plurality of testsIt is found that the weighting coefficients are 1/4,1/2 and 1/4, and the formula is as follows: sf [ n ]]＝Sf[n-1]/4+Sf[n]/2+Sf[n+1](ii)/4; further, the power spectrum Sf [ n ]]Frequency f corresponding to spectral peak_pThe peak frequency is the wave spectrum peak period T_p＝1/f_p。

For Sf [ n ] in reserved frequency band]Performing zero order and first order integration to obtain wave zero order moment m of wave power spectrum₀First moment m of sum wave₁Inverting the effective wave height H by the first and second moments_m0And the average wave period Ta, as shown in steps 203 and 204, respectively:

203 according to the formula of discrete integrationPerforming reserved signal frequency band integration on the power spectrum S (f) to obtain wave zero order moment m₀The spectral characteristic effective wave height of the power spectrum S (f) is obtained

204 according to the formula of discrete integralPerforming reserved frequency band integration on the power spectrum S (f) to obtain a first moment m of the wave₁The average wave period T of the spectral characteristics of the power spectrum S (f) is obtained_a＝m₀/m₁Wherein, S (f)_n) The power spectrum s (f) corresponding to the nth frequency is shown as a function value.

205 according to the formula theta_1n＝arctan(B₁/A₁) Obtaining the time sequence of the average direction of the wave, and calculating the spectrum peak value frequency f of the time sequence of the average direction of the wave_pCorresponding angle D_mean＝θ_1n(f_p) And D is_meanAs the average direction of the wave corresponding to the peak frequency of the wave.

206 according to the formula θ_2n＝0.5*arctan(B₂/A₂) Obtaining a wave main wave direction time sequence, and calculating the spectrum peak value frequency f of the wave main wave direction time sequence_pCorresponding angleD_m＝θ_2n(f_p) And D is_mAs the main wave direction corresponding to the peak wave frequency.

The above embodiments are only for illustrating the technical concept and features of the present invention, and the purpose thereof is to enable those skilled in the art to understand the contents of the present invention and implement the present invention accordingly, and not to limit the protection scope of the present invention accordingly. All equivalent changes or modifications made in accordance with the spirit of the present disclosure are intended to be covered by the scope of the present disclosure.