阅读说明：本技术 光纤光栅形变测量系统时间同步方法和装置 (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.)

1. The time synchronization method of the fiber bragg grating deformation measurement system for the airborne distributed POS is characterized by comprising the following steps:

judging whether the wavelength variation data measured by each fiber grating sensor channel stored by the fiber grating demodulator and the data of any IMU (inertial measurement unit) installed at a sub-node in the distributed POS (point of sale) lose packets or not, acquiring the data of the packet loss time of the wavelength variation data and the data of the IMU at the sub-node by using a cubic spline interpolation method to obtain complete data, and respectively recording the complete data as the wavelength variation data of each channel and the IMU data of the sub-node;

calculating the deflection angle data of the node where the IMU is located by inertial navigation calculation of the child node IMU data, and recording the deflection angle data as the deflection angle data;

respectively interpolating the deflection angle data and the wavelength variable quantity data of each channel into 2000Hz high-frequency data by using cubic spline interpolation;

removing abnormal values in the high-frequency deflection angle data and the high-frequency wavelength variation data of each channel by using a clustering method based on a density space, and removing partial noise in the data;

respectively removing the denoised high-frequency deflection angle data and the slowly-varying component in the high-frequency wavelength variable quantity data of each channel by using an empirical mode decomposition method to obtain vibration components of the high-frequency deflection angle data and the slowly-varying component in the high-frequency wavelength variable quantity data of each channel, and recording the vibration components as the vibration components of the high-frequency deflection angle data and the vibration components of the high-frequency wavelength variable quantities of each channel;

respectively extracting the first section and the last section of vibration data of the two groups of data of high-frequency deflection angle vibration components and high-frequency wavelength variation of each channel in the whole measurement process, and respectively normalizing the maximum value of the amplitude of the two sections of vibration data of the two groups of data;

determining the sampling constant delay and the actual sampling frequency of the fiber bragg grating deformation measurement system through waveform matching;

and determining a more accurate time tag of the high-frequency wavelength variation data of the fiber grating deformation measurement system, thereby realizing the time synchronization of the fiber grating deformation measurement system and the sub-node IMU.

2. The method according to claim 1, wherein the determining whether the wavelength variation data measured by each fiber grating sensor channel stored in the fiber grating demodulator and the IMU data installed at the subnode in the distributed POS are lost, and obtaining data at the time of packet loss by using a cubic spline interpolation method to obtain complete data, and the obtaining of the complete data as the wavelength variation data of each channel and the IMU data of the subnode respectively includes:

judging whether the sub-node IMU data is lost or not according to the difference value of the time labels of the sub-node IMU data, and if the difference value of the time labels of the two adjacent sub-node IMU data is more than 1.5 times of the sampling interval set by the sensor, considering that the packet is lost;

judging whether the data of the fiber grating deformation measurement system is lost or not according to the difference value of the time labels of the wavelength variation data of the fiber grating sensor stored by the fiber grating demodulator, and if the difference value of the time labels of the wavelength variation data of two adjacent times is more than 1.5 times of the sampling interval set by the sensor, considering that the data is lost;

if packet loss exists in the sub-node IMU data and the fiber bragg grating sensor wavelength variation data, obtaining data at the packet loss moment by using a cubic spline interpolation method respectively.

3. The method for time synchronization of a fiber bragg grating deformation measurement system for an airborne distributed POS according to claim 1, wherein the calculating of the bending angle data of the node where the subnode IMU is located by inertial navigation solution of the subnode IMU data, which is recorded as the bending angle data, comprises:

determining an initial value required by attitude calculation;

reading the angular rate output by a gyroscope and the specific force output by an accelerometer in the IMU data of the child nodes, and updating quaternion to update the attitude matrix;

from the attitude matrixCalculating an attitude angle and ensuring that the attitude angle is between-2 pi and 2 pi;

and calculating the speed and the position of the sub-node IMU according to a specific force equation.

4. The time synchronization method for the FBG strain measurement system for the airborne distributed POS according to claim 1, wherein the interpolation of the deflection angle data and the wavelength variation data of each channel into 2000Hz high frequency data by cubic spline interpolation comprises:

respectively calculating cubic spline interpolation functions of the deflection angle data and the wavelength variation data of each channel by using a three-bending moment method, and further obtaining a solving formula of the high-frequency data;

and respectively obtaining 2000Hz high-frequency deflection angle data and 2000Hz high-frequency wavelength variable data of each channel through a cubic spline interpolation function, and respectively recording the data as the high-frequency deflection angle data and the high-frequency wavelength variable data of each channel.

5. The time synchronization method for the FBG deformation measurement system for the airborne distributed POS according to claim 1, characterized in that a clustering method based on density space is used to remove abnormal values in the high-frequency deflection angle data and the high-frequency wavelength variation data of each channel and remove part of noise in the data;

by utilizing a concept based on density clustering, if the number of objects contained in a certain area in a clustering space is not less than a given threshold value, clustering the objects into a cluster;

connecting adjacent regions with high enough density, and eliminating isolated points which do not belong to the clusters, and defining the isolated points as noise or abnormal data points.

6. The method for synchronizing time of a FBG (fiber bragg grating) deformation measurement system for the airborne distributed POS according to claim 1, wherein the step of removing the denoised high-frequency deflection angle data and the slowly varying component in the high-frequency wavelength variation data of each channel by using an empirical mode decomposition method to obtain the vibration components of the denoised high-frequency deflection angle data and the slowly varying component in the high-frequency wavelength variation data of each channel comprises the following steps:

dividing the denoised high-frequency deflection angle data and high-frequency wavelength variable quantity data into n sections x according to time_i,0(t) (i ═ 1,2,3 … n), for 100 seconds each, processed with subsequent steps, piece by piece;

data segment x separated from the previous step_i(t, k) (i is 1,2,3 … n), where i refers to the processing of the ith time period, k refers to the number of times the subsequent steps are performed, and k has an initial value ofZero; determination of x by derivation_iAll extreme points of (t, k) and forming a lower envelope e for the extreme points by cubic spline interpolation_mini(t, k) forming an upper envelope e for the maxima_maxi(t,k)；

Calculating the mean value of the upper and lower envelopes in the previous step

Judgment of d_i(t, k) whether it is a connotative modal component; if d is_i(t, k) is the connotative modal component, then the data segment x is passed_i(t, k) removing the component and obtaining a residual signal x_i(t,k+1)＝x_i(t,k)-d_i(t, k); if not, repeating the previous steps until d_i(t, k) is the IMF component sequence, calculate x_i(t, k) removing the contained modal component d_iResidual signal x after (t, k)_i(t,k+1)＝x_i(t,k)-d_i(t,k)；

For x_i(t, k +1) repeating the steps until the IMF component sequence can not be extracted to obtain x_i(t, m); wherein m is the number of the included modal components corresponding to the ith section of data separated in the previous step;

for each data segment x_i(t, k) (i ═ 1,2,3 … n) and corresponding IMF sequences d_i(t, k) (k is 1,2,3, …, m), and d_i(t, k) (k ═ 1,2,3, …, m) is ordered by frequency from high to low, and d is assigned_i(t, k) (k is 1,2,3, …, m), the IMF components of lower frequency in several orders are merged, so as to obtain the main vibration component y with zero mean value of each data segment_i(t) (i is 1,2,3 … n), and removing the slowly varying components of the high-frequency wavelength variation data and the high-frequency deflection angle data.

7. The method for synchronizing time of a FBG (fiber bragg grating) deformation measurement system for the airborne distributed POS according to claim 1, wherein the steps of respectively extracting the first section and the last section of vibration data of the two groups of data of the high-frequency deflection angle vibration component and the high-frequency wavelength variation of each channel in the whole measurement process and respectively normalizing the maximum values of the amplitudes of the two sections of vibration data of the two groups of data comprise:

the wavelength variation vibration component and the deflection angle vibration component are respectively segmented according to the time interval of 50ms and divided into p-segment signals y_i(t) (i ═ 1,2,3, …, p), i.e. the real-time frequency and amplitude are updated at 50ms intervals. For each segment of signals in the p segments of signals, sequentially executing the subsequent steps;

the wavelet function is a same group of function sequences obtained by the same wavelet mother function through stretching and translation; the wavelet mother function is a function which changes in a limited time range and has a zero average value;

computing

wherein the integration interval is 50 ms;

