光纤光栅形变测量系统时间同步方法和装置
阅读说明：本技术 光纤光栅形变测量系统时间同步方法和装置 (Time synchronization method and device for fiber bragg grating deformation measurement system ) 是由 宫晓琳 孙一弘 刘刚 房建成 田珂珂 符倚伦 丁孝双 于 2020-06-05 设计创作，主要内容包括：本公开提供了机载分布式POS用光纤光栅形变测量系统时间同步方法,基于分布式POS测量的挠曲角和光纤光栅形变测量系统测量的柔性基线形变量的内在相关性,对两者数据基于密度空间的聚类方法的去噪、高频插值、经验模态分解法获取振动分量、小波变换确定振动分量的实时振动频率及幅值等过程,并通过波形匹配确定光纤光栅形变测量系统数据的实际采样频率及采样常值延迟,确定光纤光栅形变测量系统高频波长变化量数据的更加准确的时间标签,获取分布式POS传递对准时刻所需的光纤光栅传感器波长变化量数据。该方法提高光纤光栅传感器数据的采样时间精度,提高其与分布式POS的时间同步精度。本公开还提出机载分布式POS用光纤光栅形变测量系统时间同步装置。(The invention provides a time synchronization method of a fiber grating deformation measurement system for an airborne distributed POS, which is characterized in that based on the internal correlation of a deflection angle measured by the distributed POS and a flexible baseline deformation measured by the fiber grating deformation measurement system, the processes of denoising, high-frequency interpolation and empirical mode decomposition of the two data based on a density space clustering method are used for obtaining a vibration component, wavelet transformation is used for determining the real-time vibration frequency and amplitude of the vibration component, the actual sampling frequency and sampling constant delay of the fiber grating deformation measurement system data are determined through waveform matching, a more accurate time label of the high-frequency wavelength variation data of the fiber grating deformation measurement system is determined, and the fiber grating sensor wavelength variation data required by the distributed POS transmission alignment moment are obtained. The method improves the sampling time precision of the fiber bragg grating sensor data and improves the time synchronization precision of the fiber bragg grating sensor data and the distributed POS. The disclosure also provides a time synchronization device of the fiber bragg grating deformation measurement system for the airborne distributed POS.)
技术领域
本公开涉及航空航天技术领域，具体而言，涉及光纤光栅形变测量系统时间同步方法和装置，具体涉及机载分布式POS用光纤光栅形变测量系统时间同步方法和装置。
背景技术
目前，多任务遥感载荷是机载对地观测的发展方向之一。如机载分布式阵列天线合成孔径雷达(Synthetic Aperture Radar，SAR)和柔性多基线干涉 SAR等。对于装备多任务遥感载荷的高性能航空遥感系统，需要同时对位于机翼等柔性基线上的所有遥感载荷点的运动参数实现高精度测量。
位置姿态测量系统(Position and Orientation System，POS)是目前机载对地观测遥感载荷获取位置、速度、姿态等运动参数的主要手段。对于装备了多个或多种观测载荷的机载对地观测系统，由于多个或多种观测载荷安装在飞机的不同位置，飞机弹性变形导致载荷间的空间关系随时间复杂变化，必须建立机载分布式时空基准系统，即机载分布式POS。
机载分布式POS一般由一个高精度主POS和多个子惯性测量单元 (InertialMeasurement Unit，IMU)组成。主POS一般安装在载机机舱内；子IMU一般分布安装在机体的不同位置(包括机翼)，受重量、体积等限制，子IMU往往选择中低精度的IMU，需要依靠主POS的高精度位置、速度、姿态等运动参数对其进行传递对准以实现所在处运动信息的精确测量。由于飞机机翼和机体存在弹性变形，主、子系统间的杆臂柔性时变，因此确定各载荷之间瞬息万变的相对运动关系成为通过传递对准获得子节点高精度运动参数的前提。光纤光栅传感器是形变测量的一个常用手段，广泛应用于结构形变测量及故障诊断中。光纤光栅形变测量系统主要包括光纤光栅传感器和解调仪两部分。
由于机翼等柔性基线的形变幅值可达数厘米甚至数十厘米且频率可高达数十赫兹，1毫秒的时间同步误差就可导致毫米级的形变测量误差。因此，光纤光栅传感器用于测量分布式POS传递对准所需的柔性基线运动信息时，首先要进行二者的高精度时间同步。然而，目前光纤光栅形变测量主要用于更新频率较低的场合，对于与其他设备高频协同同步工作以及高频时间同步相关研究尚未见报道。目前光纤光栅形变测量系统无法实现高频时间同步的问题主要包括系统存在数据时间延迟和实际采样频率不等于设定值两方面。光纤光栅形变测量系统数据采集存在时间延迟的原因为：分布式POS 高精度传递对准需要精确测量各节点间的结构形变，而要获取多个方向的变形信息就需要在分布式POS分布的载机结构上布设多条光纤光栅传感器。出于成本等因素考虑，光纤光栅解调仪通过光开关切换在每个采样区间逐次读取各条传感器数据，从而导致不同光纤光栅通道数据采集的常值延迟。除此之外，传感器数据采集本身存在的采集设备控制过程和传输延迟最终也会叠加到数据采集时间上。而光纤光栅形变测量系统实际采样频率不等于设定值的原因为：光纤光栅解调仪数据采集系统存在时钟误差(计算机时钟误差或晶振误差等)。上述数据时间延迟和实际采样频率不等于设定值这两方面问题导致光纤光栅形变测量系统数据与分布式POS系统数据时间不同步，无法精确地获得分布式POS系统传递对准所需时刻的各节点相对运动数据。
发明内容
为了解决现有技术中的技术问题，本公开实施例提供了机载分布式POS 用光纤光栅形变测量系统时间同步方法和装置，能够基于波形识别和信号处理等技术，可将光纤光栅传感器数据与分布式POS数据进行时间同步，从而通过光纤光栅传感器获取分布式POS传递对准所需时刻的柔性基线运动信息。
第一方面，本公开实施例提供了机载分布式POS用光纤光栅形变测量系统时间同步方法，所述方法包括：判断光纤光栅解调仪存储的每个光纤光栅传感器通道测量的波长变化量数据和分布式POS中任意一个安装于子节点处的IMU的数据是否丢包，并使用三次样条插值法获取二者丢包时刻的数据，得到完整数据，分别记为各通道波长变化量数据和子节点IMU数据；通过对子节点IMU数据的惯性导航解算，计算出该IMU所在节点的挠曲角数据，记为挠曲角数据；使用三次样条插值将挠曲角数据和各通道波长变化量数据分别插值为2000Hz的高频数据；使用基于密度空间的聚类方法去除高频挠曲角数据和各通道高频波长变化量数据中的异常值，去除数据中的部分噪声；使用经验模态分解法分别去掉去噪后的高频挠曲角数据和各通道高频波长变化量数据中的缓变分量，得到二者的振动分量，记为高频挠曲角数据振动分量和各通道高频波长变化量振动分量；分别提取高频挠曲角振动分量和各通道高频波长变化量这两组数据在全程测量过程中各自测量数据的第一段和最后一段振动数据，并分别对这两组数据的两段振动数据的幅值最大值进行归一化；通过波形匹配确定光纤光栅形变测量系统的采样常值延迟及实际采样频率；确定光纤光栅形变测量系统高频波长变化量数据的更加准确的时间标签，从而实现光纤光栅形变测量系统和子节点IMU的时间同步。
第二方面，本公开实施例提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现上述的方法的步骤。
第三方面，本公开实施例提供了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现上述的方法的步骤。
第四方面，本公开实施例提供了机载分布式POS用光纤光栅形变测量系统时间同步装置，所述装置包括：判断模块，用于判断光纤光栅解调仪存储的每个光纤光栅传感器通道测量的波长变化量数据和分布式POS中任意一个安装于子节点处的IMU的数据是否丢包，并使用三次样条插值法获取二者丢包时刻的数据，得到完整数据，分别记为各通道波长变化量数据和子节点IMU数据；解算模块，用于通过对子节点IMU数据的惯性导航解算，计算出该IMU所在节点的挠曲角数据，记为挠曲角数据；插值模块，用于使用三次样条插值将挠曲角数据和各通道波长变化量数据分别插值为 2000Hz的高频数据；噪声去除模块，用于使用基于密度空间的聚类方法去除高频挠曲角数据和各通道高频波长变化量数据中的异常值，去除数据中的部分噪声；振动分量获取模块，用于使用经验模态分解法分别去掉去噪后的高频挠曲角数据和各通道高频波长变化量数据中的缓变分量，得到二者的振动分量，记为高频挠曲角数据振动分量和各通道高频波长变化量振动分量；提取模块，用于分别提取高频挠曲角振动分量和各通道高频波长变化量这两组数据在全程测量过程中各自测量数据的第一段和最后一段振动数据，并分别对这两组数据的两段振动数据的幅值最大值进行归一化；确定模块，用于通过波形匹配确定光纤光栅形变测量系统的采样常值延迟及实际采样频率；时间同步模块，用于确定光纤光栅形变测量系统高频波长变化量数据的更加准确的时间标签，从而实现光纤光栅形变测量系统和子节点IMU的时间同步。
本发明提供的机载分布式POS用光纤光栅形变测量系统时间同步方法和装置，利用分布式POS中子节点惯性测量单元导航解算的挠曲角与光纤光栅形变测量系统波长变化量之间的关系通过波形匹配等方式实现光纤光栅形变测量系统与分布式POS的时间同步。首先使用三次样条插值补齐子节点惯性测量单元数据和各通道光纤光栅传感器的波长变化量数据中丢包的数据，之后对该子节点惯性测量单元数据进行导航解算获取该子节点的形变挠曲角数据并将光纤光栅传感器数据和挠曲角数据插值为高频。在对高频挠曲角数据和高频波长变化量数据去噪后使用经验模态分解法分别获取其振动分量以及对应的瞬时振动频率和幅值；最后基于高频挠曲角数据和高频波长变化量数据的振动分量的瞬时振动频率和幅值通过波形匹配确定光纤光栅形变测量系统的采样常值延迟及实际采样频率。由此最终确定光纤光栅形变测量系统高频波长变化量数据的更加准确的时间标签，从而实现时间同步，为获取分布式POS数据采样时刻各节点的高精度相对运动信息提供前提，进而辅助分布式POS提供机载多任务遥感载荷的柔性基线上多个节点的高精度运动信息，辅助机载多任务遥感载荷对地观测获取更高精度的成像信息。
附图说明
为了更清楚地说明本公开实施例的技术方案，下面对实施例描述中所需要使用的附图作简单地介绍：
图1为本发明一个实施例中的机载分布式POS用光纤光栅形变测量系统时间同步方法的步骤流程示意图；
图2为本发明另一实施例中的机载分布式POS用光纤光栅形变测量系统时间同步方法的步骤流程示意图；
图3为本发明一个实施例中的机载分布式POS用光纤光栅形变测量系统时间同步装置的结构示意图；
图4为本发明一个实施例中的机载分布式POS用光纤光栅形变测量系统时间同步装置的硬件框图；
图5为本发明一个实施例中的计算机可读存储介质的示意图。
具体实施方式
下面结合附图和实施例对本申请进行进一步的详细介绍。
在下述介绍中，术语“第一”、“第二”仅为用于描述的目的，而不能理解为指示或暗示相对重要性。下述介绍提供了本公开的多个实施例，不同实施例之间可以替换或者合并组合，因此本申请也可认为包含所记载的相同和/ 或不同实施例的所有可能组合。因而，如果一个实施例包含特征A、B、C，另一个实施例包含特征B、D，那么本申请也应视为包括含有A、B、C、D 的一个或多个所有其他可能的组合的实施例，尽管该实施例可能并未在以下内容中有明确的文字记载。
为了使本发明的目的、技术方案及优点更加清楚明白，以下通过实施例，并结合附图，对本发明机载分布式POS用光纤光栅形变测量系统时间同步方法和装置的具体实施方式进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
如图1所示，为一个实施例中的机载分布式POS用光纤光栅形变测量系统时间同步方法的流程示意图，具体包括以下步骤：
步骤11，判断光纤光栅解调仪存储的每个光纤光栅传感器通道测量的波长变化量数据和分布式POS中任意一个安装于子节点处的IMU的数据是否丢包，并使用三次样条插值法获取二者丢包时刻的数据，得到完整数据，分别记为各通道波长变化量数据和子节点IMU数据。
步骤12，通过对子节点IMU数据的惯性导航解算，计算出该IMU所在节点的挠曲角数据，记为挠曲角数据。
步骤13，使用三次样条插值将挠曲角数据和各通道波长变化量数据分别插值为2000Hz的高频数据。
步骤14，使用基于密度空间的聚类方法去除高频挠曲角数据和各通道高频波长变化量数据中的异常值，去除数据中的部分噪声。
步骤15，使用经验模态分解法分别去掉去噪后的高频挠曲角数据和各通道高频波长变化量数据中的缓变分量，得到二者的振动分量，记为高频挠曲角数据振动分量和各通道高频波长变化量振动分量。
步骤16，分别提取高频挠曲角振动分量和各通道高频波长变化量这两组数据在全程测量过程中各自测量数据的第一段和最后一段振动数据，并分别对这两组数据的两段振动数据的幅值最大值进行归一化。
步骤17，通过波形匹配确定光纤光栅形变测量系统的采样常值延迟及实际采样频率。
步骤18，确定光纤光栅形变测量系统高频波长变化量数据的更加准确的时间标签，从而实现光纤光栅形变测量系统和子节点IMU的时间同步。
在本实施例中，利用分布式POS中子节点惯性测量单元导航解算的挠曲角与光纤光栅形变测量系统波长变化量之间的关系通过波形匹配等方式实现光纤光栅形变测量系统与分布式POS的时间同步。首先使用三次样条插值补齐子节点惯性测量单元数据和各通道光纤光栅传感器的波长变化量数据中丢包的数据，之后对该子节点惯性测量单元数据进行导航解算获取该子节点的形变挠曲角数据并将光纤光栅传感器数据和挠曲角数据插值为高频。在对高频挠曲角数据和高频波长变化量数据去噪后使用经验模态分解法分别获取其振动分量以及对应的瞬时振动频率和幅值；最后基于高频挠曲角数据和高频波长变化量数据的振动分量的瞬时振动频率和幅值通过波形匹配确定光纤光栅形变测量系统的采样常值延迟及实际采样频率。由此最终确定光纤光栅形变测量系统高频波长变化量数据的更加准确的时间标签，从而实现时间同步，为获取分布式POS数据采样时刻各节点的高精度相对运动信息提供前提，进而辅助分布式POS提供机载多任务遥感载荷的柔性基线上多个节点的高精度运动信息，辅助机载多任务遥感载荷对地观测获取更高精度的成像信息。
为了更加清晰与准确地理解与应用本公开所涉及的机载分布式POS用光纤光栅形变测量系统时间同步方法，进行以下示例。需要说明的是，本公开所保护的范围不限于以下示例。
如图2所示，为本发明另一实施例中的机载分布式POS用光纤光栅形变测量系统时间同步方法的步骤流程示意图。
具体的，如图2所示，本发明的具体方法实施如下：
1、判断光纤光栅解调仪存储的每个光纤光栅传感器通道测量的波长变化量数据和分布式POS中任意一个安装于子节点处的IMU的数据是否丢包，并使用三次样条插值法获取二者丢包时刻的数据，得到完整数据，分别记为各通道波长变化量数据和子节点IMU数据。具体实施方式为：
通过子节点IMU数据的时间标签的差值判断子节点IMU数据是否丢包，若相邻两次子节点IMU数据的时间标签的差值大于1.5倍的传感器设定的采样间隔，则认为丢包。通过光纤光栅解调仪存储的光纤光栅传感器波长变化量数据的时间标签的差值判断光纤光栅形变测量系统数据是否丢包，若相邻两次波长变化量数据的时间标签的差值大于1.5倍的传感器设定的采样间隔，则认为丢包。若子节点IMU数据和光纤光栅传感器波长变化量数据存在丢包，分别使用三次样条插值法获取二者丢包时刻的数据。三次样条插值法的步骤如下：
使用三弯矩法分别计算挠曲角数据和各通道波长变化量数据的三次样条插值函数，即可由三次样条插值函数计算得到高频数据。三次样条插值函数的定义说明如下：
对于待插值的数据段的节点x0,x1,x2,…,xn，其对应的每个节点的插值前的数据值为y0,y1,y2,…,yn，其中n为待插值的数据段的数据个数。如果函数f(x) 在节点x0,x1,x2,…,xn处的函数值为：
f(xj)＝yj,j＝0,1,2,…,n
并且关于这个节点集的三次样条函数s(x)满足插值条件：
s(xj)＝yj,j＝0,1,2,…,n
则称这个三次样条函数s(x)为三次样条插值函数。在求解三次样条插值函数后，可根据该函数获取高频插值数据。
此处选择使用三弯矩法计算三次样条插值函数s(x)，具体计算步骤如下：
三次样条插值函数s(x)在第k个小区间[xk,xk+1]上的表达式为：
s(x)＝sk1(x-xk)3+sk2(x-xk)2+sk3(x-xk)1+sk4
其中，sk1、sk2、sk3和sk4分别为三次样条插值函数s(x)的三次项、二次项、一次项和常数项系数，通过下式确定sk1、sk2、sk3和sk4即可计算出三次样条插值函数s(x)：
其中：lk为步长，即插值前的数据对应的采样时间间隔。Mk为对应第k个数据的三次样条插值函数s(x)的二阶导数，二阶导数Mk由下式确定：
其中，μk和λk为由步长计算得到的比例系数。
上述方程为欠定方程组，需根据边界条件再确定两个方程组进而求解出二阶导数Mk。对于开始和结束时刻均为静止状态的测量过程，认为两端点的二阶导数值为零，由此得到下式两个方程进而求解出二阶导数Mk。
其中，M0为待插值数据第一个采样点对应的三次样条插值函数的二阶导数的值，Mn为待插值数据最后一个采样点对应的三次样条插值函数的二阶导数的值。
由此计算得到三次样条插值函数s(x)，进而可计算得到插值后的补齐了丢包数据的高频光纤光栅波长变化量数据和挠曲角数据，记为各通道波长变化量数据和子节点IMU数据。
2、通过对子节点IMU数据的惯性导航解算，计算出该IMU所在节点的挠曲角数据，记为挠曲角数据。具体实施方式如下：
1)确定IMU姿态解算所需的初始值，包括：
a)地球地理信息：地球自转角速度ωie、地球长半径Re、地球椭圆度e、初始重力加速度g0，子节点IMU处的初始纬度la(1)、初始高度h(1)，当地卯酉圈主曲率半径Rxt(1)和子午圈主曲率半径Ryt(1)，WGS84坐标系下考虑地理位置的初始时刻重力加速度g(1)计算所需参数gk1和gk2，公式如下：
Rxt(1)＝Re·(1+(e·sin(la(1)))2)
Ryt(1)＝Re·(1-(e·(2-3sin(la(1)))2))
其中，gk1和gk2为WGS84坐标系下考虑地理位置的重力加速度计算公式中所需参数，为常值；
b)子节点IMU处的位置、姿态和速度初始值以及传感器数据采样周期：初始纬度la(1)、初始经度lon(1)、初始高度h(1)、东向初始速度
c)姿态四元数Q(k)＝[q0(k) q1(k) q2(k) q3(k)]T的初始值 Q(1)＝[q0(1) q1(1)q2(1) q3(1)]T及地理坐标系到载体坐标系的姿态矩阵的初始值：
其中，k表示第k时刻，此处k＝1。
2)读取子节点IMU数据中陀螺仪输出的角速率和加速度计输出的比力，进行四元数更新进而更新姿态矩阵，具体步骤如下：
a)计算坐标系间的转换角速度，并计算姿态角变化值：
需要计算的坐标系间的转换角速度有如下几项：
①计算地球自转角速度，即地球相对地心惯性坐标系(i系)的转动角速度在地理坐标系(t系)下的分解
其中，la(k)为k时刻的纬度，lon(k)为k时刻的经度。
②计算地理坐标系(t系)相对地球坐标系(e系)的转动角速度在地理坐标系下的分解
其中，la(k)为k时刻的纬度，lon(k)为k时刻的经度，为k时刻的东向速度，Rxt(k)为k时刻的当地卯酉圈主曲率半径，Ryt(k)为k时刻的当地子午圈主曲率半径。
③计算地理坐标系(t系)相对地心惯性坐标系(i系)的转动角速度在惯性坐标系下分解
④计算地理坐标系(t系)相对地心惯性坐标系(i系)的转动角速度在载体坐标系下分解
其中，Ti,j(i＝1,2,3；j＝1,2,3)为k时刻从地理坐标系(t系)到IMU本体坐标系 (b系)的转换矩阵中第i行第j列的元素。
⑤计算载体坐标系(b系)相对地理坐标系(t系)的转动角速度在载体坐标系(b系)下分解
由此可以分别计算出IMU安装处航向角ψ(k)、俯仰角θ(k)和横滚角γ(k) 的角度变化值δψ(k)、δθ(k)、δγ(k)：
其中，T为子节点IMU数据的采样周期。
b)使用角增量法更新四元数，计算公式如下：
四元数微分方程可写成矩阵形式如下：
使用定时采样增量法对上述四元数微分方程式解算如下：由于IMU采样时间间隔是相同的，因此上述微分方程式是关于四元数 Q(k)＝[q0(k) q1(k) q2(k) q3(k)]T的齐次线性方程，解为：
其中，Q(k)、Q(k+1)分别为k和k+1时刻的姿态四元数，
将上式进行泰勒级数展开并取前四阶近似，由于：
ΔΘ3(k)＝ΔΘ2(k)·ΔΘ(k)＝-Δθ2(k)ΔΘ(k)
ΔΘ4(k)＝ΔΘ2(k)·ΔΘ2(k)＝Δθ4(k)I
因此可得更新后的四元数表达式如下：
其中，
使用更新后的四元数Q(k+1)更新姿态矩阵
3)由姿态矩阵计算姿态角，并保证其在-2π到2π之间。具体步骤如下：
更新后的IMU的姿态矩阵的转置矩阵为由计算子IMU的航向角ψ(k+1)、俯仰角θ(k+1)和横滚角γ(k+1)，将记为
其中Tl'm(k+1)为矩阵中第l行、第m列的元素，l＝1,2,3，m＝1,2,3；则子IMU航向角ψ(k+1)、俯仰角θ(k+1)和横滚角γ(k+1)的值计算如下：
θ(k+1)＝arcsin(T3'2(k+1))
由于航向角ψ(k+1)、俯仰角θ(k+1)和横滚角γ(k+1)的取值范围分别定义为[0,2π]、
θ(k+1)＝θ(k+1)
4)根据比力方程计算子节点IMU测量中心的速度和位置
首先，将地心惯性坐标系(i系)下的比力转换为地理坐标系(t系)下的比力分量
其中：
然后，计算k时刻当地卯酉圈主曲率半径Rxt(k)和子午圈主曲率半径 Ryt(k)，并由此计算当地重力加速度g(k)；计算公式如下：
Rxt(k)＝Re·(1+(e·sin(la(k)))2)
Ryt(k)＝Re·(1-(e·(2-3sin(la(k)))2))
最后根据比力方程计算k+1时刻的地理坐标系下的速度Vt(k+1)、纬度 la(k+1)、经度lon(k+1)和高度h(k+1)。
计算k+1时刻的速度的公式如下：
其中：
计算k+1时刻的位置的公式如下：
纬度：
经度：
高度：
不断重复步骤2)至4)，直至子节点IMU的n1个数据均完成导航解算，得到整个测量过程所有时刻的挠曲角数据，记为挠曲角数据。
3、使用三次样条插值将挠曲角数据和各通道波长变化量数据分别插值为2000Hz的高频数据。具体实施过程如下：
使用三弯矩法分别计算挠曲角数据和各通道波长变化量数据的三次样条插值函数，进而得到高频数据的求解公式。三次样条插值函数的定义说明如下：
对于待插值的数据段的节点x0,x1,x2,…,xn，其对应的每个节点的插值前的数据值为y0,y1,y2,…,yn，其中n为待插值的数据段的数据个数。如果函数f(x) 在节点x0,x1,x2,…,xn处的函数值为：
f(xj)＝yj,j＝0,1,2,…,n
并且关于这个节点集的三次样条函数s(x)满足插值条件：
s(xj)＝yj,j＝0,1,2,…,n
则称这个三次样条函数s(x)为三次样条插值函数。在求解三次样条插值函数后，即可根据该函数获取高频插值数据。
这里选择使用三弯矩法计算三次样条插值函数s(x)，具体步骤如下：
三次样条插值函数s(x)在第k个小区间[xk,xk+1]上的表达式为：
s(x)＝sk1(x-xk)3+sk2(x-xk)2+sk3(x-xk)1+sk4
其中，sk1、sk2、sk3和sk4分别为三次样条插值函数s(x)的三次项、二次项、一次项和常数项系数，sk1、sk2、sk3和sk4由下式确定：
其中：lk为步长，即插值前的数据对应的采样时间间隔。Mk为对应第k个数据的三次样条插值函数s(x)的二阶导数，二阶导数Mk由下式确定：
其中，μk和λk为由步长计算得到的比例系数。
上述方程为欠定方程组，需根据边界条件再确定两个方程组进而求解出二阶导数Mk。对于开始和结束时刻均为静止状态的测量过程，认为两端点的二阶导数值为零，由此得到下式两个方程进而求解出二阶导数Mk。
其中，M0为待插值数据第一个采样点对应的三次样条插值函数的二阶导数的值，Mn为待插值数据最后一个采样点对应的三次样条插值函数的二阶导数的值。
通过上述三次样条插值分别得到2000Hz的高频挠曲角数据和各通道 2000Hz的高频波长变化量数据，分别记为高频挠曲角数据和各通道高频波长变化量数据。
4、使用基于密度空间的聚类方法去除高频挠曲角数据和各通道高频波长变化量数据中的异常值，即去除数据中的部分噪声。具体实施方式如下：
对于待去噪的数据集合Mr(即高频挠曲角数据和高频波长变化量数据) 的各个元素mri(i＝1,2,…,n)，计算DBSCAN聚类法的扫描半径E和最小包含点数MinPts。之后利用E及MinPts判断每个点是否属于某个簇，若不属于则该点为孤立点(噪声)，将其剔除。其中，扫描半径E及最小包含点数MinPts 的确定方法如下：
给定数据集P＝{p(i)；i＝0,1,...n}，对于任意点p(i)，计算点p(i)到整个数据集合的子集S＝{p(1),p(2),...,p(i-1),p(i+1),...,p(n)}中所有点之间的距离，距离按照从小到大的顺序排序，假设排序后的距离集合为 D＝{d(1),d(2),...,d(k-1),d(k),d(k+1),...,d(n)}，将d(k)称为p(i)对应的k- 距离值，将k-距离突变处对应的k-距离值作为扫描半径E，将k-距离突变处对应的k值作为最小包含点数MinPts。
5、使用经验模态分解法分别去掉去噪后的高频挠曲角数据和各通道高频波长变化量数据中的缓变分量，得到二者均值近似为零的振动分量，记为高频挠曲角数据振动分量和各通道高频波长变化量振动分量。具体实施方式如下：
1)将去噪后的高频挠曲角数据和高频波长变化量数据按照时间分为n段 xi,0(t)(i＝1,2,3…n)，每段100秒，逐段使用步骤2)至步骤6)进行处理。
2)对于步骤1)分出的数据段xi(t,k)(i＝1,2,3…n)，这里的i指的是处理第 i段时间段，k指的是执行2)至4)步骤的次数，k初值为零。通过求导确定xi(t,k)的所有极值点，并用三次样条插值法对极小值点形成下包络 emini(t,k)，对极大值形成上包络emaxi(t,k)。三次样条插值法的实施步骤如下：
使用三弯矩法分别计算挠曲角数据和各通道波长变化量数据的三次样条插值函数，进而得到高频数据的求解公式。三次样条插值函数的定义说明如下：
对于待插值的数据段的节点x0,x1,x2,…,xn，其对应的每个节点的插值前的数据值为y0,y1,y2,…,yn，其中n为待插值的数据段的数据个数。如果函数f(x) 在节点x0,x1,x2,…,xn处的函数值为：
f(xj)＝yj,j＝0,1,2,…,n
并且关于这个节点集的三次样条函数s(x)满足插值条件：
s(xj)＝yj,j＝0,1,2,…,n
则称这个三次样条函数s(x)为三次样条插值函数。在求解三次样条插值函数后，可根据该函数获取高频插值数据。
这里使用三弯矩法计算三次样条插值函数s(x)，具体步骤如下：
三次样条插值函数s(x)在第k个小区间[xk,xk+1]上的表达式为：
s(x)＝sk1(x-xk)3+sk2(x-xk)2+sk3(x-xk)1+sk4
其中，sk1、sk2、sk3和sk4分别为三次样条插值函数s(x)的三次项、二次项、一次项和常数项系数，sk1、sk2、sk3和sk4由下式确定：
其中：lk为步长，即插值前的数据对应的采样时间间隔。Mk为对应第k个数据的三次样条插值函数s(x)的二阶导数，二阶导数Mk由下式确定：
其中，μk和λk为由步长计算得到的比例系数。
上述方程为欠定方程组，需根据边界条件再确定两个方程组进而求解出二阶导数Mk。对于开始和结束时刻均为静止状态的测量过程，认为两端点的二阶导数值为零，由此得到下式两个方程进而求解出二阶导数Mk。
其中，M0为待插值数据第一个采样点对应的三次样条插值函数的二阶导数的值，Mn为待插值数据最后一个采样点对应的三次样条插值函数的二阶导数的值。
3)计算步骤2)中上下包络的均值
4)判断di(t,k)是否为内涵模态分量(Intrinsic Mode Functions,IMF)，判断方法为判断di(t,k)是否满足两个条件：①在整个数据段内，极值点的个数和过零点的个数必须相等或相差最多不能超过一个；②在任意时刻，由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的平均值为零。若di(t,k)是内涵模态分量，则从数据段xi(t,k)中去掉该分量并得到剩余信号 xi(t,k+1)＝xi(t,k)-di(t,k)；若不是，则重复步骤2)至4)直至di(t,k)为IMF 分量序列，计算xi(t,k)去掉内含模态分量di(t,k)之后的剩余信号 xi(t,k+1)＝xi(t,k)-di(t,k)。
5)对xi(t,k+1)重复步骤2)至4)，直至无法提取出IMF分量序列，得
5)对xi(t,k+1)重复步骤2)至4)，直至无法提取出IMF分量序列，得到xi(t,m)。其中，m为对应于步骤1)分出的第i段数据的内含模态分量的个数。
6)对每个数据段xi(t,k)(i＝1,2,3…n)及对应的IMF序列：
di(t,k)(k＝1,2,3,…,m)，将di(t,k)(k＝1,2,3,…,m)按照频率由高到低排序。
并将di(t,k)(k＝1,2,3,…,m)中频率较低的几阶IMF分量合并，即可获取每个数据段的均值为零的主要振动分量yi(t)(i＝1,2,3…n)，剔除高频波长变化量数据和高频挠曲角数据的缓变分量。
6、分别提取步骤(5)得到的高频挠曲角振动分量和各通道高频波长变化量这两组数据在全程测量过程中各自测量数据的第一段和最后一段振动数据，并分别对这两组数据的两段振动数据的幅值最大值进行归一化。具体实施方式如下：
1)将波长变化量振动分量以及挠曲角振动分量分别按照50ms的时间间隔进行分段，划分为p段信号yi(t)(i＝1,2,3,…,p)，即实时频率和幅值以50ms 为一个时间间隔更新。对于p段信号中的每段信号，均依次执行步骤2)至6)。
2)选择一个小波函数基其中ai,k为尺度因子，bi,k为伸缩因子，将ai,k和bi,k作为初始值。这里的小波函数基是由同一小波母函数经过伸缩和平移得到的同一组函数序列。小波母函数是在有限时间范围内变化且平均值为零的函数。
3)计算
小波系数的计算公式如下：
其中，积分区间为50ms。
4)判断CWTfi,k(ai,k,bi,k)是否达到极大值，若不满足，改变平移因子bi和尺度因子ai值，重复步骤4)。直至得到对于yi,k(t)而言CWTfi,k(ai,k,bi,k)达到极大值的尺度因子ai,k和伸缩因子bi,k。
5)根据步骤5)得到的CWTfi,k(ai,k,bi,k)达到极大值的尺度因子ai,k和伸缩因子bi,k计算yi,k(t)的瞬时振幅ei,k(t)和瞬时频率ωi,k(t)；计算方法如下：
CWTfi,k(ai,k,bi,k)可分解为实部SRi,k(t)和虚部SIi,k(t)两部分；信号的瞬时振幅可由下式确定：
信号的瞬时频率可由下式确定：
6)将波长变化量振动分量及挠曲角的振动分量每一数据段 yi(t)(i＝1,2,3,…,p)的瞬时频率和瞬时幅值按照时间顺序连接，得到整个波长变化量振动分量及挠曲角的振动分量的瞬时频率和瞬时幅值。根据瞬时幅值的变化范围寻找开始测量后的第一段振动和结束测量前的最后一段振动v1(t) 和v2(t)。
7、基于具体实施方式6得到的高频挠曲角振动分量和各通道波长变化量的各自两段振动数据，通过波形匹配确定光纤光栅形变测量系统的采样常值延迟及实际采样频率。具体实施方式如下：
首先确定各通道的传感器数据的实际采样频率，各通道实际采样频率相同。步骤如下：根据其中一个传感器通道的波长变化量振动分量和挠曲角振动分量的开始之后和结束之前的两个振动段v1(t)和v2(t)，通过两次振动的波峰波谷之间的时间差计算光纤光栅传感器和IMU的采样间隔的比例，进而以IMU数据采样时间为基准计算光纤光栅传感器的实际采样频率。公式如下：
其中，T1IMU为IMU在测量开始之后的第一个振动段的振幅最大的时刻对应的数据位置，T2IMU为IMU在测量结束之前的最后一个振动段的振幅最大的时刻对应的数据位置，T1IMU-T2IMU即IMU在测量全程的第一个振动段的振幅最大的时刻和最后一个振动段的振幅最大的时刻之间的IMU数据的数据个数。T1FBG为光纤光栅形变测量系统在测量开始之后的第一个振动段的振幅最大的时刻对应的数据位置，T2FBG为光纤光栅形变测量系统在测量结束之前的最后一个振动段的振幅最大的时刻对应的数据位置，T1FBG-T2FBG即光纤光栅形变测量系统在测量全程的第一个振动段的振幅最大的时刻和最后一个振动段的振幅最大的时刻之间的光纤光栅波长变化量数据的数据个数。
之后将各通道的波长变化量振动分量和挠曲角振动分量的第一个振动时刻的波峰波谷的时间差值作为各通道光纤光栅传感器数据相对于IMU数据的常值延迟Tdelay(j)j＝1,2,3,…,m，其中m为光纤光栅传感器的通道个数，j 表示求得的常值延迟是光纤光栅形变测量系统第j通道的光纤光栅传感器相对于IMU数据的常值延迟。
8、确定光纤光栅形变测量系统高频波长变化量数据的更加准确的时间标签，从而实现光纤光栅形变测量系统和IMU的时间同步，具体实施方式如下：
根据具体实施方式7得到的光纤光栅形变测量系统各通道的采样常值延迟及实际采样频率，结合IMU数据的时间标签给具体实施方式3得到的高频波长变化量数据打上更加准确的时间标签，之后根据所需频率降频并将对应时刻POS时间标签赋给光纤光栅，从而实现光纤光栅形变测量系统和IMU 的时间同步。
光纤光栅传感器波长变化量更加准确的时间标签计算公式如下：
TFBG(1,j)＝TIMU(1)+Tdelay(j)j＝1,2,3,…,m
其中，TFBG(k,j)为光纤光栅形变测量系统第j通道的光纤光栅传感器在第k时刻的光纤光栅波长变化量数据的时间标签，TIMU(1)为与第1个光纤光栅波长变化量数据对应的IMU数据的时间标签，fFBG为7中求得的光纤光栅形变测量系统的实际采样频率，Tdelay(j)为7中求得的光纤光栅形变测量系统第j通道的光纤光栅传感器数据与IMU数据之间的常值时间延迟。
由此，可实现光纤光栅传感器数据与分布式POS数据的更高精度时间同步，从而为通过光纤光栅传感器获取分布式POS传递对准所需时刻的柔性基线运动信息提供前提，进而辅助分布式POS提供机载多任务遥感载荷的柔性基线上多个节点的高精度运动信息，辅助机载多任务遥感载荷对地观测获取更高精度的成像信息。
综上所述，针对机载分布式POS与光纤光栅传感器的高精度时间同步问题，基于分布式POS测量的挠曲角和光纤光栅形变测量系统测量的柔性基线形变量的内在相关性，对两者数据进行基于密度空间的聚类方法的去噪、经验模态分解法和波形匹配等过程，确定光纤光栅形变测量系统数据的实际采样频率及采样常值延迟，进而获取分布式POS传递对准时刻所需的光纤光栅传感器波长变化量数据。该方法能够提高光纤光栅传感器数据的采样时间精度，进而提高了其与分布式POS的时间同步精度，为通过光纤光栅形变测量系统获取多个POS测量节点间的形变量进而获取相对运动信息提供了前提及应用可能。同时为光纤光栅传感器与其他设备进行高数据更新率及高时间同步精度的协同工作提供了数据处理思路。
基于同一发明构思，还提供了机载分布式POS用光纤光栅形变测量系统时间同步装置。由于此装置解决问题的原理与前述机载分布式POS用光纤光栅形变测量系统时间同步方法相似，因此，该装置的实施可以按照前述方法的具体步骤实现，重复之处不再赘述。
如图3所示，为一个实施例中的机载分布式POS用光纤光栅形变测量系统时间同步装置的结构示意图。该机载分布式POS用光纤光栅形变测量系统时间同步装置10包括：判断模块100、解算模块200、插值模块300、噪声去除模块400、振动分量获取模块500、提取模块600、确定模块700 和时间同步模块800。
其中，判断模块100用于判断光纤光栅解调仪存储的每个光纤光栅传感器通道测量的波长变化量数据和分布式POS中任意一个安装于子节点处的 IMU的数据是否丢包，并使用三次样条插值法获取二者丢包时刻的数据，得到完整数据，分别记为各通道波长变化量数据和子节点IMU数据；解算模块200用于通过对子节点IMU数据的惯性导航解算，计算出该IMU所在节点的挠曲角数据，记为挠曲角数据；插值模块300用于使用三次样条插值将挠曲角数据和各通道波长变化量数据分别插值为2000Hz的高频数据；噪声去除模块400用于使用基于密度空间的聚类方法去除高频挠曲角数据和各通道高频波长变化量数据中的异常值，去除数据中的部分噪声；振动分量获取模块500用于使用经验模态分解法分别去掉去噪后的高频挠曲角数据和各通道高频波长变化量数据中的缓变分量，得到二者的振动分量，记为高频挠曲角数据振动分量和各通道高频波长变化量振动分量；提取模块600用于分别提取高频挠曲角振动分量和各通道高频波长变化量这两组数据在全程测量过程中各自测量数据的第一段和最后一段振动数据，并分别对这两组数据的两段振动数据的幅值最大值进行归一化；确定模块700用于通过波形匹配确定光纤光栅形变测量系统的采样常值延迟及实际采样频率；时间同步模块 800用于确定光纤光栅形变测量系统高频波长变化量数据的更加准确的时间标签，从而实现光纤光栅形变测量系统和子节点IMU的时间同步。
图4是图示根据本公开的实施例的机载分布式POS用光纤光栅形变测量系统时间同步装置的硬件框图。如图4所示，根据本公开实施例的机载分布式POS用光纤光栅形变测量系统时间同步装置40包括存储器401和处理器402。机载分布式POS用光纤光栅形变测量系统时间同步装置40中的各组件通过总线系统和/或其它形式的连接机构(未示出)互连。
存储器401用于存储非暂时性计算机可读指令。具体地，存储器401可以包括一个或多个计算机程序产品，计算机程序产品可以包括各种形式的计算机可读存储介质，例如易失性存储器和/或非易失性存储器。易失性存储器例如可以包括随机存取存储器(RAM)和/或高速缓冲存储器(cache)等。非易失性存储器例如可以包括只读存储器(ROM)、硬盘、闪存等。
处理器402可以是中央处理单元(CPU)或者具有数据处理能力和/或指令执行能力的其它形式的处理单元，并且可以控制机载分布式POS用光纤光栅形变测量系统时间同步装置40中的其它组件以执行期望的功能。在本公开的一个实施例中，所述处理器402用于运行存储器401中存储的计算机可读指令，使得机载分布式POS用光纤光栅形变测量系统时间同步装置40 执行上述机载分布式POS用光纤光栅形变测量系统时间同步方法。机载分布式POS用光纤光栅形变测量系统时间同步装置与上述机载分布式POS用光纤光栅形变测量系统时间同步方法描述的实施例相同，在此将省略其重复描述。
图5是图示根据本公开的实施例的计算机可读存储介质的示意图。如图 5所示，根据本公开实施例的计算机可读存储介质500其上存储有非暂时性计算机可读指令501。当所述非暂时性计算机可读指令501由处理器运行时，执行参照上述描述的根据本公开实施例的机载分布式POS用光纤光栅形变测量系统时间同步方法。
以上，根据本公开实施例的机载分布式POS用光纤光栅形变测量系统时间同步方法和装置，以及计算机可读存储介质，能够基于波形识别和信号处理等技术，可将光纤光栅传感器数据与分布式POS数据进行时间同步，从而通过光纤光栅传感器获取分布式POS传递对准所需时刻的柔性基线运动信息的有益效果。
以上结合具体实施例描述了本公开的基本原理，但是，需要指出的是，在本公开中提及的优点、优势、效果等仅是示例而非限制，不能认为这些优点、优势、效果等是本公开的各个实施例必须具备的。另外，上述公开的具体细节仅是为了示例的作用和便于理解的作用，而非限制，上述细节并不限制本公开为必须采用上述具体的细节来实现。
本公开中涉及的器件、装置、设备、系统的方框图仅作为例示性的例子并且不意图要求或暗示必须按照方框图示出的方式进行连接、布置、配置。如本领域技术人员将认识到的，可以按任意方式连接、布置、配置这些器件、装置、设备、系统。诸如“包括”、“包含”、“具有”等等的词语是开放性词汇，指“包括但不限于”，且可与其互换使用。这里所使用的词汇“或”和“和”指词汇“和/或”，且可与其互换使用，除非上下文明确指示不是如此。这里所使用的词汇“诸如”指词组“诸如但不限于”，且可与其互换使用。
另外，如在此使用的，在以“至少一个”开始的项的列举中使用的“或”指示分离的列举，以便例如“A、B或C的至少一个”的列举意味着A或B或C，或AB或AC或BC，或ABC(即A和B和C)。此外，措辞“示例的”不意味着描述的例子是优选的或者比其他例子更好。
还需要指出的是，在本公开的系统和方法中，各部件或各步骤是可以分解和/或重新组合的。这些分解和/或重新组合应视为本公开的等效方案。
可以不脱离由所附权利要求定义的教导的技术而进行对在此所述的技术的各种改变、替换和更改。此外，本公开的权利要求的范围不限于以上所述的处理、机器、制造、事件的组成、手段、方法和动作的具体方面。可以利用与在此所述的相应方面进行基本相同的功能或者实现基本相同的结果的当前存在的或者稍后要开发的处理、机器、制造、事件的组成、手段、方法或动作。因而，所附权利要求包括在其范围内的这样的处理、机器、制造、事件的组成、手段、方法或动作。
提供所公开的方面的以上描述以使本领域的任何技术人员能够做出或者使用本公开。对这些方面的各种修改对于本领域技术人员而言是非常显而易见的，并且在此定义的一般原理可以应用于其他方面而不脱离本公开的范围。因此，本公开不意图被限制到在此示出的方面，而是按照与在此公开的原理和新颖的特征一致的最宽范围。
为了例示和描述的目的已经给出了以上描述。此外，此描述不意图将本公开的实施例限制到在此公开的形式。尽管以上已经讨论了多个示例方面和实施例，但是本领域技术人员将认识到其某些变型、修改、改变、添加和子组合。
- 上一篇：一种医用注射器针头装配设备
- 下一篇：闲置光纤资源在线监测系统