具体实施方式

：

为使本发明的目的、技术方案和优点更加清楚、明白，下面结合附图和技术方案对本发明作进一步详细的说明。

本发明实施例，参见图1所示，提供一种水下动态重力数据采集半实物仿真方法，用于在实验室环境下模拟水下重力辅助惯性导航航迹，包含如下内容：

S101、依据初始设定参数获取水下潜器预设航迹，并依据预设航迹确定其在每个采样点上的航迹数据，该航迹数据包含：位置、速度和姿态信息；

S102、根据航迹数据通过计算角增量和速度增量来模拟水下动态重力测量数据，该水下动态重力测量数据包含惯性角增量和速度增量及比力；

S103、针对水下动态重力数据，通过扣除有害加速度以及该点正常重力影响来获取水下采样点重力异常值；依据重力异常值及各采样点实际重力异常值的差值进行统计分析，评估半实物仿真精度。

通过对水下动态重力数据采集半实物仿真来实现对潜器沿水下航迹的重力数据测量采集过程的物理呈现，可为进行水下重力辅助惯性导航关键技术的深入研究提供与实际测量环境相似的可靠便捷的实验条件，从而节省大量的人力、物力和时间成本，并对仿真进行评估，提升模拟效果和精确度。

进一步地，本发明实施例中，通过设定水下潜器初始姿态、速度、位置以及运动加速度和角速度参数，基于航迹推算方法获取潜器的预设航迹；并按照惯导数据采样率，从预设航迹确定潜器在每个采样点上的姿态、速度、位置。

记当地水平坐标系为n系，载体坐标系为b系，俯仰、横滚、方位欧拉角向量为 A＝[θγψ]^T、欧拉角速率向量ω＝[ω_θω_γω_ψ]^T、b系加速度a^b＝[0a_y 0]^T，则水下航行运动航迹设置满足如下微分方程组：

其中

假设初值A(t₀)＝[θ₀γ₀ψ₀]^T，p(t₀)＝[B₀L₀ H₀]^T，R_M为子午圈曲率半径，R_N为卯酉圈曲率半径，B为大地纬度，H为大地高，为由当前时刻欧拉角实时计算求得的姿态矩阵。

通过分段设置机动动作，即航迹分段输入参数ω_θ、ω_γ、ω_ψ和a_y中的一个或多个，采用一阶欧拉法，依据惯导采样率设定步长参数求解(1)式时变微分方程，解得航迹参数A，vⁿ和 p，完成水下航迹预设。

惯性角增量与速度增量输出仿真，通过捷联惯导反演算法原理，求取角增量与速度增量信息，模拟惯性系统的测量输出。

记t_m时刻的姿态阵为t_m-1时刻的n系到t_m时刻的n系的转换阵为时刻的导航系相对惯性系的转动角速度在导航系下的投影为t_m-1时刻到t_m时刻的b系转动的等效旋转矢量时间段T＝t_m-t_m-1内比力引起的速度增量为时间段 T＝t_m-t_m-1内有害加速度引起的速度增量为时刻地球自转角速度在n系下的投影为时刻n系相对e系的角速度在n系下的投影为时刻载体所在位置处的重力值为完成潜器航迹设计后，航迹上各采样点的姿态、速度和位置参数均为已知，因此和等量均已知或可计算，由这些量即可求解惯性传感器的增量采样信息。

首先，角增量计算公式为：

令初值Δθ₀＝0。

速度增量计算公式为：

式中

令初值Δv₀＝0。依据(2)、(3)式完成角增量Δθ_m和速度增量Δv_m的仿真，生成惯性输出角增量与速度增量。

水下动态重力测量比力输出仿真，顾及海水层质量影响，将海洋重力异常延拓至水下采样点航迹处；以该采样点实际重力值替代正常重力值代入捷联惯导反演算法获取角增量与速度增量信息，然后确定比力，该比力值可视为水下动态重力数据测量输出。

假定海水层的密度为δ₀，则海水层质量影响为

T_z＝2×2πGδ₀h (6)

