CN113188566A - Airborne distributed POS data fusion method - Google Patents

Abstract

The present disclosure provides an airborne distributed POS data fusion method, which specifically includes: calculating redundant position and attitude information where a child node is located; calculating a measurement vector of the child node, and using a measurement bootstrapping strategy to generate a bootstrapping amount measurement set; use the Metropolis‑Hastings sampling criterion to sample the elements of the bootstrap measurement set to obtain a valid measurement set; calculate the weight and measurement noise matrix of each valid measurement vector in the valid measurement set; use sequential filtering Data fusion is performed on the child node, and the motion parameters of the child node are corrected by the fusion result; the above steps are respectively repeated for other child nodes without data fusion. The method utilizes the motion information of all sub-nodes, reduces the adverse effect of measurement noise uncertainty on the accuracy of transfer alignment estimation, and uses the data fusion result for motion parameter correction, thus improving the motion parameters of each sub-node. estimation accuracy.

Description

An Airborne Distributed POS Data Fusion Method

technical field

The invention relates to the technical field of airborne multi-task remote sensing load node information fusion, in particular to a data fusion method of an airborne distributed position and attitude measurement system (Position and Orientation System, POS), which can be used to upgrade the airborne distributed POS sub-nodes The estimation accuracy of the motion parameters.

Background technique

The airborne position and orientation system (Position and Orientation System, POS) is the main means for obtaining the position, speed and attitude information of the airborne aircraft. Using the above information, the correction of remote sensing data can be realized and the imaging quality can be improved. Therefore, POS is often used in aerial remote sensing imaging, aerial photogrammetry and other fields. For example, remote sensing payloads such as POS-assisted Synthetic Aperture Rada (SAR), lidar, multispectral scanners, digital cameras, and large-field infrared scanners are used to meet the needs of high-precision motion imaging.

With the development of flight platform technology, the same carrier began to install multiple or multiple remote sensing equipment to achieve simultaneous observation of multiple observation windows, typically SAR, visible light camera, imaging spectrometer and lidar working at the same time, as well as multiple antennas Array antenna SAR working simultaneously. This multi-mission remote sensing mode integrating multiple or multiple remote sensing payloads has become an important development direction of aerial remote sensing systems. Since multiple observation loads are installed in different positions of the carrier aircraft, the use of a traditional single POS must not meet the needs of measuring the position and attitude information of multiple nodes. Therefore, in actual use, multiple Inertial Measurement Units (IMU) will be installed on the carrier to form a distributed position and attitude measurement system, namely distributed POS, to realize the measurement of multi-node motion parameters.

Distributed POS usually consists of one or a small number of high-precision main POS and multiple medium/low-precision sub-IMUs. These two parts are respectively installed near each remote sensing load on both sides of the airframe or wing to form a distributed multi-sensor system. The location where the main POS is located is called the main node, and the location where each sub-IMU is located is called the child node.

In order to realize the accurate measurement of the motion information at each sub-node in the distributed POS, the sub-nodes need to correct the result of their own strapdown solution according to the high-precision position, speed, attitude and other motion information provided by the main POS, and use it as The final motion parameters of the point, a process called transfer alignment. Due to the differences in the body deformation, lever arm and other factors at each sub-node, after completing the one-to-one transfer alignment from the main node to each sub-node, the accuracy of the motion parameters at each sub-node is different, ranging from high to low. In particular, the flexural deformation of the sub-nodes far from the center of the body is the most complicated, and the transfer alignment accuracy is low. If each child node can use the motion information of other nodes for data fusion, the accuracy of the motion parameters of each child node will be further improved.

At present, the data fusion algorithms for multi-sensor systems can be divided into two types: the first is centralized fusion, which fuses the measurement data of all sensors at one time, and is proved to be the smallest in information loss and the best in the world at the same time. Excellent, but for multi-sensor systems such as airborne distributed POS, the multiplication and inversion of high-dimensional matrices are involved at the time of fusion. If this fusion method is used in airborne distributed POS, it will have a large computational load and real-time performance. The disadvantage of low and low fault tolerance. The second is distributed fusion, in which the raw data information of multiple sensors is filtered through multiple filters, and then processed centrally by the fusion center. For example, the patent No. CN201810153913.5 adopts the distributed fusion method. The inverse of the covariance matrix obtained by the sub-IMU transfer and alignment is used as the weight matrix of data fusion, and the information fusion of position, velocity and attitude is given respectively. specific method. However, this method needs to calculate the weight matrix with the total number of child nodes for data fusion. The calculation of the weight matrix is complicated, and it involves the inversion of the matrix. There is a problem that the matrix falls into singularity, and the total calculation amount is also large. At the same time, the method does not consider the correlation of the deflection motion between the sub-nodes, which will adversely affect the accuracy of the data fusion results.

Commonly used structures in distributed fusion are federated filtering and sequential filtering. Both distributed fusion structures are based on Kalman filtering. Among them, federated filtering is a kind of decentralized filtering method. Its structure includes several sub-filters and a main filter. At the same time, the principle of information distribution is used to realize the distribution of information between the sub-filters and the main filter. Federated filtering has the advantages of flexible design and good fault tolerance, and has been widely used in integrated navigation systems. However, federated filtering is aimed at the fusion of sensor data with complementary measurement functions installed at a certain point, and only one sub-IMU is installed in each sub-node in the airborne distributed POS, so federated filtering cannot be formed for any sub-node. The common reference system and several subsystems required for use cannot constitute a federated filtering structure, so federated filtering is not suitable for data fusion of airborne distributed POS. Sequential filtering is to perform Kalman filtering on multiple sets of measurement data in turn, and use the final result as the data fusion result. For the airborne distributed POS, the deformation data measured by the fiber grating sensor can be used to convert the motion parameter information of each sub-node to any other sub-node, that is, any sub-node can get the rest of the sub-nodes to convert to the location where they are located. The redundant motion parameter information is used to form multiple measurement vectors, which provides the possibility to use sequential filtering to filter and estimate multiple measurement vectors of each child node, and the sequential filtering has the advantages of less computation and simple structure. The advantages.

SUMMARY OF THE INVENTION

The technical solution of the present invention is: to overcome the deficiencies of the prior art, an airborne distributed POS data fusion method is proposed, and the estimation accuracy of the motion parameters of each sub-node of the airborne distributed POS is improved.

The technical solution of the present invention is: an airborne distributed POS data fusion method. The specific steps are as follows:

( ₁ ) Calculate the redundant position and attitude information at the child node at time tk;

(2) Calculate the measurement _vector of the child node at time tk, and use the measurement bootstrapping strategy to generate a bootstrapping measurement set;

(3) Use the Metropolis-Hastings sampling criterion to sample the elements in the bootstrap measurement set to obtain the effective measurement set of the child nodes;

(4) Calculate the weight and measurement noise matrix of each valid measurement vector in the valid measurement set;

(5) use sequential filtering to perform data fusion on the child node, and use the fusion result for the correction of the motion parameter of the child node at time t _k ;

(6) Repeat steps (1) to (5) to other child nodes that do not perform data fusion and motion parameter correction, until all child nodes complete the data fusion at time t _k and the correction of motion parameters;

(7) t _k =t _k+1 , perform steps (1) to (6) until the data fusion of all child nodes at all times and the correction of motion parameters are completed.

The specific steps for obtaining redundant position and attitude information at the child nodes at time t _k in the above step (1) are as follows:

1) Sub-nodes for strapdown solution

The definition of the relevant reference coordinate system includes: e is the earth coordinate system; the navigation coordinate system at the main node and the child nodes are the northeast geographic coordinate system, the navigation coordinate system of the main node is represented by n, and the navigation coordinate system of the Jth child node and calculation The navigation coordinate system is represented by n _J and n' _J respectively, J=1,2,...,N, N is the number of child nodes; the origin of the carrier coordinate system is the center of gravity of the carrier, the x-axis is to the right along the horizontal axis of the carrier, and the y-axis Forward along the vertical axis of the carrier, the z-axis is upward along the vertical axis of the carrier, the coordinate system is fixed on the carrier, and is called the front-right upper carrier coordinate system, and b _m and b _J represent the carrier coordinates of the main node and the J-th child node respectively. system, J=1,2,...,N.

其中，L^J、λ^J和H^J分别表示第J个子节点的纬度、经度和高度；ψ^J、θ^J和γ^J分别表示第J个子节点的航向角、俯仰角和横滚角；

和

分别表示第J个子节点在n_J系下的东向速度、北向速度和天向速度。Each sub-node will get the position [L ^J λ ^J H ^J ] ^T , the attitude [ψ ^J θ ^J γ ^J ] ^T and the velocity at the time t _k after the strapdown solution.

Among them, L ^J , λ ^J and H ^J represent the latitude, longitude and altitude of the J-th child node, respectively; ψ ^J , θ ^J and γ ^J represent the heading angle, pitch angle and roll angle of the J-th child node, respectively;

and

respectively represent the easting velocity, northing velocity and sky velocity of the Jth child node under the n _J system.

2) Acquisition of redundant position and attitude information of child nodes

a) Acquisition of redundant location information of child nodes

Through the fiber grating sensor installed on the wing, the strain information at any sub-node at time t _k can be obtained, and then its displacement relative to its initial position and any two sub-nodes J, J' (J, J'=1, 2,...,N; J≠J′) relative position information ΔP _J←J′ , ΔP _J←J′ represents the difference in latitude, longitude and height between any two child nodes J and J′.

At time t _k , the position obtained by the Jth child node after strapdown solution is [L ^J λ ^J H ^J ] ^T , J=1,2,...,N, and the other N-1 child nodes after strapdown solution The position of is expressed as [L ^J′ λ ^J′ H ^J′ ] ^T (J′=1,2,…,N,J′≠J). The redundant location information at the Jth child node can be expressed as:

[L ^J←J′ λ ^J←J′ H ^J←J′ ] ^T = [L ^J′ λ ^J′ H ^J′ ] ^T +ΔP _J←J′ (J,J′＝1,2,…,N ; J≠J′)

其余N-1个位置信息[L^J←J′ λ^J←J′ H^J←J′]^T(J,J′＝1,2,…,N；J≠J′)重新记为

Thus, plus the position information obtained by the strapdown solution, there are N calculated values of latitude, longitude and altitude at the Jth child node, namely [L ^J λ ^J H ^J ] ^T and [L ^J←J′ λ ^{J ←J′} H ^J←J′ ] ^T (J, J′=1,2,…,N,J′≠J). Rewrite the position information [L ^J λ ^J H ^J ] ^T obtained by the direct strapdown solution of the current child node J as

The remaining N-1 pieces of position information [L ^J←J′ λ ^J←J′ H ^J←J′ ] ^T (J,J′=1,2,…,N; J≠J′) are rewritten as

b) Acquisition of redundant attitude information of child nodes

第J个子节点处的载体坐标系(b_J系)与该点导航坐标系(n_J系)之间的转换矩阵(姿态矩阵)表示为

其它N-1个子节点处的姿态矩阵表示为

第J个子节点处的冗余姿态矩阵可表示为：Using the angular displacement information of each sub-node at time t _k measured by the fiber grating sensor, the carrier coordinate system of any two sub-nodes J, J' (J, J' = 1, 2, ..., N; J≠J') can be constructed transformation matrix between

The transformation matrix (attitude matrix) between the carrier coordinate system (b _J system) at the Jth child node and the point navigation coordinate system (n _J system) is expressed as

The attitude matrix at the other N-1 child nodes is expressed as

The redundant pose matrix at the Jth child node can be expressed as:

in,

At time _tk , there are N pose matrices at the Jth child node. Using the N attitude matrices, N attitude angles can be calculated. The specific calculation method of the attitude angle is as follows:

记为：

Will

Record as:

中第x行、第y列的元素，x＝1,2,3，y＝1,2,3；则第J个子节点的姿态角，包括航向角ψ^J、俯仰角θ^J和横滚角γ^J的主值，即

和

分别为：

Among them, T _Jxy is a matrix

The elements in the xth row and the yth column, x=1,2,3, y=1,2,3; then the attitude angle of the Jth child node, including the heading angle ψ ^J , the pitch angle θ ^J and the roll angle The principal value of γ ^J , i.e.

and

They are:

[-π，+π]；那么，ψ^J、θ^J和γ^J的真值可由下式确定：The value ranges of the heading angle ψ ^J , the pitch angle θ ^J and the roll angle γ ^J are defined as [0, 2π],

[-π, +π]; then, the true values of ψ ^J , θ ^J and γ ^J can be determined by:

对应的N-1个冗余姿态角。将由子节点J自身的姿态矩阵

解算出的姿态角重新记为

由其余N-1个姿态矩阵

解算出的姿态角记为

According to the above calculation method of attitude angle, N-1 redundant attitude matrices at the Jth child node can also be calculated.

The corresponding N-1 redundant attitude angles. Will be composed of the attitude matrix of the child node J itself

The calculated attitude angle is re-recorded as

By the remaining N-1 pose matrices

The calculated attitude angle is recorded as

In the above step (2), the measurement vector of the child node at time t _k is calculated, and the measurement bootstrapping strategy is used to generate the specific steps of the bootstrapping measurement set:

1) Calculate the measurement vector of the child node

其表达式如下：At time t _k , the latitude, longitude and height of master node m after lever arm compensation are respectively recorded as L ^m , λ ^m , H ^m , and the heading angle, pitch angle and roll angle of master node m are respectively recorded as ψ ^m , θ ^m , γ ^m . For any child node J (J=1,2,...,N), N measurement vectors at time t _k can be obtained

Its expression is as follows:

分别代表t_k时刻子节点J的第i个纬度、经度、高度与主节点m经杆臂补偿后的纬度、经度、高度的差值；

分别代表t_k时刻子节点J的第i个航向角、俯仰角、横滚角与主节点m对应的航向角、俯仰角、横滚角的差值。in,

respectively represent the difference between the ith latitude, longitude and height of the child node J at time _tk and the latitude, longitude and height of the main node m after compensation by the lever arm;

respectively represent the difference between the ith heading angle, pitch angle, and roll angle of the child node J at time _tk and the heading angle, pitch angle, and roll angle corresponding to the main node m.

2) Use the measurement bootstrapping strategy to generate a bootstrapping measurement set

的基础上，通过增加扰动的方式生成自举量测

具体步骤如下：in the measurement vector

On the basis of , the bootstrap measurements are generated by adding perturbations

Specific steps are as follows:

表示以原始量测向量

为基础产生的第c个自举量测，c＝1,2,…,l，l表示自举量测

的总个数，l＝5；H^J是系统的量测矩阵，x^J是系统的状态向量，v^J是量测噪声，且v^J满足是零均值的高斯白噪声，其协方差为

和v^J具有相同的统计特性，即

满足是零均值的高斯白噪声，其协方差

in,

represents the original measurement vector

The c-th bootstrap measurement generated based on c=1,2,...,l, where l represents the bootstrap measurement

The total number of , l=5; H ^J is the measurement matrix of the system, x ^J is the state vector of the system, v ^J is the measurement noise, and v ^J is Gaussian white noise with zero mean, and its covariance is

has the same statistical properties as v ^J , namely

Satisfy is white Gaussian noise with zero mean, and its covariance

定义自举量测集合：

其中

代表原始量测向量，即

Through the above method, N×l bootstrap measurement data can be obtained

Define the bootstrap measurement set:

in

represents the original measurement vector, i.e.

In the above step (3), the Metropolis-Hastings sampling criterion is used to sample the elements in the bootstrap measurement set to obtain the effective measurement set of the child nodes. The specific steps are as follows:

1) Calculate the credibility and acceptance probability

中随意选择两个集合

再从这两个集合中分别随机抽取一个元素

和

并计算对应的可信度

和

可信度的计算公式如下所示：According to the Metropolis-Hastings sampling principle, from N bootstrap measurement sets

Randomly choose two sets from

Then randomly select an element from each of these two sets

and

and calculate the corresponding reliability

and

The formula for calculating reliability is as follows:

是量测预测均值，

代表量测噪声v^J协方差矩阵

的行列式。in

is the mean of the measurement forecast,

represents the measurement noise v ^J covariance matrix

determinant of .

和

的基础上，按照下式计算接受概率：getting credibility

and

On the basis of , the acceptance probability is calculated according to the following formula:

的取值为

和1之间的最小值。the acceptance probability

value of

and the smallest value between 1.

2) Build an effective measurement set

First, a random number χ that satisfies the uniform distribution U(0,1) is generated. The method for determining the effective measurement set Ω ^J is as follows:

和

对应的接受概率大于等于生成的随机数χ，则将

作为有效量测集合Ω^J中的元素；反之，将

作为有效量测集合Ω^J中的元素。上述过程即为确定有效量测向量的采样过程。In the formula, if the randomly selected elements

and

The corresponding acceptance probability is greater than or equal to the generated random number χ, then the

as an element in the effective measurement set Ω ^J ; otherwise, the

as an element in the effective measurement set Ω ^J. The above process is the sampling process for determining the effective measurement vector.

Repeat L sampling, L=2N, and define L effective measurement vectors in the effective measurement set as z ^J (1), z ^J (2),...z ^J (L), these effective measurement vectors are composed of Valid measurement set Ω ^J .

The specific steps of calculating the weight of each effective measurement vector and the measurement noise matrix in the effective measurement set in the above step (4) are as follows:

1) Calculate the consistency distance and consistency matrix

其计算方法如下：Define the confidence distance between any two measurement vectors z ^J (ξ) and z ^J (ζ) in the effective measurement set Ω ^J as

Its calculation method is as follows:

和一致性矩阵Ψ^J，其计算方法如下：Calculate the consistency distance on this basis

and the consistency matrix Ψ ^J , which is calculated as follows:

表示Ω^J中任意两个量测向量之间置信距离的最大值。in,

represents the maximum confidence distance between any two measurement vectors in Ω ^J.

2) Calculate the weight of the effective measurement vector and the measurement noise matrix

令权向量

此时，

可以用来表示相应的有效量测向量被有效量测集合Ω^J中所有元素的总体支持程度。After obtaining the consistency matrix Ψ ^J , calculate the maximum modulus eigenvalue λ of the matrix and the corresponding eigenvector β, and unite β to get

Let the weight vector

at this time,

It can be used to represent the overall support degree of all elements in the effective measurement set Ω ^J for the corresponding effective measurement vector.

其计算方法如下：According to the basic principle of variance calculation in mathematical statistics, the measurement noise matrix corresponding to each element in the effective measurement set is calculated.

Its calculation method is as follows:

In the above step (5), sequential filtering is used to perform data fusion on the child node, and the fusion result is used to correct the motion parameters of the child node at time t _k . The specific steps are as follows:

1) Establish an airborne distributed POS delivery alignment model

速度误差

位置误差(δL^J、δλ^J、δH^J)、陀螺常值漂移

和加计常值偏置

N为子系统的个数。其中，

分别为第J个子节点的东向、北向、天向失准角；

分别为第J个子节点在n_J系下的东向、北向和天向速度误差；δL^J、δλ^J、δH^J分别为第J个子节点的纬度误差、经度误差、高度误差；

和

分别为第J个子节点处子IMU在其载体坐标系下(b_J系)陀螺常值漂移和加计常值偏置；

分别为第J个子节点处子IMU在b_J系x轴、y轴和z轴上的陀螺仪随机常值漂移；

分别为第J个子节点处子IMU在b_J系x轴、y轴和z轴上的加速度计常值偏置。因此系统状态向量有如下形式：The transfer alignment model adopts the matching method of "position + attitude", and the system state vector is a 15-dimensional state vector, including the platform misalignment angle of the child node J

speed error

Position error (δL ^J , δλ ^J , δH ^J ), gyro constant drift

and add-up constant offset

N is the number of subsystems. in,

are the east, north and sky misalignment angles of the J-th child node;

are the easting, northing and sky velocity errors of the J-th child node under the n _J system, respectively; δL ^J , δλ ^J , and δH ^J are the latitude error, longitude error, and altitude error of the J-th child node;

and

are respectively the gyro constant drift and the accumulating constant offset of the child IMU at the Jth child node in its carrier coordinate system (b _J system);

are the random constant drift of the gyroscope on the x-axis, y-axis and z-axis of the b _J system of the child IMU at the Jth child node, respectively;

are the constant accelerometer offsets of the child IMU at the Jth child node on the x-axis, y-axis and z-axis of the b _J system, respectively. Therefore, the system state vector has the following form:

a) Establishment of the equation of state

The state equation of the system has the following form:

为零均值高斯白噪声，其方差阵

由IMU中陀螺仪和加速度计的噪声水平决定。F^J为状态转移矩阵，其具体表达形式如下：where system noise

White Gaussian noise with zero mean, its variance matrix

Determined by the noise level of the gyroscope and accelerometer in the IMU. ^FJ is the state transition matrix, and its specific expression is as follows:

in,

和

分别为子节点J处子IMU在导航坐标系(n_J系)中的东向比力、北向比力和天向比力分量；

和

分别为子节点J在导航坐标系(n_J系)下的东向速度、北向速度和天向速度；

和

分别为子节点J沿子午圈和卯酉圈的主曲率半径；T为滤波周期；G^J为子节点J的系统噪声矩阵，其具体表达形式如下：Among them, ω _ie represents the angular velocity of the earth's rotation;

and

are the eastward specific force, northward specific force and skyward specific force components of the child IMU at the child node J in the navigation coordinate system (n _J system), respectively;

and

are the easting speed, northing speed and sky speed of the child node J in the navigation coordinate system (n _J system) respectively;

and

are respectively the principal radius of curvature of the sub-node J along the meridian and 卯unitary circles; T is the filter period; G ^J is the system noise matrix of the sub-node J, and its specific expression is as follows:

b) Establishment of measurement equations

The specific expressions of the L effective measurement vectors z ^J (τ) (τ=1,2,...,L) of the child node J are:

分别表示子节点J的第τ个纬度、经度、高度计算值与主节点m经杆臂补偿后的纬度、经度、高度的差值；

和

分别表示子节点J第τ个航向角、俯仰角、横滚角计算值与主节点m对应的航向角、俯仰角、横滚角的差值；in,

respectively represent the difference between the calculated value of the τth latitude, longitude and height of the child node J and the latitude, longitude and height of the main node m after compensation by the lever arm;

and

respectively represent the difference between the calculated value of the τth heading angle, pitch angle and roll angle of the child node J and the heading angle, pitch angle and roll angle corresponding to the main node m;

The expression of the measurement equation:

The measurement matrix corresponding to z ^J (τ)

in,

且令

表示矩阵T_a中第s行、第t列的元素，即T_a可以表示为：T _a is the attitude matrix of the main node, i.e.

and make

Represents the elements of the s-th row and the t-th column in the matrix T _a , that is, T _a can be expressed as:

2) Use sequential filtering for data fusion

和传递对准模型可以得到当前t_k时刻的一步预测估计值

计算公式如下所示：First, the time update is performed, and for the Jth child node, according to the filter estimated value of the previous time t _k-1

and the transfer alignment model can get a one-step prediction estimate at the current time t _k

The calculation formula is as follows:

Among them, Γ ^J is the transition matrix of the discrete system obtained by discretizing the state transition matrix F ^J of the continuous system.

Similarly, according to the estimated mean square error matrix P ^J (k-1|k-1) at the previous moment and the transfer alignment model, the one-step forecast mean square error matrix of the current moment t _k can be obtained. The calculation formula is as follows:

Among them, Ξ ^J is the noise-driven array of the discrete system obtained by discretizing the noise-driven array G ^J of the continuous system.

第J个子节点的数据融合公式如下：Then use the effective measurement vector obtained in step (3) to perform measurement update, and record the L measurement vectors of the _Jth child node at the current time tk as

The data fusion formula of the Jth child node is as follows:

进行量测更新时，有：when using

When performing a measurement update, there are:

和

分别为

所对应的量测矩阵和量测噪声矩阵；

为使用

进行量测更新时计算得到的滤波增益矩阵；

表示使用

进行量测更新时计算得到的状态估计值；

为

对应的估计均方误差阵；I为与

维度相同的单位阵。in,

and

respectively

The corresponding measurement matrix and measurement noise matrix;

for use

The filter gain matrix calculated when the measurement update is performed;

indicate use

The state estimate calculated when the measurement update is performed;

for

The corresponding estimated mean square error matrix; I is the

A unit matrix of the same dimensions.

依次进行量测更新时，有：when using

When performing measurement updates in sequence, there are:

和

分别为

对应的量测矩阵和量测噪声矩阵；

为使用

进行量测更新时计算得到的滤波增益矩阵；

表示使用

进行量测更新时计算得到的状态估计值；

为

对应的估计均方误差阵。in,

and

respectively

Corresponding measurement matrix and measurement noise matrix;

for use

The filter gain matrix calculated when the measurement update is performed;

indicate use

The state estimate calculated when the measurement update is performed;

for

The corresponding estimated mean squared error matrix.

作为当前t_k时刻第J个子节点处最终的数据融合结果，即

其中，

包含t_k时刻第J个子节点数据融合后的纬度误差δL^J(k)、经度误差δλ^J(k)和高度误差δH^J(k)，还包含数据融合后在n_J系下的东向速度误差

北向速度误差

和天向速度误差

以及数据融合后的东向失准角

北向失准角

和天向失准角

Execute the above data fusion formula, when the above process runs to τ=L, the obtained

As the final data fusion result at the Jth child node at the current _tk time, that is,

in,

It includes the latitude error δL ^J (k), longitude error δλ ^J (k) and altitude error δH ^J (k) after data fusion of the Jth child node at time t _k , and also includes the eastward velocity under the n _J system after data fusion error

northbound speed error

and sky speed error

and the easting misalignment angle after data fusion

North misalignment angle

and sky misalignment

3) Motion parameter correction

Use the above fusion result to correct the result of the strapdown solution of the _Jth child node at time tk. The calibration steps are as follows:

a) Position correction

L ^J (k)| _new =L ^J (k)-δL ^J (k),λ ^J (k)| _new =λ ^J (k)-δλ ^J (k), H ^J (k)| _new =H ^J (k)-δH ^J (k)

Among them, L ^J (k), λ ^J (k), H ^J (k) represent the latitude, longitude and altitude obtained by the strapdown solution of the J-th child node at time t _k , respectively; L ^J (k) | _new , λ ^J (k)| _new and H ^J ( _k )| _new respectively represent the corrected latitude, longitude and altitude of the Jth child node at time tk.

b) Speed correction

分别代表t_k时刻第J个子节点捷联解算得到的东向速度、北向速度和天向速度；

分别代表t_k时刻第J个子节点校正后的东向速度、北向速度和天向速度。in,

respectively represent the easting velocity, northing velocity and sky velocity obtained by the strapdown solution of the _Jth child node at time tk;

respectively represent the corrected easting speed, northing speed and sky speed of the _Jth child node at time tk.

c) Attitude correction

Calculate the transformation matrix between the navigation coordinate system n _J of the J-th child node at time t _k and the calculated navigation coordinate system n _J ′

为：

其中，

为t_k时刻J个子节点进行捷联解算后得到的姿态矩阵；Corrected transformation matrix

for:

in,

is the attitude _matrix obtained by the strapdown solution for J child nodes at time tk;

计算t_k时刻第J个子节点的航向角ψ^J(k)|_new、俯仰角θ^J(k)|_new和横滚角γ^J(k)|_new。Using the modified transformation matrix

Calculate the heading angle ψ ^J (k)| _new , the pitch angle θ ^J (k)| _new and the roll angle γ ^J (k)| _new of the J-th child node at time t _k .

Repeat steps (1) to (5) for other sub-nodes that do not perform data fusion and motion parameter correction in the above step (6), until all sub-nodes complete the data fusion and motion parameter correction at time t _k , and the specific steps are as follows:

1) Let the number of the child nodes that have completed steps (1) to (5) at time t _k be k _N , and the initial value of k _N is 1;

2) If k _N <N, and N is the number of child nodes, it means that there are still child nodes that have not completed data fusion and motion parameter correction, k _N =k _N +1, repeat the steps for the child nodes numbered k _N (1) to (5);

3) If k _N =N, stop step (6), indicating that all child nodes have completed data fusion and motion parameter correction at time t _k .

The advantages of the present invention compared with the prior art are:

Focusing on the goal of improving the overall accuracy of airborne distributed POS multi-node motion parameter measurement, the measurement bootstrapping strategy and sequential filtering are combined. The motion parameter information of all child nodes of the former is used to obtain measurement vectors that are closer to the true value of each child node, and then use these measurement vectors to perform sequential filtering, and the obtained data fusion results are used for the correction of child node motion parameters. . The method makes full use of the motion parameter information of all sub-nodes, overcomes the adverse effect of measurement noise uncertainty on the estimation accuracy of transfer alignment, and finally improves the estimation accuracy of motion parameters of sub-nodes.

Description of drawings

Fig. 1 is the flow chart of the present invention;

FIG. 2 is a structural diagram of a data fusion method based on measurement bootstrapping and sequential filtering of the present invention.

Detailed ways

As shown in Figure 1, the concrete method of the present invention is implemented as follows:

1. Calculate the redundant position and attitude information of the child nodes at time tk, and the specific calculation steps are as _follows :

(1) Sub-nodes for strapdown solution

The definition of the relevant reference coordinate system includes: e is the earth coordinate system; the navigation coordinate system at the main node and the child nodes are the northeast geographic coordinate system, the navigation coordinate system of the main node is represented by n, and the navigation coordinate system of the Jth child node and calculation The navigation coordinate system is represented by n _J and n' _J respectively, J=1,2,...,N, N is the number of child nodes; the origin of the carrier coordinate system is the center of gravity of the carrier, the x-axis is to the right along the horizontal axis of the carrier, and the y-axis Forward along the vertical axis of the carrier, and the z-axis is upward along the vertical axis of the carrier. The coordinate system is fixed on the carrier and is called the front-right upper carrier coordinate system, and b _m and b _J represent the carrier coordinates of the main node and the J-th child node respectively. system, J=1,2,...,N.

其中，L^J、λ^J和H^J分别表示第J个子节点的纬度、经度和高度；ψ^J、θ^J和γ^J分别表示第J个子节点的航向角、俯仰角和横滚角；

和

分别表示第J个子节点在n_J系下的东向速度、北向速度和天向速度。Each sub-node will get the position [L ^J λ ^J H ^J ] ^T , the attitude [ψ ^J θ ^J γ ^J ] ^T and the velocity at the time t _k after the strapdown solution.

Among them, L ^J , λ ^J and H ^J represent the latitude, longitude and altitude of the J-th child node, respectively; ψ ^J , θ ^J and γ ^J represent the heading angle, pitch angle and roll angle of the J-th child node, respectively;

and

respectively represent the easting velocity, northing velocity and sky velocity of the Jth child node under the n _J system.

(2) Acquisition of redundant position and attitude information of child nodes

a) Acquisition of redundant location information of child nodes

Through the fiber grating sensor installed on the wing, the strain information at any sub-node at time t _k can be obtained, and then its displacement relative to its initial position and any two sub-nodes J, J' (J, J'=1, 2,...,N; J≠J′) relative position information ΔP _J←J′ , ΔP _J←J′ represents the difference in latitude, longitude and height between any two child nodes J and J′.

At time t _k , the position obtained by the Jth child node after strapdown solution is [L ^J λ ^J H ^J ] ^T , J=1,2,...,N, and the other N-1 child nodes after strapdown solution The position of is expressed as [L ^J′ λ ^J′ H ^J′ ] ^T (J′=1,2,…,N,J′≠J). The redundant location information at the Jth child node can be expressed as:

[L ^J←J′ λ ^J←J′ H ^J←J′ ] ^T = [L ^J′ λ ^J′ H ^J′ ] ^T +ΔP _J←J′ (J,J′＝1,2,…,N ; J≠J′)

其余N-1个位置信息[L^J←J′ λ^J←J′ H^J←J′]^T(J,J′＝1,2,…,N；J≠J′)重新记为

Thus, plus the position information obtained by the strapdown solution, there are N calculated values of latitude, longitude and altitude at the Jth child node, namely [L ^J λ ^J H ^J ] ^T and [L ^J←J′ λ ^{J ←J′} H ^J←J′ ] ^T (J, J′=1,2,…,N,J′≠J). Rewrite the position information [L ^J λ ^J H ^J ] ^T obtained by the direct strapdown solution of the current child node J as

The remaining N-1 pieces of position information [L ^J←J′ λ ^J←J′ H ^J←J′ ] ^T (J,J′=1,2,…,N; J≠J′) are rewritten as

b) Acquisition of redundant attitude information of child nodes

第J个子节点处的载体坐标系(b_J系)与该点导航坐标系(n_J系)之间的转换矩阵(姿态矩阵)表示为

其它N-1个子节点处的姿态矩阵表示为

第J个子节点处的冗余姿态矩阵可表示为：Using the angular displacement information of each sub-node at time t _k measured by the fiber grating sensor, the carrier coordinate system of any two sub-nodes J, J' (J, J' = 1, 2, ..., N; J≠J') can be constructed transformation matrix between

The transformation matrix (attitude matrix) between the carrier coordinate system (b _J system) at the Jth child node and the point navigation coordinate system (n _J system) is expressed as

The attitude matrix at the other N-1 child nodes is expressed as

The redundant pose matrix at the Jth child node can be expressed as:

in,

At time _tk , there are N pose matrices at the Jth child node. Using the N attitude matrices, N attitude information, that is, N attitude angles, can be calculated. The specific calculation method of the attitude angle is as follows:

记为：Will

Record as:

中第x行、第y列的元素，x＝1,2,3，y＝1,2,3；则第J个子节点的姿态角，包括航向角ψ^J、俯仰角θ^J和横滚角γ^J的主值，即

和

分别为：Among them, T _Jxy is a matrix

The elements in the xth row and the yth column, x=1,2,3, y=1,2,3; then the attitude angle of the Jth child node, including the heading angle ψ ^J , the pitch angle θ ^J and the roll angle The principal value of γ ^J , i.e.

and

They are:

[-π，+π]；那么，ψ^J、θ^J和γ^J的真值可由下式确定：The value ranges of the heading angle ψ ^J , the pitch angle θ ^J and the roll angle γ ^J are defined as [0, 2π],

[-π, +π]; then, the true values of ψ ^J , θ ^J and γ ^J can be determined by:

对应的N-1个冗余姿态角。将由子节点J自身的姿态矩阵

解算出的姿态角重新记为

由其余N-1个姿态矩阵

解算出的姿态角记为

According to the above attitude angle calculation method, N-1 redundant attitude matrices at the Jth child node can also be calculated.

The corresponding N-1 redundant attitude angles. Will be composed of the attitude matrix of the child node J itself

The calculated attitude angle is re-recorded as

By the remaining N-1 pose matrices

The calculated attitude angle is recorded as

2. Calculate the measurement _vector of the child node at time tk, and use the measurement bootstrapping strategy to generate a bootstrapping measurement set. The specific steps are as follows:

(1) Calculate the measurement vector of the child node

其表达式如下：At time t _k , the latitude, longitude and height of master node m after lever arm compensation are recorded as L ^m , λ ^m , and H ^m respectively, and the heading angle, pitch angle and roll angle of master node m are respectively recorded as ψ ^m , θ ^m , γ ^m . For any child node J (J=1,2,...,N), N measurement vectors at time t _k can be obtained

Its expression is as follows:

分别代表t_k时刻子节点J的第i个纬度、经度、高度与主节点m经杆臂补偿后的纬度、经度、高度的差值；

分别代表t_k时刻子节点J的第i个航向角、俯仰角、横滚角与主节点m对应的航向角、俯仰角、横滚角的差值。in,

respectively represent the difference between the ith latitude, longitude and height of the child node J at time _tk and the latitude, longitude and height of the main node m after compensation by the lever arm;

respectively represent the difference between the ith heading angle, pitch angle, and roll angle of the child node J at time _tk and the heading angle, pitch angle, and roll angle corresponding to the main node m.

(2) Use the measurement bootstrapping strategy to generate a bootstrapping measurement set

的基础上，通过增加扰动的方式生成自举量测

具体步骤如下：in the measurement vector

On the basis of , the bootstrap measurements are generated by adding perturbations

Specific steps are as follows:

表示以原始量测向量

为基础产生的第c个自举量测，c＝1,2,…,l，l表示自举量测

的总个数，l＝5；H^J是系统的量测矩阵，x^J是系统的状态向量，v^J是量测噪声，且v^J满足是零均值的高斯白噪声，其协方差为

和v^J具有相同的统计特性，即

满足是零均值的高斯白噪声，其协方差

in,

represents the original measurement vector

The c-th bootstrap measurement generated based on c=1,2,...,l, where l represents the bootstrap measurement

The total number of , l=5; H ^J is the measurement matrix of the system, x ^J is the state vector of the system, v ^J is the measurement noise, and v ^J is Gaussian white noise with zero mean, and its covariance is

has the same statistical properties as v ^J , namely

Satisfy is white Gaussian noise with zero mean, and its covariance

定义自举量测集合：

其中

代表原始量测向量，即

Through the above method, N×l bootstrap measurement data can be obtained

Define the bootstrap measurement set:

in

represents the original measurement vector, i.e.

3. Use the Metropolis-Hastings sampling criterion to sample the elements in the bootstrap measurement set obtained in step 2 to obtain the effective measurement set of the child nodes. The specific steps are as follows:

(1) Calculate the reliability and acceptance probability

中随意选择两个集合

再从这两个集合中分别随机抽取一个元素

和

并计算对应的可信度

和

可信度的计算公式如下所示：According to the Metropolis-Hastings sampling principle, from N bootstrap measurement sets

Randomly choose two sets from

Then randomly select an element from each of these two sets

and

and calculate the corresponding reliability

and

The formula for calculating reliability is as follows:

是量测预测均值，

代表量测噪声v^J协方差矩阵

的行列式。in

is the mean of the measurement forecast,

represents the measurement noise v ^J covariance matrix

determinant of .

和

的基础上，按照下式计算接受概率：getting credibility

and

On the basis of , the acceptance probability is calculated according to the following formula:

的取值为

和1之间的最小值。the acceptance probability

value of

and the smallest value between 1.

(2) Build an effective measurement set

First, a random number χ that satisfies the uniform distribution U(0,1) is generated. The method for determining the effective measurement set Ω ^J is as follows:

和

对应的接受概率大于等于生成的随机数χ，则将

作为有效量测集合Ω^J中的元素；反之，将

作为有效量测集合Ω^J中的元素。上述过程即为确定有效量测向量的采样过程。In the formula, if the randomly selected elements

and

The corresponding acceptance probability is greater than or equal to the generated random number χ, then the

as an element in the effective measurement set Ω ^J ; otherwise, the

as an element in the effective measurement set Ω ^J. The above process is the sampling process for determining the effective measurement vector.

Repeat L sampling, L=2N, and define L effective measurement vectors in the effective measurement set as z ^J (1), z ^J (2),...z ^J (L), these effective measurement vectors are composed of Valid measurement set Ω ^J .

4. For the effective measurement set obtained in step 3, calculate the weight of each effective measurement vector and the measurement noise matrix, and the specific steps are as follows:

(1) Calculate the consistency distance and consistency matrix

其计算方法如下：Define the confidence distance between any two measurement vectors z ^J (ξ) and z ^J (ζ) in the effective measurement set Ω ^J as

Its calculation method is as follows:

和一致性矩阵Ψ^J，其计算方法如下：Calculate the consistency distance on this basis

and the consistency matrix Ψ ^J , which is calculated as follows:

表示Ω^J中任意两个量测向量之间置信距离的最大值。in,

represents the maximum confidence distance between any two measurement vectors in Ω ^J.

(2) Calculate the weight of the effective measurement vector and the measurement noise matrix

令权向量

此时，

可用来表示相应的有效量测向量被有效量测集合Ω^J中所有元素的总体支持程度。After obtaining the consistency matrix Ψ ^J , calculate the maximum modulus eigenvalue λ of the matrix and the corresponding eigenvector β, and unite the eigenvector β to get

Let the weight vector

at this time,

It can be used to represent the overall support degree of all elements in the effective measurement set Ω ^J for the corresponding effective measurement vector.

其计算方法如下：According to the basic principle of variance calculation in mathematical statistics, the measurement noise matrix corresponding to each element in the effective measurement set is calculated.

Its calculation method is as follows:

5. Use sequential filtering to perform data fusion on the child node, and use the fusion result to correct the motion parameters of the child node at time tk. The specific steps are as _follows :

(1) Establish an airborne distributed POS delivery alignment model

速度误差

位置误差(δL^J、δλ^J、δH^J)、陀螺常值漂移

和加计常值偏置

N为子系统的个数。其中，

分别为第J个子节点的东向、北向、天向失准角；

分别为第J个子节点在n_J系下的东向、北向和天向速度误差；δL^J、δλ^J、δH^J分别为第J个子节点的纬度误差、经度误差、高度误差；

和

分别为第J个子节点处子IMU在其载体坐标系下(b_J系)陀螺常值漂移和加计常值偏置；

分别为第J个子节点处子IMU在b_J系x轴、y轴和z轴上的陀螺仪随机常值漂移；

分别为第J个子节点处子IMU在b_J系x轴、y轴和z轴上的加速度计常值偏置。因此系统状态向量有如下形式：The transfer alignment model adopts the matching method of "position + attitude", and the system state vector is a 15-dimensional state vector, including the platform misalignment angle of the child node J

speed error

Position error (δL ^J , δλ ^J , δH ^J ), gyro constant drift

and add-up constant offset

N is the number of subsystems. in,

are the east, north and sky misalignment angles of the J-th child node;

are the easting, northing and sky velocity errors of the J-th child node under the n _J system, respectively; δL ^J , δλ ^J , and δH ^J are the latitude error, longitude error, and altitude error of the J-th child node;

and

are respectively the gyro constant drift and the accumulating constant offset of the child IMU at the Jth child node in its carrier coordinate system (b _J system);

are the random constant drift of the gyroscope on the x-axis, y-axis and z-axis of the b _J system of the child IMU at the Jth child node, respectively;

are the constant accelerometer offsets of the child IMU at the Jth child node on the x-axis, y-axis and z-axis of the b _J system, respectively. Therefore, the system state vector has the following form:

a) Establishment of the equation of state

The state equation of the system has the following form:

为零均值高斯白噪声，其方差阵

由IMU中陀螺仪和加速度计的噪声水平决定。where system noise

White Gaussian noise with zero mean, its variance matrix

Determined by the noise level of the gyroscope and accelerometer in the IMU.

In the formula, F ^J is the state transition matrix, and its specific expression is as follows:

in,

和

分别为子节点J处子IMU在导航坐标系(n_J系)中的东向比力、北向比力和天向比力分量；

和

分别为子节点J在导航坐标系(n_J系)下的东向速度、北向速度和天向速度；

和

分别为子节点J沿子午圈和卯酉圈的主曲率半径；T为滤波周期。G^J为子节点J的系统噪声矩阵，其具体表达形式如下：Among them, ω _ie represents the angular velocity of the earth's rotation;

and

are the eastward specific force, northward specific force and skyward specific force components of the child IMU at the child node J in the navigation coordinate system (n _J system), respectively;

and

are the easting speed, northing speed and sky speed of the child node J in the navigation coordinate system (n _J system) respectively;

and

are the principal curvature radii of the child node J along the meridian and 卯unitary circles, respectively; T is the filtering period. G ^J is the system noise matrix of the child node J, and its specific expression is as follows:

b) Establishment of measurement equations

The specific expressions of the L effective measurement vectors z ^J (τ) (τ=1,2,...,L) of the child node J are:

分别表示子节点J的第τ个纬度、经度、高度计算值与主节点m经杆臂补偿后的纬度、经度、高度的差值；

和

分别表示子节点J第τ个航向角、俯仰角、横滚角计算值与主节点m对应的航向角、俯仰角、横滚角的差值；in,

respectively represent the difference between the calculated value of the τth latitude, longitude and height of the child node J and the latitude, longitude and height of the main node m after compensation by the lever arm;

and

respectively represent the difference between the calculated value of the τth heading angle, pitch angle and roll angle of the child node J and the heading angle, pitch angle and roll angle corresponding to the main node m;

The measurement equation has the following expression:

有如下形式：The measurement matrix corresponding to z ^J (τ)

Has the following form:

in

且令

表示矩阵T_a中第s行、第t列的元素，即T_a可以表示为：T _a is the attitude matrix of the main node, i.e.

and make

Represents the elements of the s-th row and the t-th column in the matrix T _a , that is, T _a can be expressed as:

Since any sub-node of each sub-node in the airborne distributed POS can obtain the redundant motion parameter information of the other sub-nodes, the redundant measurement vector can be obtained according to the method of step 2. Therefore, for sub-node J , its measurement vector can be expressed as:

对应的量测矩阵可表示为

measurement vector

The corresponding measurement matrix can be expressed as

(2) Data fusion using sequential filtering

和传递对准模型可以得到当前t_k时刻的一步预测估计值

计算公式如下所示：First, the time update is performed, and for the Jth child node, according to the filter estimated value of the previous time t _k-1

and the transfer alignment model can get a one-step prediction estimate at the current time t _k

The calculation formula is as follows:

Among them, Γ ^J is the transition matrix of the discrete system obtained by discretizing the state transition matrix F ^J of the continuous system.

Similarly, according to the estimated mean square error matrix P ^J (k-1|k-1) at the previous moment and the transfer alignment model, the one-step forecast mean square error matrix of the current moment t _k can be obtained. The calculation formula is as follows:

Among them, Ξ ^J is the noise-driven array of the discrete system obtained by discretizing the continuous-system noise-driven array G ^J.

第J个子节点的数据融合公式如下：Then use the effective measurement vector obtained in step (3) to perform measurement update, and record the L measurement vectors of the _Jth child node at the current time tk as

The data fusion formula of the Jth child node is as follows:

进行量测更新时，有：when using

When performing a measurement update, there are:

和

分别为

所对应的量测矩阵和量测噪声矩阵；

为使用

进行量测更新时计算得到的滤波增益矩阵；

表示使用

进行量测更新时计算得到的状态估计值；

为

对应的估计均方误差阵；I为与

维度相同的单位阵。in,

and

respectively

The corresponding measurement matrix and measurement noise matrix;

for use

The filter gain matrix calculated when the measurement update is performed;

indicate use

The state estimate calculated when the measurement update is performed;

for

The corresponding estimated mean square error matrix; I is the

A unit matrix of the same dimensions.

依次进行量测更新时，有：when using

When performing measurement updates in sequence, there are:

和

分别为

对应的量测矩阵和量测噪声矩阵；

为使用

进行量测更新时计算得到的滤波增益矩阵；

表示使用

进行量测更新时计算得到的状态估计值；

为

对应的估计均方误差阵。in,

and

respectively

Corresponding measurement matrix and measurement noise matrix;

for use

The filter gain matrix calculated when the measurement update is performed;

indicate use

The state estimate calculated when the measurement update is performed;

for

The corresponding estimated mean squared error matrix.

作为当前t_k时刻第J个子节点处最终的数据融合结果，即

其中，

包含t_k时刻第J个子节点数据融合后的纬度误差δL^J(k)、经度误差δλ^J(k)和高度误差δH^J(k)，还包含数据融合后在n_J系下的东向速度误差

北向速度误差

和天向速度误差

以及数据融合后的东向失准角

北向失准角

和天向失准角

Execute the above data fusion formula, when the above process runs to τ=L, the obtained

As the final data fusion result at the Jth child node at the current _tk time, that is,

in,

It includes the latitude error δL ^J (k), longitude error δλ ^J (k) and altitude error δH ^J (k) of the Jth child node after data fusion at time t _k , and also includes the eastward velocity under the n _J system after data fusion error

northbound speed error

and sky speed error

and the easting misalignment angle after data fusion

North misalignment angle

and sky misalignment

(3) Motion parameter correction

Use the above fusion result to correct the result of the strapdown solution of the _Jth child node at time tk. The calibration steps are as follows:

a) Position correction

Among them, L ^J (k), λ ^J (k), H ^J (k) represent the latitude, longitude and altitude obtained by the strapdown solution of the J-th child node at time t _k , respectively; L ^J (k) | _new , λ ^J (k)| _new and H ^J ( _k )| _new respectively represent the corrected latitude, longitude and altitude of the Jth child node at time tk.

b) Speed correction

分别代表t_k时刻第J个子节点捷联解算得到的东向速度、北向速度和天向速度；

分别代表t_k时刻第J个子节点校正后的东向速度、北向速度和天向速度。in,

respectively represent the easting velocity, northing velocity and sky velocity obtained by the strapdown solution of the _Jth child node at time tk;

respectively represent the corrected easting speed, northing speed and sky speed of the _Jth child node at time tk.

c) Attitude correction

Calculate the transformation matrix between the navigation coordinate system n _J of the J-th child node at time t _k and the calculated navigation coordinate system n _J ′

为：Corrected transformation matrix

for:

为t_k时刻J个子节点进行捷联解算后得到的姿态矩阵；in,

is the attitude _matrix obtained by the strapdown solution for J child nodes at time tk;

计算t_k时刻第J个子节点的航向角ψ^J(k)|_new、俯仰角θ^J(k)|_new和横滚角γ^J(k)|_new。Using the modified transformation matrix

Calculate the heading angle ψ ^J (k)| _new , the pitch angle θ ^J (k)| _new and the roll angle γ ^J (k)| _new of the J-th child node at time t _k .

6. Repeat steps 1 to 5 for other child nodes that have not undergone data fusion and motion parameter correction until all child nodes complete data fusion and motion parameter correction at time tk. The specific steps are as _follows :

(1) Assume that the number of child nodes that have completed steps 1 to 5 at time t _k is k _N , and the initial value of k _N is 1;

(2) If k _N <N, and N is the number of child nodes, it means that there are still child nodes that have not completed data fusion and motion parameter correction, k _N =k _N +1, repeat for the child nodes numbered k _N the first 5 steps;

(3) If k _N =N, stop step 6, indicating that all child nodes have completed data fusion and motion parameter correction at time t _k .

7. t _k =t _k+1 , and perform steps 1 to 6 until the data fusion of all child nodes at all times and the correction of motion parameters are completed.

Contents that are not described in detail in the specification of the present invention belong to the prior art known to those skilled in the art. For those skilled in the art, according to the idea of the present invention, there will be changes in the specific embodiments and application scope. In conclusion, the contents of this specification should not be construed as limiting the present invention.

Claims (8)

1. A data fusion method for an airborne distributed position and attitude measurement system comprises the following specific steps:

1.1 calculating t_kRedundant position and attitude information at the time child node;

1.2 calculating t_kMeasuring vectors of the sub-nodes at the moment, and generating a bootstrap measurement set by using a measurement bootstrap strategy;

1.3 sampling elements in the bootstrap measurement set by using a Metropolis-Hastings sampling criterion to obtain an effective measurement set of the child nodes;

1.4 calculating the weight and the measurement noise matrix of each effective measurement vector in the effective measurement set;

1.5 data fusion for the child node using sequential filtering, applying the fusion result to the child node t_kCorrecting the motion parameters at the moment;

1.6 repeating steps 1.1 to 1.5 for other sub-nodes without data fusion and motion parameter correction until all sub-nodes finish t_kData fusion and correction of motion parameters at the moment;

1.7 t_k＝t_k+1and executing the steps 1.1 to 1.6 until the data fusion and the correction of the motion parameters of all the child nodes at all the time are completed.

2. The data fusion method of the airborne distributed position and attitude measurement system according to claim 1, characterized in that: in the step 1.1, t is calculated_kThe specific steps of the redundant position and attitude information at the time child node are as follows:

2.1 child node strapdown resolution

The definition of the relevant reference coordinate system includes: e is a terrestrial coordinate system; the navigation coordinate systems of the main node and the sub-nodes are northeast geographic coordinate systems, the navigation coordinate system of the main node is represented by n, and the navigation coordinate system and the calculation navigation coordinate system of the J-th sub-node are respectively represented by n_J and n′_JJ is 1,2, …, N is the number of child nodes; the origin of the carrier coordinate system is the center of gravity of the carrier, the x-axis is right along the transverse axis of the carrier, the y-axis is forward along the longitudinal axis of the carrier, the z-axis is upward along the vertical axis of the carrier, and the coordinate system is fixed on the carrier and is called as a right-front upper carrier coordinate system, and b is used_m、b_JA carrier coordinate system representing the main POS and the jth sub-node, J ═ 1,2, …, N;

each child node is solved through strapdown to obtain t_kThe position [ L ] of each time^J λ^J H^J]^TAttitude [ psi^J θ^Jγ^J]^TAnd velocity

wherein ,L^J、λ^J and H^JRespectively representing the latitude, longitude and altitude of the J-th sub-node; psi^J、θ^J and γ^JRespectively representing a course angle, a pitch angle and a roll angle of the J-th sub-node;

and

respectively representing the east speed, the north speed and the sky speed of the J-th child node;

2.2 acquisition of redundant position and attitude information of child nodes

a) Acquisition of child node redundant location information

T can be obtained by a fiber grating sensor arranged on the wing_kStrain information of any sub-node at any moment, and further displacement of the strain information relative to the initial position of the strain information and any two sub-nodes JAnd J' relative position information DeltaP_J←J′，ΔP_J←J′Represents the difference in latitude, longitude, and altitude between any two child nodes J, J ', J, J' is 1,2, …, N; j is not equal to J';

at t_kAt the moment, the J-th child node is subjected to strapdown solution to obtain a position [ L^J λ^J H^J]^TJ-1, 2, …, N, and the other N-1 child nodes are designated as [ L^J′ λ^J′ H^J′]^TJ '═ 1,2, …, N, J' ≠ J; the redundant location information at the jth child node may be expressed as:

[L^J←J′ λ^J←J′ H^J←J′]^T＝[L^J′ λ^J′ H^J′]^T+ΔP_J←J′(J,J′＝1,2,…,N；J≠J′)

therefore, the J-th sub-node has N calculated values of latitude, longitude and altitude, namely [ L ] in total, in addition to the position information obtained by self strapdown calculation^J λ^J H^J]^T and [L^J←J′ λ^J←J′ H^J←J′]^TJ, J '═ 1,2, …, N, J' ≠ J; position information [ L ] obtained by directly resolving the current child node J in a strapdown manner^J λ^J H^J]^TIs newly recorded as

Remaining N-1 location information [ L^J←J′λ^J←J′ H^J←J′]^TIs newly recorded as

b) Acquisition of child node redundant attitude information

T measured by fiber grating sensor_kThe angular displacement information of each child node at the moment can construct a conversion matrix between carrier coordinate systems of any two child nodes J, J

Carrier coordinate system b at jth sub-node_JIs calculated with the point and navigates the coordinate system n_JThe transformation matrix between the systems is expressed as

The transformation matrix is also called an attitude matrix, and the attitude matrices at the other N-1 sub-nodes are expressed as

The redundant attitude matrix at the jth child node can be represented as:

wherein ,

at t_kAt the moment, N attitude matrixes exist at the J-th sub-node; by utilizing the N attitude matrixes, N attitude information, namely N attitude angles, can be calculated; the specific calculation method of the attitude angle is as follows:

will be provided with

Is recorded as:

wherein ,T_JxyIs a matrix

The x-th row and the y-th column in the formula (I), wherein x is 1,2,3, and y is 1,2, 3; the attitude angle of the J-th sub-node, including the heading angle psi^JAngle of pitch theta^JAnd roll angle γ^JPrincipal value of, i.e.

And

respectively as follows:

course angle psi^JAngle of pitch theta^JAnd roll angle γ^JAre respectively defined as [0, 2 pi ]]、

[-π，+π](ii) a Then, ψ^J、θ^J and γ^JThe true value of (c) can be determined by:

according to the calculation method of the attitude angle, N-1 redundant attitude matrixes at the J-th sub-node can be calculated

Corresponding N-1 redundant pose information, J ═ 1,2, …, N; j is not equal to J'; the attitude matrix of the child node J itself

The calculated attitude angle is recorded again

From the rest N-1 attitude matrices

The attitude angle calculated is recorded as

3. The data fusion method of the airborne distributed position and attitude measurement system according to claim 2, characterized in that: in the step 1.2, t is calculated_kThe method comprises the following specific steps of measuring vectors of the sub-nodes at the moment and generating a bootstrap measurement set by using a measurement bootstrap strategy:

3.1 calculating the measurement vector of the child node

t_kThe latitude, longitude and altitude of the time master node m after lever arm compensation are respectively recorded as L^m、λ^m、H^mThe heading angle, pitch angle and roll angle of the master node m are respectively recorded as psi^m、θ^m、γ^m(ii) a For any sub-node J, J ═ 1,2, …, N, t can be obtained_kN measured vector at time

The expression is as follows:

wherein ,

respectively represents t_kThe ith latitude, longitude and altitude of the time sub-node J are different from the latitude, longitude and altitude of the master node m after lever arm compensation;

respectively represents t_kThe difference value of the ith course angle, the pitch angle and the roll angle of the sub-node J at the moment and the course angle, the pitch angle and the roll angle corresponding to the main node m;

3.2 generating a bootstrap measurement set by using a measurement bootstrap strategy

In the measurement vector

Based on the measurement result, the bootstrap measurement is generated by increasing the disturbance

The method comprises the following specific steps:

wherein ,

representing the original measurement vector

The c-th bootstrap measurement generated on the basis, c ═ 1,2, …, l, l denote bootstrap measurements

L is 5; h^JIs a measurement matrix of the system, x^JIs the state vector of the system, v^JIs measurement noise, and v^JWhite Gaussian noise satisfying zero mean with a covariance of

and v^JHaving the same statistical properties, i.e.

White Gaussian noise satisfying zero mean, its covariance

By the method, NxL bootstrap measurement data can be obtained

Defining a bootstrap measurement set:

wherein

Representing the original measurement vector, i.e.

4. The data fusion method of the airborne distributed position and attitude measurement system according to claim 3, characterized in that: in step 1.3, the Metropolis-Hastings sampling criterion is used to sample elements in the bootstrap measurement set to obtain an effective measurement set of child nodes, and the specific steps are as follows:

4.1 calculate confidence and acceptance probability

From N bootstrap measurements, a set was collected according to Metropolis-Hastings sampling principle

In the random selection of two sets

Then randomly extracting an element from the two sets

And

and calculating corresponding credibility

And

the calculation formula of the reliability is as follows:

wherein

Is a measure of the mean value of the prediction,

representative of the measurement noise v^JCovariance matrix

Determinant of (4);

in obtaining confidence level

And

on the basis of (1), calculating the acceptance probability according to the following formula:

i.e. probability of acceptance

Is taken as

And a minimum value between 1;

4.2 building efficient metrology collections

Firstly, a random number χ satisfying uniform distribution U (0,1) is generated, and an effective measurement set Ω^JThe determination method of (2) is as follows:

in the formula, if the elements are randomly extracted

And

if the corresponding acceptance probability is greater than or equal to the generated random number χ, then it will be

As an effective measurement set omega^JThe elements of (1); otherwise, will

As an effective measurement set omega^JThe elements of (1); the above process is a sampling process for determining an effective measurement vector;

repeating the sampling for L times, wherein L is 2N, and L effective measurement vectors in the effective measurement set are respectively defined as z^J(1),z^J(2),…z^J(L) these effective measurement vectors constitute an effective measurement set omega^J。

5. The data fusion method of the airborne distributed position and attitude measurement system according to claim 4, characterized in that: the specific steps of calculating the weight of each effective measurement vector in the effective measurement set and measuring the noise matrix in step 1.4 are as follows:

5.1 calculate the consistency distance and consistency matrix

Defining an effective metrology set omega^JAny two of the measurement vectors z^J(ξ) and z^JThe confidence distance between (ζ) is

The calculation method is as follows:

calculating a consistency distance based thereon

And the consistency matrix Ψ^JThe calculation method is as follows:

wherein ,

represents omega^JThe maximum value of the confidence distance between any two measurement vectors;

5.2 calculating weights of effective measurement vectors and measurement noise matrix

Obtain the consistency matrix Ψ^JThen, the maximum module eigenvalue lambda and the corresponding eigenvector beta of the matrix are calculated, and the eigenvector beta is unitized to obtain

Order weight vector

At this time, the process of the present invention,

can be used to represent the corresponding valid measurement vector is valid measurement set omega^JThe overall degree of support of all elements in;

according to the basic principle of calculating the variance in the mathematical statistics, the measurement noise matrix corresponding to each element in the effective measurement set is calculated

The calculation method is as follows:

6. the data fusion method of the airborne distributed position and attitude measurement system according to claims 1 and 4, characterized in that: in the step 1.5, sequential filtering is used for carrying out data fusion on the child nodes, and the fusion result is used for the child nodes t_kThe correction of the moment motion parameters comprises the following specific steps:

6.1 building airborne distributed POS transfer alignment model

The transfer alignment model adopts a matching mode of 'position + attitude', and the system state vector is a 15-dimensional state vector and comprises a platform misalignment angle of a child node J

Error in velocity

Position error δ L^J、δλ^J、δH^JConstant drift of gyro

And adding a constant offset

N is the number of subsystems; wherein,

east, north and sky misalignment angles of the J-th child node respectively;

respectively, J-th sub-node is at n_JEast, north and sky speed errors under the system; delta L^J、δλ^J、δH^JRespectively a latitude error, a longitude error and an altitude error of the J-th child node;

and

respectively, the sub IMU at the J-th sub-node is in the carrier coordinate system b_JThe gyro constant drift and the addition constant bias are tied down;

at the J-th child node, respectively, the child IMU is at b_JThe gyroscope on the x axis, the y axis and the z axis is subjected to random constant drift;

at the J-th child node, respectively, the child IMU is at b_JAccelerometer constant offsets in the x, y and z axes; the system state vector thus has the form:

a) establishment of equation of state

The system state equation has the form:

wherein the system noise

Is zero mean white Gaussian noise, its variance matrix

Determined by the noise levels of the gyroscopes and accelerometers in the IMU;

in the formula ,F^JThe state transition matrix is expressed in the following specific form:

wherein ,

，

wherein ,ω_ieRepresenting the rotational angular velocity of the earth;

and

respectively, the sub IMU at the sub node J in the navigation coordinate system n_JEast, north and sky specific force components in the system;

and

respectively a child node J in a navigation coordinate system n_JAn east-direction speed, a north-direction speed and a sky-direction speed under the system;

and

the main curvature radius of the child node J along the meridian circle and the prime circle is respectively; t is a filtering period; g^JThe system noise matrix of the child node J is expressed in the following specific form:

b) establishment of measurement equation

L valid measurement vectors z of child node J^J(τ), τ ═ 1,2, …, and the specific expression for L is:

wherein ,

respectively representing the difference values of the tau-th latitude, longitude and altitude calculation value of the child node J and the latitude, longitude and altitude of the master node m after lever arm compensation;

and

respectively representing differences of the tau-th course angle, the pitch angle and the roll angle calculated value of the child node J and the course angle, the pitch angle and the roll angle corresponding to the master node m;

the measurement equation has the following expression form:

z^J(τ) corresponding measurement matrix

Has the following forms:

wherein ,

T_aas a master node's attitude matrix, i.e.

And order

Representation matrix T_aThe elements in the s-th row and T-th column, i.e. T_aCan be expressed as:

6.2 data fusion Using sequential Filtering

Firstly, time updating is carried out, and for the J-th child node, the time is updated according to the last time t_k-1Filtered estimate of

And transferring the alignment model to obtain the current t_kOne-step predictive estimate of time of day

The calculation formula is as follows:

wherein ,Γ^JIs to transfer the state of the continuous system to the matrix F^JObtaining a transfer matrix of a discrete system after discretization;

in the same way, according to the estimated mean square error array P of the last moment^J(k-1| k-1) and the transfer alignment model may be derived for the current t_kOne step predictive mean square of timeError matrix, the calculation formula is as follows:

wherein ,Ξ^JIs to drive the noise of the continuous system to the array G^JDiscretizing to obtain a noise driving array of a discrete system;

then, the effective measurement vector obtained in the step 1.3 is used for measurement updating, and the J-th sub-node is located at the current t_kThe L measured direction at that time are respectively recorded as

The data fusion formula of the J-th child node is as follows:

when in use

When the measurement is updated, the following steps are carried out:

wherein ,

and

are respectively as

The corresponding measurement matrix and the measurement noise matrix;

to use

Calculating a filter gain matrix during measurement updating;

indicating use of

Calculating to obtain a state estimation value during measurement updating; p₁ ^J(k | k) is

A corresponding estimated mean square error matrix; i is and

unit arrays with the same dimension;

when in use

When carrying out measurement update in proper order, have:

wherein ,

and

are respectively as

The corresponding measurement matrix and the measurement noise matrix;

to use

Calculating a filter gain matrix during measurement updating;

indicating use of

Calculating to obtain a state estimation value during measurement updating;

is composed of

A corresponding estimated mean square error matrix;

executing the data fusion formula, and when the flow runs to the condition that tau is equal to L, obtaining the result

As the current t_kThe final data fusion result at the J-th child node at time instant, i.e.

wherein ,

containing t_kLatitude error delta L after data fusion of J-th sub-node at moment^J(k) Longitude error δ λ^J(k) And height error deltaH^J(k) And also includes n after data fusion_JEast-down velocity error under tie

Error in north direction velocity

And speed error in the sky

And east misalignment angle after data fusion

Angle of north misalignment

And angle of vertical misalignment

6.3 motion parameter correction

Using the above fusion result to t_kAnd correcting the result of the J-th sub-node strapdown calculation at the moment, wherein the correction steps are as follows:

a) position correction

wherein ,L^J(k)、λ^J(k)、H^J(k) Respectively represents t_kThe J-th sub-node at the moment is obtained by the quick-link solutionLatitude, longitude and altitude to; l is^J(k)|_new、λ^J(k)|_new、H^J(k)|_newRespectively represents t_kLatitude, longitude and altitude after J-th child node correction at the moment;

b) velocity correction

wherein ,

respectively represents t_kAt the moment, the J-th child node is subjected to strapdown calculation to obtain an east-direction speed, a north-direction speed and a sky-direction speed;

respectively represents t_kThe east speed, the north speed and the sky speed after the J-th sub-node correction at the moment;

c) attitude correction

Calculating t_kNavigation coordinate system n of J-th sub-node at moment_JAnd calculating a navigation coordinate system n_J' conversion matrix between

Modified transformation matrix

Comprises the following steps:

wherein ,

is t_kPerforming strapdown calculation on the J-th sub-node at the moment to obtain an attitude matrix;

using modified transformation matrices

Calculating t_kCourse angle psi of sub-node J at time^J(k)|_newAngle of pitch theta^J(k)|_newAnd roll angle γ^J(k)|_new。

7. The data fusion method of the airborne distributed position and attitude measurement system according to claim 1, characterized in that: in the step 1.6, the steps 1.1 to 1.5 are repeated for other child nodes which are not subjected to data fusion and motion parameter correction until all the child nodes finish t_kThe method comprises the following steps of time data fusion and motion parameter correction:

7.1 hypothesis t_kThe child node number k at which steps 1.1 to 1.5 have been completed_N，k_NIs 1;

7.2 if k_NN is the number of the child nodes, which indicates that the child nodes still do not finish data fusion and correction of motion parameters, k_N＝k_N+1, pair number k_NThe child node of (1) repeats steps 1.1 to 1.5;

7.3 if k_NN, then step 1.6 is stopped, indicating that all child nodes have completed t_kAnd (4) data fusion of time and correction of motion parameters.

8. The data fusion method of the airborne distributed position and attitude measurement system according to claim 1, characterized in that: in the step 1.7, t is enabled_k＝t_k+1And executing the steps 1.1 to 1.6 until the data fusion and the correction of the motion parameters of all the child nodes at all the time are completed.

Images