CN113188566B - Airborne distributed POS data fusion method - Google Patents

Airborne distributed POS data fusion method Download PDF

Info

Publication number
CN113188566B
CN113188566B CN202110309944.7A CN202110309944A CN113188566B CN 113188566 B CN113188566 B CN 113188566B CN 202110309944 A CN202110309944 A CN 202110309944A CN 113188566 B CN113188566 B CN 113188566B
Authority
CN
China
Prior art keywords
measurement
node
matrix
sub
child node
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110309944.7A
Other languages
Chinese (zh)
Other versions
CN113188566A (en
Inventor
宫晓琳
丁孝双
孙一弘
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN202110309944.7A priority Critical patent/CN113188566B/en
Publication of CN113188566A publication Critical patent/CN113188566A/en
Application granted granted Critical
Publication of CN113188566B publication Critical patent/CN113188566B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C23/00—Combined instruments indicating more than one navigational value, e.g. for aircraft; Combined measuring devices for measuring two or more variables of movement, e.g. distance, speed or acceleration
    • G01C23/005—Flight directors

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Navigation (AREA)

Abstract

The disclosure provides an airborne distributed POS data fusion method, which specifically comprises the following steps: calculating redundant position and posture information of the position where the child node is located; calculating a measurement vector of the sub-node, and generating a bootstrap measurement set by using a measurement bootstrap strategy; sampling the elements of the bootstrap measurement set by using a Metropolis-Hastings sampling criterion to obtain an effective measurement set; calculating the weight and the measurement noise matrix of each effective measurement vector in the effective measurement set; performing data fusion on the child nodes by using sequential filtering, and correcting the motion parameters of the child nodes by using fusion results; and repeating the steps respectively for other child nodes which do not perform data fusion. The method utilizes the motion information of all the sub-nodes, reduces the adverse effect of measurement noise uncertainty on transmission alignment estimation precision, and uses the data fusion result for correcting the motion parameters, thereby improving the estimation precision of the motion parameters of all the sub-nodes.

Description

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 for improving the estimation accuracy of airborne distributed POS sub-node motion parameters.
Background
An on-board position and attitude measurement system (Position and Orientation System, POS) is the primary means for acquiring vehicle position, velocity and attitude information. By using the information, remote sensing data can be corrected, and imaging quality is improved. Therefore, POS is often used in the fields of aerial remote sensing imaging, aerial photogrammetry, and the like. For example, the POS is used for assisting remote sensing loads such as synthetic aperture radar (Synthetic Aperture Rada, SAR), laser radar, multispectral scanner, digital camera, large-field infrared scanner and the like so as to meet the requirement of high-precision motion imaging.
With the development of the flight platform technology, the same carrier is provided with a plurality of or a plurality of remote sensing devices to realize synchronous observation of a plurality of observation windows, and the SAR, the visible light camera, the imaging spectrometer and the laser radar typically work simultaneously, and the array antenna SAR with a plurality of antennas working simultaneously. Such a multi-tasking remote sensing mode integrating multiple or multiple remote sensing loads has become an important development direction for aerial remote sensing systems. Because a plurality of observation loads are installed at different positions of the carrier, the requirement of measuring the position and posture information of multiple nodes cannot be met by using a traditional single POS. Therefore, in actual use, a plurality of inertial measurement units (Inertial Measurement Unit, IMU) are installed on the carrier to form a distributed position and posture measurement system, i.e. a distributed POS, so as to realize measurement of multi-node motion parameters.
The distributed POS is usually composed of one or a small number of high-precision main POS and a plurality of medium/low-precision sub-IMUs, and the two parts are respectively arranged near remote sensing loads on two sides of a machine body or a wing to form a distributed multi-sensor system. The location of the main POS is called a main node, and the locations of the sub-IMUs are called sub-nodes.
In order to realize accurate measurement of motion information at each sub-node in the distributed POS, the sub-node needs to correct the self strapdown calculation result according to the motion information such as high-precision position, speed, gesture and the like provided by the main POS, and uses the result as the final motion parameter of the point, which is called transfer alignment. Because of the differences of factors such as machine body deformation, lever arms and the like at each sub-node, after one-to-one transmission alignment of the main node to each sub-node is completed, the precision of the motion parameters at each sub-node is different and has high or low. Particularly, the deflection deformation condition of the sub-node far from the center of the machine body is the most complex, and the transmission alignment precision is lower. If each child node can use the motion information of the other nodes to perform data fusion, the motion parameter precision of each child node is further improved.
The current data fusion algorithms for multi-sensor systems can be divided into two categories: the first is centralized fusion, which completes the fusion calculation of the measurement data of all sensors at one time, and proves to be the least information loss and the globally optimal at the same time, but for multi-sensor systems such as airborne distributed POS, the multiplication and inversion of the high-dimensional matrix are involved at the fusion time, and if the fusion method is adopted in the airborne distributed POS, the disadvantages of large calculation load, low real-time performance and low fault tolerance are caused. The second is distributed fusion, which filters the original data information of a plurality of sensors respectively through a plurality of filters and then performs centralized processing through a fusion center. For example, the patent with the patent number of CN201810153913.5 adopts a distributed fusion method, and the inverse of the covariance matrix obtained by sub IMU transfer alignment is used as a weight matrix for data fusion, thereby respectively providing a specific method for fusing position, speed and gesture information. However, the method needs to calculate a weight matrix with the total number of sub nodes for data fusion, the calculation of the weight matrix is complex, the inversion of the matrix is involved, the problem that the matrix falls into singular exists, and the total calculation amount is large. Meanwhile, the method does not consider the relevance of flexural motion among the child nodes, and the accuracy of the data fusion result is also adversely affected.
Common structures in distributed fusion are federal filtering and sequential filtering. Both distributed fusion structures are based on kalman filtering. The federal filtering belongs to one of decentralized filtering methods, and the structure of the federal filtering comprises a plurality of sub-filters and a main filter, and meanwhile, information distribution between the sub-filters and the main filter is realized by applying an information distribution principle. The federal filtering has the advantages of flexible design and good fault tolerance, and can be widely applied to the integrated navigation system. However, the federal filtering aims at the fusion of sensor data with complementarity of a plurality of measurement functions installed at a certain point, and each sub-node in the airborne distributed POS is only provided with one sub-IMU, so that a public reference system and a plurality of subsystems required during the use of the federal filtering cannot be formed for any sub-node, namely, a federal filtering structure cannot be formed, and therefore, the federal filtering is not suitable for the data fusion of the airborne distributed POS. The sequential filtering is to sequentially perform Kalman filtering on multiple sets of measurement data, and the final result is used as a data fusion result. For the airborne distributed POS, the deformation data measured by the fiber bragg grating sensor can be utilized to convert the motion parameter information of each sub-node to any other sub-node, namely, any sub-node can obtain the redundant motion parameter information converted to the position by other sub-nodes, so that a plurality of measurement vectors are formed, the possibility is provided for carrying out filtering estimation on the plurality of measurement vectors of each sub-node by adopting sequential filtering, and the sequential filtering has the advantages of smaller calculated amount and simple structure.
Disclosure of Invention
The technical solution of the invention is as follows: the method for fusing the airborne distributed POS data overcomes the defects of the prior art, and improves the estimation precision of the motion parameters of all sub-nodes of the airborne distributed POS.
The technical scheme of the invention is as follows: an airborne distributed POS data fusion method. The method comprises the following specific steps:
(1) Calculating t k Redundant position and posture information at the time sub-node;
(2) Calculating t k Measuring vectors of the time child nodes and generating a bootstrap measuring set by utilizing a measurement bootstrap strategy;
(3) Sampling elements in the bootstrap measurement set by using a Metropolis-Hastings sampling criterion to obtain an effective measurement set of the child node;
(4) Calculating the weight and the measurement noise matrix of each effective measurement vector in the effective measurement set;
(5) Data fusion is carried out on child nodes by using sequential filtering, and fusion results are used for the child node t k Correcting moment motion parameters;
(6) Repeating the steps (1) to (5) for other child nodes which do not perform data fusion and motion parameter correction until all child nodes complete t k Data fusion at moment and correction of motion parameters;
(7)t k ＝t k+1 and (3) executing the steps (1) to (6) until the data fusion and the correction of the motion parameters of all the child nodes at all the moments are completed.
Obtaining t in the step (1) k The redundant position and gesture information at the time sub-node comprises the following specific steps:
1) Sub-node strapdown calculation
The definition of the relevant reference coordinate system includes: e is the earth 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 of the J sub node and the calculated navigation coordinate system are respectively represented by n J and n′J J=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 rightward 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, the coordinate system is fixed on the carrier and is called the right front upper carrier coordinate system, and b is used m 、b J The carrier coordinate systems representing the master node and the J-th child node, j=1, 2, …, N, respectively.
Each child node obtains t after strapdown calculation k The respective positions [ L ] of the moments J λ J H J ] T Posture [ phi ] J θ J γ J ] T Sum speed of wherein ,LJ 、λ J and HJ Respectively representing the latitude, longitude and altitude of the J-th child node; psi phi type J 、θ J and γJ Respectively representing the course angle, the pitch angle and the roll angle of the J-th sub-node; /> and />Respectively represent the J-th sub-node at n J The east, north and sky speeds of the tie down.
2) Sub-node redundant position and attitude information acquisition
a) Acquisition of child node redundant location information
By means of fiber bragg grating sensors mounted on the wing, t can be obtained k Strain information at any sub-node at any moment is obtained to obtain displacement relative to the initial positionAnd relative position information Δp between any two child nodes J, J ' (J, J ' =1, 2, …, N; j+.j ') J←J′ ，ΔP J←J′ Representing the difference in latitude, longitude, and altitude between any two child nodes J, J'.
At t k At moment, the position obtained by strapdown calculation of the J-th child node is [ L ] J λ J H J ] T J=1, 2, …, N, the positions of the other N-1 child nodes after strapdown resolution are denoted as L J′ λ J′ H J′ ] T (J '=1, 2, …, N, J' noteqj). The redundant location information at the J-th 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′)
thus, adding the position information obtained by self strapdown calculation, the J-th child node has N latitude, longitude and altitude calculated values, namely [ L ] J λ J H J ] T and [LJ←J′ λ J←J′ H J←J′ ] T (J, J '=1, 2, …, N, J' noteqj). Position information [ L ] obtained by direct strapdown calculation of current child node J J λ J H J ] T Re-recorded asThe remaining N-1 position information [ L J←J′ λ J←J′ H J←J′ ] T (J, J '=1, 2, …, N; j+.j') is recommenced as
b) Sub-node redundant attitude information acquisition
T measured by fiber bragg grating sensor k Angular displacement information of each child node at moment can be constructedThe transformation matrix between the carrier coordinate systems of two child nodes J, J ' (J, J ' =1, 2, …, N; j+.j ') is intendedCarrier coordinate system at the J-th child node (b J System) and the point navigation coordinate system (n J The system) is expressed as a transformation matrix (posture matrix)The gesture matrix at the other N-1 child nodes is denoted +.>The redundant pose matrix at the J-th child node may be expressed as:
wherein ,
at t k At time, there are N gesture matrices at the J-th child node. With the N gesture matrices, N gesture angles can be calculated. The specific calculation method of the attitude angle is as follows:
will beThe method is characterized by comprising the following steps: />
wherein ,TJxy Is a matrixX=1, 2,3, y=1, 2,3; the pose angle of the J-th child node includes the heading angle ψ J Pitch angle theta J And roll angle gamma J Main value of (i.e.)> and />The method comprises the following steps of:
heading angle psi J Pitch angle theta J And roll angle gamma J The value ranges of (2) are respectively defined as [0,2 pi ]]、 [-π，+π]The method comprises the steps of carrying out a first treatment on the surface of the Then, ψ J 、θ J and γJ The true value of (2) can be determined by:
according to the attitude angle calculating method, N-1 redundant attitude matrixes at the J-th sub-node can be calculatedCorresponding to N-1 redundant attitude angles. The attitude matrix to be defined by the child node J itself The calculated attitude angle is remembered as +.>From the remaining N-1 gesture matrices +.>The calculated attitude angle is recorded as +.>
Calculating t in the step (2) k The measuring vector of the time child node and the bootstrap measuring strategy are utilized to generate the bootstrap measuring set, which comprises the following specific steps:
1) Calculating a measurement vector for a sub-node
t k The latitude, longitude and altitude of the moment main node m after the compensation of the lever arm are respectively recorded as L m 、λ m 、H m The course angle, pitch angle and roll angle of the main node m are respectively recorded as psi m 、θ m 、γ m . For any child junction J (j=1, 2, …, N), t can be derived k N measurement vectors of timeThe expression is as follows:
wherein ,respectively represent t k The difference value between the ith latitude, longitude and altitude of the moment sub-node J and the latitude, longitude and altitude of the main node m after the compensation of the lever arm; />Respectively represent t k And the difference value of the i-th course angle, pitch angle and roll angle of the moment sub-node J and the course angle, pitch angle and roll angle corresponding to the main node m.
2) Generating bootstrap measurement sets using measurement bootstrap strategy
In measuring vectorOn the basis of (1) generating bootstrap measures by adding disturbances +.>The method comprises the following specific steps:
wherein ,representing the original measurement vector->Based on the c-th bootstrap measurement, c=1, 2, …, l, l represents the bootstrap measurement +. >I=5; h J Is the measurement matrix of the system, x J Is a state vector of the system, v J Is the measurement noise and v J White gaussian noise satisfying zero mean with covariance +.> and vJ Has the same statistical properties, i.e.)>Gaussian white noise satisfying zero mean, covariance +>
By the method, N×l bootstrap measurement data can be obtainedDefining a bootstrap measurement set: /> wherein />Representing the original measurement vector, i.e
In the step (3), elements in the bootstrap measurement set are sampled by using a Metropolis-Hastings sampling rule to obtain an effective measurement set of the child node, which comprises the following specific steps:
1) Calculating confidence and acceptance probability
From N bootstrap measurement sets according to Metropolis-Hastings sampling principleTwo sets are selected at will->Randomly extracting an element from the two sets and />And calculates the corresponding confidence level +.> and />The calculation formula of the credibility is as follows:
wherein Is the measurement prediction mean value->Representing measurement noise v J Covariance matrix->Is a determinant of (2).
In obtaining the credibility and />On the basis of (a), the acceptance probability is calculated according to the following formula:
i.e. probability of acceptanceThe value of (2) is +.>And a minimum value between 1.
2) Constructing an active measurement set
Firstly, a random number χ satisfying the uniform distribution U (0, 1) is generated, and an effective measurement set Ω is obtained J The determination method of (2) is as follows:
in the formula, if randomly extracted elements and />The corresponding acceptance probability is larger than or equal to the generated random number χ, thenAs an effective measurement set Ω J Elements of (a) and (b); conversely, will->As an effective measurement set Ω J Is a component of the group. The above process is the sampling process for determining the effective measurement vector.
Repeating the sampling for L times, wherein L=2N, and defining L effective measurement vectors in the effective measurement set as z respectively J (1),z J (2),…z J (L) the effective measurement vectors form an effective measurement set Ω J 。
The specific steps of calculating the weight and the measurement noise matrix of each effective measurement vector in the effective measurement set in the step (4) are as follows:
1) Calculating a coherence distance and a coherence matrix
Defining an active measurement set Ω J Any two measurement vectors z J(ξ) and zJ Confidence distance between (ζ) isThe calculation method is as follows:
on the basis of which the consistency distance is calculatedAnd a consistency matrix ψ J The calculation method is as follows:
wherein ,representing omega J The maximum value of the confidence distance between any two measurement vectors.
2) Calculating weights of effective measurement vectors and measurement noise matrix
Obtaining a consistency matrix ψ J Then, the maximum modulus eigenvalue lambda and the corresponding eigenvector beta of the matrix are calculated, and beta is unitized to obtainLet weight vector->At this time, a->Can be used to represent the corresponding valid measurement vector is effectively measured by the set Ω J The overall degree of support for all elements in a system.
According to the basic principle of variance calculation in mathematical statistics, calculating to obtain a measurement noise matrix corresponding to each element in the effective measurement setThe calculation method is as follows:
in the step (5), the sub-node is subjected to data fusion by using sequential filtering, and the fusion result is used for the sub-node t k The correction of the moment motion parameters comprises the following specific steps:
1) Establishing an onboard distributed POS transfer alignment model
The transfer alignment model adopts a matching mode of 'position and gesture', and the system state vector is a 15-dimensional state vector and comprises a platform misalignment angle of a sub-node JSpeed error->Position error (δL) J 、δλ J 、δH J ) Gyro constant drift->And add constant bias ∈ ->N is the number of subsystems. Wherein (1)> The east, north and sky misalignment angles of the J-th sub-node are respectively;respectively the J-th sub-node is at n J The tangential, north and sky speed errors; δL (delta L) J 、δλ J 、δH J Latitude error, longitude error and altitude error of the J-th child node respectively; / >Andthe sub-IMUs at the J-th sub-node respectively (b) under its carrier coordinate system J System) gyro constant drift and add constant bias; />sub-IMU at the J-th sub-node at b J Random constant drift of gyroscopes on the x-axis, y-axis and z-axis; />sub-IMU at the J-th sub-node at b J Constant offsets of accelerometers in x-axis, y-axis and z-axis are tied. The system state vector is thus of the form:
a) Establishment of the equation of state
The system state equation has the following form:
in which the system noiseZero mean Gaussian white noise with variance matrix +.>Determined by the noise level of the gyroscopes and accelerometers in the IMU. F (F) J The state transition matrix is specifically expressed as follows:
wherein ,
wherein ,ωie Indicating the rotational angular velocity of the earth; and />The sub-IMUs at sub-nodes J are respectively in a navigation coordinate system (n J Tie) east to north to force and natural to force components; /> and />Respectively, the child nodes J are arranged in a navigation coordinate system (n J Tie) east, north and sky speeds; /> and />The main curvature radius of the subnode J along the meridian circle and the mortise unitary circle is respectively; t is a filtering period; g J The system noise matrix is a subnode J, and the specific expression form is as follows:
b) Establishment of measurement equation
L effective measurement vectors z of sub-node J J (τ) (τ=1, 2, …, L) is expressed as:
wherein ,respectively representing the difference value of the tau latitude, longitude and altitude calculated value of the child node J and the latitude, longitude and altitude of the main node m after the compensation of the lever arm; /> and />Respectively representing the difference value of the tau-th course angle, pitch angle and roll angle calculated value of the subnode J and the course angle, pitch angle and roll angle corresponding to the main node m;
expression of the measurement equation:
z J (tau) corresponding measurement matrix
wherein ,
T a is the gesture matrix of the main node, i.eAnd let->Representation matrix T a Elements of row s, column T, i.e. T a Can be expressed as:
2) Data fusion using sequential filtering
Firstly, time updating is carried out, and aiming at the J-th child node, according to the last time t k-1 Filtered estimates of (2)And the transfer alignment model can obtain the current t k One-step predictive estimate of time instant +.>The calculation formula is as follows:
wherein ,ΓJ Is to transfer the state of the continuous system into a matrix F J And discretizing to obtain a transfer matrix of the discrete system.
Similarly, the mean square error matrix P is estimated according to the last moment J (k-1|k-1) and the transfer alignment model can obtain the current t k The calculation formula of the one-step prediction mean square error matrix at the moment is as follows:
wherein ,ΞJ Is to drive the noise of continuous system into array G J Discretizing to obtain the final productThe noise of the dispersion system drives the array.
Then, the effective measurement vector obtained in the step (3) is utilized to carry out measurement update, and the J-th sub-node is at the current t k The L measurement vectors at the moment are respectively recorded asThe data fusion formula for the J-th child node is as follows:
when in useWhen the measurement is updated, there are:
wherein , and />Respectively->The corresponding measuring matrix and measuring noise matrix; />For use +.>A filtering gain matrix is calculated when measurement updating is carried out; />Indicate use +.>Calculating a state estimation value when measuring and updating; />Is->A corresponding estimated mean square error matrix; i is AND->A unit array of the same dimensions.
When in useWhen the measurement is updated in sequence, there are:
wherein , and />Respectively->The corresponding measuring matrix and measuring noise matrix; />For use +.>A filtering gain matrix is calculated when measurement updating is carried out; />Indicate use +.>Calculating a state estimation value when measuring and updating; />Is->A corresponding estimated mean square error matrix.
Executing the data fusion formula, and when the flow runs to tau=L, obtainingAs the current t k The final data fusion result at the J-th child node at time instant, i.e. +. > wherein ,/>Comprises t k Latitude error delta L after data fusion of J-th child node at moment J (k) Longitude error delta lambda J (k) And a height error δH J (k) Also comprises the following step n after data fusion J East speed error of the system->North speed error->And the error of the tangential velocity->And east misalignment angle after data fusion +.>North misalignment angle->And the angle of misalignment of the sky->
3) Motion parameter correction
The fusion result is used for t k And correcting the result of the J-th child node strapdown calculation at the moment. The correction 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)
wherein ,LJ (k)、λ J (k)、H J (k) Respectively represent t k The J-th child node strapdown at the moment calculates the latitude, longitude and altitude; l (L) J (k)| new 、λ J (k)| new 、H J (k)| new Respectively represent t k Latitude, longitude, and altitude after the moment J child node correction.
b) Speed correction
wherein ,respectively represent t k The east speed, the north speed and the sky speed obtained by the moment J child node strapdown calculation; />Respectively represent t k The corrected east, north and sky speeds of the moment J child node.
c) Posture correction
Calculating t k Navigation coordinate system n of J-th sub-node at moment J And calculating a navigation coordinate system n J ' conversion matrix
Modified transformation matrixThe method comprises the following steps: /> wherein ,/>At t k The attitude matrix is obtained after strapdown calculation is carried out on the J child nodes at the moment;
Using the modified conversion matrixCalculating t k Heading angle psi of moment J child node J (k)| new Pitch angle theta J (k)| new And roll angle gamma J (k)| new 。
Repeating steps (1) to (5) for other sub-nodes which do not perform data fusion and motion parameter correction in step (6) until all sub-nodes complete t k The method comprises the following specific steps of data fusion at moment and correction of motion parameters:
1) Let t be k The child node number k at which steps (1) to (5) have been completed N ，k N The initial value is 1;
2) If k N N is the number of the child nodes, the fact that the child nodes still do not finish data fusion and correction of motion parameters is indicated, and k is calculated N ＝k N +1, pair number k N Repeating steps (1) to (5) for child nodes of (a);
3) If k N If N, stop step (6) to indicate that all child nodes have completed t k Data fusion at the moment and correction of motion parameters.
Compared with the prior art, the invention has the advantages that:
the measurement bootstrap strategy and the sequential filtering are combined around the aim of measuring the overall accuracy of the airborne distributed POS multi-node motion parameters. The motion parameter information of all the sub-nodes is obtained by obtaining measurement vectors of which each sub-node is closer to a true value, then sequentially filtering the measurement vectors, and using the obtained data fusion result for correcting the motion parameters of the sub-nodes. The method fully utilizes the motion parameter information of all the sub-nodes, overcomes the adverse effect of measurement noise uncertainty on transmission alignment estimation accuracy, and finally improves the estimation accuracy of the motion parameters of the sub-nodes.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a block diagram of a data fusion method based on metrology bootstrapping and sequential filtering in accordance with the present invention.
Detailed Description
As shown in fig. 1, the specific method of the present invention is implemented as follows:
1. calculating t k Redundant position and posture information at the time sub-node is calculated as follows:
(1) Sub-node strapdown calculation
The definition of the relevant reference coordinate system includes: e is the earth coordinate system; navigation at a master node and child nodesThe coordinate systems are the geographic coordinate system of northeast and north days, the navigation coordinate system of the main node is represented by n, and the navigation coordinate system of the J-th sub node and the calculated navigation coordinate system are respectively represented by n J and n′J J=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 rightward 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, the coordinate system is fixed on the carrier and is called the right front upper carrier coordinate system, and b is used m 、b J The carrier coordinate systems representing the master node and the J-th child node, j=1, 2, …, N, respectively.
Each child node obtains t after strapdown calculation k The respective positions [ L ] of the moments J λ J H J ] T Posture [ phi ] J θ J γ J ] T Sum speed of wherein ,LJ 、λ J and HJ Respectively representing the latitude, longitude and altitude of the J-th child node; psi phi type J 、θ J and γJ Respectively representing the course angle, the pitch angle and the roll angle of the J-th sub-node; /> and />Respectively represent the J-th sub-node at n J The east, north and sky speeds of the tie down.
(2) Sub-node redundant position and attitude information acquisition
a) Acquisition of child node redundant location information
By means of fiber bragg grating sensors mounted on the wing, t can be obtained k The strain information at any sub-node at a moment in time is then used to obtain the displacement of the sub-node relative to its initial position and the relative position information Δp between any two sub-nodes J, J ' (J, J ' =1, 2, …, N; j+.j ') J←J′ ，ΔP J←J′ Representing between any two child nodes J, JLatitude, longitude, and altitude.
At t k At moment, the position obtained by strapdown calculation of the J-th child node is [ L ] J λ J H J ] T J=1, 2, …, N, the positions of the other N-1 child nodes after strapdown resolution are denoted as L J′ λ J′ H J′ ] T (J '=1, 2, …, N, J' noteqj). The redundant location information at the J-th 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′)
thus, adding the position information obtained by self strapdown calculation, the J-th child node has N latitude, longitude and altitude calculated values, namely [ L ] J λ J H J ] T and [LJ←J′ λ J←J′ H J←J′ ] T (J, J '=1, 2, …, N, J' noteqj). Position information [ L ] obtained by direct strapdown calculation of current child node J J λ J H J ] T Re-recorded asThe remaining N-1 position information [ L J←J′ λ J←J′ H J←J′ ] T (J, J '=1, 2, …, N; j+.j') is recommenced as
b) Sub-node redundant attitude information acquisition
T measured by fiber bragg grating sensor k The angular displacement information of each sub-node at the moment can construct a transformation matrix between the carrier coordinate systems of any two sub-nodes J, J ' (J, J ' =1, 2, …, N; J +.J ')Carrier coordinate system at the J-th child node (b J System) and the point navigation coordinate system (n J The system) is expressed as a transformation matrix (posture matrix)The gesture matrix at the other N-1 child nodes is denoted +.>The redundant pose matrix at the J-th child node may be expressed as:
wherein ,
at t k At time, there are N gesture matrices at the J-th child node. N pieces of gesture information, namely N gesture angles, can be calculated by using the N gesture matrixes. The specific calculation method of the attitude angle is as follows:
will beThe method is characterized by comprising the following steps:
wherein ,TJxy Is a matrixX=1, 2,3, y=1, 2,3; the pose angle of the J-th child node includes the heading angle ψ J Pitch angle theta J And roll angle gamma J Main value of (i.e.)> and />The method comprises the following steps of:
heading angle psi J Pitch angle theta J And roll angle gamma J The value ranges of (2) are respectively defined as [0,2 pi ]]、 [-π，+π]The method comprises the steps of carrying out a first treatment on the surface of the Then, ψ J 、θ J and γJ The true value of (2) can be determined by:
according to the attitude angle calculating method, N-1 redundant attitude matrixes at the J-th sub-node can be calculatedCorresponding to N-1 redundant attitude angles. The posture matrix to be defined by child node J itself +.>The calculated attitude angle is remembered as +.>From the remaining N-1 gesture matrices +.>The calculated attitude angle is recorded as +.>
2. Calculating t k The measurement vector of the time child node and the bootstrap measurement strategy are utilized to generate a bootstrap measurement set, and the specific steps are as follows:
(1) Calculating a measurement vector for a sub-node
t k The latitude, longitude and altitude of the moment main node m after the compensation of the lever arm are respectively recorded as L m 、λ m 、H m The course angle, pitch angle and roll angle of the main node m are respectively recorded as psi m 、θ m 、γ m . For any child junction J (j=1, 2, …, N), t can be derived k N measurement vectors of timeThe expression is as follows:
wherein ,respectively represent t k The difference value between the ith latitude, longitude and altitude of the moment sub-node J and the latitude, longitude and altitude of the main node m after the compensation of the lever arm; />Respectively represent t k And the difference value of the i-th course angle, pitch angle and roll angle of the moment sub-node J and the course angle, pitch angle and roll angle corresponding to the main node m.
(2) Generating bootstrap measurement sets using measurement bootstrap strategy
In measuring vectorOn the basis of (1) generating bootstrap measures by adding disturbances +.>The method comprises the following specific steps:
wherein ,representing the original measurement vector->Based on the c-th bootstrap measurement, c=1, 2, …, l, l represents the bootstrap measurement +.>I=5; h J Is the measurement matrix of the system, x J Is a state vector of the system, v J Is the measurement noise and v J White gaussian noise satisfying zero mean with covariance +.> and vJ Has the same statistical properties, i.e.)>Gaussian white noise satisfying zero mean, covariance +>
By the method, N×l bootstrap measurement data can be obtainedDefining a bootstrap measurement set: /> wherein />Representing the original measurement vector, i.e
3. Sampling elements in the bootstrap measurement set obtained in the step 2 by using a Metropolis-Hastings sampling criterion to obtain an effective measurement set of the child node, wherein the method comprises the following specific steps of:
(1) Calculating confidence and acceptance probability
From N bootstrap measurement sets according to Metropolis-Hastings sampling principleTwo sets are selected at will->Randomly extracting an element from the two sets and />And calculates the corresponding confidence level +.> and />The calculation formula of the credibility is as follows:
wherein Is the measurement prediction mean value->Representing measurement noise v J Covariance matrix->Is a determinant of (2).
In obtaining the credibility and />On the basis of (a), the acceptance probability is calculated according to the following formula:
i.e. probability of acceptanceThe value of (2) is +.>And a minimum value between 1.
(2) Constructing an active measurement set
Firstly, a random number χ satisfying the uniform distribution U (0, 1) is generated, and an effective measurement set Ω is obtained J The determination method of (2) is as follows:
in the formula, if randomly extracted elements and />The corresponding acceptance probability is larger than or equal to the generated random number χ, thenAs an effective measurement set Ω J Elements of (a) and (b); conversely, will->As an effective measurement set Ω J Is a component of the group. The above process is the sampling process for determining the effective measurement vector.
Repeating the sampling for L times, wherein L=2N, and defining L effective measurement vectors in the effective measurement set as z respectively J (1),z J (2),…z J (L) the effective measurement vectors form an effective measurement set Ω J 。
4. Aiming at the effective measurement set obtained in the step 3, the weight and the measurement noise matrix of each effective measurement vector are calculated, and the specific steps are as follows:
(1) Calculating a coherence distance and a coherence matrix
Defining an active measurement set Ω J Any two measurement vectors z J(ξ) and zJ Confidence distance between (ζ) is The calculation method is as follows:
on the basis of which the consistency distance is calculatedAnd a consistency matrix ψ J The calculation method is as follows:
wherein ,representing omega J The maximum value of the confidence distance between any two measurement vectors.
(2) Calculating weights of effective measurement vectors and measurement noise matrix
Obtaining a consistency matrix ψ J Then, the maximum modulus eigenvalue lambda and the corresponding eigenvector beta of the matrix are calculated, and the eigenvector beta is unitized to obtainLet weight vector->At this time, a->Can be used to represent the corresponding valid measurement vector is effectively measured by the set Ω J The overall degree of support for all elements in a system.
According to the basic principle of variance calculation in mathematical statistics, calculating to obtain a measurement noise matrix corresponding to each element in the effective measurement setThe calculation method is as follows:
5. data fusion is carried out on child nodes by using sequential filtering, and fusion results are used for the child node t k The correction of the moment motion parameters comprises the following specific steps:
(1) Establishing an onboard distributed POS transfer alignment model
The transfer alignment model adopts a matching mode of 'position and gesture', and the system state vector is a 15-dimensional state vector and comprises a platform misalignment angle of a sub-node JSpeed error->Position error (δL) J 、δλ J 、δH J ) Gyro constant drift- >And add constant bias ∈ ->N is the number of subsystems. Wherein (1)> The east, north and sky misalignment angles of the J-th sub-node are respectively;respectively the J-th sub-node is at n J The tangential, north and sky speed errors; δL (delta L) J 、δλ J 、δH J Latitude error, longitude error and altitude error of the J-th child node respectively; />Andthe sub-IMUs at the J-th sub-node respectively (b) under its carrier coordinate system J System) gyro constant drift and add constant bias; />sub-IMU at the J-th sub-node at b J Random constant drift of gyroscopes on the x-axis, y-axis and z-axis; />sub-IMU at the J-th sub-node at b J Constant offsets of accelerometers in x-axis, y-axis and z-axis are tied. The system state vector is thus of the form: />
a) Establishment of the equation of state
The system state equation has the following form:
in which the system noiseZero mean Gaussian white noise with variance matrix +.>Determined by the noise level of the gyroscopes and accelerometers in the IMU.
in the formula,FJ The state transition matrix is specifically expressed as follows:
wherein ,
/>
wherein ,ωie Indicating the rotational angular velocity of the earth; and />The sub-IMUs at sub-nodes J are respectively in a navigation coordinate system (n J Tie) east to north to force and natural to force components; / > and />Respectively, the child nodes J are arranged in a navigation coordinate system (n J Tie) east, north and sky speeds; /> and />The main curvature radius of the subnode J along the meridian circle and the mortise unitary circle is respectively; t is the filter period. G J The system noise matrix is a subnode J, and the specific expression form is as follows:
b) Establishment of measurement equation
L effective measurement vectors z of sub-node J J (τ) (τ=1, 2, …, L) is expressed as:
/>
wherein ,respectively representing the difference value of the tau latitude, longitude and altitude calculated value of the child node J and the latitude, longitude and altitude of the main node m after the compensation of the lever arm; /> and />Respectively representing the difference value of the tau-th course angle, pitch angle and roll angle calculated value of the subnode J and the course angle, pitch angle and roll angle corresponding to the main node m;
the measurement equation has the following expression:
z J (tau) corresponding measurement matrixThe method comprises the following steps:
wherein
T a Is the gesture matrix of the main node, i.eAnd let->Representation matrix T a Elements of row s, column T, i.e. T a Can be expressed as:
because any one of the sub-nodes in the airborne distributed POS can obtain redundant motion parameter information of the position where the rest of the sub-nodes are converted, a redundant measurement vector can be obtained according to the method of the step 2, and therefore, the measurement vector can be expressed as follows for the sub-node J:
Measuring vectorThe corresponding measurement matrix can be expressed as +.>
(2) Data fusion using sequential filtering
Firstly, time updating is carried out, and aiming at the J-th child node, according to the last time t k-1 Filtered estimates of (2)And the transfer alignment model can obtain the current t k One-step predictive estimate of time instant +.>The calculation formula is as follows:
wherein ,ΓJ Is to transfer the state of the continuous system into a matrix F J Discrete articleAnd (5) obtaining a transfer matrix of the discrete system after the conversion.
Similarly, the mean square error matrix P is estimated according to the last moment J (k-1|k-1) and the transfer alignment model can obtain the current t k The calculation formula of the one-step prediction mean square error matrix at the moment is as follows:
wherein ,ΞJ Is to drive the continuous system noise into the array G J And discretizing to obtain a noise driving array of the discrete system.
Then, the effective measurement vector obtained in the step (3) is utilized to carry out measurement update, and the J-th sub-node is at the current t k The L measurement vectors at the moment are respectively recorded asThe data fusion formula for the J-th child node is as follows:
when in useWhen the measurement is updated, there are:
wherein , and />Respectively->The corresponding measuring matrix and measuring noise matrix; />For use +.>A filtering gain matrix is calculated when measurement updating is carried out; / >Indicate use +.>Calculating a state estimation value when measuring and updating; />Is->A corresponding estimated mean square error matrix; i is AND->A unit array of the same dimensions.
When in useWhen the measurement is updated in sequence, there are:
wherein , and />Respectively->The corresponding measuring matrix and measuring noise matrix; />For use +.>A filtering gain matrix is calculated when measurement updating is carried out; />Indicate use +.>Calculating a state estimation value when measuring and updating; />Is->A corresponding estimated mean square error matrix.
Executing the data fusion formula, and when the flow runs to tau=L, obtainingAs the current t k The final data fusion result at the J-th child node at time instant, i.e. +.> wherein ,/>Comprises t k Latitude error delta L after data fusion of J-th child node at moment J (k) Longitude error delta lambda J (k) And a height error δH J (k) Also comprises the following step n after data fusion J East speed error of the system->North speed error->And the error of the tangential velocity->And east misalignment angle after data fusion +.>North misalignment angle->And the angle of misalignment of the sky->
(3) Motion parameter correction
The fusion result is used for t k And correcting the result of the J-th child node strapdown calculation at the moment. The correction steps are as follows:
a) Position correction
/>
wherein ,LJ (k)、λ J (k)、H J (k) Respectively represent t k The J-th child node strapdown at the moment calculates the latitude, longitude and altitude; l (L) J (k)| new 、λ J (k)| new 、H J (k)| new Respectively represent t k Latitude, longitude, and altitude after the moment J child node correction.
b) Speed correction
wherein ,respectively represent t k The east speed, the north speed and the sky speed obtained by the moment J child node strapdown calculation; />Respectively represent t k The corrected east, north and sky speeds of the moment J child node.
c) Posture correction
Calculating t k Navigation coordinate system n of J-th sub-node at moment J And calculating a navigation coordinate system n J ' conversion matrix
Modified transformation matrixThe method comprises the following steps:
wherein ,at t k The attitude matrix is obtained after strapdown calculation is carried out on the J child nodes at the moment;
using the modified conversion matrixCalculating t k Heading angle psi of moment J child node J (k)| new Pitch angle theta J (k)| new And roll angle gamma J (k)| new 。
6. Repeating steps 1 to 5 for other sub-nodes which do not perform data fusion and motion parameter correction until all sub-nodes complete t k The method comprises the following specific steps of data fusion at moment and correction of motion parameters:
(1) Let t be k Step 1 to step 1 are completed at the moment5 with child node number k N ，k N An initial value of 1;
(2) If k N N is the number of the child nodes, the fact that the child nodes still do not finish data fusion and correction of motion parameters is indicated, and k is calculated N ＝k N +1, pair number k N Repeating the first 5 steps;
(3) If k N If N, stop step 6, indicating that all child nodes have completed t k Data fusion at the moment and correction of motion parameters.
7、t k ＝t k+1 Steps 1 to 6 are performed until the data fusion and the correction of the motion parameters of all the child nodes at all times are completed.
What is not described in detail in the present specification belongs to the prior art known to those skilled in the art. Variations in the detailed description and the application scope will occur to those skilled in the art upon consideration of the teachings of the present invention. In view of the foregoing, this description should not be construed as limiting the invention.