式中G是万有引力常数，h为潜器潜水深度。扣除海水层质量影响后，利用地球重力场模型 (C_nm,S_nm)计算重力异常Δg的垂直梯度变化，从而将海洋面重力异常向下延拓至潜器潜水深度。

式中，M为地球总质量，R表示地球平均半径，系地心纬度、经度，表示缔合勒让德函数。

依据WGS84正常重力公式计算采样点正常重力值：

将延拓后的海洋实际重力值代入(5)式，然后由(3)式计算速度增量，并对时间求微分可得比力：

注意式中Δv_m由实际重力值代入(4)式求得。在此基础上加入重力测量白噪声，生成水下动态重力数据输出比力值。

惯性测量数据解算。由生成的惯性测量数据输出仿真，依据捷联惯导算法原理，确定潜器的姿态、速度、位置。

捷联惯导更新算法是一种递推算法，即根据前一时刻潜器的姿态、速度、位置，以及更新周期内加速度计与陀螺仪的采样输出，推算后一时刻潜器的姿态、速度、位置。

姿态更新，由矩阵链乘法则：

式中，和分别表示t_m-1和t_m时刻的姿态矩阵。

陀螺在[t_m-1,t_m]，[t_m-2,t_m-1]时间段内进行角增量等间隔采样为和Δθ_m-1，令T＝t_m-t_m-1，则有Δθ_m

(11)

式(10)～式(13)即为捷联惯导数值递推姿态更新原理。

速度更新中，

其中，和分别为t_m-1和t_m时刻的速度。有害加速度引起的速度增量计算按 (5)式，比力积分速度增量计算公式为：

式中Δθ_m和Δθ_m-1分别是在时间段[t_m-1,t_m]，[t_m-2,t_m-1]内角增量输出，Δv_m和Δv_m-1分别是在时间段[t_m-1,t_m]，[t_m-2,t_m-1]内速度增量输出。为t_m-1时刻的姿态矩阵。

位置更新，采用梯形积分法对式(1)的位置微分式进行离散化，得位置更新算法：

式中，M_pv(m-1/2)可采用线性外推，对矩阵整体M_pv进行外推；也可对矩阵元素中的位置变量B，H外推，再构造矩阵M_pv。

水下动态重力测量数据处理，由惯性输出仿真通过惯性测量数据解算实时提供潜器位置、姿态、速度及加速度，从水下动态重力测量数据输出中扣除科里奥利加速度与运动加速度影响，生成水下采样点的重力异常值。

n系下的水下动态重力测量矢量方程描述为：

式中，gⁿ为重力加速度。为每个采样点上动态重力测量数据的比力输出值；为科里奥利加速度。

本设计侧重于水下标量重力测量，为避免姿态误差的传递累积，根据重力矢量模不变量原理，将基于捷联惯导的动态重力测量方程改写为：

式中，

将重力仿真生成值经数据处理后获取的重力异常值与各采样点实际重力异常值的差值进行统计分析，评估水下动态重力数据采集半实物仿真系统的仿真精度。

进一步，基于上述的方法，本发明实施例还提供一种水下动态重力数据采集半实物仿真系统，用于在实验室环境下模拟水下重力辅助惯性导航航迹上的重力测量数据的实时采集，，包含：预设模块、模拟模块和评估模块，其中，

预设模块，用于依据初始设定参数获取水下潜器预设航迹，并依据预设航迹确定其在每个采样点上的航迹数据，该航迹数据包含：位置、速度和姿态信息；

模拟模块，用于根据航迹数据通过计算角增量和速度增量来模拟水下动态重力数据，该水下动态重力测量数据包含惯性角增量和速度增量及比力；

评估模块，用于针对水下动态测量重力数据，通过扣除有害加速度以及对应采样点处正常重力影响来获取水下采样点重力异常值；依据重力异常值及各采样点实际重力异常值的差值进行统计分析，评估半实物仿真精度。

为验证本案方案有效性，下面结合仿真数据做进一步解释说明：

设置水下航迹参数如下，航迹起始点经纬度pos0＝(110.9°7.5°)；潜水深度100m；初始姿态att0＝(0°0°0°)；初始速度v0＝(0 0 0)预设；机动动作以0.0172m/s²的加速度，匀加速300s，速度达到10节后，匀速直线运动1h，然后左转90°，匀速运动1h，然后右转90°，匀速运动1h，航迹结束。整个仿真时长为3h。在上述航迹参数设置的基础上，采样率设置为100HZ，得到图3～7所示的各采样点的姿态、速度、位置，如图3a到图7a分别所示的航迹归算俯仰角、横滚角、航向角、推算速度和推算位置；利用惯导反演算法生成惯性输出的角增量与速度增量，在不顾及陀螺以及加速度计的常值漂移和随机游走误差的影响下进行惯性导航解算，得到潜器的姿态、速度、位置，如图3b～图7b分别所示的无误差惯导解算俯仰角、横滚角、航向角、推算速度和推算位置。理论上解算得到的各采样点的姿态、速度、位置应与设计航迹上各采样点的姿态、速度、位置完全相等，但受模型误差和计算机字长的限制，两者之间存在有微小差异，如果姿态角误差在角秒级、速度误差在0.01m/s以下，则认为仿真得到的惯性输出信息符合精度要求。各采样点的姿态、速度、位置误差如图3c～图7c分别所示俯仰角误差、横滚角误差、航向角误差、推算速度误差和推算位置误差，误差统计列于表1到表3。

表1姿态误差统计信息

表2速度误差统计信息 表3位置误差统计信息

由表1-表3可知，在仿真时长3h时间段内，姿态角的最大误差为0.29″，速度误差的最大值为0.002m/s，位置误差最大值为8.5m。由于潜器位置误差通常以海里为单位，位置误差在量级上完全能够满足水下导航对数据源仿真精度的要求，同时也验证了仿真惯性元器件信息的正确性。

将水下航迹生成采样点的实际重力值、姿态、速度、位置代入捷联惯导反演算法，确定各采样点处重力仪的三个加速度计的输出仿真比力值，如图8所示，(a)、(b)、(c)分别表示重力仪x、y、z方向加速度计比力输出。获得重力仪的实际测量值-比力后，截除潜器机动转向时间段的测量数据，然后利用惯导解算得到的各采样点的速度、位置、对重力测量值进行科里奥利加速度与运动加速度改正，即可确定该点处的重力异常。理论上，仿真测量获得的重力异常值应等于设计航迹上各采样点的实际重力异常值，但由于惯导仿真误差与计算误差的影响，仿真测量重力异常值与实际海洋重力异常值存在差异，其差异值如图9～11所示，误差统计列于表4，最大误差为0.026mGal，均方误差为0.0096mGal，实验结果说明所设计的水下动态重力数据采集半实物仿真系统符合捷联式动态重力测量原理，可满足实验室水下重力测量的仿真要求。

表4重力仿真精度评估/单位：mGal

除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对步骤、数字表达式和数值并不限制本发明的范围。

基于上述的系统，本发明实施例还提供一种服务器，包括：一个或多个处理器；存储装置，用于存储一个或多个程序，当所述一个或多个程序被所述一个或多个处理器执行，使得所述一个或多个处理器实现上述的系统。

基于上述的系统，本发明实施例还提供一种计算机可读介质，其上存储有计算机程序，其中，该程序被处理器执行时实现上述的系统。

本发明实施例所提供的装置，其实现原理及产生的技术效果和前述系统实施例相同，为简要描述，装置实施例部分未提及之处，可参考前述系统实施例中相应内容。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统和装置的具体工作过程，可以参考前述系统实施例中的对应过程，在此不再赘述。

在这里示出和描述的所有示例中，任何具体值应被解释为仅仅是示例性的，而不是作为限制，因此，示例性实施例的其他示例可以具有不同的值。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

附图中的流程图和框图显示了根据本发明的多个实施例的系统、系统和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/ 或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

在本申请所提供的几个实施例中，应该理解到，所揭示的系统、装置和系统，可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。

所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个处理器可执行的非易失的计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备 (可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述系统的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭示的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。