# US 10,768,312 B1 · Integrity analysis method based on kinematic-to-kinematic relative positioning scenario

> Agents: this data is available as JSON from the patentgrep API and MCP server — see https://patentgrep.com/llms.txt

- Number: US 10,768,312 B1
- Document ID: US10768312B1
- Status: Granted
- Kind: B1
- Application number: 16/713,322
- Filed: 2019-12-13
- Granted: 2020-09-08
- Earliest claimed: 2019-11-07
- Assignees: Beihang University
- Applicants: Beihang University
- Inventors: Zhipeng Wang; Yanbo Zhu; Xiaopeng Hou; Kun Fang
- Web page: https://patentgrep.com/patent/US10768312B1

## CPC

- G01S 19/20: Integrity monitoring, fault detection or fault isolation of space segment
- G01S 19/44: Carrier phase ambiguity resolution; Floating ambiguity; LAMBDA [Least-squares AMBiguity Decorrelation Adjustment] method
- G01S 19/51: Relative positioning
- G01S: RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES

## Abstract

The present invention provides an integrity analysis method based on a kinematic-to-kinematic relative positioning scenario, including the following steps: a) establishing a kinematic-to-kinematic relative positioning model, and inputting navigation data; b) calculating a float solution of an integer ambiguity; c) detecting and correcting cycle slips based on a total electron content rate; d) calculating a probability of correct fix and a probability of incorrect fix for the integer ambiguity; e) determining a fault to be detected and a satellite fault probability; 0 calculating a standard deviation δ_(v|CF) and a position domain deviation b_m; and g) calculating an integrity risk value of a carrier phase. The present invention provides an integer ambiguity calculation algorithm for a kinematic-to-kinematic positioning system in the case of a long baseline, to calculate carrier phase integrity.

## Claims

### Claim 1

