Sea wave spectrum calculation method based on inclined attitude

Abstract

The invention discloses a sea wave spectrum calculation method based on a tilt attitude method, and the method comprises the following steps: obtaining N groups of horizontal system displacements, east tilt rates and north tilt rates within a preset time; respectively defining a horizontal system displacement DZ, an east slope rate IE and a north slope rate IN; calculating a co-cross spectrum and a self-cross spectrum between the variables to obtain a corresponding sequence; and calculating the first M coefficients of Fourier decomposition of the directional spectrum according to the sequences C11, Q12, Q13, C22, C33 and C23, M being greater than or equal to 2, synthesizing the M coefficients into the directional spectrum, and smoothing the directional spectrum. According to the method, the direction spectrum, the average wave direction and the main wave direction are obtained through the cross spectrum among the three time sequences of the east inclination rate, the north inclination rate and the horizontal system vertical displacement sequence, and compared with a traditional method, accurate measurement of the wave direction spectrum can be achieved under the condition that obvious displacement cannot be generated due to the fact that an anchor system is straightened.

Description

Technical Field

[0001] The present invention relates to the field of ocean observation applications, and in particular to a method for calculating ocean wave spectra based on tilt attitude. Background Art

[0002] When studying ocean waves, waves are usually considered to be the superposition of multiple groups of random sine waves. Each group of sine waves has a different frequency but a height mean of zero, so the average value of the wave surface height obtained by superposition is zero, that is, it does not change with time. Generally speaking, the height variance of waves changes very slowly within a very short observation time interval (10-30 minutes), and the changes in waves during this period can be regarded as a quasi-stationary random process. Therefore, using spectral methods to describe waves as random processes has become a very important research approach. The wave spectrum includes the power spectrum S(f) and the directional spectrum S(f,θ). The wave power spectrum S(f) is the density division of the reaction energy, which gives the distribution of the energy provided by the waves of different frequencies relative to the frequency. Generally speaking, for the wave components, the energy is mainly provided by the component waves concentrated in a narrow frequency band, while the frequencies at both ends of the power spectrum are very small or very large, and the energy provided is very small. This wave spectrum feature is reflected in the objective manifestation of the ocean with the naked eye: the waves in a certain sea area are uneven in height and length, but the wave heights of small frequencies (long-period gravity waves) and high frequencies (short-period capillary waves, clutter waves, etc.) are very small, and there are some large waves with obvious periodicity that occupy a major position. The directional spectrum S (f, θ) is an energy spectrum that reflects the internal directional structure of the waves, which can give the distribution of the energy of each component wave in different directions relative to the frequency. For a given wave frequency, the directional spectrum can give the energy in different directional intervals, that is, the distribution of the energy of the component waves of a given frequency relative to the wave direction.

[0003] At home and abroad, wave measurement systems (inertial wave sensors, pressure wave sensors, etc.) are mainly installed in buoys to obtain wave characteristic parameters, and the directional spectrum, average wave direction and main wave direction are obtained by cross-spectral calculation between the three time series of eastward displacement, northward displacement and vertical displacement (D_Z). However, when the buoy is continuously pulled by the ocean current and is in tension without obvious displacement movement with the wave, errors will occur in the directional spectrum calculation. Summary of the invention

[0004] In view of the above problems, the present invention proposes a wave spectrum calculation method based on tilt posture, which mainly solves the problem that when the wave buoy is continuously pulled by the ocean current and is in a taut state without obvious displacement movement with the waves, errors may occur in the directional spectrum calculation.

[0005] In order to solve the above technical problems, the technical solution of the present invention is as follows:

[0006] A method for calculating a wave spectrum based on an inclined attitude comprises the following steps:

[0007] Obtain N groups of horizontal system acceleration, pitch angle, roll angle and heading angle within a preset time, N≥2; perform tilt rate attitude solution on the pitch angle, roll angle and heading angle to obtain the east tilt rate I_E and the north tilt rate I_N; perform two frequency domain numerical integrations on the Z-axis acceleration of the horizontal system to obtain the horizontal system vertical displacement sequence D_Z, calculate the self-intersection spectrum of the horizontal system vertical displacement sequence D_Z, and obtain the sequence C 11 ; Calculate the orthogonal spectrum of the horizontal system vertical displacement sequence D_Z and the eastward inclination rate I_E to obtain the sequence Q 12 Calculate the orthogonal spectrum of the horizontal system vertical displacement sequence D_Z and the north direction inclination rate I_N to obtain the sequence Q 13 Calculate the self-intersection spectrum of the eastward inclination I_E to obtain the sequence C 22 ; Calculate the self-intersection spectrum of the north inclination I_N to obtain the sequence C 33 Calculate the cross spectrum of the east inclination I_E and the north inclination I_N to obtain the sequence C 23 According to sequence C 11 , Q 12 , Q 13 , C 22 , C 33 and C 23 The first M coefficients of the Fourier decomposition of the directional spectrum S(f,θ) are calculated, where M≥2, the M coefficients are synthesized into the directional spectrum S(f,θ), and the directional spectrum S(f,θ) is smoothed.

[0008] In some embodiments, a fast Fourier transform is performed on N groups of horizontal vertical displacement sequences D_Z to obtain N groups of real parts and imaginary parts, and then the sum of the squares of the real parts and the imaginary parts is calculated and divided by the number of sampling points N to obtain a rough power spectrum of the horizontal vertical displacement sequence D_Z, and the rough power spectrum is smoothed to obtain a smoothed power spectrum S(f).

[0009] In some embodiments, the spectrum peak frequency f of the power spectrum S(f) of the horizontal system vertical displacement sequence D_Z is calculated. p , get the wave peak period T p The power spectrum S(f) adaptively adjusts the filter bandwidth to 1 / (T p +2) <f p <1 / (T p -2), only retain the component wave components within K seconds before and after the main wave period, 3≤K≤17.

[0010] In some embodiments, according to the discrete integral formula The power spectrum S(f) is integrated in the reserved signal frequency band to obtain the zero-order moment m0 ​​of the wave, and the spectral characteristic effective wave height of the power spectrum S(f) is obtained.

[0011] In some embodiments, according to the discrete integral formula The power spectrum S(f) is integrated in the reserved frequency band to obtain the first-order moment m1 of the wave, and the spectral characteristic average wave period T of the power spectrum S(f) is obtained. a =m0 / m1, where S(fn) represents the power spectrum S(f) function value corresponding to the nth frequency.

[0012] In some embodiments, according to the formula θ 1n =arctan(B1 / A1), get the wave average direction time series, calculate the spectrum peak frequency f of the wave average direction time series p The corresponding angle D mean =θ 1n (f p ), and D mean The average direction of the waves at the peak frequency of the waves.

[0013] In some embodiments, according to the formula θ 2n =0.5*arctan(B2 / A2), get the main wave direction time series of the wave, calculate the spectrum peak frequency f of the main wave direction time series of the wave p The corresponding angle D m =θ 2n (f p ), and D m The main direction of the wave corresponding to the peak frequency of the wave.

[0014] The beneficial effects of the present invention are as follows: the present method adopts the cross-spectrum between three time series, namely the east direction tilt rate, the north direction tilt rate and the horizontal system vertical displacement sequence, to obtain the directional spectrum, the average wave direction and the main wave direction. Compared with the traditional method, the present method can solve the problem that when the wave buoy is continuously pulled by the ocean current and is in a taut state without obvious displacement movement with the wave, errors will occur in the directional spectrum calculation. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 The present invention is a flowchart of a method for calculating a wave spectrum based on tilted attitude disclosed in an embodiment of the present invention. DETAILED DESCRIPTION

[0016] In order to make the purpose, technical solution and advantages of the present invention clearer and more specific, the content of the present invention is further described in detail below in conjunction with the accompanying drawings and specific implementation methods. It is understood that the specific embodiments described herein are only used to explain the present invention, rather than to limit the present invention. It should also be noted that, for the convenience of description, only the parts related to the present invention are shown in the accompanying drawings, rather than all the contents.

[0017] like Figure 1 As shown, this embodiment proposes a wave spectrum calculation method based on tilt attitude, which can be used for onboard calculation to improve the performance of the wave measurement system, and mainly includes the following steps:

[0018] 101, obtain N groups of horizontal system acceleration, east direction tilt rate, and north direction tilt rate within a preset time, where N ≥ 2.

[0019] In this embodiment, the 3-axis acceleration, 3-axis angular acceleration, and 3-axis magnetometer of the carrier system are collected at the set sampling frequency, and the complementary fusion filtering algorithm based on the 9-axis parameters is simultaneously performed to perform real-time attitude solution to obtain the 3-axis attitude angles (pitch, roll, yaw), and the horizontal system acceleration A_Z is obtained, and A_Z is integrated twice in the frequency domain to obtain the horizontal system displacement D_Z.

[0020] The sampling frequency set in advance is 2 / 4 / 8Hz, and the sampling time length is 4 / 8 / 16 minutes. When the sampling time is up, the wave calculation program will be entered.

[0021] 102. The tilt rate is calculated by combining the three-axis attitude angles (pitch, roll, yaw). The east tilt rate I_E and the north tilt rate I_N are calculated according to the following formulas:

[0022] I_E=sin(y)sin(p) / cos(p)-cos(y)sin(r) / cos(p)cos(r)

[0023] I_N=cos(y)sin(p) / cos(p)+sin(y)sin(r) / cos(p)cos(r)

[0024] In this embodiment, the east-direction inclination rate I_E is subjected to FFT transformation, and the real part / imaginary part are respectively stored in a predetermined storage space; the north-direction inclination rate I_N is subjected to FFT transformation, and the real part / imaginary part are respectively stored in a predetermined storage space;

[0025] Calculate the co-concordance spectrum C between each variable xy , self-cross spectrum Q xy , and obtain the corresponding sequence. The calculation formulas are

[0026] C xy=(Re[x]*Re[y]-Im[x]*Im[y]) / N, Q xy =(Im[x]*Re[y]+Re[x]*Im[y]) / N,

[0027] The specific steps are as follows: 103-108, in no particular order:

[0028] 103. Perform two frequency domain numerical integrations on the Z-axis acceleration of the horizontal system to obtain the horizontal system vertical displacement sequence D_Z. Calculate the self-intersection spectrum of the horizontal system vertical displacement sequence D_Z to obtain the sequence C 11 ;

[0029] 104, calculate the orthogonal spectrum of the horizontal system vertical displacement sequence D_Z and the eastward inclination rate I_E, and obtain the sequence Q 12 ;

[0030] 105, calculate the orthogonal spectrum of the horizontal system vertical displacement sequence D_Z and the north direction inclination rate I_N, and obtain the sequence Q 13 ;

[0031] 106, calculate the self-intersection spectrum of the eastward inclination I_E and obtain the sequence C 22 ;

[0032] 107, calculate the self-intersection spectrum of the north inclination I_N and obtain the sequence C 33 ;

[0033] 108, calculate the co-confluence spectrum of the eastward inclination I_E and the northward inclination I_N, and obtain the sequence C 23 ;

[0034] 109, according to sequence C 11 , Q 12 , Q 13 , C 22 , C 33 and C 23 Calculate the first M coefficients of the Fourier decomposition of the directional spectrum S(f,θ), M ≥ 2, synthesize the M coefficients into the directional spectrum S(f,θ), and smooth the directional spectrum S(f,θ).

[0035] In this embodiment, M is 5, so the first five coefficients are A0, A1, A2, B1, B2, and the calculation formulas are:

[0036]

[0037]

[0038] The coefficient k is the number of waves obtained in each measurement, which can be obtained by performing zero-crossing statistics in the wave time domain in the early stage of the program. According to the Fourier series decomposition formula: S(f,θ)=(1 / 2+A1cosθ+B1sinθ+A2cos2θ+B2sin2θ) / π, the five Fourier coefficients are substituted respectively to calculate the wave directional spectrum S(f,θ). After obtaining the original directional spectrum, the directional spectrum needs to be smoothed. After many tests, it is found that the best effect is achieved when the weight coefficient is 1 / 4, 1 / 2, and 1 / 4. The smoothing formula is:

[0039] In the onboard calculation process of the directional spectrum S(f,θ), it is a two-dimensional variable function. If the frequency components are 4096 groups (4Hz, the frequency interval is 4 / 4096) and the directional components are 360 ​​groups (the 360-degree range is divided into 1-degree intervals) for calculation, 5760K of memory is required, which is far beyond the 256K memory of STM32F427. Moreover, the calculation takes a long time and cannot meet the real-time calculation requirements. Combined with the physical characteristics of ocean waves, the calculation grid resolution is appropriately reduced during the onboard calculation process, the frequency component f is discretized into 128 groups (4Hz, the frequency interval is 4 / 128), and the directional component θ is discretized into 72 groups (the 360-degree range is divided into 5-degree intervals), and finally 9216 groups of S(f,θ) values ​​are obtained, which only requires 36K of memory.

[0040] This method uses the cross-spectrum between three time series: the eastward tilt rate, the northward tilt rate, and the horizontal system vertical displacement sequence to obtain the directional spectrum, average wave direction and main wave direction. Compared with the traditional method, it can solve the problem that when the wave buoy is continuously pulled by the ocean current and is in a taut state without obvious displacement movement with the wave, errors in the directional spectrum calculation will occur.

[0041] Furthermore, the wave characteristic parameter T is inverted by the power spectrum S(f) and directional spectrum S(f,θ). p ,T a , H m0 ,D mean ,D m ,as follows:

[0042] 201, perform fast Fourier transform on N groups of horizontal vertical displacement sequences D_Z to obtain N groups of real parts and imaginary parts, then calculate the sum of the squares of the real parts and the imaginary parts and divide it by the number of sampling points N to obtain a rough power spectrum of the horizontal vertical displacement sequence D_Z, and smooth the rough power spectrum to obtain a smoothed power spectrum S(f).

[0043] 202, calculate the spectrum peak frequency f of the power spectrum S(f) of the horizontal system vertical displacement sequence D_Z p , get the wave peak period T p , the power spectrum S(f) adaptively adjusts the filter bandwidth to 1 / (Tp +2) <f p <1 / (T p -2), only retain the component wave components within K seconds before and after the main wave period, 3≤K≤17.

[0044] In this embodiment, K is 10.

[0045] In this embodiment, a three-point sliding average is performed on the rough power spectrum obtained in step 201. After multiple tests, it is found that the best effect is achieved when the weight coefficients are 1 / 4, 1 / 2, and 1 / 4. The formula is: Sf[n] = Sf[n-1] / 4 + Sf[n] / 2 + Sf[n+1] / 4; further, the frequency f corresponding to the peak value of the power spectrum Sf[n] p is the peak frequency, then the wave spectrum peak period T p =1 / f p .

[0046] Perform zero-order and first-order integration on Sf[n] in the reserved frequency band to obtain the zero-order moment m0 ​​and the first-order moment m1 of the wave power spectrum, respectively. The effective wave height H is inverted by the zero-order moment and the first-order moment. m0 and the average wave period Ta, the formulas are shown in steps 203 and 204 respectively:

[0047] 203, according to the discrete integral formula The power spectrum S(f) is integrated in the reserved signal frequency band to obtain the zero-order moment m0 ​​of the wave, and the spectral characteristic effective wave height of the power spectrum S(f) is obtained.

[0048] 204, according to the discrete integral formula The power spectrum S(f) is integrated in the reserved frequency band to obtain the first-order moment m1 of the wave, and the spectral characteristic average wave period T of the power spectrum S(f) is obtained. a =m0 / m1, where S(f n ) represents the power spectrum s(f) function value corresponding to the nth frequency.

[0049] 205, according to the formula θ 1n =arctan(B1 / A1), get the wave average direction time series, calculate the spectrum peak frequency f of the wave average direction time series p The corresponding angle D mean =θ 1n (f p ), and D mean The average direction of the waves at the peak frequency of the waves.

[0050] 206, according to the formula θ 2n=0.5*arctan(B2 / A2), get the main wave direction time series of the wave, calculate the spectrum peak frequency f of the main wave direction time series of the wave p The corresponding angle D m =θ 2n (f p ), and D m The main direction of the wave corresponding to the peak frequency of the wave.

[0051] The above embodiments are only for illustrating the technical concept and features of the present invention, and their purpose is to enable ordinary technicians in the field to understand the content of the present invention and implement it accordingly, and they cannot be used to limit the protection scope of the present invention. Any equivalent changes or modifications made based on the essence of the content of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for calculating sea wave spectrum based on tilt attitude, characterized in that: The following steps are involved: Obtain N groups of horizontal system acceleration, eastward tilt rate, and northward tilt rate within a preset time, N ≥ 2; They are defined as the horizontal system acceleration A_Z, the eastward inclination I_E, and the northward inclination I_N respectively; Perform two frequency domain numerical integrations on the Z-axis acceleration of the horizontal system to obtain the horizontal system vertical displacement sequence D_Z, calculate the self-intersection spectrum of the horizontal system vertical displacement sequence D_Z, and obtain the sequence C 11 ; Calculate the orthogonal spectrum of the horizontal system vertical displacement sequence D_Z and the eastward inclination rate I_E to obtain the sequence Q 12 ; Calculate the orthogonal spectrum of the horizontal system vertical displacement sequence D_Z and the north direction tilted by I_N, and get the sequence Q 13 ; Calculate the self-intersection spectrum of the eastward inclination I_E and obtain the sequence C 22 ; Calculate the self-intersection spectrum of the north inclination I_N and obtain the sequence C 33 ; Calculate the co-confluence spectrum of the eastward inclination I_E and the northward inclination I_N to obtain the sequence C 23 ; According to sequence C 11 , Q 12 , Q 13 , C 22 , C 33 and C 23 The first M coefficients of the Fourier decomposition of the directional spectrum S(f,θ) are calculated, where M≥2, the M coefficients are synthesized into the directional spectrum S(f,θ), and the directional spectrum S(f,θ) is smoothed.

2. The method for calculating the wave spectrum based on the tilted attitude according to claim 1, characterized in that: Perform fast Fourier transform on N groups of horizontal vertical displacement sequences D_Z to obtain N groups of real parts and imaginary parts, then calculate the square sum of the real part and the imaginary part and divide it by the number of sampling points N to obtain the rough power spectrum of the horizontal vertical displacement sequence D_Z, and smooth the rough power spectrum to obtain the smoothed power spectrum S(f).

3. The method for calculating the wave spectrum based on the tilted attitude according to claim 2, characterized in that: Calculate the spectrum peak frequency f of the power spectrum S(f) of the horizontal system vertical displacement sequence D_Z p , get the wave peak period T p The power spectrum S(f) adaptively adjusts the filter bandwidth to 1 / (T p +2) <f p <1 / (T p -2), only retain the component wave components within K seconds before and after the main wave period, 3≤K≤17.

4. The method for calculating the wave spectrum based on the tilted attitude according to claim 3, characterized in that: According to the discrete integral formula The power spectrum S(f) is integrated in the reserved signal frequency band to obtain the zero-order moment m0 ​​of the wave, and the spectral characteristic effective wave height of the power spectrum S(f) is obtained.

5. The method for calculating the wave spectrum based on the tilted attitude according to claim 4, characterized in that: According to the discrete integral formula The power spectrum S(f) is integrated in the reserved frequency band to obtain the first-order moment m1 of the wave, and the spectral characteristic average wave period T of the power spectrum S(f) is obtained. a =m0 / m1, where S(f n ) represents the power spectrum s(f) function value corresponding to the nth frequency.

6. The method for calculating the wave spectrum based on the tilted attitude according to claim 5, characterized in that: According to the formula θ 1n =arctan(B1 / A1), get the wave average direction time series, calculate the spectrum peak frequency f of the wave average direction time series p The corresponding angle D mean =θ 1n (f p ), and D mean As the average wave direction corresponding to the wave peak frequency, B1 and A1 are the first-order Fourier coefficients of the directional spectrum.

7. The method for calculating the wave spectrum based on the tilted attitude according to claim 6, characterized in that: According to the formula θ 2n =0.5*arctan(B2 / A2), get the main wave direction time series of the wave, calculate the spectrum peak frequency f of the main wave direction time series of the wave p The corresponding angle D m =θ 2n (f p ), and D m As the main direction of the wave corresponding to the wave peak frequency, B2 and A2 are the second-order Fourier coefficients of the directional spectrum.

Images