determining CWTf_i,k(a_i,k,b_i,k) If not, changing the translation factor b_iAnd a scale factor a_iValue, repeat this step until y is obtained_i,k(t) CWTf for_i,k(a_i,k,b_i,k) Scale factor a up to maximum_i,kAnd scaling factor b_i,k；

CWTf obtained by the operation of the previous step_i,k(a_i,k,b_i,k) Scale factor a up to maximum_i,kAnd scaling factor b_i,kCalculating y_i,kInstantaneous amplitude e of (t)_i,k(t) and instantaneous frequency ω_i,k(t)；

The vibration component of the wavelength variation and the vibration component of the deflection angle are transmitted to each data segment y_i(t) (i ═ 1,2,3, …, p) instantaneous frequencies and instantaneous amplitudes are connected in time sequence, and instantaneous frequencies and instantaneous amplitudes of vibration components of the entire wavelength variation and the bending angle are obtained;

searching the first section of vibration after starting measurement and the last section of vibration v before finishing measurement according to the variation range of the instantaneous amplitude₁(t) and v₂(t)。

8. The method for time synchronization of a fiber grating deformation measurement system for an airborne distributed POS according to claim 1, wherein the determining the sampling constant delay and the actual sampling frequency of the fiber grating deformation measurement system by waveform matching includes:

determining an actual sampling frequency of the sensor data of each channel; wherein the actual sampling frequency of each channel is the same.

9. The method for time synchronization of a FBG (fiber bragg grating) deformation measurement system for the airborne distributed POS according to claim 8, wherein the determining the actual sampling frequency of the sensor data of each channel comprises: first and last vibration sections v of vibration component according to wavelength variation and deflection angle vibration component of one sensor channel₁(t) and v₂(t)；

And calculating the proportion of sampling intervals of the fiber grating sensor and the IMU through the time difference between the wave crests and the wave troughs of the two-time vibration, and calculating the actual sampling frequency of the fiber grating sensor by taking the IMU data sampling time as a reference.

10. The time synchronization method for the FBG strain measurement system according to claim 9, wherein the determining the sampling constant delay of the sensor data of each channel comprises: the first to separate the wavelength variation vibration component and the deflection angle vibration component of each channelThe time difference value of the wave crest and the wave trough at one vibration moment is taken as the constant delay T of the data of the fiber bragg grating sensor of each channel relative to the IMU data_delay(j)j＝1,2,3,…,m；

Wherein m is the number of channels of the fiber grating sensor, and j represents the obtained constant delay, which is the constant delay of the fiber grating sensor of the jth channel of the fiber grating deformation measurement system relative to the IMU data.

11. The method for synchronizing time of a fiber bragg grating deformation measurement system for an airborne distributed POS according to claim 1, wherein determining a more accurate time tag for high frequency wavelength variation data of the fiber bragg grating deformation measurement system, thereby achieving time synchronization of the fiber bragg grating deformation measurement system and the sub-node IMU includes:

according to the sampling constant delay and the actual sampling frequency of each channel of the fiber bragg grating deformation measurement system obtained in the previous step, a time tag with higher precision is marked on the high-frequency wavelength variable data obtained in the previous step by combining the time tag of the IMU data;

and reducing the frequency according to the required frequency and assigning the POS time label at the corresponding moment to the fiber bragg grating, thereby realizing the time synchronization of the fiber bragg grating deformation measurement system and the IMU.

12. Airborne distributed is fiber grating deformation measurement system time synchronizer for POS, its characterized in that, the device includes:

the judging module is used for judging whether the wavelength variation data measured by each fiber grating sensor channel stored in the fiber grating demodulator and the data of any IMU (inertial measurement unit) installed at a sub-node in the distributed POS (point of sale) lose packets or not, acquiring the data of the packet loss time of the wavelength variation data and the data of the IMU at the sub-node by using a cubic spline interpolation method to obtain complete data, and respectively recording the complete data as the wavelength variation data of each channel and the IMU data of the sub-node;

the resolving module is used for calculating the bending angle data of the node where the IMU is located through inertial navigation resolving of the IMU data of the sub-node, and recording the bending angle data as the bending angle data;

the interpolation module is used for respectively interpolating the deflection angle data and the wavelength variation data of each channel into 2000Hz high-frequency data by using cubic spline interpolation;

the noise removal module is used for removing abnormal values in the high-frequency deflection angle data and the high-frequency wavelength variation data of each channel by using a clustering method based on a density space, and removing partial noise in the data;

the vibration component acquisition module is used for respectively removing the denoised high-frequency deflection angle data and the slowly-varying component in the high-frequency wavelength variable quantity data of each channel by using an empirical mode decomposition method to obtain vibration components of the high-frequency deflection angle data and the slowly-varying component in the high-frequency wavelength variable quantity data of each channel, and the vibration components are recorded as the vibration components of the high-frequency deflection angle data and the vibration components of the high-frequency wavelength variable quantities of each;

the extraction module is used for respectively extracting the first section and the last section of vibration data of the two groups of data of the high-frequency deflection angle vibration component and the high-frequency wavelength variation of each channel in the whole measurement process, and respectively normalizing the maximum value of the amplitude of the two sections of vibration data of the two groups of data;

the determining module is used for determining the sampling constant delay and the actual sampling frequency of the fiber bragg grating deformation measuring system through waveform matching;

and the time synchronization module is used for determining a more accurate time label of the high-frequency wavelength variation data of the fiber grating deformation measurement system, so that the time synchronization of the fiber grating deformation measurement system and the sub-node IMU is realized.

Technical Field

The disclosure relates to the technical field of aerospace, in particular to a time synchronization method and device for a fiber bragg grating deformation measurement system, and in particular relates to a time synchronization method and device for a fiber bragg grating deformation measurement system for an airborne distributed POS.

Background

At present, the multitask remote sensing load is one of the development directions of airborne earth observation. Such as airborne distributed array antenna Synthetic Aperture Radar (SAR) and flexible multi-baseline interference SAR. For a high-performance aerial remote sensing system equipped with a multitask remote sensing load, high-precision measurement of motion parameters of all remote sensing load points on flexible baselines such as wings needs to be realized at the same time.

A Position and Orientation measurement System (POS) is a main means for acquiring motion parameters such as Position, speed and Orientation of an airborne earth-observation remote sensing load at present. For an airborne earth observation system equipped with a plurality of or a plurality of kinds of observation loads, as the plurality of or a plurality of kinds of observation loads are arranged at different positions of an airplane, the space relation among the loads is changed along with time in a complex way due to the elastic deformation of the airplane, and an airborne distributed space-time reference system, namely an airborne distributed POS, has to be established.

The airborne distributed POS generally consists of a high-precision main POS and a plurality of sub Inertial Measurement Units (IMUs). The host POS is typically mounted within the aircraft cabin; the sub-IMUs are generally distributed and installed at different positions (including wings) of the machine body and limited by weight, volume and the like, the sub-IMUs often select IMUs with medium and low precision, and the sub-IMUs need to be transmitted and aligned by means of motion parameters such as high-precision positions, speeds, postures and the like of the main POS to achieve accurate measurement of motion information of the positions. Because the wings and the body of the airplane are elastically deformed and the lever arms between the main system and the sub system are flexible and time-varying, the determination of the relative motion relation of instantaneous change between all loads becomes the premise of obtaining the high-precision motion parameters of the sub nodes through transfer alignment. The fiber grating sensor is a common means for deformation measurement, and is widely applied to structural deformation measurement and fault diagnosis. The fiber grating deformation measuring system mainly comprises a fiber grating sensor and a demodulator.

Since the deformation amplitude of flexible baselines such as wings can reach several centimeters or even tens of centimeters and the frequency can reach tens of hertz, a 1 millisecond time synchronization error can cause a millimeter-scale deformation measurement error. Therefore, when the fiber grating sensor is used for measuring the flexible baseline motion information required by the distributed POS to transfer alignment, high-precision time synchronization of the two is firstly carried out. However, currently, fiber bragg grating deformation measurement is mainly used in the occasions with low updating frequency, and no report is found on the research related to high-frequency cooperative synchronization work and high-frequency time synchronization with other devices. The problems that the existing fiber bragg grating deformation measurement system cannot realize high-frequency time synchronization mainly include two aspects of data time delay and unequal actual sampling frequency to a set value of the system. The reason that the time delay exists in the data acquisition of the fiber bragg grating deformation measurement system is as follows: the distributed POS high-precision transmission alignment needs to accurately measure the structural deformation among nodes, and a plurality of fiber bragg grating sensors need to be arranged on a distributed carrier structure of the distributed POS to obtain deformation information in a plurality of directions. In consideration of cost and other factors, the fiber grating demodulator successively reads data of each sensor in each sampling interval through switching of the optical switch, so that constant delay of data acquisition of different fiber grating channels is caused. In addition, the control process and transmission delay of the acquisition equipment in the sensor data acquisition process are finally superposed on the data acquisition time. The reason why the actual sampling frequency of the fiber bragg grating deformation measurement system is not equal to the set value is as follows: the data acquisition system of the fiber grating demodulator has clock errors (computer clock errors or crystal oscillator errors and the like). The two problems of data time delay and unequal actual sampling frequency to a set value cause that the time of the fiber bragg grating deformation measurement system data is asynchronous with that of the distributed POS system data, and the relative motion data of each node at the time required by the transmission and alignment of the distributed POS system cannot be accurately obtained.

Disclosure of Invention

In order to solve technical problems in the prior art, the embodiment of the disclosure provides a time synchronization method and a time synchronization device for a fiber bragg grating deformation measurement system for an airborne distributed POS, which can perform time synchronization on fiber bragg grating sensor data and distributed POS data based on techniques such as waveform recognition and signal processing, so that flexible baseline movement information of a distributed POS at a time required for transfer alignment is acquired through a fiber bragg grating sensor.

In a first aspect, an embodiment of the present disclosure provides a time synchronization method for a fiber bragg grating deformation measurement system for an airborne distributed POS, where the method includes: judging whether the wavelength variation data measured by each fiber grating sensor channel stored by the fiber grating demodulator and the data of any IMU (inertial measurement unit) installed at a sub-node in the distributed POS (point of sale) lose packets or not, acquiring the data of the packet loss time of the wavelength variation data and the data of the IMU at the sub-node by using a cubic spline interpolation method to obtain complete data, and respectively recording the complete data as the wavelength variation data of each channel and the IMU data of the sub-node; calculating the deflection angle data of the node where the IMU is located by inertial navigation calculation of the child node IMU data, and recording the deflection angle data as the deflection angle data; respectively interpolating the deflection angle data and the wavelength variable quantity data of each channel into 2000Hz high-frequency data by using cubic spline interpolation; removing abnormal values in the high-frequency deflection angle data and the high-frequency wavelength variation data of each channel by using a clustering method based on a density space, and removing partial noise in the data; respectively removing the denoised high-frequency deflection angle data and the slowly-varying component in the high-frequency wavelength variable quantity data of each channel by using an empirical mode decomposition method to obtain vibration components of the high-frequency deflection angle data and the slowly-varying component in the high-frequency wavelength variable quantity data of each channel, and recording the vibration components as the vibration components of the high-frequency deflection angle data and the vibration components of the high-frequency wavelength variable quantities of each channel; respectively extracting the first section and the last section of vibration data of the two groups of data of high-frequency deflection angle vibration components and high-frequency wavelength variation of each channel in the whole measurement process, and respectively normalizing the maximum value of the amplitude of the two sections of vibration data of the two groups of data; determining the sampling constant delay and the actual sampling frequency of the fiber bragg grating deformation measurement system through waveform matching; and determining a more accurate time tag of the high-frequency wavelength variation data of the fiber grating deformation measurement system, thereby realizing the time synchronization of the fiber grating deformation measurement system and the sub-node IMU.

In a second aspect, the disclosed embodiments provide a computer-readable storage medium, on which a computer program is stored, which when executed by a processor, implements the steps of the method described above.

In a third aspect, the disclosed embodiments provide a computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, the processor implementing the steps of the method described above when executing the program.

In a fourth aspect, an embodiment of the present disclosure provides a time synchronization apparatus for a fiber bragg grating deformation measurement system for an airborne distributed POS, where the apparatus includes: the judging module is used for judging whether the wavelength variation data measured by each fiber grating sensor channel stored in the fiber grating demodulator and the data of any IMU (inertial measurement unit) installed at a sub-node in the distributed POS (point of sale) lose packets or not, acquiring the data of the packet loss time of the wavelength variation data and the data of the IMU at the sub-node by using a cubic spline interpolation method to obtain complete data, and respectively recording the complete data as the wavelength variation data of each channel and the IMU data of the sub-node; the resolving module is used for calculating the bending angle data of the node where the IMU is located through inertial navigation resolving of the IMU data of the sub-node, and recording the bending angle data as the bending angle data; the interpolation module is used for respectively interpolating the deflection angle data and the wavelength variation data of each channel into 2000Hz high-frequency data by using cubic spline interpolation; the noise removal module is used for removing abnormal values in the high-frequency deflection angle data and the high-frequency wavelength variation data of each channel by using a clustering method based on a density space, and removing partial noise in the data; the vibration component acquisition module is used for respectively removing the denoised high-frequency deflection angle data and the slowly-varying component in the high-frequency wavelength variable quantity data of each channel by using an empirical mode decomposition method to obtain vibration components of the high-frequency deflection angle data and the slowly-varying component in the high-frequency wavelength variable quantity data of each channel, and the vibration components are recorded as the vibration components of the high-frequency deflection angle data and the vibration components of the high-frequency wavelength variable quantities of each; the extraction module is used for respectively extracting the first section and the last section of vibration data of the two groups of data of the high-frequency deflection angle vibration component and the high-frequency wavelength variation of each channel in the whole measurement process, and respectively normalizing the maximum value of the amplitude of the two sections of vibration data of the two groups of data; the determining module is used for determining the sampling constant delay and the actual sampling frequency of the fiber bragg grating deformation measuring system through waveform matching; and the time synchronization module is used for determining a more accurate time label of the high-frequency wavelength variation data of the fiber grating deformation measurement system, so that the time synchronization of the fiber grating deformation measurement system and the sub-node IMU is realized.

According to the time synchronization method and device for the fiber bragg grating deformation measurement system for the airborne distributed POS, provided by the invention, the time synchronization of the fiber bragg grating deformation measurement system and the distributed POS is realized by utilizing the relation between the deflection angle calculated by the navigation of the sub-node inertia measurement unit in the distributed POS and the wavelength variation of the fiber bragg grating deformation measurement system through the modes of waveform matching and the like. The method comprises the steps of firstly, completing data of a sub-node inertia measurement unit and data of packet loss in wavelength variable quantity data of each channel fiber grating sensor by using cubic spline interpolation, then carrying out navigation calculation on the data of the sub-node inertia measurement unit to obtain deformation deflection angle data of the sub-node, and interpolating the data of the fiber grating sensor and the deflection angle data into high frequency. After denoising the high-frequency deflection angle data and the high-frequency wavelength variable quantity data, respectively acquiring vibration components, corresponding instantaneous vibration frequency and amplitude by using an empirical mode decomposition method; and finally, determining the sampling constant delay and the actual sampling frequency of the fiber grating deformation measurement system through waveform matching based on the instantaneous vibration frequency and amplitude of the vibration component of the high-frequency deflection angle data and the high-frequency wavelength variation data. Therefore, a more accurate time label of the high-frequency wavelength variation data of the fiber grating deformation measurement system is finally determined, time synchronization is achieved, the premise is provided for obtaining high-precision relative motion information of each node at the sampling moment of distributed POS data, the distributed POS is assisted to provide high-precision motion information of a plurality of nodes on a flexible base line of an airborne multi-task remote sensing load, and the airborne multi-task remote sensing load is assisted to observe the ground to obtain higher-precision imaging information.

Drawings

In order to more clearly illustrate the technical solutions of the embodiments of the present disclosure, the drawings needed to be used in the description of the embodiments are briefly introduced as follows:

FIG. 1 is a schematic flow chart illustrating steps of a time synchronization method for a FBG deformation measurement system for airborne distributed POS according to an embodiment of the present invention;

FIG. 2 is a schematic flowchart illustrating steps of a time synchronization method for a FBG deformation measurement system for airborne distributed POS according to another embodiment of the present invention;

FIG. 3 is a schematic structural diagram of a time synchronization apparatus of a FBG deformation measurement system for airborne distributed POS according to an embodiment of the present invention;

FIG. 4 is a hardware block diagram of a time synchronization apparatus of a FBG deformation measurement system for onboard distributed POS according to an embodiment of the present invention;

FIG. 5 is a schematic diagram of a computer-readable storage medium in one embodiment of the invention.

Detailed Description

The present application will now be described in further detail with reference to the accompanying drawings and examples.

In the following description, the terms "first" and "second" are used for descriptive purposes only and are not intended to indicate or imply relative importance. The following description provides embodiments of the disclosure, which may be combined or substituted for various embodiments, and this application is therefore intended to cover all possible combinations of the same and/or different embodiments described. Thus, if one embodiment includes feature A, B, C and another embodiment includes feature B, D, then this application should also be considered to include an embodiment that includes one or more of all other possible combinations of A, B, C, D, even though this embodiment may not be explicitly recited in text below.

In order to make the objects, technical solutions and advantages of the present invention more clearly understood, the following describes in detail specific embodiments of a method and an apparatus for time synchronization of a fiber bragg grating deformation measurement system for an airborne distributed POS according to the present invention by using an embodiment and with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

Fig. 1 is a schematic flow chart of a time synchronization method of a fiber bragg grating deformation measurement system for an airborne distributed POS in an embodiment, which specifically includes the following steps:

and step 11, judging whether the wavelength variation data measured by each fiber grating sensor channel stored in the fiber grating demodulator and the data of any IMU (inertial measurement unit) installed at the subnode in the distributed POS (point of sale) lose packets or not, acquiring the data of the packet loss time of the wavelength variation data and the data of the IMU at the subnode by using a cubic spline interpolation method to obtain complete data, and respectively recording the complete data as the wavelength variation data of each channel and the IMU data of the subnode.

And step 12, calculating the bending angle data of the node where the IMU is located by inertial navigation calculation of the child node IMU data, and recording the bending angle data as the bending angle data.

And step 13, interpolating the deflection angle data and the wavelength variation data of each channel into 2000Hz high-frequency data by using cubic spline interpolation.

And step 14, removing abnormal values in the high-frequency deflection angle data and the high-frequency wavelength variation data of each channel by using a clustering method based on a density space, and removing partial noise in the data.

And step 15, respectively removing the denoised high-frequency deflection angle data and the slowly-varying component in the high-frequency wavelength variation data of each channel by using an empirical mode decomposition method to obtain vibration components of the denoised high-frequency deflection angle data and the denoised slowly-varying component in the high-frequency wavelength variation data of each channel, and recording the vibration components as the vibration components of the high-frequency deflection angle data and the vibration components of the high-frequency wavelength variation of each channel.

And step 16, respectively extracting the first section and the last section of vibration data of the two groups of data of the high-frequency deflection angle vibration component and the high-frequency wavelength variation of each channel in the whole measurement process, and respectively normalizing the maximum value of the amplitude of the two sections of vibration data of the two groups of data.

And step 17, determining the sampling constant delay and the actual sampling frequency of the fiber bragg grating deformation measurement system through waveform matching.

And step 18, determining a more accurate time tag of the high-frequency wavelength variation data of the fiber grating deformation measurement system, thereby realizing the time synchronization of the fiber grating deformation measurement system and the sub-node IMU.

In this embodiment, the time synchronization between the fiber grating deformation measurement system and the distributed POS is realized by using the relationship between the deflection angle calculated by the navigation of the sub-node inertia measurement unit in the distributed POS and the wavelength variation of the fiber grating deformation measurement system through waveform matching and the like. The method comprises the steps of firstly, completing data of a sub-node inertia measurement unit and data of packet loss in wavelength variable quantity data of each channel fiber grating sensor by using cubic spline interpolation, then carrying out navigation calculation on the data of the sub-node inertia measurement unit to obtain deformation deflection angle data of the sub-node, and interpolating the data of the fiber grating sensor and the deflection angle data into high frequency. After denoising the high-frequency deflection angle data and the high-frequency wavelength variable quantity data, respectively acquiring vibration components, corresponding instantaneous vibration frequency and amplitude by using an empirical mode decomposition method; and finally, determining the sampling constant delay and the actual sampling frequency of the fiber grating deformation measurement system through waveform matching based on the instantaneous vibration frequency and amplitude of the vibration component of the high-frequency deflection angle data and the high-frequency wavelength variation data. Therefore, a more accurate time label of the high-frequency wavelength variation data of the fiber grating deformation measurement system is finally determined, time synchronization is achieved, the premise is provided for obtaining high-precision relative motion information of each node at the sampling moment of distributed POS data, the distributed POS is assisted to provide high-precision motion information of a plurality of nodes on a flexible base line of an airborne multi-task remote sensing load, and the airborne multi-task remote sensing load is assisted to observe the ground to obtain higher-precision imaging information.

The following example is performed to more clearly and accurately understand and apply the time synchronization method of the fiber bragg grating deformation measurement system for airborne distributed POS according to the present disclosure. It should be noted that the protection scope of the present disclosure is not limited to the following examples.

Fig. 2 is a schematic flow chart illustrating steps of a time synchronization method for a fiber bragg grating deformation measurement system for onboard distributed POS according to another embodiment of the present invention.

Specifically, as shown in fig. 2, the specific method of the present invention is implemented as follows:

1. judging whether the wavelength variation data measured by each fiber grating sensor channel stored in the fiber grating demodulator and the data of any IMU (inertial measurement unit) installed at the sub-node in the distributed POS (point of sale) lose packets or not, acquiring the data of the packet loss time of the wavelength variation data and the data of the IMU at the sub-node by using a cubic spline interpolation method to obtain complete data, and recording the complete data as the wavelength variation data of each channel and the IMU data of the sub-node respectively. The specific implementation mode is as follows:

and judging whether the sub-node IMU data is lost or not according to the difference value of the time labels of the sub-node IMU data, and if the difference value of the time labels of the two adjacent sub-node IMU data is more than 1.5 times of the sampling interval set by the sensor, considering that the packet is lost. And judging whether the data of the fiber grating deformation measurement system is lost or not according to the difference value of the time labels of the wavelength variation data of the fiber grating sensor stored by the fiber grating demodulator, and if the difference value of the time labels of the wavelength variation data of two adjacent times is more than 1.5 times of the sampling interval set by the sensor, determining that the data is lost. If packet loss exists in the sub-node IMU data and the fiber bragg grating sensor wavelength variation data, obtaining data at the packet loss moment by using a cubic spline interpolation method respectively. The procedure of cubic spline interpolation is as follows:

and respectively calculating cubic spline interpolation functions of the deflection angle data and the wavelength variation data of each channel by using a three-bending moment method, namely calculating high-frequency data by using the cubic spline interpolation functions. The definition of the cubic spline interpolation function is explained as follows:

node x for a data segment to be interpolated₀,x₁,x₂,…,x_nThe data value before interpolation of each corresponding node is y₀,y₁,y₂,…,y_nAnd n is the data number of the data segment to be interpolated. If the function f (x) is at node x₀,x₁,x₂,…,x_nThe function value of (b) is:

f(x_j)＝y_j,j＝0,1,2,…,n

and the cubic spline function s (x) for this set of nodes satisfies the interpolation condition:

s(x_j)＝y_j,j＝0,1,2,…,n

this cubic spline function s (x) is called the cubic spline interpolation function. After solving the cubic spline interpolation function, high frequency interpolation data can be obtained according to the function.

The cubic spline interpolation function s (x) is calculated by using a three-bending moment method, and the specific calculation steps are as follows:

cubic spline interpolation function s (x) in the k-th cell [ x ]_k,x_k+1]The expression above is:

s(x)＝s_k1(x-x_k)³+s_k2(x-x_k)²+s_k3(x-x_k)¹+s_k4

wherein s is_k1、s_k2、s_k3And s_k4Coefficient of cubic term, quadratic term, primary term and constant term of cubic spline interpolation function s (x), respectively, and s is determined by the following formula_k1、s_k2、s_k3And s_k4The cubic spline interpolation function s (x) can be calculated:

wherein: l_kThe step size is the sampling time interval corresponding to the data before the plug value. M_kThe second derivative of the cubic spline interpolation function s (x) corresponding to the kth data, the second derivative M_kIs determined by the following formula:

wherein, mu_kAnd λ_kIs the scaling factor calculated from the step size.

The above equation is an underdetermined equation set, and two equation sets are determined according to the boundary condition to solve the second derivative M_k. For the measurement process with the start time and the end time both in a static state, the second derivative value of two end points is considered to be zero, so that two equations with the following formula are obtained, and the second derivative M is solved_k。

Wherein M is₀The value of the second derivative of the cubic spline interpolation function corresponding to the first sampling point of the data to be interpolated, M_nAnd the value of the second derivative of the cubic spline interpolation function corresponding to the last sampling point of the data to be interpolated.

And calculating to obtain a cubic spline interpolation function s (x), and calculating to obtain interpolated high-frequency fiber grating wavelength variation data and interpolated deflection angle data which are used for supplementing the packet loss data and are recorded as the wavelength variation data of each channel and the IMU data of the sub-nodes.

2. And calculating the bending angle data of the node where the IMU is located by inertial navigation calculation of the sub-node IMU data, and recording the bending angle data as the bending angle data. The specific implementation mode is as follows:

1) determining initial values required for IMU attitude solution, including:

a) the geographic information of the earth: angular velocity of rotation omega of the earth_ieEarth major radius Re, earth ellipticity e, initial gravitational acceleration g₀Initial latitude la (1), initial height h (1) and main curvature radius R of local mortise and unitary ring at subnode IMU_xt(1) And the meridian principal radius of curvature R_yt(1) In WGS84, the gravity acceleration g (1) at the initial time of the geographical position is considered to calculate the required parameter g_k1And g_k2The formula is as follows:

R_xt(1)＝Re·(1+(e·sin(la(1)))²)

R_yt(1)＝Re·(1-(e·(2-3sin(la(1)))²))

wherein, g_k1And g_k2The parameters required in the gravity acceleration calculation formula of the geographical position are considered under the WGS84 coordinate system and are constant values;

b) initial values of position, attitude and speed at the position of the child node IMU and a sensor data sampling period: initial latitude la (1), initial longitude lon (1), initial altitude h (1)East initial velocity

Initial velocity in north direction

Initial speed of the sky

Initial heading angle psi (1), initial pitch angle theta (1) and initial roll angle gamma (1), gyroscope and accelerometer sampling period T in the sub-node IMU.

c) Attitude quaternion q (k) ═ q₀(k) q₁(k) q₂(k) q₃(k)]^TIs equal to [ Q ] Q (1)₀(1) q₁(1)q₂(1) q₃(1)]^TAnd attitude matrix from the geographic coordinate system to the carrier coordinate systemInitial value of (a):

where k denotes the kth time, where k is 1.

2) Reading the angular rate output by a gyroscope and the specific force output by an accelerometer in the IMU data of the child nodes, and updating quaternion and further updating an attitude matrix, wherein the method specifically comprises the following steps:

a) calculating the conversion angular speed between the coordinate systems, and calculating the attitude angle change value:

the conversion angular velocity between the coordinate systems to be calculated has the following terms:

firstly, calculating the rotational angular velocity of the earth, namely the decomposition of the rotational angular velocity of the earth relative to an inertial coordinate system (i system) of the earth center under a geographic coordinate system (t system)

Where la (k) is the latitude at time k, and lon (k) is the longitude at time k.

Secondly, calculating the decomposition of the rotation angular speed of the geographic coordinate system (t system) relative to the earth coordinate system (e system) under the geographic coordinate system

Where la (k) is the latitude at time k, lon (k) is the longitude at time k,east velocity at time k, R_xt(k) Local prime radius of curvature at time k, R_yt(k) Is the local meridian principal radius of curvature at time k.

Thirdly, calculating the rotation angular velocity of the geographic coordinate system (t system) relative to the earth center inertial coordinate system (i system) and decomposing the rotation angular velocity under the inertial coordinate system

Fourthly, calculating the rotation angular velocity of the geographic coordinate system (t system) relative to the earth center inertial coordinate system (i system) and decomposing the rotation angular velocity under the carrier coordinate system

Wherein, T_i,jAnd (i is 1,2, 3; j is 1,2,3) is an element in the ith row and the jth column in the transformation matrix from the geographic coordinate system (t system) to the IMU body coordinate system (b system) at the moment k.

Calculating the rotation angular speed of the carrier coordinate system (system b) relative to the geographic coordinate system (system t) and decomposing the rotation angular speed under the carrier coordinate system (system b)

Thus, angle change values psi (k), theta (k) and gamma (k) of a heading angle psi (k), a pitch angle theta (k) and a roll angle gamma (k) at the IMU installation position can be respectively calculated:

wherein T is the sampling period of the IMU data of the child nodes.

b) And updating the quaternion by using an angle increment method, wherein the calculation formula is as follows:

the quaternion differential equation can be written in matrix form as follows:

the quaternion differential equation is solved by using a timing sampling increment method as follows: since the IMU sampling interval is phaseSimilarly, the above differential equation is related to the quaternion q (k) ═ q₀(k) q₁(k) q₂(k) q₃(k)]^TThe solution of the homogeneous linear equation of (a) is:

wherein Q (k) and Q (k +1) are attitude quaternions at the time of k and k +1, respectively,

for the value of the calculated decomposition of the rotation angular velocity of the carrier coordinate system (system b) relative to the geographic coordinate system (system t) at the time k in the carrier coordinate system (system b), dt integration interval is one IMU data sampling period, and delta theta is [ theta ═ theta [ theta ] ]_xθ_yθ_z]^TIs composed of

Starting from the moment k, the integral in one IMU data sampling period is obtained, namely the decomposition of the angular increment caused by the rotating angular speed of the carrier coordinate system (system b) relative to the geographic coordinate system (system t) in one IMU data sampling period is carried out in the carrier coordinate system (system b).

The above equation is subjected to taylor series expansion and the first four orders of approximation are taken, since:

ΔΘ³(k)＝ΔΘ²(k)·ΔΘ(k)＝-Δθ²(k)ΔΘ(k)

ΔΘ⁴(k)＝ΔΘ²(k)·ΔΘ²(k)＝Δθ⁴(k)I

the quaternion expression thus updated may be as follows:

wherein the content of the first and second substances,

and I is a third-order unit array.

Updating the attitude matrix using the updated quaternion Q (k +1)

3) From the attitude matrixThe attitude angle is calculated and guaranteed to be between-2 pi and 2 pi. The method comprises the following specific steps:

updated IMU's attitude matrixIs transposed matrix ofByCalculating heading angle psi (k +1), pitch angle theta (k +1) and roll angle gamma (k +1) of the sub IMU, and calculating the heading angle phi (k +1), the pitch angle theta (k +1) and the roll angle gamma (k +1) of the sub IMUIs marked as

Wherein T is_l'_m(k +1) is a matrixThe elements in the l-th row and the m-th column are 1,2,3, and m is 1,2, 3; the values for the sub-IMU heading angle ψ (k +1), pitch angle θ (k +1), and roll angle γ (k +1) are calculated as follows:

θ(k+1)＝arcsin(T₃'₂(k+1))

the ranges of heading angle psi (k +1), pitch angle theta (k +1) and roll angle gamma (k +1) are respectively defined as 0,2 pi]、

[-π,π](ii) a The heading angle ψ (k +1), the pitch angle θ (k +1), and the roll angle γ (k +1) are thus corrected as follows:

θ(k+1)＝θ(k+1)

4) calculating the speed and position of the sub-node IMU measurement center according to the specific force equation

Firstly, the specific force in the geocentric inertial coordinate system (i system) is converted into the specific force component in the geographic coordinate system (t system)

Wherein:

the projection of the acceleration of the sub-node IMU relative to the earth center inertial coordinate system (i system) under the sub-node IMU body coordinate system (b system).

Is a sub-sectionThe acceleration of the point IMU relative to the earth's center inertial frame (system i) is projected under the geographic frame (system t) at the sub-node IMU.

Then, calculating the main curvature radius R of the local prime circle at the time k_xt(k) And the meridian principal radius of curvature R_yt(k) And calculating therefrom a local gravitational acceleration g (k); the calculation formula is as follows:

R_xt(k)＝Re·(1+(e·sin(la(k)))²)

R_yt(k)＝Re·(1-(e·(2-3sin(la(k)))²))

finally, calculating the speed V under the geographic coordinate system at the moment k +1 according to a specific force equation^t(k +1), latitude la (k +1), longitude lon (k +1), and altitude h (k + 1).

The formula for calculating the velocity at time k +1 is as follows:

wherein:

is the projection of the acceleration of the sub-node IMU relative to the earth's center inertial frame (frame i) under the geographic frame (frame t) at the sub-node IMU.

As coordinates of the earthThe projection of the rotational angular velocity of the system (e-system) relative to the earth's center inertial frame (i-system) under the geographic frame (t-system) at the child node IMU.

Is a projection of the rotational angular velocity of the geographic coordinate system (t system) at the child node IMU relative to the earth coordinate system (e system) under the geographic coordinate system (t system) at the child node IMU.

The formula for calculating the position at time k +1 is as follows:

latitude:

longitude:

height:

wherein: r_xt(k) The main curvature radius of the local unitary fourth of twelve earthly branches at the time k, R_yt(k) Is the radius of the meridian principal curvature. k is 1,2,3 … n₁，n₁The number of data in the whole measuring process of the child node IMU which is supplemented with the lost data is increased;

continuously repeating the steps 2) to 4) until n of the child node IMU₁And (4) completing navigation calculation of all the data to obtain deflection angle data of all the moments in the whole measuring process, and recording the deflection angle data as the deflection angle data.

3. The deflection angle data and the wavelength variation data of each channel are interpolated into 2000Hz high frequency data, respectively, using cubic spline interpolation. The specific implementation process is as follows:

and respectively calculating cubic spline interpolation functions of the deflection angle data and the wavelength variation data of each channel by using a three-bending moment method, and further obtaining a solving formula of the high-frequency data. The definition of the cubic spline interpolation function is explained as follows:

node x for a data segment to be interpolated₀,x₁,x₂,…,x_nThe data value before interpolation of each corresponding node isy₀,y₁,y₂,…,y_nAnd n is the data number of the data segment to be interpolated. If the function f (x) is at node x₀,x₁,x₂,…,x_nThe function value of (b) is:

f(x_j)＝y_j,j＝0,1,2,…,n

and the cubic spline function s (x) for this set of nodes satisfies the interpolation condition:

s(x_j)＝y_j,j＝0,1,2,…,n

this cubic spline function s (x) is called the cubic spline interpolation function. After the cubic spline interpolation function is solved, high-frequency interpolation data can be obtained according to the function.

The cubic spline interpolation function s (x) is calculated by using a three-bending moment method, and the specific steps are as follows:

cubic spline interpolation function s (x) in the k-th cell [ x ]_k,x_k+1]The expression above is:

s(x)＝s_k1(x-x_k)³+s_k2(x-x_k)²+s_k3(x-x_k)¹+s_k4

wherein s is_k1、s_k2、s_k3And s_k4Coefficient of cubic term, quadratic term, primary term and constant term of cubic spline interpolation function s (x), respectively_k1、s_k2、s_k3And s_k4Is determined by the following formula:

wherein: l_kThe step size is the sampling time interval corresponding to the data before the plug value. M_kThe second derivative of the cubic spline interpolation function s (x) corresponding to the kth data, the second derivative M_kIs determined by the following formula:

wherein, mu_kAnd λ_kIs the scaling factor calculated from the step size.

The above equation is an underdetermined equation set, and two equation sets are determined according to the boundary condition to solve the second derivative M_k. For the measurement process with the start time and the end time both in a static state, the second derivative value of two end points is considered to be zero, so that two equations with the following formula are obtained, and the second derivative M is solved_k。

Wherein M is₀The value of the second derivative of the cubic spline interpolation function corresponding to the first sampling point of the data to be interpolated, M_nAnd the value of the second derivative of the cubic spline interpolation function corresponding to the last sampling point of the data to be interpolated.

The 2000Hz high-frequency deflection angle data and the 2000Hz high-frequency wavelength variation data of each channel are respectively obtained through the cubic spline interpolation and are respectively recorded as the high-frequency deflection angle data and the high-frequency wavelength variation data of each channel.

4. And removing abnormal values in the high-frequency deflection angle data and the high-frequency wavelength variation data of each channel by using a clustering method based on a density space, namely removing part of noise in the data. The specific implementation mode is as follows:

for a data set M to be denoised_r(i.e., high-frequency deflection angle data and high-frequency wavelength variation data) of each element m_ri(i ═ 1,2, …, n), the scan radius E and the minimum contained point number MinPts of the DBSCAN clustering method were calculated. And then, E and MinPts are used for judging whether each point belongs to a certain cluster, if not, the point is an isolated point (noise), and the isolated point is removed. The method for determining the scanning radius E and the minimum contained point number MinPts is as follows:

given a dataset P ═ { P (i); i 0, 1.. n }, for any point p (i), calculating the distances between all points in the subset of the point p (i) to the whole data set S { p (1), p (2),. p, p (i-1), p (i +1),. p, p (n)) }, the distances being sorted in order from small to large, assuming that the sorted distance set is D { D (1), D (2),. D, D (k-1), D (k), D +1),. D, D (n) }, D (k) is referred to as p (i) corresponding k-distance value, taking the k-distance value corresponding to the k-distance mutation as the scanning radius E, and taking the k value corresponding to the k-distance mutation as the minimum point number of MinPts.

5. And respectively removing the denoised high-frequency deflection angle data and the slowly-varying component in the high-frequency wavelength variable quantity data of each channel by using an empirical mode decomposition method to obtain a vibration component with the average value of the denoised high-frequency deflection angle data and the slowly-varying component in the high-frequency wavelength variable quantity data of each channel, and recording the vibration component as the vibration component of the high-frequency deflection angle data and the vibration component of the high-frequency wavelength variable quantity of. The specific implementation mode is as follows:

1) dividing the denoised high-frequency deflection angle data and high-frequency wavelength variable quantity data into n sections x according to time_i,0(t) (i ═ 1,2,3 … n), processed 100 seconds each, using steps 2) to 6) piece by piece.

2) For the data segment x separated in step 1)_i(t, k) (i ═ 1,2,3 … n), where i refers to the number of times that the ith time period is processed, k refers to the number of times that steps 2) to 4) are performed, and k has an initial value of zero. Determination of x by derivation_iAll extreme points of (t, k) and forming a lower envelope e for the extreme points by cubic spline interpolation_mini(t, k) forming an upper envelope e for the maxima_maxi(t, k). The implementation steps of the cubic spline interpolation method are as follows:

and respectively calculating cubic spline interpolation functions of the deflection angle data and the wavelength variation data of each channel by using a three-bending moment method, and further obtaining a solving formula of the high-frequency data. The definition of the cubic spline interpolation function is explained as follows:

node x for a data segment to be interpolated₀,x₁,x₂,…,x_nThe data value before interpolation of each corresponding node is y₀,y₁,y₂,…,y_nAnd n is the data number of the data segment to be interpolated. If the function f (x) is at node x₀,x₁,x₂,…,x_nThe function value of (b) is:

f(x_j)＝y_j,j＝0,1,2,…,n

and the cubic spline function s (x) for this set of nodes satisfies the interpolation condition:

s(x_j)＝y_j,j＝0,1,2,…,n

this cubic spline function s (x) is called the cubic spline interpolation function. After solving the cubic spline interpolation function, high frequency interpolation data can be obtained according to the function.

The cubic spline interpolation function s (x) is calculated by using a three-bending moment method, and the specific steps are as follows:

cubic spline interpolation function s (x) in the k-th cell [ x ]_k,x_k+1]The expression above is:

s(x)＝s_k1(x-x_k)³+s_k2(x-x_k)²+s_k3(x-x_k)¹+s_k4

wherein s is_k1、s_k2、s_k3And s_k4Coefficient of cubic term, quadratic term, primary term and constant term of cubic spline interpolation function s (x), respectively_k1、s_k2、s_k3And s_k4Is determined by the following formula:

wherein: l_kThe step size is the sampling time interval corresponding to the data before the plug value. M_kThe second derivative of the cubic spline interpolation function s (x) corresponding to the kth data, the second derivative M_kIs determined by the following formula:

wherein, mu_kAnd λ_kIs the scaling factor calculated from the step size.

The above equation is an underdetermined equation set, and two equation sets are determined according to the boundary condition to solve the second derivative M_k. For the measurement process with the start time and the end time both in a static state, the second derivative value of two end points is considered to be zero, so that two equations with the following formula are obtained, and the second derivative M is solved_k。

Wherein M is₀The value of the second derivative of the cubic spline interpolation function corresponding to the first sampling point of the data to be interpolated, M_nAnd the value of the second derivative of the cubic spline interpolation function corresponding to the last sampling point of the data to be interpolated.

3) Calculating the mean value of the upper envelope and the lower envelope in the step 2)

From data segment x_iSubtracting the mean value m from (t, k)_i(t,k)：d_i(t,k)＝x_i(t,k)-m_i(t,k)

4) Judgment of d_i(t, k) whether the (t, k) is an Intrinsic Mode Functions (IMF) or not, and the determination method is d_i(t, k) whether two conditions are satisfied: firstly, in the whole data segment, the number of extreme points and the number of zero-crossing points must be equal or the difference cannot exceed one at most; ② at any time, the average value of the upper envelope formed by the local maximum value points and the lower envelope formed by the local minimum value points is zero. If d is_i(t, k) is the connotative modal component, then the data segment x is passed_i(t, k) removing the component and obtaining a residual signal x_i(t,k+1)＝x_i(t,k)-d_i(t, k); if not, repeating steps 2) to 4) until d_i(t, k) is the IMF component sequence, calculate x_i(t, k) removing the contained modal component d_iResidual signal x after (t, k)_i(t,k+1)＝x_i(t,k)-d_i(t,k)。

5) For x_i(t, k +1) repeating the steps 2) to 4) until the IMF component sequence can not be extracted, thus obtaining

5) For x_i(t, k +1) repeating the steps 2) to 4) until the IMF component sequence can not be extracted, and obtaining x_i(t, m). Wherein m is the number of the included modal components corresponding to the ith section of data separated in the step 1).

6) For each data segment x_i(t, k) (i ═ 1,2,3 … n) and the corresponding IMF sequences:

d_i(t, k) (k is 1,2,3, …, m), and d_i(t, k) (k ═ 1,2,3, …, m) is ordered by frequency from high to low.

And d is_i(t, k) (k is 1,2,3, …, m), the IMF components of lower frequency in several orders are merged, so as to obtain the main vibration component y with zero mean value of each data segment_i(t) (i is 1,2,3 … n), and removing the slowly varying components of the high-frequency wavelength variation data and the high-frequency deflection angle data.

6. And (4) respectively extracting the first section and the last section of vibration data of the two groups of data of the high-frequency bending angle vibration component and the high-frequency wavelength variation of each channel obtained in the step (5) in the whole measurement process, and respectively normalizing the maximum values of the amplitudes of the two sections of vibration data of the two groups of data. The specific implementation mode is as follows:

1) the wavelength variation vibration component and the deflection angle vibration component are respectively segmented according to the time interval of 50ms and divided into p-segment signals y_i(t) (i ═ 1,2,3, …, p), i.e. the real-time frequency and amplitude are updated at 50ms intervals. Steps 2) to 6) are performed in sequence for each of the p segments of signals.

2) Selecting a wavelet function basisWherein a is_i,kIs a scale factor, b_i,kFor the scaling factor, willa_i,kAnd b_i,kAs an initial value. The wavelet function base is the same group of function sequences obtained by the same wavelet mother function through stretching and shifting. The wavelet mother function is a function that varies over a finite time range and has an average value of zero.

3) Computing

And y_i,k(t) degree of similarity between them, here using wavelet coefficients CWTf_i,k(a_i,k,b_i,k)(i＝1,2,3n；k＝1,2,3 … p), the larger the wavelet coefficient, the higher the degree of similarity.

The wavelet coefficient is calculated as follows:

wherein the integration interval is 50 ms.

4) Determining CWTf_i,k(a_i,k,b_i,k) If not, changing the translation factor b_iAnd a scale factor a_iValue, repeat step 4). Until y is obtained_i,k(t) CWTf for_i,k(a_i,k,b_i,k) Scale factor a up to maximum_i,kAnd scaling factor b_i,k。

5) CWTf obtained according to step 5)_i,k(a_i,k,b_i,k) Scale factor a up to maximum_i,kAnd scaling factor b_i,kCalculating y_i,kInstantaneous amplitude e of (t)_i,k(t) and instantaneous frequency ω_i,k(t); the calculation method is as follows:

CWTf_i,k(a_i,k,b_i,k) Decomposable into real part SR_i,k(t) and imaginary part SI_i,k(t) two parts; the instantaneous amplitude of the signal can be determined by:

the instantaneous frequency of the signal can be determined by:

6) the vibration component of the wavelength variation and the vibration component of the deflection angle are transmitted to each data segment y_i(t) (i ═ 1,2,3, …, p) are connected in chronological order to obtain the instantaneous frequency and instantaneous amplitude of the vibration component of the entire wavelength variation and the deflection angle. Searching for the first section vibration after starting measurement according to the variation range of the instantaneous amplitudeAnd the last vibration v before the end of the measurement₁(t) and v₂(t)。

7. Based on the high-frequency bending angle vibration component obtained in embodiment 6 and the respective two-segment vibration data of the wavelength variation of each channel, the sampling constant delay and the actual sampling frequency of the fiber grating deformation measurement system are determined by waveform matching. The specific implementation mode is as follows:

the actual sampling frequency of the sensor data of each channel is determined first, and the actual sampling frequency of each channel is the same. The method comprises the following steps: two vibration sections v after the start and before the end of the vibration component and the deflection angle vibration component according to the amount of wavelength change of one sensor channel₁(t) and v₂And (t) calculating the proportion of the sampling intervals of the fiber grating sensor and the IMU through the time difference between the wave crests and the wave troughs of the two-time vibration, and further calculating the actual sampling frequency of the fiber grating sensor by taking the IMU data sampling time as a reference. The formula is as follows:

wherein, T_1IMUFor the data position, T, corresponding to the IMU at the moment of maximum amplitude of the first vibration segment after the start of the measurement_2IMUData position corresponding to the moment in time when the amplitude of the last vibration segment of the IMU is maximum before the end of the measurement, T_1IMU-T_2IMUNamely, the data number of the IMU data between the moment when the amplitude of the first vibration segment of the whole measurement range is maximum and the moment when the amplitude of the last vibration segment is maximum. T is_1FBGThe data position T corresponding to the time when the amplitude of the first vibration section of the fiber bragg grating deformation measurement system is maximum after the measurement is started_2FBGThe data position T corresponding to the time when the amplitude of the last vibration section of the fiber bragg grating deformation measurement system is maximum before the measurement is finished_1FBG-T_2FBGNamely the data number of the fiber grating wavelength variation data of the fiber grating deformation measuring system between the moment when the amplitude of the first vibration section is maximum and the moment when the amplitude of the last vibration section is maximum in the whole measuring process.

Then, the time difference between the wave crest and the wave trough of the wavelength variation vibration component and the first vibration moment of the deflection angle vibration component of each channel is taken as the constant delay T of the fiber bragg grating sensor data of each channel relative to the IMU data_delay(j) j is 1,2,3, …, m, where m is the number of channels of the fiber grating sensor, and j indicates that the obtained constant delay is the constant delay of the fiber grating sensor of the jth channel of the fiber grating deformation measurement system relative to the IMU data.

8. Determining a more accurate time tag of the high-frequency wavelength variation data of the fiber grating deformation measurement system, thereby realizing the time synchronization of the fiber grating deformation measurement system and the IMU, wherein the specific implementation mode is as follows:

according to the sampling constant delay and the actual sampling frequency of each channel of the fiber grating deformation measurement system obtained in the embodiment 7, a more accurate time tag is marked on the high-frequency wavelength variation data obtained in the embodiment 3 by combining the time tag of the IMU data, and then the frequency is reduced according to the required frequency and the POS time tag at the corresponding moment is given to the fiber grating, so that the time synchronization of the fiber grating deformation measurement system and the IMU is realized.

The more accurate time tag calculation formula of the wavelength variation of the fiber grating sensor is as follows:

T_FBG(1,j)＝T_IMU(1)+T_delay(j)j＝1,2,3,…,m

wherein, T_FBG(k, j) is a time label of fiber grating wavelength variation data of a fiber grating sensor of a j channel of the fiber grating deformation measurement system at the k moment, T_IMU(1) Is a time stamp of IMU data corresponding to the 1 st FBG wavelength variation data, f_FBGFor the actual sampling frequency, T, of the fiber grating deformation measuring system calculated in 7_delay(j) And 7, obtaining the constant time delay between the data of the fiber grating sensor of the j channel of the fiber grating deformation measurement system and the IMU data.

Therefore, higher-precision time synchronization of the data of the fiber grating sensor and the data of the distributed POS can be realized, thereby providing a premise for acquiring flexible baseline motion information of the distributed POS at a moment required for transferring and aligning through the fiber grating sensor, further assisting the distributed POS to provide high-precision motion information of a plurality of nodes on a flexible baseline of the airborne multi-task remote sensing load, and assisting the airborne multi-task remote sensing load to observe the ground to acquire higher-precision imaging information.

In summary, for the problem of high-precision time synchronization between the airborne distributed POS and the fiber grating sensor, based on the internal correlation between the deflection angle measured by the distributed POS and the flexible baseline deformation measured by the fiber grating deformation measurement system, the two data are subjected to the processes of denoising, empirical mode decomposition, waveform matching and the like based on the density space clustering method, the actual sampling frequency and the sampling constant delay of the fiber grating deformation measurement system data are determined, and then the fiber grating sensor wavelength variation data required by the distributed POS at the transmission alignment time are obtained. The method can improve the sampling time precision of the fiber bragg grating sensor data, further improve the time synchronization precision of the fiber bragg grating sensor data and the distributed POS, and provide a premise and application possibility for obtaining deformation quantities among a plurality of POS measuring nodes and further obtaining relative motion information through a fiber bragg grating deformation measuring system. Meanwhile, a data processing idea is provided for cooperative work of the fiber grating sensor and other equipment with high data updating rate and high time synchronization precision.

Based on the same inventive concept, the invention also provides a time synchronization device of the fiber bragg grating deformation measurement system for the airborne distributed POS. Because the principle of solving the problems of the device is similar to the time synchronization method of the fiber bragg grating deformation measurement system for the airborne distributed POS, the implementation of the device can be realized according to the specific steps of the method, and repeated parts are not repeated.

Fig. 3 is a schematic structural diagram of a time synchronization apparatus of a fiber bragg grating deformation measurement system for onboard distributed POS in an embodiment. The time synchronization device 10 of the fiber bragg grating deformation measurement system for the airborne distributed POS includes: the device comprises a judgment module 100, a calculation module 200, an interpolation module 300, a noise removal module 400, a vibration component acquisition module 500, an extraction module 600, a determination module 700 and a time synchronization module 800.

The judging module 100 is configured to judge whether wavelength variation data measured by each fiber grating sensor channel stored in the fiber grating demodulator and data of any IMU installed at a subnode in the distributed POS are lost, and obtain data at the time of packet loss by using a cubic spline interpolation method to obtain complete data, which are respectively recorded as wavelength variation data of each channel and IMU data of the subnode; the calculation module 200 is configured to calculate, through inertial navigation calculation on the sub-node IMU data, deflection angle data of a node where the IMU is located, and record the deflection angle data as the deflection angle data; the interpolation module 300 is configured to use cubic spline interpolation to interpolate the deflection angle data and the wavelength variation data of each channel into 2000Hz high-frequency data; the noise removing module 400 is configured to remove abnormal values in the high-frequency deflection angle data and the high-frequency wavelength variation data of each channel by using a clustering method based on a density space, and remove part of noise in the data; the vibration component obtaining module 500 is configured to use an empirical mode decomposition method to remove the denoised high-frequency deflection angle data and the slowly varying component in the high-frequency wavelength variation data of each channel, to obtain vibration components of the denoised high-frequency deflection angle data and the slowly varying component in the high-frequency wavelength variation data of each channel, and to record the vibration components as the vibration components of the high-frequency deflection angle data and the vibration components of the high-frequency wavelength variation of each channel; the extracting module 600 is configured to respectively extract the first section and the last section of vibration data of the two sets of data, i.e., the high-frequency deflection angle vibration component and the high-frequency wavelength variation of each channel, in the whole measurement process, and respectively normalize the maximum amplitude values of the two sections of vibration data of the two sets of data; the determining module 700 is configured to determine a sampling constant delay and an actual sampling frequency of the fiber grating deformation measurement system through waveform matching; the time synchronization module 800 is configured to determine a more accurate time tag of the high-frequency wavelength variation data of the fiber grating deformation measurement system, so as to implement time synchronization between the fiber grating deformation measurement system and the sub-node IMU.

Fig. 4 is a hardware block diagram illustrating a time synchronization apparatus of a fiber bragg grating deformation measurement system for an onboard distributed POS according to an embodiment of the present disclosure. As shown in fig. 4, the time synchronization apparatus 40 for the fiber bragg grating deformation measurement system for onboard distributed POS according to the embodiment of the present disclosure includes a memory 401 and a processor 402. The components of the fiber grating deformation measurement system time synchronizer 40 for the onboard distributed POS are interconnected via a bus system and/or other form of connection mechanism (not shown).

The memory 401 is used to store non-transitory computer readable instructions. In particular, memory 401 may include one or more computer program products that may include various forms of computer-readable storage media, such as volatile memory and/or non-volatile memory. Volatile memory can include, for example, Random Access Memory (RAM), cache memory (or the like). The non-volatile memory may include, for example, Read Only Memory (ROM), a hard disk, flash memory, and the like.

The processor 402 may be a Central Processing Unit (CPU) or other form of processing unit having data processing capabilities and/or instruction execution capabilities, and may control other components in the fiber grating deformation measurement system time synchronizer 40 for onboard distributed POS to perform desired functions. In an embodiment of the present disclosure, the processor 402 is configured to execute the computer readable instructions stored in the memory 401, so that the time synchronization apparatus 40 for the fiber grating deformation measurement system for onboard distributed POS performs the above-mentioned method for time synchronization of the fiber grating deformation measurement system for onboard distributed POS. The time synchronization device of the fiber bragg grating deformation measurement system for the onboard distributed POS is the same as the above-described embodiment of the time synchronization method of the fiber bragg grating deformation measurement system for the onboard distributed POS, and a repeated description thereof will be omitted.

Fig. 5 is a schematic diagram illustrating a computer-readable storage medium according to an embodiment of the present disclosure. As shown in fig. 5, a computer-readable storage medium 500 according to an embodiment of the disclosure has non-transitory computer-readable instructions 501 stored thereon. The non-transitory computer readable instructions 501, when executed by a processor, perform the method for time synchronization of a fiber grating deformation measurement system for onboard distributed POS according to the embodiments of the present disclosure described above.

In the above, according to the time synchronization method and apparatus for the fiber bragg grating deformation measurement system for the airborne distributed POS and the computer-readable storage medium of the embodiment of the present disclosure, based on the techniques such as waveform recognition and signal processing, time synchronization can be performed on the fiber bragg grating sensor data and the distributed POS data, so that the beneficial effect of obtaining the flexible baseline movement information of the distributed POS at the time required for transmission alignment through the fiber bragg grating sensor is achieved.

The foregoing describes the general principles of the present disclosure in conjunction with specific embodiments, however, it is noted that the advantages, effects, etc. mentioned in the present disclosure are merely examples and are not limiting, and they should not be considered essential to the various embodiments of the present disclosure. Furthermore, the foregoing disclosure of specific details is for the purpose of illustration and description and is not intended to be limiting, since the disclosure is not intended to be limited to the specific details so described.

The block diagrams of devices, apparatuses, systems referred to in this disclosure are only given as illustrative examples and are not intended to require or imply that the connections, arrangements, configurations, etc. must be made in the manner shown in the block diagrams. These devices, apparatuses, devices, systems may be connected, arranged, configured in any manner, as will be appreciated by those skilled in the art. Words such as "including," "comprising," "having," and the like are open-ended words that mean "including, but not limited to," and are used interchangeably therewith. The words "or" and "as used herein mean, and are used interchangeably with, the word" and/or, "unless the context clearly dictates otherwise. The word "such as" is used herein to mean, and is used interchangeably with, the phrase "such as but not limited to".

Also, as used herein, "or" as used in a list of items beginning with "at least one" indicates a separate list, such that, for example, a list of "A, B or at least one of C" means A or B or C, or AB or AC or BC, or ABC (i.e., A and B and C). Furthermore, the word "exemplary" does not mean that the described example is preferred or better than other examples.

It is also noted that in the systems and methods of the present disclosure, components or steps may be decomposed and/or re-combined. These decompositions and/or recombinations are to be considered equivalents of the present disclosure.

Various changes, substitutions and alterations to the techniques described herein may be made without departing from the techniques of the teachings as defined by the appended claims. Moreover, the scope of the claims of the present disclosure is not limited to the particular aspects of the process, machine, manufacture, composition of matter, means, methods and acts described above. Processes, machines, manufacture, compositions of matter, means, methods, or acts, presently existing or later to be developed that perform substantially the same function or achieve substantially the same result as the corresponding aspects described herein may be utilized. Accordingly, the appended claims are intended to include within their scope such processes, machines, manufacture, compositions of matter, means, methods, or acts.

The previous description of the disclosed aspects is provided to enable any person skilled in the art to make or use the present disclosure. Various modifications to these aspects will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other aspects without departing from the scope of the disclosure. Thus, the present disclosure is not intended to be limited to the aspects shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

The foregoing description has been presented for purposes of illustration and description. Furthermore, this description is not intended to limit embodiments of the disclosure to the form disclosed herein. While a number of example aspects and embodiments have been discussed above, those of skill in the art will recognize certain variations, modifications, alterations, additions and sub-combinations thereof.