An integrity analysis method based on a kinematic-to-kinematic relative positioning scenario, wherein, the analysis method comprises the following steps:
a) establishing a kinematic-to-kinematic relative positioning model for a navigation system, and inputting navigation data;
b) calculating a float solution of an integer ambiguity;
c) detecting and correcting cycle slips based on a total electron content rate, comprising:
step 1: performing linear combination of wide-lane phase minus narrow-lane pseudorange and detecting an ionospheric total electron content rate for the foregoing solution data respectively:
[
Δ
⁢
N
1
P
⁡
(
k
)
-
Δ
⁢
N
2
P
⁡
(
k
)
]
=
N
WL
P
⁡
(
k
-
1
)
-
N
WL
P
(
⁢
k
)
=
[
⁢
Φ
1
P
⁡
(
k
-
1
)
-
Φ
2
P
⁡
(
k
-
1
)
-
f
1
⁢
P
1
P
⁡
(
k
-
1
)
+
f
2
⁢
P
2
P
⁡
(
k
-
1
)
λ
WL
⁡
(
f
1
+
f
2
)
]
-
[
⁢
Φ
1
P
⁡
(
k
)
-
Φ
2
P
⁡
(
k
)
-
⁢
f
1
⁢
P
1
P
⁡
(
k
)
+
f
2
⁢
P
2
P
⁡
(
k
)
λ
WL
⁡
(
f
1
+
f
2
)
]
;
⁢
[
λ
1
⁢
Δ
⁢
N
1
P
⁡
(
k
)
-
λ
2
⁢
Δ
⁢
N
2
P
⁡
(
k
)
]
=
4
⁢
0
.
3
×
1
⁢
0
16
⁢
(
γ
-
1
)
⁢
Δ
⁢
⁢
t
·
TECR
Φ
⁡
(
k
)
f
1
2
-
λ
1
⁡
[
Φ
1
P
⁡
(
k
)
-
Φ
1
P
⁡
(
k
-
1
)
]
+
λ
2
⁡
[
Φ
2
P
⁡
(
k
)
-
Φ
2
P
⁡
(
k
-
1
)
]
;
where f1, f2 represent two frequencies of a received double-frequency satellite signal; λ1, λ2 represent wavelengths at corresponding frequencies;
γ
=
f
1
2
f
2
2
;
λ
WL
=
c
(
f
1
-
f
2
)
;
TECRΦ(k) represents an ionospheric total electron content rate at an epoch K; Φ1P(k) and Φ1P(k−1) represent carrier phase measurement values between two adjacent epochs at the frequency f1; P1P(k) and P1P(k−1) represent pseudorange measurement values between two adjacent epochs at the frequency f1; ΔN1P(k) and ΔN2P(k) represent cycle slips that are detected at the frequencies f1 and f2 using a cycle as a unit respectively; c represents a light velocity; NWLP(k) and NWLP(k−1) represent integer ambiguities between two adjacent epochs;
step 2: calculating influence of a pseudo random error on cycle slip detection;
SPIEmax=max{|δP1P(k)|,|δ2P(k)|},δP1P(k),δP2P(k)∈(−∞,+∞);
where δP1P(k) and δP2P(k) represent pseudorange measurements at f1 and f2 at the epoch K respectively, and SPIEmax represents a maximum value of a pseudo random error;
step 3: calculating cycle slip search areas;
SSΔN1P(k):ΔN1P(k)−δΔN1P(k)≤ΔN1P(k)cand≤ΔN1P(k)+δΔN1P(k);
SSΔN2P(k):ΔN2P(k)−δΔN2P(k)≤ΔN2P(k)cand≤ΔN2P(k)+δΔN2P(k);
where ΔN1P(k)cand represents a correct cycle slip search area at L1, and ΔN2P(k)cand represents a correct cycle slip search area at L2;
step 4: calculating a wide-lane cycle slip variation by using a cycle slip search rule 1; and calculating a total electron content rate threshold by using a cycle slip search rule 2;
TECRΦ(k)=H{λ1[Φ1P(k)−Φ1P(k−1)]+λ2[Φ2P(k−1)−Φ2P(k)]+[λ1ΔN1P(k)cand−λ2ΔN2P(k)cand]};
wherein
H
=
f
1
2
4
⁢
0
.
3
×
1
⁢
0
16
⁢
(
γ
-
1
)
⁢
Δ
⁢
t
;
d) calculating a probability of correct fix PCF and a probability of incorrect fix PIF for the integer ambiguity;
calculating a probability of a candidate ambiguity value Z:
P
⁡
(
a
⋓
=
z
)
=
∏
i
=
1
m
⁢
[
Φ
⁡
(
1
-
2
⁢
I
i
T
⁡
(
a
-
z
)
2
⁢
σ
i
/
I
)
+
Φ
⁡
(
1
+
2
⁢
I
i
T
⁡
(
a
-
z
)
2
⁢
σ
i
/
I
)
-
1
]
;
where ă represents an integer ambiguity estimated using a bootstrap method; a represents a correctly fixed integer ambiguity; σi/I represents a standard deviation of ambiguities, and the standard deviation depends on a set I={1, 2, . . . , i−1}; Ii represents the ith column vector of a unit-lower-triangular matrix L−T generated from decomposition of LDLT of a fuzzy covariance matrix; m represents the number of fixed integer ambiguities, and
Φ
⁡
(
x
)
=
∫
-
∞
x
⁢
(
1
/
2
⁢
π
)
⁢
exp
⁡
(
-
1
2
⁢
v
2
)
⁢
dv
;
a−z represents a difference obtained by subtracting an integer ambiguity calculated using a correlation method from an integer ambiguity calculated using a bootstrap method;
when the correctly fixed integer ambiguity a is equal to the candidate ambiguity value Z, the probability of correct fix PCF is calculated by:
P
⁡
(
a
⋓
=
a
)
=
∏
i
=
1
m
⁢
[
2
⁢
Φ
⁡
(
1
2
⁢
σ
i
/
I
)
-
1
]
;
the probability of incorrect fix PIF is calculated by using an exclusion relation between a correct fix event and an incorrect fix event;
e) determining a fault to be detected and a satellite fault probability;
f) calculating a standard deviation δv|CF and a position domain deviation bm;
wherein the following models are established: y=Gx+ε and
{
x
^
=
Sy
S
=
(
G
T
⁢
WG
)
-
1
⁢
GW
;
and an ambiguity float solution-relative position vector matrix Ŝ and a related estimated error covariance matrix is obtained according to the foregoing equations:
S
^
=
[
x
^
N
^
]
;
P
=
[
P
x
^
P
x
^
⁢
N
^
P
N
^
⁢
x
^
P
N
^
]
;
where G represents a coefficient matrix; ϵ represents a total satellite error; {circumflex over (x)} represents an estimated value in a solution state; W represents a relation matrix between a pseudorange and an error; {circumflex over (N)} represents an estimated value of an integer ambiguity; S represents a state vector; P represents a covariance matrix of position and float ambiguity state estimated error;
step 1: calculating the kth integer ambiguity element Nz(k);
Nz(k)=round([01×3Zk]S(k-1));
where Zk represents the kth row of a transition matrix Z;
Z1:k=[01×(k-1) 1 01×(n-k)]; S(k-1) represents an updated state vector S based on fixed ambiguities 1→K;
step 2: calculating the standard deviation δv|CF;
continuously updating position vector estimation and remaining ambiguity real values S;
K=P(k-1)[01×3Zk]T([01×3Zk]TP(k-1)[01×3Zk]T)−1;
S(k)=S(k-1)+K(NZ(k)−[01×3Zk]TS(k-1));
P(k)=(I−K[01×3Zk])P(k-1);
where K represents a gain; P represents a covariance matrix of position and float ambiguity state estimated error; P(k) represents a covariance matrix of an updated state vector after the kth ambiguity is fixed; and it is assumed that δv|CF=√{square root over (P(3,3)(m))} after m iterations;
step 3: calculating the position domain deviation bm;
K=P[0m×3Z1:m]T([0m×3Z1:m]P[0m×3Z1:m]T)−1;
S(m)=Ŝ+K(NZ(1:m)−[0m×3Z1:m]TŜ);
where Z1:m represents a transition matrix Z of a row vector 1 to m; Z1:m=[Im×m 0m×(n-m)]; NZ(1:m) represents a fixed solution for fix of 1 to m integer ambiguities;
the position domain deviation bm caused by a relative position error of the mth incorrect candidate pair (a−z)m is:
bm=SCF(m)−SIF(m)=K(a−z)m;
where am represents a correctly fixed ambiguity vector; zm represents an incorrectly fixed ambiguity vectors; SCF(m) represents a correct state obtained after fix of m ambiguities; SIF(m) represents an incorrect state obtained after fix of m ambiguities;
SCF(m)=Ŝ+K(am−[0m×3Z1:m]Ŝ);
SIF(m)=Ŝ+K(zm−[0m×3Z1:m]Ŝ);
g) calculating an integrity risk value IH0 of a carrier phase for the GNSS;
I
H
⁢
0
=
1
-
(
1
-
P
⁢
{

x
^
v
-
x
v

>
VAL
|
C
⁢
F
}
)
⁢
PCF
-
∑
n
=
1
k
⁢
(
1
-
P
⁢
{

x
^
v
-
x
v

>
VAL
|
IF
n
}
)
⁢
PIF
n
;
where PVALCF=P{|{circumflex over (x)}v−xv|>VAL|CF} represents a probability that a vertical error is greater than an alarm threshold under the condition of correct fix; PCF represents a probability of correct fix; PIFn represents a probability of incorrect fix; PVALIF=P{|{circumflex over (x)}v−xv|>VAL|IFn} represents a probability that a vertical error is greater than an alarm threshold under the condition of incorrect fix; and
h) applying the integrity risk value IH0 in the navigation system to improve position precision of the navigation system on the kinematic-to-kinematic relative positioning scenario by comparing the integrity risk value IH0 to a zero-fault integrity value IH0req, wherein a geometric mechanism of the navigation system is considered to be available when the integrity risk value IH0 is less than the zero-fault integrity value.