Claims (8)

1. The data fusion method of the airborne distributed position and posture measurement system comprises the following specific steps:
1.1 calculating t k Redundant position and posture information at the time sub-node;
1.2 calculation of t k Measuring vectors of the time child nodes and generating a bootstrap measuring set by utilizing 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 node;
1.4, calculating the weight and the measurement noise matrix of each effective measurement vector in the effective measurement set;
1.5 data fusion of child node using sequential filtering, fusion result is used for the child node t k Correcting moment motion parameters;
1.6 repeating steps 1.1 to 1.5 for other child nodes not subjected to data fusion and motion parameter correction until all child nodes finish t k Data fusion at moment and correction of motion parameters;
1.7 t k ＝t k+1 steps 1.1 to 1.6 are performed until the data fusion and the correction of the motion parameters of all the child nodes at all times are completed.
2. The method for data fusion of an on-board distributed position and orientation measurement system according to claim 1, wherein: calculating t in step 1.1 k The redundant position and posture information at the time sub-node comprises the following specific steps:
2.1 sub-nodes for strapdown resolution
The definition of the relevant reference coordinate system includes: e is the earth 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 of the J sub node and the calculated navigation coordinate system are respectively represented by n J and n′J J=1, 2, N being the number of child nodes; the origin of the carrier coordinate system is the center of gravity of the carrier, the x axis is rightward 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, the coordinate system is fixed on the carrier and is called the right front upper carrier coordinate system, and b is used m 、b J The carrier coordinate systems representing the main POS and the J-th child node, j=1, 2, N;
each child node obtains t after strapdown calculation k The respective positions [ L ] of the moments J λ J H J ] T Posture [ phi ] J θ J γ J ] T Sum speed of wherein ,LJ 、λ J and HJ Respectively representing the latitude, longitude and altitude of the J-th child node; psi phi type J 、θ J and γJ Respectively representing the course angle, the pitch angle and the 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
Acquiring t through fiber bragg grating sensors mounted on wings k Strain information at any child node at any time is obtained, and displacement of the child node relative to the initial position of the child node is obtained, and relative position information delta P between any two child nodes J, J J←J′ ，ΔP J←J′ Representing the difference in latitude, longitude, altitude between any two child nodes J, J ', J, J' =1, 2, N; j is not equal to J';
at t k At moment, the position obtained by strapdown calculation of the J-th child node is [ L ] J λ J H J ] T J=1, 2, N, the positions of the other N-1 child nodes after strapdown solution are denoted as L J′ λ J′ H J′ ] T J '=1, 2, N, J' noteqj; the redundant location information at the J-th 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′)
Thus, adding the position information obtained by self strapdown calculation, the J-th child node has N latitude, longitude and altitude calculated values, namely [ L ] J λ J H J ] T and [LJ←J′ λ J←J′ H J←J′ ] T J, J '=1, 2, …, N, J' noteqj; position information [ L ] obtained by direct strapdown calculation of current child node J J λ J H J ] T Re-recorded asThe remaining N-1 position information [ L J←J′ λ J←J′ H J←J′ ] T Re-recorded as
b) Sub-node redundant attitude information acquisition
T measured by fiber bragg grating sensor k Angular displacement information of each child node at moment, and constructing transformation matrix between carrier coordinate systems of any two child nodes J, J'Carrier coordinate System b at the J-th child node J Calculating a navigation coordinate system n with the point J The transformation matrix between the lines is denoted +.> The transformation matrix is also called the pose matrix, and the pose matrix at the other N-1 sub-nodes is denoted +.>The redundant pose matrix at the J-th child node may be expressed as:
wherein ,
at t k At the moment, N gesture matrixes exist at the J-th sub-node; n pieces of gesture information, namely N gesture angles, can be calculated by utilizing the N gesture matrixes; the specific calculation method of the attitude angle is as follows:
will beThe method is characterized by comprising the following steps:
wherein ,TJxy Is a matrixX=1, 2,3, y=1, 2,3; the pose angle of the J-th child node includes the heading angle ψ J Pitch angle theta J And roll angle gamma J Main value of (i.e.)> and />The method comprises the following steps of:
heading angle psi J Pitch angle theta J And roll angle gamma J The value ranges of (2) are respectively defined as [0,2 pi ]]、 [-π，+π]The method comprises the steps of carrying out a first treatment on the surface of the Then, ψ J 、θ J and γJ The true value of (2) can be determined by:
according to the attitude angle calculating method, N-1 redundant attitude matrixes at the J-th sub-node are calculatedCorresponding N-1 redundant pose information, J' =1, 2, …, N; j is not equal to J'; the posture matrix to be defined by child node J itself +.>The calculated attitude angle is remembered as +.>From the remaining N-1 gesture matrices +.>The calculated attitude angle is recorded as
3. The method for data fusion of an on-board distributed position and orientation measurement system according to claim 2, wherein: calculating t in step 1.2 k The measuring vector of the time child node and the bootstrap measuring strategy are utilized to generate the bootstrap measuring set, which comprises the following specific steps:
3.1 calculating the measurement vector of the child node
t k The latitude, longitude and altitude of the moment main node m after the compensation of the lever arm are respectively recorded as L m 、λ m 、H m The course angle, pitch angle and roll angle of the main node m are respectively recorded as psi m 、θ m 、γ m The method comprises the steps of carrying out a first treatment on the surface of the For any child junction J, j=1, 2, …, N, t is obtained k N measurement vectors at a time instantThe expression is as follows:
wherein ,respectively represent t k The difference value between the ith latitude, longitude and altitude of the moment sub-node J and the latitude, longitude and altitude of the main node m after the compensation of the lever arm; />Respectively represent t k The difference value of the i-th course angle, pitch angle and roll angle of the moment sub-node J and the course angle, pitch angle and roll angle corresponding to the main node m;
3.2 generating bootstrap measurement sets using measurement bootstrap strategy
In measuring vectorOn the basis of (1) generating bootstrap measures by adding disturbances +.>The method comprises the following specific steps:
wherein ,representing the original measurement vector->Based on the generated c-th bootstrap measurement, c=1,2, …, l, l represents bootstrap quantity +.>I=5; h J Is the measurement matrix of the system, x J Is a state vector of the system, v J Is the measurement noise and v J White gaussian noise satisfying zero mean with covariance +.> and vJ Has the same statistical properties, i.e.)>Gaussian white noise satisfying zero mean, covariance +>
By the method, N multiplied by l bootstrap measurement data are obtainedDefining a bootstrap measurement set: /> wherein />Representing the original measurement vector, i.e.)>
4. The method for data fusion of an on-board distributed position and orientation measurement system according to claim 3, wherein: in the step 1.3, elements in the bootstrap measurement set are sampled by using a Metropolis-Hastings sampling rule to obtain an effective measurement set of the child node, and the specific steps are as follows:
4.1 calculating confidence and probability of acceptance
From N bootstrap measurement sets according to Metropolis-Hastings sampling principleTwo sets are selected at will->Then randomly extracting an element from the two sets>Andand calculates the corresponding confidence level +.> and />The calculation formula of the credibility is as follows:
where z is the measurement prediction mean value,representing measurement noise v J Covariance matrix->Is a determinant of (2);
in obtaining the credibility and />On the basis of (a), the acceptance probability is calculated according to the following formula:
i.e. probability of acceptanceThe value of (2) is +.>And a minimum value between 1;
4.2 construction of an active measurement set
Firstly, a random number χ satisfying the uniform distribution U (0, 1) is generated, and an effective measurement set Ω is obtained J The determination method of (2) is as follows:
in the formula, if randomly extracted elements and />The corresponding acceptance probability is greater than or equal to the generated random number χ, then +.>As an effective measurement set Ω J Elements of (a) and (b); conversely, will->As an effective measurement set Ω J Elements of (a) and (b); the above-mentioned passThe process is the sampling process for determining the effective measurement vector;
repeating the sampling for L times, wherein L=2N, and defining L effective measurement vectors in the effective measurement set as z respectively J (1),z J (2),…z J (L) the effective measurement vectors form an effective measurement set Ω J 。
5. The method for data fusion of an on-board distributed position and orientation measurement system according to claim 4, wherein: the specific steps of calculating the weight and the measurement noise matrix of each effective measurement vector in the effective measurement set in the step 1.4 are as follows:
5.1 calculating the coherence distance and the coherence matrix
Defining an active measurement set Ω J Any two measurement vectors z J(ξ) and zJ Confidence distance between (ζ) isThe calculation method is as follows:
on the basis of which the consistency distance is calculatedAnd a consistency matrix ψ J The calculation method is as follows:
wherein ,representing omega J A maximum value of the confidence distance between any two measurement vectors;
5.2 calculating the weight of the effective measurement vector and the measurement noise matrix
Obtaining a consistency matrix ψ J Then, the maximum modulus eigenvalue lambda and the corresponding eigenvector beta of the matrix are calculated, and the eigenvector beta is unitized to obtainLet weight vector->At this time, a->Representing the corresponding valid measurement vector is valid to the set Ω J The overall degree of support for all elements in (a);
according to the basic principle of variance calculation in mathematical statistics, calculating to obtain a measurement noise matrix corresponding to each element in the effective measurement setThe calculation method is as follows:
6. The method for data fusion of an on-board distributed position and orientation measurement system according to claim 4, wherein: in the step 1.5, the data fusion is carried out on the child node by using sequential filtering, and the fusion result is used for the child node t k The correction of the moment motion parameters comprises the following specific steps:
6.1 building an airborne distributed POS transfer alignment model
The transfer alignment model adopts a matching mode of 'position and gesture', and the system state vector is a 15-dimensional state vector and comprises a platform misalignment angle of a sub-node JSpeed error->Position error δL J 、δλ J 、δH J Gyro constant drift->And add constant bias ∈ ->N is the number of subsystems; wherein (1)>The east, north and sky misalignment angles of the J-th sub-node are respectively;respectively the J-th sub-node is at n J The tangential, north and sky speed errors; δL (delta L) J 、δλ J 、δH J Latitude error, longitude error and altitude error of the J-th child node respectively; />Andsub-IMUs at the J-th sub-node respectively in their carrier coordinate system b J Tying constant drift of the lower gyroscope and adding constant bias; />sub-IMU at the J-th sub-node at b J Random constant drift of gyroscopes on the x-axis, y-axis and z-axis; />Respectively the J-th sub-nodesub-IMU at b J Constant offsets of accelerometers in x-axis, y-axis and z-axis are tied; the system state vector is thus of the form:
a) Establishment of the equation of state
The system state equation has the following form:
in which the system noiseZero mean Gaussian white noise, variance matrix thereofDetermined by the noise levels of the gyroscopes and accelerometers in the IMU;
in the formula,FJ The state transition matrix is specifically expressed as follows:
wherein ,
wherein ,ωie Indicating the rotational angular velocity of the earth; and />Sub IMU at sub-node J in navigation coordinate system n J East to specific force, north to specific force and natural to specific force components in the system; /> and />Respectively sub-nodes J are in a navigation coordinate system n J The tangential speed, the north speed and the sky speed of the tie down; /> and />The main curvature radius of the subnode J along the meridian circle and the mortise unitary circle is respectively; t is a filtering period; g J The system noise matrix is a subnode J, and the specific expression form is as follows:
b) Establishment of measurement equation
L effective measurement vectors z of sub-node J J (τ), τ=1, 2, …, L has the following specific expression:
wherein ,respectively representing the difference value of the tau latitude, longitude and altitude calculated value of the child node J and the latitude, longitude and altitude of the main node m after the compensation of the lever arm; / > and />Respectively representing the difference value of the tau-th course angle, pitch angle and roll angle calculated value of the subnode J and the course angle, pitch angle and roll angle corresponding to the main node m;
the measurement equation has the following expression:
z J (tau) corresponds toIs of the measurement matrix of (a)The method comprises the following steps:
wherein ,
T a is the gesture matrix of the main node, i.eAnd let->Representation matrix T a Elements of row s, column T, i.e. T a Expressed as:
6.2 data fusion Using sequential filtering
Firstly, time updating is carried out, and aiming at the J-th child node, according to the last time t k-1 Filtered estimates of (2)And transmitting the alignment model to obtain the current t k One-step predictive estimate of time instant +.>The calculation formula is as follows:
wherein ,ΓJ Is to transfer the state of the continuous system into a matrix F J A transfer matrix of a discrete system is obtained after discretization;
similarly, the mean square error matrix P is estimated according to the last moment J (k-1|k-1) and passing the alignment model to get the current t k The calculation formula of the one-step prediction mean square error matrix at the moment is as follows:
wherein ,ΞJ Is to drive the noise of continuous system into array G J The noise driving matrix of the discrete system is obtained after discretization;
then, the effective measurement vector obtained in the step 1.3 is utilized to carry out measurement update, and the J-th sub-node is at the current t k The L measurement vectors at the moment are respectively recorded asThe data fusion formula for the J-th child node is as follows:
when in useWhen the measurement is updated, there are:
wherein , and />Respectively->The corresponding measuring matrix and measuring noise matrix; />For use ofA filtering gain matrix is calculated when measurement updating is carried out; />Indicate use +.>Calculating a state estimation value when measuring and updating; p (P) 1 J (k|k) is->A corresponding estimated mean square error matrix; i is AND->A unit array with the same dimension;
when in useWhen the measurement is updated in sequence, there are:
wherein , and />Respectively->The corresponding measuring matrix and measuring noise matrix; />For use ofA filtering gain matrix is calculated when measurement updating is carried out; />Indicate use +.>Calculating a state estimation value when measuring and updating; />Is->A corresponding estimated mean square error matrix;
executing the data fusion formula, and when the flow runs to tau=L, obtainingAs the current t k The final data fusion result at the J-th child node at time instant, i.e. +.> wherein ,/>Comprises t k Latitude error delta L after data fusion of J-th child node at moment J (k) Longitude error delta lambda J (k) And a height error δH J (k) Also comprises the following step n after data fusion J East speed error of the system->North speed error->And the error of the tangential velocity->And east misalignment angle after data fusion +.>North misalignment angle->And the angle of misalignment of the sky->
6.3 motion parameter correction
The fusion result is used for t k And correcting the result of the time J child node strapdown calculation, wherein the correcting step is as follows:
a) Position correction
wherein ,LJ (k)、λ J (k)、H J (k) Respectively represent t k The J-th child node strapdown at the moment calculates the latitude, longitude and altitude; l (L) J (k)| new 、λ J (k)| new 、H J (k)| new Respectively represent t k Latitude, longitude and altitude corrected by the J th child node at the moment;
b) Speed correction
wherein ,respectively represent t k The east speed, the north speed and the sky speed obtained by the moment J child node strapdown calculation; />Respectively represent t k The east speed, the north speed and the sky speed corrected by the J-th child node at the moment;
c) Posture correction
Calculating t k Navigation coordinate system n of J-th sub-node at moment J And calculating a navigation coordinate system n J ' conversion matrix
Modified transformation matrixThe method comprises the following steps:
wherein ,at t k The J-th child node at the moment performs strapdown calculation to obtain a gesture matrix;
using the modified conversion matrixCalculating t k Heading angle psi of moment J child node J (k)| new Pitch angle theta J (k)| new And roll angle gamma J (k)| new 。
7. The method for data fusion of an on-board distributed position and orientation measurement system according to claim 1, wherein: repeating steps 1.1 to 1.5 for other child nodes which do not perform data fusion and motion parameter correction in step 1.6 until all child nodes complete t k The method comprises the following specific steps of:
7.1 assume t k The child node number k at which steps 1.1 to 1.5 have been completed N ，k N An initial value of 1;
7.2 if k N If N is less than N, N is the number of the child nodes, the fact that the child nodes still have incomplete data is indicatedFusion and correction of motion parameters, k N ＝k N +1, pair number k N Repeating steps 1.1 to 1.5;
7.3 if k N If N, stop step 1.6, indicating that all child nodes have completed t k Data fusion at the moment and correction of motion parameters.
8. The method for data fusion of an on-board distributed position and orientation measurement system according to claim 1, wherein: in the step 1.7, t is set k ＝t k+1 And executing 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 moments are completed.
CN202110309944.7A 2021-03-23 2021-03-23 Airborne distributed POS data fusion method Active CN113188566B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110309944.7A CN113188566B (en) 2021-03-23 2021-03-23 Airborne distributed POS data fusion method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110309944.7A CN113188566B (en) 2021-03-23 2021-03-23 Airborne distributed POS data fusion method

