基于惯性测量的海浪谱计算方法

1.本发明涉及海洋观测应用领域，尤其涉及一种基于惯性测量的海浪谱计算方法。


背景技术：

2.研究海洋波浪，通常将波浪视为多组随机的正弦波叠加而成，每一组正弦波的频率不同但高度均值为零，从而叠加得到波面高度的平均值为零，即不随时间发生变化。通常来说，波浪在很短的观测时间间隔内(10
‑
30分钟)，高度方差的变化很缓慢，可将海浪在此期间变化视为准平稳随机过程。因此，利用谱方法以随机过程来描述海浪成为一种非常重要研究途径。海浪谱包括功率谱s(f)和方向谱s(f,θ)。海浪功率谱s(f)是反应能量的密度分度，给出不同频率组成波提供的能量相对于频率的分布。一般来说，对于海浪组成波中，能量主要由集中在一狭窄频率带内的组成波提供，而功率谱两端的频率很小或很大者，提供的能量很小。这种海浪谱特征反映到海洋的肉眼客观的表现为：某海区海浪高低长短不齐，但其中小频率(长周期重力波)和高频率(短周期毛细波、碎杂波等)的波高很小，而其中会有一些周期性明显的大浪占据着主要地位。方向谱s(f,θ)为反应海浪内部方向结构的能谱，可给出不同方向上各组成波的能量相对于频率的分布。就给定波浪频率而言，方向谱可给出不同方向间隔内的能量，也就是给定频率的组成波能量相对于波向的分布。
3.国内外主要将波浪测量系统(惯性波浪传感器、压力式波浪传感器等)安装在浮标内用于获取波浪特征参数，通过横滚角(roll)、俯仰角(pitch)、垂直方向位移(d_z)这3个时间序列间的交叉谱计算来获取方向谱、平均波向和主波向(称为prz方法)。但当浮标尺寸形状较大时，容易造成浮标晃动倾斜角较小甚至无倾斜的状况，导致方向谱计算出现误差。


技术实现要素：

4.针对上述问题，本发明提出一种基于惯性测量的海浪谱计算方法，主要解决大尺寸浮标晃动倾斜角较小或无倾斜导致方向谱计算存在误差的问题。
5.为解决上述技术问题，本发明的技术方案如下：
6.一种基于惯性测量的海浪谱计算方法，包括以下步骤：
7.获取预设时间内的n组水平系加速度，n≥2；对所述水平系的x轴和y轴加速度进行矢量分解，得到地理坐标系正东方向和正北方向的加速度分量，分别定义为东方向加速度a_e以及北方向加速度a_n；对所述水平系的z轴加速度进行2次频域数值积分得到水平系垂直位移序列d_z，计算所述水平系垂直位移序列d_z的自交谱，得到序列c
11
；计算所述水平系垂直位移序列d_z与所述东方向加速度a_e的正交谱，得到序列q
12
；计算所述水平系垂直位移序列d_z与所述北方向加速度a_n的正交谱，得到序列q
13
；计算所述东方向加速度a_e的自交谱，得到序列c
22
；计算所述北方向加速度a_n的自交谱，得到序列c
33
；计算所述东方向加速度a_e和所述北方向加速度a_n的协交谱，得到序列c
23
；根据序列c
11
、q
12
、q
13
、c
22
、c
33
和c
23
计算方向谱s(f,θ)的傅里叶分解的前m个系数，m≥2，将m个所述系数合成所述方向谱s(f,θ)，并对所述方向谱s(f,θ)进行平滑。
8.在一些实施方式中，对n组水平系垂直位移序列d_z进行快速傅氏变换，获取n组实部和虚部，然后求取实部和虚部的平方和并除以采样点数n，得到水平系垂直位移序列d_z的粗糙功率谱，并对所述粗糙功率谱进行平滑，得到平滑后的功率谱s(f)。
9.在一些实施方式中，计算所述水平系垂直位移序列d_z的功率谱s(f)的谱峰值频率f
p
，得到波浪峰值周期t
p
，所述功率谱s(f)自适应调整滤波器带宽为1/(t
p
+5)<f
p
<1/(t
p
‑
5)，仅保留主波周期前后k秒范围内的组成波成分，5≤k≤15。
10.在一些实施方式中，按离散积分公式对所述功率谱s(f)进行保留信号频段积分得到波浪零阶矩m0，求取所述功率谱s(f)的谱特征有效波高
11.在一些实施方式中，按离散积分公式对所述功率谱s(f)进行保留频段积分得到波浪一阶矩m1，求取所述功率谱s(f)的谱特征平均波周期t
a
＝m0/m1，其中，s(fn)表示第n个频率所对应的功率谱s(f)函数值。
12.在一些实施方式中，按公式θ
1n
＝arctan(b1/a1)，得到波浪平均方向时间序列，计算所述波浪平均方向时间序列的谱峰值频率f
p
所对应的角度d
mean
＝θ
1n
(f
p
)，并将d
mean
作为波浪峰值频率处所对应的波浪平均方向。
13.在一些实施方式中，按公式θ
2n
＝0.5*arctan(b2/a2)，得到波浪主波向时间序列，计算所述波浪主波向时间序列的谱峰值频率f
p
所对应的角度d
m
＝θ
2n
(f
p
)，并将d
m
作为波浪峰值频率处所对应的波浪主方向。
14.本发明的有益效果为：本方法采用东方向加速度、北方向加速度、水平系垂直位移序列3个时间序列间的交叉谱来获取方向谱、平均波向和主波向，相对于传统方法，可解决因浮标尺寸形状较大而造成浮标晃动倾斜角较小甚至无倾斜状况而造成的方向谱计算误差。
附图说明
15.图1为本发明实施例公开的基于惯性测量的海浪谱计算方法的流程示意图。
具体实施方式
16.为使本发明的目的、技术方案及优点更加清楚、明确，下面结合附图和具体实施方式对本发明的内容做进一步详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释本发明，而非对本发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本发明相关的部分而非全部内容。
17.如图1所示,本实施例提出了一种基于惯性测量的海浪谱计算方法，可用于板载计算，以提高波浪测量系统的性能，主要包括以下步骤：
18.101，获取预设时间内的n组水平系加速度，n≥2。
19.在本实施例中，按照设定的采样频率进行载体系3轴加速度、3轴角加速度、3轴磁力计采集，并同步进行基于9轴参数的互补融合滤波算法进行实时姿态解算获取3轴姿态角(pitch、roll、yaw)，获取水平系3轴加速度(a_x、a_y、a_z)。提前设置的采样频率2/4/8hz可
选、采样时间长度4/8/16分钟可选，当采样时间到，进入波浪计算程序。
20.102，结合航向角yaw对水平系的x轴(a_x)和y轴(a_y)加速度进行矢量分解，得到地理坐标系正东方向和正北方向的加速度分量，分别定义为东方向加速度a_e以及北方向加速度a_n；
21.在本实施例中，对东方向加速度a_e做fft变换，实部/虚部分别存储在预定的存储空间；对北方向加速度a_n做fft变换，实部/虚部分别存储在预定的存储空间；
22.计算各变量之间的协交谱c
xy
、自交谱q
xy
，获得对应的序列，计算公式分别为
23.c
xy
＝(re[x]*re[y]
‑
im[x]*im[y])/n，q
xy
＝(im[x]*re[y]+re[x]*im[y])/n，具体如步骤103
‑
108，无先后顺序之分：
[0024]
103，对水平系的z轴加速度进行2次频域数值积分得到水平系垂直位移序列d_z，计算水平系垂直位移序列d_z的自交谱，得到序列c
11
；
[0025]
104，计算水平系垂直位移序列d_z与东方向加速度a_e的正交谱，得到序列q
12
；
[0026]
105，计算水平系垂直位移序列d_z与北方向加速度a_n的正交谱，得到序列q
13
；
[0027]
106，计算东方向加速度a_e的自交谱，得到序列c
22
；
[0028]
107，计算北方向加速度a_n的自交谱，得到序列c
33
；
[0029]
108，计算东方向加速度a_e和北方向加速度a_n的协交谱，得到序列c
23
；
[0030]
109，根据序列c
11
、q
12
、q
13
、c
22
、c
33
和c
23
计算方向谱s(f,θ)的傅里叶分解的前m个系数，m≥2，将m个系数合成方向谱s(f,θ)，并对方向谱s(f,θ)进行平滑。
[0031]
在本实施例中，m取5，因此前5个系数为a0，a1，a2，b1，b2，计算公式分别为：a0＝c
11
；a1＝q
12
；b1＝q
13
；a2＝c
22
+c
33
；b2＝2c
23
，其中系数k为每次测量得到的波的个数，可由程序前期做波浪时域跨零统计时得到。根据傅里叶级数分解公式：s(f,θ)＝a0+a
1 cosθ+b
1 sinθ+a
2 cos2θ+b
2 sin2θ，分别代入5项傅里叶系数，计算得到波浪方向谱s(f,θ)。因仅采用前5项来近似代替方向谱函数的傅里叶级数，会对真实结果产生一定的误差。因此，还需对方向谱进行平滑，多次试验后发现，权重系数取2/5、3/4、1/5时效果最好，平滑公式为：
[0032][0033]
方向谱s(f,θ)的板载计算过程中，因是二维变量函数。若按照频率分量为4096组(4hz，频率间隔为4/4096)，方向分量为360组(360度范围按照1度间隔划分)来代入计算，需要5760k的内存，远远超出stm32f427的256k内存。而且，计算耗时较长，不能满足实时计算需求。结合海浪物理特性，在板载计算过程中适当降低计算网格分辨率，将频率分量f离散为128组(4hz，频率间隔为4/128)，方向分量θ离散为72组(360度范围按照5度间隔划分)，最终得到9216组s(f,θ)值，只需要36k内存。
[0034]
本方法采用东方向加速度、北方向加速度、水平系垂直位移序列3个时间序列间的交叉谱来获取方向谱、平均波向和主波向，相对于传统方法，可解决因浮标尺寸形状较大而造成浮标晃动倾斜角较小甚至无倾斜状况而造成的方向谱计算误差。具体地，其原理为：通过两个倾斜角度(横滚、俯仰)与垂直方向位移进行方向谱计算的浮标在小型波浪载体浮标上的效果较好，因为小型浮标的随波性明显，其随波在两个角度上均可以发生明显的晃动，因此可以很好的识别方向。但在一些大型的浮标上，其稳定性较好，随波不一定会发生两个轴上的明显晃动，因此在大型浮标上提供晃动角度的方法计算方向谱效果不好(很大的波
浪才会引起大型浮标的很小的倾斜角度，大型浮标对倾斜角度变化不是很敏感)。但不管在大型/小型浮标上，水平方向的两个加速度的灵敏度要高于角度，及时一个很小的波浪冲击过来，在大型浮标的东和北两个方向上都能测量到明显的加速度变化，因此通过东/北两个方向的加速度和垂直方向的位移的协谱计算波浪谱，具有更好的平台适应性。
[0035]
进一步地，通过功率谱s(f)和方向谱s(f,θ)反演了波浪特征参数t
p
,t
a
，h
m0
,d
mean
,d
m
，如下：
[0036]
201，对n组水平系垂直位移序列d_z进行快速傅氏变换，获取n组实部和虚部，然后求取实部和虚部的平方和并除以采样点数n，得到水平系垂直位移序列d_z的粗糙功率谱，并对粗糙功率谱进行平滑，得到平滑后的功率谱s(f)。
[0037]
202，计算水平系垂直位移序列d_z的功率谱s(f)的谱峰值频率f
p
，得到波浪峰值周期t
p
，功率谱s(f)自适应调整滤波器带宽为1/(t
p
+5)<f
p
<1/(t
p
‑
5)，仅保留主波周期前后k秒范围内的组成波成分，5≤k≤15。
[0038]
在本实施例中，k取10。
[0039]
在本实施例中，对步骤201得到的粗糙功率谱进行3点滑动平均，多次试验后发现，权重系数取1/4,1/2,1/4时效果最好，公式为：sf[n]＝sf[n
‑
1]/4+sf[n]/2+sf[n+1]/4；进一步，功率谱sf[n]谱峰值所对应的频率f
p
为峰值频率，则波浪谱峰值周期t
p
＝1/f
p
。
[0040]
对保留频段内的sf[n]进行零阶次和一阶次积分，分别得到波浪功率谱的波浪零阶矩m0和波浪一阶矩m1，通过零阶矩和一阶矩反演有效波高h
m0
和平均波周期ta，公式分别如步骤203、204所示：
[0041]
203，按离散积分公式对功率谱s(f)进行保留信号频段积分得到波浪零阶矩m0，求取功率谱s(f)的谱特征有效波高
[0042]
204，按离散积分公式对功率谱s(f)进行保留频段积分得到波浪一阶矩m1，求取功率谱s(f)的谱特征平均波周期t
a
＝m0/m1，其中，s(f
n
)表示第n个频率所对应的功率谱s(f)函数值。
[0043]
205，按公式θ
1n
＝arctan(b1/a1)，得到波浪平均方向时间序列，计算波浪平均方向时间序列的谱峰值频率f
p
所对应的角度d
mean
＝θ
1n
(f
p
)，并将d
mean
作为波浪峰值频率处所对应的波浪平均方向。
[0044]
206，按公式θ
2n
＝0.5*arctan(b2/a2)，得到波浪主波向时间序列，计算波浪主波向时间序列的谱峰值频率f
p
所对应的角度d
m
＝θ
2n
(f
p
)，并将d
m
作为波浪峰值频率处所对应的波浪主方向。
[0045]
上述实施例只是为了说明本发明的技术构思及特点，其目的是在于让本领域内的普通技术人员能够了解本发明的内容并据以实施，并不能以此限制本发明的保护范围。凡是根据本发明内容的实质所做出的等效的变化或修饰，都应涵盖在本发明的保护范围内。