### Claim 2 (depends on 1)

The analysis method as described in claim 1, wherein, the calculating a float solution of an integer ambiguity comprises the following steps:
establishing a function model y=f(x,β) corresponding to m satellite-collect points (x1,y1,z1), (x2,y2,z2), . . . , (xm,ym,zm), and calculating a minimum value β of a residual sum of squares, to minimize the residual sum of squares;
wherein S(β)=Σi=1m ri(β)2; and a residual ri=yi−f(xi,β) is calculated;
the minimum value is calculated by taking a derivative of the foregoing equation:
∂
s
∂
β
j
=
2
⁢
∑
i
⁢
r
i
⁢
∂
r
i
∂
β
j
=
0
⁢
⁢
(
j
=
1
,
…
⁢
,
n
)
.

### Claim 3 (depends on 1)

The analysis method as described in claim 1, wherein, a maximum allowable number of satellite faults is calculated according to Nfault,max=max{r∈1, . . . , Nsat|Pmultiple(r+1, Pevent,1, . . . , Pevent,Nsat+Nconst)≤PFAVLT_THRES};
where Pevent,i=PIFi represents a prior probability of an independent error event i; Nfault,max represents a probability of a multi-fault modes consisting of r independent events; Nsat represents the number of satellites; Nconst represents the number of constellations; Pmultiple represents a probability that multiple satellite are faulty; and PFAVLT-THRES represents an alarm threshold.

### Claim 4 (depends on 3)

The analysis method as described in claim 3, wherein, Nfault,max of a combination corresponding to all subsets is determined according to the Nfault,max, and is used to determine the satellite fault probability.

### Claim 5 (depends on 1)

The analysis method as described in claim 1, wherein, a conditional variance σi/I2, of an ambiguity i is defined, the conditional variance depends on previously fixed integer ambiguities of the set I={1, 2, . . . , i−1}, and an item (i, i) in a diagonal matrix D is generated from decomposition of a float solution estimated error covariance matrix LDLT of cycle ambiguities.

### Claim 6 (depends on 1)

The analysis method as described in claim 1, wherein, by calculating the influence of the pseudo random error on cycle slip detection, a maximum value δΔN1P(k)max of cycle slips detected at the frequency f1 using a cycle as a unit, a maximum value δΔN2P(k)max of cycle slips detected at the frequency f2 using a cycle as a unit, and a difference δ[ΔN1P(k)−ΔN2P(k)]max between the maximum values of cycle slips detected at the frequencies f1 and f2 using a cycle as a unit is calculated:
δ
⁡
[
Δ
⁢
⁢
N
1
P
⁡
(
k
)
-
Δ
⁢
N
2
P
⁡
(
k
)
]
=
int
⁡
[
f
1
⁢
δ
⁢
P
1
P
⁡
(
k
)
+
f
2
⁢
δ
⁢
P
2
P
⁡
(
k
)
λ
WL
⁡
(
f
1
+
f
2
)
]
;
δ
⁢
Δ
⁢
⁢
N
1
P
⁡
(
k
)
=
int
⁡
[
-
λ
2
(
λ
1
-
λ
2
)
⁢
f
1
⁢
δ
⁢
⁢
P
1
P
⁡
(
k
)
+
f
2
⁢
δ
⁢
⁢
P
2
P
⁡
(
k
)
λ
WL
⁡
(
f
1
+
f
2
)
]
;
δΔ
⁢
⁢
N
2
P
⁡
(
k
)
=
int
⁡
[
-
λ
1
(
λ
1
-
λ
2
)
⁢
f
1
⁢
δ
⁢
⁢
P
1
P
⁡
(
k
)
+
f
2
⁢
δ
⁢
⁢
P
2
P
⁡
(
k
)
λ
WL
⁡
(
f
1
+
f
2
)
]
.

### Claim 7 (depends on 1)

The analysis method as described in claim 1, wherein, the wide-lane cycle slip variation is calculated by using the maximum value SPIEmax of the pseudo random error.

### Claim 8 (depends on 1)

The analysis method as described in claim 1, wherein, a threshold model of the total electron content rate is further used to reduce remaining cycle slip candidate pairs, and is further used to revise carrier phase measurement values at the frequencies f1 and f2 for candidate pairs ΔN1P(k)cand at the frequency f1 and candidate pairs ΔN2P(k)cand at the frequency f2 by the following equations:
=Φ1P(k)+ΔN1P(k)cand;
=Φ2P(k)+ΔN2P(k)cand.

### Claim 9 (depends on 1)

The analysis method as described in claim 1, wherein, when the probability of correct fix PCF and the probability of incorrect fix PIF are calculated for the integer ambiguity, a conditional variance of an ambiguity i is defined as σi/I2, the conditional variance depends on previously fixed integer ambiguities of the set I={1, 2, . . . , i−1}, and an item (i, i) in a diagonal matrix D is generated from decomposition of a float solution estimated error covariance matrix LDLT of cycle ambiguities.

### Claim 10 (depends on 1)

The analysis method as described in claim 1, wherein, the probability of correct fix PCF is calculated, and for the calculation of the probability of incorrect fix PIF, a value of (a−z) is used, and a candidate value of interest in an incorrect fix event is calculated.

### Claim 11 (depends on 1)

The analysis method as described in claim 1, further comprising, in step h):
when the integrity risk value IH0 is greater than the zero-fault integrity value IH0req, determining that the float solution of a next ambiguity is to be fixed, without traversing all of the candidate ambiguity values.