Publications (2)

Publication Number Publication Date
CN113188566A CN113188566A (en) 2021-07-30
CN113188566B true CN113188566B (en) 2023-09-29

Family

ID=76973665

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110309944.7A Active CN113188566B (en) 2021-03-23 2021-03-23 Airborne distributed POS data fusion method

Country Status (1)

Country Link
CN (1) CN113188566B (en)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104655152A (en) * 2015-02-11 2015-05-27 北京航空航天大学 Onboard distributed type POS real-time transmission alignment method based on federal filtering
CN107728182A (en) * 2017-09-18 2018-02-23 北京航空航天大学 Flexible more base line measurement method and apparatus based on camera auxiliary
CN108387227A (en) * 2018-02-22 2018-08-10 北京航空航天大学 The multinode information fusion method and system of airborne distribution POS
CN108458709A (en) * 2018-02-22 2018-08-28 北京航空航天大学 The airborne distributed POS data fusion method and device of view-based access control model subsidiary
WO2020233290A1 (en) * 2019-05-17 2020-11-26 东南大学 Dual-filter-based transfer alignment method under dynamic deformation
CN112525191A (en) * 2021-02-08 2021-03-19 北京航空航天大学 Airborne distributed POS transfer alignment method based on relative strapdown calculation

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104655152A (en) * 2015-02-11 2015-05-27 北京航空航天大学 Onboard distributed type POS real-time transmission alignment method based on federal filtering
CN107728182A (en) * 2017-09-18 2018-02-23 北京航空航天大学 Flexible more base line measurement method and apparatus based on camera auxiliary
CN108387227A (en) * 2018-02-22 2018-08-10 北京航空航天大学 The multinode information fusion method and system of airborne distribution POS
CN108458709A (en) * 2018-02-22 2018-08-28 北京航空航天大学 The airborne distributed POS data fusion method and device of view-based access control model subsidiary
WO2020233290A1 (en) * 2019-05-17 2020-11-26 东南大学 Dual-filter-based transfer alignment method under dynamic deformation
CN112525191A (en) * 2021-02-08 2021-03-19 北京航空航天大学 Airborne distributed POS transfer alignment method based on relative strapdown calculation

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
Multi-Node Transfer Alignment Based on Mechanics Modeling for Airborne DPOS;Gong Xiaolin et al.;《IEEE SENSORS JOURNAL》;第18卷(第2期);669-679 *
Multisensor Time Synchronization Error Modeling and Compensation Method for Distributed POS;Li Jianli et al.;《IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT》;第65卷(第11期);2637-2645 *
基于Metropolis-Hastings采样的多传感器集合卡尔曼滤波算法;胡振涛等;《电子学报》;第45卷(第4期);868-873 *
机载分布式POS传递对准建模与仿真;房建成等;《中国惯性技术学报》;第20卷(第4期);379-385 *
联邦模糊自适应卡尔曼滤波在速度+姿态传递对准中的应用;宋丽君;秦永元;;测控技术(第05期);135-138 *

