本发明涉及海洋波浪观测领域，具体涉及一种捷联惯导式测波方法及系统。
背景技术：
在海洋波浪观测领域，加速度式波浪观测设备因其不受时空范围影响、观测周期范围广等优势，自问世30余年来一直是海洋波浪观测的主要手段。长期以来，人们对加速度测波技术(加速度-位移积分，信号滤波等)和波浪数据处理方法(时域特征统计，频域谱分析等)倾注了极大的关注。对基于高性能嵌入式技术的加速度波浪传感器的研制，则是随着电子技术的发展始于最近20余年。加速度测波技术中，一个重要的环节就是通过惯性测量单元imu(3轴角加速度计、3轴加速度计、3轴磁力计或电子罗盘)获取载体平台(浮标、船舶、海洋平台)在自然地理坐标系下的的加速度、姿态角，进而通过二次积分、数字滤波技术、跨零法统计获取波浪的特征参数。目前，波浪传感器中通过imu单元获取自然地理坐标系下加速度和姿态角的方法，主要分为平台惯导式和捷联惯导式这两类。
平台惯导系统机电结构复杂，精度要求较高，对系统的加工工艺、装配精度都需要很高的要求，成本更高，往往一个惯性平台的加工费用就高达整个惯导系统40％以上，且系统稳定性较差，不利于在小体积的波浪传感器上使用。捷联惯导式通过固定在载体上的“数字惯导平台”代替机电式惯性平台，大大减少了系统的体积和制造成本，提高系统的稳定性，但对姿态解算算法、嵌入式计算机处理性能、惯性敏感元件的精度提出了更高要求。因此，基于捷联惯导式波浪测量技术会是下一代波浪传感器的发展方向。捷联惯导式测波技术能够有效的消除海洋浮动平台自身在动态测量环境中带来的不利影响(升沉、纵摇、横摇等)，可稳定获取载体平台在自然地理坐标系下姿态角、航向角、加速度等参数，提高波浪特征参数的计算精度。
根据惯性导航系统的基本原理，运动载体的速度及位移都是基于载体坐标系下加速度2次积分获取位移获取，在积分过程中误差及初始速度不确定性会随着积分在速度和位移中不断累加，需采用合理的数值积分策略和滤波技术方能得到可靠的位移精度。以往的研究多采用硬件积分和硬件滤波器的方法实现，ribiero的研究表明硬件积分器对正弦稳态信号积分效果较理想，而对于随机性较强的波浪信号积分则误差很大，相对来说数值积分可取得较高的精度并且更加灵活。硬件滤波器的滤波性能是固定的，对于周期随机性较强的波浪来说，不具有自适应性。不同的数值积分处理策略产生的积分结果差别很大，探索合理的加速度-位移积分策略是目前研究的热点，其中包括加速度信号滤波方法选取(fir,iir)，不同采样频率下积分方案选取(梯形，simpson等)，位移趋势项剔除(滑动滤波、最小二乘、经验模态分解等)。
在数据处理平台方面，因为受限于嵌入式微处理的性能，前期的波浪传感器多是按照指定的采样频率(2hz,4hz或者8hz)进行imu单元采样，完成一段时间的采样后(10min，20min或60min)，再进行积分运算得到相同频率的姿态角、位移数组，然后按照跨零法进行该观测时间段内的波浪特征统计。随着嵌入式处理器性能的快速发展，按照较高的频率边采样边实时计算输出姿态、速度、位移等信息已变为可行，这样可提高仪器的实时观测能力及标准化。
基于以上几点，目前基于稳定惯性平台式波浪传感器存在的仪器成本高、系统可靠性低、缺乏实时观测计算能力、体积和质量较大不利于安装使用等问题。因此，从提高波浪传感器观测精度，提高仪器性能，降低仪器成本、体积和质量，促进仪器标准化和创新自主知识产权的角度出发，捷联惯导式波浪观测技术将是下一代主流。
技术实现要素：
本发明的目的在于针对上述现有技术中存在的问题，提出了一种捷联惯导式测波方法及系统，实现实时采样及计算输出姿态、速度和位移信息，并提高了仪器的实时观测能力及标准化，提高了仪器性能，并降低了仪器成本、体积和质量。
为达到上述发明的目的，本发明通过以下技术方案实现：
本发明还公开一种捷联惯导式测波方法，包括步骤如下：
步骤1，获取载体随波浪在三维空间的运动轨迹数据和运动姿态数据；
步骤2，对运动轨迹数据和运动姿态数据进行捷联惯导姿态解算，获取载体的欧拉姿态角和姿态矩阵，并根据姿态矩阵获得载体在自然地理坐标系下的加速度；
步骤3，对自然坐标系下的加速度，进行二次积分，获取载体的位移量；
步骤4，根据载体的位移量，进行跨零法波浪统计，将最大出现率的方向角作为主波向。
进一步，所述步骤1还包括对运动轨迹数据和运动姿态数据的原始数据进行平均值滤波和误差补偿的操作。所述捷联惯导姿态解算采用互补滤波和四元数法进行姿态解算。
进一步，所述步骤2还包括关于姿态解算的判断过程，包括：如果是首次姿态解算，则首先根据加速度计和陀螺仪的数据求出初始欧拉姿态角，进而初始化四元数；若不是首次姿态解算，则直接才用互补算法，用加速度对陀螺仪进行修正获取准确稳定的陀螺仪数值，带入四元数微分方程；有了初始化的四元数量，再结合数据数据采样周期及修正后的陀螺仪数据，即可对四元数进行更新，获取姿态矩阵
进一步，所述二次积分是加速度、速度和位移的数值积分，包括如下步骤：
步骤31，分别在三轴方向上继续对加速度做积分运算，自然地理坐标系下的某一轴加速度为an，经过中值滤波处理防止加速度传感器的脉冲干扰，得到an1；
步骤32，通过高通滤波器去除加速度传感器的零点漂移趋势，得到an2；
步骤33，去除直流重力加速度分量g，得到该轴真实的运动加速度an3；
步骤34，根据公式v(k)＝v(k-1)+ts·a(k)，数值离散化后进行时域一次数值积分得到速度v0；
步骤35，通过高通滤波去除速度序列的低频漂移得到速度v1；
步骤36，根据公式s(k)＝s(k-1)+ts·v(k)，数值离散化后进行时域二次数值积分得到位移s0；
步骤37，通过高通滤波器去除位移序列的低频漂移得到位移s1。
进一步，所述跨零法波浪统计是统计一段时间段内的波浪特征参数，所述波浪特征参数包括波高、波周期和波浪方向，所述波高包括平均波高havg、1/3大波高h1/3、1/10大波高h1/10和最大波高hmax；所述波周期包括平均波周期tavg、1/3大波周期t1/3、1/10大波周期t1/10和最大波周期tmax；所述波浪方向包括平均波向davg和主波向ddm。
进一步，所述波高和波周期的统计方法如下：
步骤41，取z轴位移平均值为波高零基准面；
步骤42，寻找相邻的跨零点，并标记相应编号；
步骤43，计算相邻跨零点内的波周期、波高，并存储在指定波高和波周期数组序列；
步骤44，波高和波周期序列按照大小排列：统计最大波高、平均波高、有效波高、1/3大波高、1/10大波高及相应周期。
进一步，所波浪方向的统计方法为：根据x，y轴位移序列合成得到波浪方向序列，结合磁力计校正得到真实的波浪方向；以22.5°为间隔，划分波向为16个方位，统计每个方向波向出现频率；取最大出现率的方位角作为主波向。
进一步，所述高通滤波采用emd经验模态分解方法，所述emd经验模态分解方法是将信号的不同尺度的波动逐级分开，产生一系列具有不同尺度的本征模函数。
本发明还公开一种捷联惯导式测波系统，采用上述方法，包括全姿态imu测量单元和数据处理单元，所述全姿态imu测量单元包括三轴陀螺仪、三轴加速度计和三轴磁强计；所述数据处理单元包括有：
姿态数据接收模块，用以获取载体随波浪在三维空间的运动轨迹；
姿态解算模块，用以对对运动轨迹数据和运动姿态数据进行捷联惯导姿态解算，获取载体的欧拉姿态角和姿态矩阵，并根据姿态矩阵获得载体在自然地理坐标系下的加速度；
积分模块，用以对自然坐标系下的加速度，进行二次积分，获取载体的位移量；
统计模块，用以根据载体的位移量，进行跨零法波浪统计，将最大出现率的方向角作为主波向。
本发明的一种捷联惯导式测波方法及系统，可实现系统的体积小、重量轻等特点，测波方法在动态波浪测量环境中具有较好的稳定性和可靠性，可代替传统惯性稳定平台式波浪传感器的应用场合，降低波浪观测成本，提高观测精度和布放灵活度。
附图说明
图1为本发明的一种捷联惯导式测波方法的步骤图。
图2是本发明的一种捷联惯导式测波方法的整体技术路线。
图3是本发明的一种捷联惯导式测波方法的数据流及整体姿态解算流程。
图4是本发明的一种捷联惯导式测波方法的四元数的姿态解算互补滤波算法流程。
图5是本发明的一种捷联惯导式测波方法加速度-位移数值积分流程。
图6是本发明的一种捷联惯导式测波方法的跨零法统计z轴波浪特征流程。
图7是本发明的一种捷联惯导式测波系统的结构框图。
具体实施方式
下面结合附图和实施例对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部实施例。
参看图1，为本发明实施例的一种捷联惯导式测波方法的流程图，包括步骤如下：
步骤1，获取载体随波浪在三维空间的运动轨迹数据和运动姿态数据；
步骤2，对运动轨迹数据和运动姿态数据进行捷联惯导姿态解算，获取载体的欧拉姿态角和姿态矩阵，并根据姿态矩阵获得载体在自然地理坐标系下的加速度；
步骤3，对自然坐标系下的加速度，进行二次积分，获取载体的位移量；
步骤4，根据载体的位移量，进行跨零法波浪统计，将最大出现率的方向角作为主波向。
参看图2至图6，为本发明的一种捷联惯导式测波方法整体及部分算法实现流程。
参看图1和图2，作为一个实施例，所述步骤1还包括对运动轨迹数据和运动姿态数据的原始数据进行平均值滤波和误差补偿的操作。所述捷联惯导姿态解算采用互补滤波和四元数法进行姿态解算。本发明实施例中，采用包括三轴加速度计、三轴陀螺仪和三轴磁强计的9轴全姿态imu，stm32f407芯片为数据处理器，其以用户设定的采样频率(2hz，4hz，8hz)对imu采样，从而获取三维空间的载体运动轨迹和运动姿态的原始数据。在获取原始数据后，采用互补滤波算法，用加速度计来补偿陀螺仪，获取稳定、准确的陀螺仪数据。进而，利用四元数法解算姿态矩阵和姿态欧拉角。进而，通过姿态矩阵完成坐标系转换，获取载体在自然地理坐标系下的加速度。进而，将自然坐标系下的加速度数据进行二次数值积分与高通滤波得到位移。进而，跨零法统计指定时间段内的位移序列，对z轴位移序列统计可以得到波高和波周期，对x，y序列统计并通过三角函数合成可以得到载体波向，结合磁力计得到的方位角进行校正，即可得到真实地理坐标系下的波浪方向。
下面结合图3至图6分别对各步骤做详细描述。
参看图3，所述步骤2，采用基于四元数的互补滤波算法，主要是完成坐标系的变换，将载体坐标系下的加速度变换到自然地理坐标系。首先用加速度传感器数据补偿修正陀螺仪数据，获得在载体坐标系下稳定可靠陀螺仪数据；进而对姿态进行解算，获得载体欧拉姿态角和姿态矩阵
参看图4，首先stm32f407完成对全姿态imu9150的初始化，启动对其中三轴陀螺仪和三轴加速度仪的原始数据采样，对原始数据进行平均值滤波和误差补偿，如果系统是第一次姿态解算则首先根据加速度计和陀螺仪的数据求出初始欧拉姿态角，进而初始化四元数。如果系统不是第一次姿态解算，则说明系统已经进行过了欧拉角和四元数初始化，则直接进入互补算法阶段，用加速度值对陀螺仪进行修正获取准确稳定的陀螺仪数值，带入四元数微分方程。有了初始化的四元数量，再结合数据数据采样周期及修正后的陀螺仪数据，即可对四元数进行更新，获取姿态矩阵
上述姿态解算及坐标系变换，已经获取载体在自然地理坐标系下的三轴加速度，需对加速度进行二次积分，获取位移。加速度-速度-位移的积分策略如图5所示，以z轴积分为例，自然地理坐标系下的z轴加速度为an，经过中值滤波处理防止加速度传感器的脉冲干扰得到an1，进一步通过高通滤波器祛除加速度传感器的零点漂移趋势得到an2，进一步去除直流重力加速度分量g，得到z轴真实的运动加速度an3，进一步根据公式v(k)＝v(k-1)+ts.a(k)，数值离散化后进行时域一次数值积分得到速度v0，进一步通过高通滤波祛除速度序列的低频漂移得到速度v1，进一步根据公式s(k)＝s(k-1)+ts.v(k)，数值离散化后进行时域二次数值积分得到位移s0，进一步通过高通滤波器祛除位移序列的低频漂移得到位移s1，供波浪特征统计使用。
本发明实施例中祛除趋势漂移项的高通滤波采用基于emd经验模态分解方法，emd分解是将信号的不同尺度的波动逐级分开，产生一系列具有不同尺度的本征模函数(imf)。本专利中位移序s(t)列经emd分解后，产生12组imf，根据信号采样频率和重力波周期范围(2～30秒)，选取imf1(代表3～4秒尺度的波动)，imf2(代表3.5～5秒尺度的波动)，imf3(代表5～10秒尺度的波动)，imf4(代表12～20秒尺度的波动)，imf5(代表20～30秒尺度的波动)，剔除低频趋势项imf6-imf11(30秒以上的更大尺度的波动，视为漂移趋势项)，最终位移序列s0(t)计算公式为：s(t)＝s0(t)-(imf6+…+imf11)。
emd分解方法是依据信号数据本身固有的特性自然地展开分解，无需设置先验的分解函数基，也无需预测信号类型。相比于低通滑动滤波法、小波变化法、最小二乘法均具有较强的自适应性，更适合应用在随机波浪加速度信号处理中。
对于海洋中随机性很强的波浪研究，可将其运动看作为随机过程，作为一种随机现象加以分析。应用统计学的概念和方法可以分析海浪的对外表观(波高、波周期、波向)。如图6所示，采用跨零方法进行波浪时域特征统计分析。
波高和波周期的统计方法为：首选取z轴位移平均值为波高零基准面；进一步寻找相邻的跨零点，并标记相应编号；进一步计算相邻跨零点内的波周期、波高，并存储在指定波高和波周期数组序列；进一步波高和波周期序列按照大小排列；统计平均波高havg，平均波周期tavg，1/3大波高h1/3，1/3大波周期t1/3，1/10大波高h1/10，1/10大波周期t1/10，最大波高hmax，最大波周期tmax，平均波向davg；
波向统计方法为：根据x，y轴位移序列合成得到波浪方向序列，结合磁力计校正得到真实的波浪方向；以22.5°为间隔，划分波向为16个方位，统计每个方向波向出现频率；最大出现率的方位角作为主波向ddm；使用x,y轴位移合成的波浪方向，是传感器随波浪真实的运动方向，相对于使用载体纵倾和横倾斜计算浮标倾斜方位角的方法反演波浪方向更加准确。而且，改方法仅需要加速度数据即可，不需使用陀螺仪数据，计算方法更加直观简单，容易实现。
参看图7，为本发明实施例的一种捷联惯导式测波系统，采用上述捷联惯导式测波方法，包括全姿态imu测量单元和数据处理单元，所述数据处理单元包括有：
姿态数据接收模块，用以获取载体随波浪在三维空间的运动轨迹；
姿态解算模块，用以对对运动轨迹数据和运动姿态数据进行捷联惯导姿态解算，获取载体的欧拉姿态角和姿态矩阵，并根据姿态矩阵获得载体在自然地理坐标系下的加速度；
积分模块，用以对自然坐标系下的加速度，进行二次积分，获取载体的位移量；
统计模块，用以根据载体的位移量，进行跨零法波浪统计，将最大出现率的方向角作为主波向；
所述全姿态imu测量单元，为9轴全姿态imu，包括三轴加速度计、三轴陀螺仪和三轴磁强计，stm32f407以用户设定的采样频率(2hz，4hz，8hz)对imu采样。，包括三轴陀螺仪、三轴加速度计和三轴磁强计；用来测量载体随波浪在三维空间的运动规矩和运动姿态，三轴的加速度计用来测量载体坐标系下x，y，z轴的加速度值，三轴的角加速度计用来测量俯仰角pitch、翻滚角roll、偏航角yaw，三轴的磁力计用来测量地磁场在载体坐标系下x，y，z轴的矢量，获取载体的方位角，用于校正波浪方向。
所述姿态数据接收模块、姿态解算模块、积分模块和统计模块是固化在数据处理单元的芯片中的逻辑电路，参看图6，9轴全姿态imu单元mpu9150通过iic总线与stm32f407连接，用于波浪姿态的测量；eeprom通过spi总线与stm32f407连接，预存储用户设置的参数，如：imu采样频率，波浪特征统计间隔等；实时时钟rtc通过iic总线与stm32f407连接，为系统提供时间和闹钟；stm32f407通过专用sdc接口连接存储卡，移植fat32文件系统，用于存储系统原始数据、中间变量、波浪统计特征等；通过通用i/o口连接lcd显示屏，显示实时姿态数据及波浪统计数据；通过uart串口连接上位机，与上位机进行通讯，用于用户进行相关预设值。
上述实施例仅用以说明本发明而并非限制本发明所描述的技术方案；因此，尽管本说明书参照上述的各个实施例对本发明已进行了详细的说明，但是，本领域的普通技术人员应当理解，仍然可以对本发明进行修改或者等同替换；而一切不脱离本发明的精神和范围的技术方案及其改进，其均应涵盖在本发明的权利要求范围当中。