CN115165019A - 一种在航gnss潮位测量方法 - Google Patents

Abstract

本发明公开了一种在航GNSS潮位测量方法，属于导航技术领域，用于获取潮位序列，包括：安装测量装置，计算海面大地高程序列，采用变分模态分解和小波变换联合消噪法提取基于大地高程的潮位序列，计算海底格网点的大地高程，获取测量区域GPS水准联测点的高程信息，采用二次多项式曲面拟合法求解海底格网点近似高程异常，大地高程与近似高程异常求和得到海底格网点的近似正常高，求解海底格网点的精确高程异常，采用双线性内插法求解测船处的高程异常，最终得到基于1985国家高程基准的潮位序列。对比现有技术，本发明提出的消噪方法提取的潮位均方根误差更小，更接近真实潮位值。利用声呐精确测量海底点的坐标和大地高，提高了垂直基准转换的精度。

Description

一种在航GNSS潮位测量方法

技术领域

本发明公开了一种在航GNSS潮位测量方法，属于导航技术领域。

背景技术

在利用多波束声呐进行海底地形测量的过程中，同步获得调查区的潮位数据，完成潮位改正是数据处理不可或缺的步骤。传统的潮位测量方式包括水尺验潮、数字式压力验潮仪验潮和长期验潮站验潮。水尺验潮方法是将水尺竖直固定于码头壁、海滩上，人工读取水面在水尺上的位置，记录海面在水尺零点上的高度，并结合钟表计时；数字式压力验潮仪验潮是通过测量海水的压力变化而间接推算出海面的升降变化；长期验潮站验潮是在国家设立的长期固定验潮点利用水尺或验潮仪观测潮位。

传统的验潮方式常受潮位模型误差、无验潮站、验潮仪丢失、耗费大量的人力资源、耗时等诸多不利因素的影响。随着GNSS载波相位差分测量技术的不断完善，GNSS-RTK、GNSS-PPK、网络RTK潮位测量方法在我国逐步推广应用。GNSS验潮是利用GPS载波相位技术得到船载天线高精度的高程信息，通过天线与海面的高差得到海面的高程序列，利用滤波方法提取出海面高程序列中的低频信息，再通过垂直基准转换得到GNSS在航潮位。

现有的GNSS验潮技术方案包括坐标系建立、GNSS接收机大地高数据的获取、海面大地高序列的滤波和垂直基准转换，具体步骤如下：

(1)测量开始前，以IMU为参考原点建立船体坐标系，船艏方向为X轴，垂直X轴指向右舷为Y轴，垂直XOY平面指向下为Z轴。将两台GNSS接收机的坐标在船体坐标系下表达。

(2)利用GNSS-RTK或GNSS-PPK或网络RTK技术获得GNSS接收机的经、纬度和大地高程，结合步骤(1)所述GNSS接收机在所建立的坐标系下的坐标、IMU测量的船体姿态，压力测深传感器测量的动态吃水来获得海面的大地高序列，并通过快速傅里叶变换(傅里叶变换法的流程图如图1)或小波强制消噪的方法去除海面的大地高序列中的高频噪声，再通过高程基准转换的方式获得1985国家高程基准下的潮位值，高程基准转换的传统方式为几何法或转换参数法。

目前传统的高程基准转换方式都没有考虑到利用多波束精确测量的海底格网点的坐标、大地高和平均大地高来进一步提高垂直基准转换的精度，GNSS验潮采用的快速傅里叶变换滤波方法不能处理非平稳信号和突变信号小波变换强制消噪小波基选取困难、缺乏自适应性，且两种滤波方法的滤波效果差。在海底地形复杂区进行潮位的垂直基准转换时，单一的几何法或转换参数法的转换精度低。

发明内容

本发明公开一种在航GNSS潮位测量方法，以解决现有技术中GNSS验潮方法滤波效果差、转换精度低的问题。

一种在航GNSS潮位测量方法，包括：

S1.安装测量装置于船舷，所述测量装置包括多波束声呐、IMU、压力测深传感器和两台GNSS接收机，IMU为惯性测量单元；

GNSS接收机独立接收卫星信号，获取GNSS接收机相位中心的经纬度和大地高程，声呐测量测区的水下地形，获得水下测量点的经纬度、大地高数据和深度，IMU获取IMU处的升沉值和测船的姿态值，所述姿态值包括横摇值、纵摇值和艏摇值，所述压力测深传感器测量测船运动时的动态吃水状态下IMU到水面的距离和船抛锚时IMU到水面的距离；

S2.综合S1获得的数据计算WGS84高程基准下的海面大地高程序列；

S3.采用变分模态分解和小波变换联合消噪法提取海面大地高程序列中基于大地高程的潮位序列；

S4.求解测量区基于WGS84高程基准的海底格网点的大地高程；

S5.获得测量区域附近6个或6个以上的GPS水准联测点的高程信息，所述高程信息同时包括GPS水准联测点的大地高程和1985国家高程基准下的高程；

S6.采用二次多项式曲面拟合法求解测量区域海底格网点的近似正常高；

S7.根据莫洛金斯基原理求解海底格网点的精确高程异常；

S8.根据S1中的GNSS接收机相位中心的经纬度，采用双线性内插法求解出测船处的精确高程异常，将S3中的潮位序列归算到1985国家高程基准下，得到基于1985国家高程基准的潮位序列。

优选地，S1包括以下子步骤：

S1.1.以IMU为参考原点建立理想船体坐标系，在制图软件中建好测量装置的三维模型，通过三维模型量取测量装置上各个部件之间的偏移关系，在理想船体坐标系下表达两台GNSS接收机的坐标，通过三维模型量取参考原点到GNSS接收机的垂直距离并记为L1，再通过压力测深传感器获得船体抛锚时IMU到水面的距离L2，即可获得GNSS接收机到水面的距离L3，计算公式如下：L3＝L1-L2；

S1.2.声呐测量时，在典型区域采用“斑片测试”获得换能器的横摇安装偏差角和纵摇安装偏差角；

S1.3.根据S1.2获得的横摇安装偏差角和纵摇安装偏差角对IMU测得的横摇值、纵摇值进行校正，若声呐***有横摇补偿功能，则只对IMU测得的纵摇值进行校正；

S1.4.用GNSS-PPK、GNSS-RTK、网络RTK或PPP的方式获得两台GNSS接收机相位中心精确的大地高和经纬度数据。

优选地，S2包括以下子步骤：

S2.1.根据S1.1中两台GNSS接收机在理想船体坐标系中的坐标和S1.3中校正后的纵摇值、横摇值，分别计算两台GNSS接收机的诱导升沉值H_i：H_i＝xsin P-ysin Rcos P-z(cos Rcos P-1)，H_i表示诱导升沉值，x、y、z表示GNSS接收机在理想船体坐标系中的坐标，P、R表示S1中IMU记录的纵摇值和横摇值；

S2.2.根据S2.1中的诱导升沉值H_i和S1中IMU处的升沉值，分别计算两台GNSS接收机处的升沉值h_G：h_G＝H_i+h₀，h₀表示S1中IMU处的升沉值；

S2.3.海面大地高程序列的计算过程为：h＝h_GPS-h_G-H_m+ΔH，h为海面的瞬时大地高，即海面大地高程序列；h_GPS为S1中GNSS接收机获得的相位中心的大地高程；H_m为L3；ΔH为S1中压力测深传感器测量测船运动时的动态吃水与S1.1中的L3之间的差值。

优选地，S3包括以下子步骤：

S3.1.采用变分模态分解法对S2得到的海面大地高序列逐次进行K层模态分解，K＝10n，n＝1，2，3……n，n为正整数，惩罚系数α＝80000，收敛容差tol＝1e-7；

S3.2.计算每一次对海面大地高序列进行K层模态分解得到的信号能量E：

n为采样点数，IMF_K(t)为每次模态分解的第K层信号序列；

S3.3.求解能量比值：

当能量比值γ趋于稳定时，γ在0值上下浮动，说明K为合适的模态分解数；

S3.4.对模态分解数为K时分解得到的第一层模态IMF₁再次进行变分模态分解，惩罚系数α＝80000，收敛容差tol＝1e-7，称为迭代分解1次，迭代分解1次的K值判断准则参考S3.2和S3.3，迭代分解1次得到的第一层模态记为IMF₁₁，第二层模态记为IMF₂₂，将IMF₁₁和IMF₂₂重构得到一新信号；

S3.5.对S3.4中的对IMF₁₁和IMF₂₂进行频谱分析，判断IMF₁₁和IMF₂₂的信号频谱组成；

S3.6.S2中的海面大地高程序列对应的最高频率ω为：

f_s为GNSS信号的采样频率；

假设海面高程序列的其他要素的频率为f_w，潮汐的频率为f_t，使用小波层数的层数z由下式决定:

将S3.4的IMF₁₁和IMF₂₂重构得到的新信号采用sym4小波分解z层后得到第z层低频系数A_z和高频系数D₁D₂………D_z，低频系数占有的频段为潮汐所在的频段，高频系数占有的频段为其他高频要素所在的频段，将高频系数置零之后与低频系数重构，即可实现海面高程序列的强制消噪，提取出的低频信号为基于WGS84高程基准的潮位序列h_{tide_84}。

优选地，S4包括以下子步骤：

S4.1.将声呐测量的水深数据进行安装偏差角校正、声速校正和粗差剔除；

S4.2.将S4.1处理后的水深数据采用S3中获取的h_{tide_84}进行潮位改正：

H_i＝|h_i-h_{tide_84}|，H_i为海底测量点的高程，h_i为S4.1处理后的水深数据；

S4.3.将S2中的海面大地高程序列按照500m*500m格网化，输出WGS84大地坐标系下的格网点的大地坐标，或者平面坐标和大地高h_Ω，将格网点各个的高程值加权取平均，得到的WGS84高程基准下的平均大地高，记为h_γ。

优选地，S6包括以下子步骤：

S6.1.在水准联测点上的高程异常与平面坐标之间，假定存在如下二次多项式曲面拟合数学模型：

a₀、a₁、a₂、a₃、a₄、a₅为模型待定参数，ζ_i为水准联测点的高程异常，当水准联测点多于6个时,组成形如下误差方程：V＝AX-ζ⁰，根据最小二乘原理可得A＝(X^TX)^-1X^Tζ⁰

按最小二乘原理解求出模型待定参数a₀、a₁、a₂、a₃、a₄、a₅的数值；

S6.2.将S4中的海底格网点的坐标代入到S6.1的二次多项式曲面拟合数学模型插值出海底格网点的近似高程异常ζ₀；

S6.3.测量区域海底格网点的近似正常高h_Γ为：h_Γ＝ζ₀+H_i，ζ₀表示近似高程异常,H_i为S4中的海底格网点的大地高程。

优选地，S7包括以下子步骤：

S7.1.获得S6得到的h_Γ；

S7.2.海底格网点的高程异常ζ由长波部分即近似高程异常ζ₀和短波部分ζ_r得到：ζ＝ζ₀+ζ_r；

S7.3.计算高程异常中的短波部分ζ_r：ζ_r＝T_r/r，T_r为地形起伏对海底格网点扰动位的影响；ζ_r为地形改正；

G为万有引力常数，G＝6.673×10^-8C³s^-2g^-1；ρ为地球平均质量密度，ρ＝2.67gC^-3，γ₀＝[(x-x_p)²+(y-y_p)²]^1/2,(x,y)为S4中的海底格网点的坐标，(x_p,y_p)为海底格网点中某个格网点P的坐标；

S7.4.似大地水准面的正常重力r：r＝r₀-0.3086h_Γ，r₀为参考椭球面上的正常重力值，位：10^-5…m/s²，

为WGS84大地坐标系中测点地理纬度，单位：度。

优选地，S8包括以下子步骤：

S8.1.根据S1中GNSS接收机相位中心的经纬度、S4中的海底格网点的大地高程、S7中的ζ，采用双线性插值法求取测船位置处的高程异常ζ₁；

S8.2.计算最终内插出的高程异常值f(x,y)：

f(f,y₁)、f(x,y₂)表示x方向插值出的高程异常值，ζ₁₁、ζ₂₁、ζ₁₂、ζ₂₂表示格网四个角点上的高程异常值，x₁、x₂、y₁、y₂表示四个角点在X方向和Y方向上的坐标，x、y表示测船在X方向和Y方向上的坐标；

S8.3.计算基于1985国家高程基准的潮位序列h_{tide_85}：h_{tide_85}＝h_{tide_84}+ζ₁。

本发明的主要有益效果为：用本发明新提出的变分模态分解和小波变换联合消噪法提取海面高程序列中的潮位的信噪比高，即本发明提出的消噪方法消噪效果好，且此种消噪方法提取的潮位、真实潮位之间的均方根误差小，即提取的GNSS在航潮位值更接近真实潮位值。利用声呐精确测量的海底格网点的坐标、大地高和平均大地高可提高垂直基准转换的精度。

附图说明

图1为现有技术中傅里叶变换法的流程图；

图2为测量装置的结构示意图；

图3为采用变分模态分解和小波变换联合消噪法提取海面高程序列中的潮位流程图；

图4为基于WGS84高程基准的潮位归算到1985国家高程基准的潮位流程图；

图5为一种基于GNSS的在航潮位测量方法的总流程图；

图6为求取基于WGS84高程基准的海底点高程示意图；

图7为计算海底格网点高程异常的示意图；

图8为计算测船处的高程异常的示意图；

图9海面的瞬时大地高；

图10为本发明实施例中基于变分模态分解和小波变换联合消噪法提取海面高程序列中的潮位图。

附图标记包括：1-GNSS接收机，2-接收机支撑杆，3-接收机支撑横梁，4-接收机安装板，5-第一立柱，6-立柱安装板，7-第二立柱，8-第三立柱，9-IMU，10-多波束发射换能器，11-多波束接收换能器，12-多波束安装板，13-压力测深传感器。

具体实施方式

下面结合具体实施例对本发明的具体实施方式做进一步说明：

一种在航GNSS潮位测量方法，包括：

S1.安装测量装置于船舷，所述测量装置包括多波束声呐、IMU9、压力测深传感器13和两台GNSS接收机1，IMU9为惯性测量单元；IMU9:Inertial…Measurement…Unit，惯性测量单元。

GNSS接收机1独立接收卫星信号，获取GNSS接收机1相位中心的经纬度和大地高程，声呐测量测区的水下地形，获得水下测量点的经纬度、大地高数据和深度，IMU9获取IMU9处的升沉值和测船的姿态值，所述姿态值包括横摇值、纵摇值和艏摇值，所述压力测深传感器13测量测船运动时的动态吃水状态下IMU9到水面的距离和船抛锚时IMU9到水面的距离；

S2.综合S1获得的数据计算WGS84高程基准下的海面大地高程序列；

S3.采用采用变分模态分解和小波变换联合消噪法提取海面大地高程序列中基于大地高程的潮位序列，如图10；

S4.求解测量区基于WGS84高程基准的海底格网点的大地高程；

S5.获得测量区域附近6个或6个以上的GPS水准联测点的高程信息，所述高程信息同时包括GPS水准联测点的大地高程和1985国家高程基准下的高程；

S6.采用二次多项式曲面拟合法求解测量区域海底格网点的近似正常高；

S7.根据莫洛金斯基原理求解海底格网点的高程异常；

S8.根据S1中的GNSS接收机1相位中心的经纬度，采用双线性内插法求解出测船处的精确高程异常，将S3中的潮位序列归算到1985国家高程基准下，得到基于1985国家高程基准的潮位序列。

S1包括以下子步骤：

S1.1.以IMU9为参考原点建立理想船体坐标系，在制图软件中建好测量装置的三维模型，通过三维模型量取测量装置上各个部件之间的偏移关系，在理想船体坐标系下表达两台GNSS接收机1的坐标，通过三维模型量取参考原点到GNSS接收机1的垂直距离并记为L1，再通过压力测深传感器13获得船体抛锚时IMU9到水面的距离L2，即可获得GNSS接收机1到水面的距离L3，计算公式如下：L3＝L1-L2；

S1.2.声呐测量时，在典型区域采用“斑片测试”获得换能器的横摇安装偏差角和纵摇安装偏差角；

S1.3.根据S1.2获得的横摇安装偏差角和纵摇安装偏差角对IMU9测得的横摇值、纵摇值进行校正，若声呐***有横摇补偿功能，则只对IMU9测得的纵摇值进行校正；

S1.4.用GNSS-PPK、GNSS-RTK、网络RTK或PPP的方式获得两台GNSS接收机1相位中心精确的大地高和经纬度数据。

S2包括以下子步骤：

S2.1.根据S1.1中两台GNSS接收机1在理想船体坐标系中的坐标和S1.3中校正后的纵摇值、横摇值，分别计算两台GNSS接收机1的诱导升沉值H_i：H_i＝xsin P-ysin Rcos P-z(cos Rcos P-1)，H_i表示诱导升沉值，x、y、z表示GNSS接收机1在理想船体坐标系中的坐标，P、R表示S1中IMU9记录的纵摇值和横摇值；

S2.2.根据S2.1中的诱导升沉值H_i和S1中IMU9处的升沉值，分别计算两台GNSS接收机1处的升沉值h_G：h_G＝H_i+h₀，h₀表示S1中IMU9处的升沉值；

S2.3.海面大地高程序列的计算过程为：h＝h_GPS-h_G-H_m+ΔH，h为海面的瞬时大地高，即海面大地高程序列；h_GPS为S1中GNSS接收机1获得的相位中心的大地高程；H_m为L3；ΔH为S1中压力测深传感器13测量测船运动时的动态吃水与S1.1的L2之间的差值。

S3包括以下子步骤：

S3.1.采用变分模态分解法对S2得到的海面大地高序列逐次进行K层模态分解，K＝10n，n＝1，2，3……n，n为正整数，惩罚系数α＝5000，收敛容差tol＝1e-7；

S3.2.计算每一次对海面大地高序列进行K层模态分解得到的信号能量E：

n为采样点数，IMF_K(t)为每次模态分解的第K层信号序列；

S3.3.求解能量比值：

当能量比值γ趋于稳定时，γ在0值上下浮动，说明K为合适的模态分解数；

S3.4.对模态分解数为K时分解得到的第一层模态IMF₁再次进行变分模态分解，惩罚系数α＝80000，收敛容差tol＝1e-7，称为迭代分解1次，迭代分解1次的K值判断准则参考S3.2和S3.3，迭代分解1次得到的第一层模态记为IMF₁₁，第二层模态记为IMF₂₂，将IMF₁₁和IMF₂₂重构得到一新信号；

S3.5.对S3.4中的对IMF₁₁和IMF₂₂进行频谱分析，判断IMF₁₁和IMF₂₂的信号频谱组成；

S3.6.S2中的海面大地高程序列对应的最高频率ω为：

f_s为GNSS信号的采样频率；

假设海面高程序列的其他要素的频率为f_w，潮汐的频率为f_t，使用小波层数的层数z由下式决定:

将S3.4所述的IMF₁₁和IMF₂₂重构后得到的信号采用sym4小波分解z层后得到第z层低频系数A_z和高频系数D₁D₂………D_z，低频系数占有的频段为潮汐所在的频段，高频系数占有的频段为其他高频要素所在的频段，将高频系数置零之后与低频系数重构，即可实现海面高程序列的强制消噪，提取出的低频信号为基于WGS84高程基准的潮位序列h_{tide_84}。

S4包括以下子步骤：

S4.1.将声呐测量的水深数据进行安装偏差角校正、声速校正和粗差剔除；

S4.2.将S4.1处理后的水深数据采用S3中获取的h_{tide_84}进行潮位改正：

H_i＝|h_i-h_{tide_84}|，H_i为海底测量点的高程，h_i为S4.1处理后的水深数据；

S4.3.将S2中的海面大地高程序列按照500m*500m格网化，输出WGS84大地坐标系下的格网点的大地坐标，或者平面坐标和大地高h_Ω，将格网点各个的高程值加权取平均，得到的WGS84高程基准下的平均大地高，记为h_γ。

S6包括以下子步骤：

S6.1.在水准联测点上的高程异常与平面坐标之间，假定存在如下二次多项式曲面拟合数学模型：

a₀、a₁、a₂、a₃、a₄、a₅为模型待定参数，ζ_i为水准联测点的高程异常，当水准联测点多于6个时,组成形如下误差方程：V＝AX-ζ⁰，根据最小二乘原理可得A＝(X^TX)^-1X^Tζ⁰

按最小二乘原理解求出模型待定参数a₀、a₁、a₂、a₃、a₄、a₅的数值；

S6.2.将S4中的海底格网点的坐标代入到S6.1的二次多项式曲面拟合数学模型插值出海底格网点的近似高程异常ζ₀；

S6.3.测量区域海底格网点的近似正常高h_Γ为：h_Γ＝ζ₀+H_i，ζ₀表示近似高程异常，H_i为S4中的海底格网点的大地高程。

S7包括以下子步骤：

S7.1.获得S6得到的h_Γ；

S7.2.海底格网点的高程异常ζ由长波部分即近似高程异常ζ₀和短波部分ζ_r得到：ζ＝ζ₀+ζ_r；

S7.3.计算高程异常中的短波部分ζ_r：ζ_r＝T_r/r，T_r为地形起伏对海底格网点扰动位的影响；ζ_r为地形改正；

G为万有引力常数，G＝6.673×10^-8C³s^-2g^-1；ρ为地球平均质量密度，ρ＝2.67gC^-3，γ₀＝[(x-x_p)²+(y-y_p)²]^1/2,(x,y)为S4中的海底格网点的坐标，(x_p,y_p)为海底格网点中某个格网点P的坐标；

S7.4.似大地水准面的正常重力r：r＝r₀-0.3086h_Γ，r₀为参考椭球面上的正常重力值，位：10^-5…m/s²，

为WGS84大地坐标系中测点地理纬度，单位：度。

S8包括以下子步骤：

S8.1.根据S1中GNSS接收机1相位中心的经纬度、S4中的海底格网点的大地高程、S7中的ζ，采用双线性插值法求取测船位置处的高程异常ζ₁；

S8.2.计算最终内插出的高程异常值f(x,y)：

f(x,y₁)、f(x,y₂)表示x方向插值出的高程异常值，ζ₁₁、ζ₂₁、ζ₁₂、ζ₂₂表示格网四个角点上的高程异常值，x₁、x₂、y₁、y₂表示四个角点在X方向和Y方向上的坐标，x、y表示测船在X方向和Y方向上的坐标；

S8.3.计算基于1985国家高程基准的潮位序列h_{tide_85}：h_{tide_85}＝h_{tide_84}+ζ₁。

测量装置如图1所示，通过固定在立柱的两个安装板焊接在测船船舷，测量装置由两台GNSS接收机1、接收机支撑杆2、接收机支撑横梁3、接收机安装板4、第一立柱5、立柱安装板6、第二立柱7、第三立柱8、多波束发射换能器10、多波束接收换能器11、IMU9、压力测深传感器13、多波束安装板12组成。

a.GNSS接收机1接收卫星发射的信号，采用GNSS-PPK、GNSS-RTK、网络RTK或

PPP的方式获取高精度的经纬度和大地高数据，采样频率为1…HZ；

b.多波束测深声呐由发射换能器和接收换能器组成，以“条带式”覆盖测量的方式获取近海高精度的海底地形；

c.IMU9与换能器水下安装，IMU9提供换能器的姿态值、升沉值，以IMU9为参考原点O建立船体坐标系，X轴水平向前指向船首方向(测量装置P安装时平行于船首方向)，Y轴过参考原点水平向右，Z轴垂直于XOY平面垂直向下；

d.压力测深传感器13的采样频率为5…HZ，在船抛锚时提供IMU9参考原点到水面的距离值，为换能器测得的水深做改正；在船运动时提供测船的动态吃水，用以计算海面的瞬时大地高；

图4为本方法的总流程图，图2为变分模态分解和小波变换联合消噪法法从海面的大地高中提取潮位(WGS84高程基准)的流程图，图9为其结果图，针对测量区域的潮汐特点，判断小波分解的层数z。以青岛半日潮海域为例，潮汐截止频率约为1/(3600*3)≈9.26*10^{- 5}Hz，高频噪声的变化周期最大为几秒至十几分钟，这里取噪声的最低频率为1/1200≈8.33*10^-4Hz，为了将潮汐信号与噪声分离，当小波分解的层数z最小取为10时0.5/2¹⁰≈4.88*10^-4Hz，则低频系数所占频段为0～4.88*10^-4Hz，高频系数所占频段为4.88*10^-4Hz～0.5Hz，潮汐频率9.26*10^-5Hz在此频段之内，同样z＝11、z＝12时，潮汐频率也在分解出的低频系数所占频段之内，并能将潮汐频率与噪声频率分开，对于青岛半日潮地区理论上10≤z≤12能将潮汐数据与噪声频率分开，再通过将高频系数置零与低频系数重构，就能将基于WGS84高程基准的潮位序列提取出来，提取出地潮位序列如图8所示。

图3为基于WGS84高程基准的潮位归算到1985国家高程基准的流程图，具体包括：

R1.通过多波束声呐的“条带式”覆盖测量的方式获取测区的水深数据；

R2.对步骤R1所述的水深数据进行安装偏差角校正、声速校正和粗差剔除，此步可在多波束数据处理软件如caris中完成；

R3采用提取出的基于WGS84高程基准的潮位对R2所述的校正之后的水深数据进行潮位改正，示意图如图5所示。

如图6，E表示测船的位置，某时刻E处于格网PBCD中，P、B、C、D格网点的高程异常在步骤k已求出；如图7，根据P、B、C、D格网点的高程异常采用双线性插值算法求出测船处的高程异常；

二次多项式曲面拟合法也可以替换为以下方法：多项式曲线拟合法、三次样条曲线拟合法、三点拟合法、分区最小二乘平面拟合法、三次多项式曲面拟合法、Hardy多面函数拟合法、薄板小挠度变形模型拟合法、离散点加权平均拟合法、移动的二次多项式拟合法。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (8)

1.一种在航GNSS潮位测量方法，其特征在于，包括：

S1.安装测量装置于船舷，所述测量装置包括多波束声呐、IMU、压力测深传感器和两台GNSS接收机，IMU为惯性测量单元；

GNSS接收机独立接收卫星信号，获取GNSS接收机相位中心的经纬度和大地高程，声呐测量测区的水下地形，获得水下测量点的经纬度、大地高数据和深度，IMU获取IMU处的升沉值和测船的姿态值，所述姿态值包括横摇值、纵摇值和艏摇值，所述压力测深传感器测量测船运动时的动态吃水状态下IMU到水面的距离和船抛锚时IMU到水面的距离；

S2.综合S1获得的数据计算WGS84高程基准下的海面大地高程序列；

S3.采用变分模态分解和小波变换联合消噪法提取WGS84高程基准下海面大地高程序列中的潮位序列；

S4.求解测量区WGS84高程基准下的海底格网点的大地高程；

S5.获得测量区域附近6个或6个以上的GPS水准联测点的高程信息，所述高程信息同时包括GPS水准联测点的WGS84高程基准下的大地高程和1985国家高程基准下的高程；

S6.采用二次多项式曲面拟合法求解测量区域海底格网点的近似正常高；

S7.根据莫洛金斯基原理求解海底格网点的高程异常；

S8.根据S1中的GNSS接收机相位中心的经纬度，采用双线性内插法求解出测船处的精确高程异常，将S3中的潮位序列归算到1985国家高程基准下，得到基于1985国家高程基准的潮位序列。

2.如权利要求1所述的一种在航GNSS潮位测量方法，其特征在于，S1包括以下子步骤：

S1.1.以IMU为参考原点建立理想船体坐标系，在制图软件中建好测量装置的三维模型，通过三维模型量取测量装置上各个部件之间的偏移关系，在理想船体坐标系下表达两台GNSS接收机的坐标，通过三维模型量取参考原点到GNSS接收机的垂直距离并记为L1，再通过压力测深传感器获得船体抛锚时IMU到水面的距离L2，即可获得GNSS接收机到水面的距离L3，计算公式如下：L3＝L1-L2；

S1.2.声呐测量时，在典型区域采用“斑片测试”获得换能器的横摇安装偏差角和纵摇安装偏差角；

S1.3.根据S1.2获得的横摇安装偏差角和纵摇安装偏差角对IMU测得的横摇值、纵摇值进行校正，若声呐***有横摇补偿功能，则只对IMU测得的纵摇值进行校正；

S1.4用GNSS-PPK、GNSS-RTK、网络RTK或PPP的方式获得两台GNSS接收机相位中心精确的大地高和经纬度数据。

3.如权利要求1所述的一种在航GNSS潮位测量方法，其特征在于，S2包括以下子步骤：

S2.1.根据S1.1中两台GNSS接收机在理想船体坐标系中的坐标和S1.3中校正后的纵摇值、横摇值，分别计算两台GNSS接收机的诱导升沉值H_i：H_i＝xsin P-ysin Rcos P-z(cosRcos P-1)，H_i表示诱导升沉值，x、y、z表示GNSS接收机在理想船体坐标系中的坐标，P、R表示S1中IMU记录的纵摇值和横摇值；

S2.2.根据S2.1中的诱导升沉值H_i和S1中IMU处的升沉值，分别计算两台GNSS接收机处的升沉值h_G：h_G＝H_i+h₀，h₀表示S1中IMU处的升沉值；

S2.3.海面大地高程序列的计算过程为：h＝h_GPS-h_G-H_m+ΔH，h为海面的瞬时大地高，即海面大地高程序列；h_GPS为S1中GNSS接收机获得的相位中心的大地高程；H_m为L3；ΔH为S1中压力测深传感器测量测船运动时的动态吃水与S1.1中的L3之间的差值。

4.如权利要求3所述的一种在航GNSS潮位测量方法，其特征在于，S3包括以下子步骤：

S3.1.采用变分模态分解法对S2得到的海面大地高序列逐次进行K层模态分解，K＝10n，n＝1，2，3......n，n为正整数，惩罚系数α＝80000，收敛容差tol＝1e-7；

S3.2.计算每一次对海面大地高序列进行K层模态分解得到的信号能量E：

n为采样点数，IMF_K(t)为每次模态分解的第K层信号序列；

S3.3.求解能量比值：

当能量比值γ趋于稳定时，γ在0值上下浮动，说明K为合适的模态分解数；

S3.4.对模态分解数为K时分解得到的第一层模态IMF₁再次进行变分模态分解，称为迭代分解1次，迭代分解1次的K值判断准则参考S3.2和S3.3，迭代分解1次得到的第一层模态记为IMF₁₁，第二层模态记为IMF₂₂，将IMF₁₁和IMF₂₂重构得到一新信号；

S3.5.对S3.4中的对IMF₁₁进行频谱分析，判断IMF₁₁的信号频谱组成；

S3.6.S2中的海面大地高程序列对应的最高频率ω为：

f_s为GNSS信号的采样频率；

假设海面高程序列的其他要素的频率为f_w，潮汐的频率为f_t，使用小波层数的层数z由下式决定：

将S3.4的IMF₁₁和IMF₂₂重构得到的新信号采用sym4小波分解z层后得到第z层低频系数A_z和高频系数D₁、D₂......D_z，低频系数占有的频段为潮汐所在的频段，高频系数占有的频段为其他高频要素所在的频段，将高频系数置零之后与低频系数重构，即可实现海面高程序列的强制消噪，提取出的低频信号为基于WGS84高程基准的潮位序列h_{tide_84}。

5.如权利要求4所述的一种在航GNSS潮位测量方法，其特征在于，S4包括以下子步骤：

S4.1.将声呐测量的水深数据进行安装偏差角校正、声速校正和粗差剔除；

S4.2.将S4.1处理后的水深数据采用S3中获取的h_{tide_84}进行潮位改正：

H_i＝|h_i-h_{tide_84}|，H_i为海底测量点的高程，h_i为S4.1处理后的水深数据；

S4.3.将S2中的海面大地高程序列按照500m*500m格网化，输出WGS84大地坐标系下的格网点的大地坐标，或者平面坐标和大地高h_Ω，将格网点各个的高程值加权取平均，得到的WGS84高程基准下的平均大地高，记为h_γ。

6.如权利要求5所述的一种在航GNSS潮位测量方法，其特征在于，S6包括以下子步骤：

S6.1.在水准联测点上的高程异常与平面坐标之间，假定存在如下二次多项式曲面拟合数学模型：

a₀、a₁、a₂、a₃、a₄、a₅为模型待定参数，ζ_i为水准联测点的高程异常，当水准联测点多于6个时，组成形如下误差方程：V＝AX-ζ⁰，根据最小二乘原理可得A＝(X^TX)^-1X^Tζ⁰

按最小二乘原理解求出模型待定参数a₀、a₁、a₂、a₃、a₄、a₅的数值；

S6.2.将S4中的海底格网点的坐标代入到S6.1的二次多项式曲面拟合数学模型插值出海底格网点的近似高程异常ζ₀；

S6.3.测量区域海底格网点的近似正常高h_Γ为：h_Γ＝ζ₀+H_i，ζ₀表示近似高程异常，H_i为S4中的海底格网点的大地高程。

7.如权利要求6所述的一种在航GNSS潮位测量方法，其特征在于，S7包括以下子步骤：

S7.1.获得S6得到的h_Γ；

S7.2.海底格网点的高程异常ζ由长波部分即近似高程异常ζ₀和短波部分ζ_r得到：ζ＝ζ₀+ζ_r；

S7.3.计算高程异常中的短波部分ζ_r：ζ_r＝T_r/r，T_r为地形起伏对海底格网点扰动位的影响；ζ_r为地形改正；

G为万有引力常数，G＝6.673×10^-8C³s^-2g^-1；ρ为地球平均质量密度，ρ＝2.67gC^-3，γ₀＝[(x-x_p)²+(y-y_p)²]^1/2，(x，y)为S4中的海底格网点的坐标，(x_p，y_p)为海底格网点中某个格网点P的坐标；

S7.4.似大地水准面的正常重力r：r＝r₀-0.3086h_Γ，r₀为参考椭球面上的正常重力值，单位：10^-5m/s²，

为WGS84大地坐标系中测点地理纬度，单位：度。

8.如权利要求7所述的一种在航GNSS潮位测量方法，其特征在于，S8包括以下子步骤：

S8.1.根据S1中GNSS接收机相位中心的经纬度、S4中的海底格网点的大地高程、S7中的ζ，采用双线性插值法求取测船位置处的高程异常ζ₁；

S8.2.计算最终内插出的高程异常值f(x，y)：

f(x，y₁)、f(x，y₂)表示x方向插值出的高程异常值，ζ₁₁、ζ₂₁、ζ₁₂、ζ₂₂表示格网四个角点上的高程异常值，x₁、x₂、y₁、y₂表示四个角点在X方向和Y方向上的坐标，x、y表示测船在X方向和Y方向上的坐标；

S8.3.计算基于1985国家高程基准的潮位序列h_{tide_85}：h_{tide_85}＝h_{tide_84}+ζ₁。

Images