Also Published As

Publication number Publication date
CN113188566A (en) 2021-07-30

Similar Documents

Publication Publication Date Title
CN110501024B (en) Measurement error compensation method for vehicle-mounted INS/laser radar integrated navigation system
CN108731670B (en) Inertial/visual odometer integrated navigation positioning method based on measurement model optimization
CN107728182B (en) Flexible multi-baseline measurement method and device based on camera assistance
CN112097763B (en) Underwater vehicle combined navigation method based on MEMS IMU/magnetometer/DVL combination
CN108387227B (en) Multi-node information fusion method and system of airborne distributed POS
CN101949703B (en) Strapdown inertial/satellite combined navigation filtering method
CN110702143B (en) Rapid initial alignment method for SINS strapdown inertial navigation system moving base based on lie group description
CN112629538A (en) Ship horizontal attitude measurement method based on fusion complementary filtering and Kalman filtering
CN111337020A (en) Factor graph fusion positioning method introducing robust estimation
CN108458709B (en) Airborne distributed POS data fusion method and device based on vision-aided measurement
CN108375383B (en) Multi-camera-assisted airborne distributed POS flexible baseline measurement method and device
CN106772524A (en) A kind of agricultural robot integrated navigation information fusion method based on order filtering
CN107014376A (en) A kind of posture inclination angle method of estimation suitable for the accurate operation of agricultural machinery
CN103884340B (en) A kind of information fusion air navigation aid of survey of deep space fixed point soft landing process
CN107727097B (en) Information fusion method and device based on airborne distributed position and attitude measurement system
CN110849360B (en) Distributed relative navigation method for multi-machine collaborative formation flight
CN113252038B (en) Course planning terrain auxiliary navigation method based on particle swarm optimization
CN110275193B (en) Cluster satellite collaborative navigation method based on factor graph
CN110285834B (en) Quick and autonomous readjustment method for double-inertial navigation system based on one-point position information
CN113291493B (en) Method and system for determining fusion attitude of multiple sensors of satellite
CN108303120B (en) Real-time transfer alignment method and device for airborne distributed POS
CN113063429A (en) Self-adaptive vehicle-mounted combined navigation positioning method
CN112325886A (en) Spacecraft autonomous attitude determination system based on combination of gravity gradiometer and gyroscope
CN111220151B (en) Inertia and milemeter combined navigation method considering temperature model under load system
CN110388942B (en) Vehicle-mounted posture fine alignment system based on angle and speed increment

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant