# A Vessel Positioning Algorithm Based on Satellite Automatic Identification System

Authors: Shexiang Ma; Jie Wang; Xin Meng; Junfeng Wang
Journal: Journal of Electrical and Computer Engineering
 (2017)

---

## Abstract

Vessels can obtain high precision positioning by using the global navigation satellite system (GNSS), but when the ship borne GNSS receiver fails, the existence of an alternative positioning system is important for the navigation safety of vessel. In this paper, a localization method based on the signals transmitted by satellite-based automatic identification system (AIS) is proposed for vessel in GNSS-denied environments. In the proposed method, the positioning model is a modification on the basis of time difference and frequency difference of arrival measurements by introducing an additional measurement, and the measurement is obtained through the interactive multiple model algorithm. The performance of the proposed strategy is evaluated through simulations, and the results validate the feasibility and reliability of vessel localization based on satellite-based AIS.


---

## Body

## 1. Introduction

Automatic identification system (AIS) is a self-reporting system designed to protect maritime security of vessel and improve maritime efficiency [1]. It plays an important role in ship collision avoidance and maritime supervision through a series of static and dynamic vessel information automatically broadcast, and the information includes latitude, longitude, course, and velocity [2]. The geographical location reported in AIS is derived by the shipboard GNSS receiver and typically with the high accuracy [3]. However, there is a problem that followed with the GNSS being widely used in navigation of maritime. GNSS is vulnerable to accidental interference [4]; the ship will not be able to locate once the GNSS signal is deliberately disturbed or the GNSS receiver fails. So it is necessary to develop a spare navigation system for the ship.AIS is a self-organized time division multiple access (TDMA) system, which not only can be self-reporting but also can receive AIS information [5]. Although AIS ignored the role of the satellite in its original design; it has been proven feasible to receive AIS signals by satellite [6, 7]. In the satellite-based AIS, the relative speed of satellite and ship is high, and the two are far apart; therefore, there are challenges for the correct detection of the AIS signal such as the problem of time delay, the high Doppler offset, and low signal-to-noise ratio (SNR) [8]. With the current level of AIS signal detection technology, the correct detection of AIS signal can be guaranteed with the improvement of synchronization algorithm, and the influence of high Doppler offset on carrier recovery is gradually decreasing [9]. The satellite-based AIS is already operational but focuses on the stage of “vessel transmitting, satellite receiving.” In view of the fact that a large number of AIS signals are likely to reach the satellite at the same time in this stage but satellite can still detect ship signals [8], the signals can definitely be received by the ship if the satellite can send information in the AIS operating frequency band according to the AIS protocol, because the possibility of AIS signal conflict is relatively low in case of the ship reception thanks to the characteristics of signal transmission. As the technology of satellite-based AIS advances, the potentiality of AIS for navigation becomes a concern and there is the investigation on ship localization using AIS signals received by satellite [10]. In this paper, we assume that, in advanced satellite-based AIS, vessels can receive AIS signals transmitted from satellite in addition to “vessel transmitting, satellite receiving,” and the information of satellite motion state is broadcast by the downlink AIS signal. On the basis of this vision, a ship localization method using AIS signals transmitted from satellite is proposed.Among the various measurements for positioning tasks, the time of arrival (TOA), the time difference of arrival (TDOA), and the frequency difference of arrival (FDOA) are very representative choices because of the potentials in attaining high localization accuracy [11, 12]. There is a lot of research on the application of TDOA to improve the positioning accuracy of the stationary target and to locate the moving target by using frequency measurements [13–15]. In addition, the positioning methods combining two kinds of measurements such as TDOA/FDOA and TDOA/DOA are also widely discussed [16, 17]. Except for reducing the number of signals required, the combination of time and frequency measurements is attractive for improvement of positioning accuracy [18, 19]. In these methods, however, only the information extracted from the received signals is used for positioning. Taking into account the fact that AIS can obtain the ship velocity and heading by connecting external sensors, in this paper, a modification positioning model based on TDOA/FDOA is proposed by introducing an additional measurement based on the interactive multiple model (IMM) algorithm [20, 21]. The method of TDOA/FDOA and the IMM algorithm are used separately for locating or tracking the target in general; given the characteristic of AIS signal carrying information, they are combined together in the study. Besides, for the purpose of making the estimated result more suitable for the vessel status, a new probability updating method for IMM is designed in this work.The solution of TDOA/FDOA measurement equation is complicated because of the high nonlinearity [22]. Taylor-series technique can linearize the equations but positioning result is easy to be affected by initial value setting [23]. The method of grid searching achieves the accuracy improvement with the sacrifice of computation [24]. In this paper, the localization results are obtained by Gauss-Newton iteration under the least squares criterion, and the solution of grid rough searching is chosen as the starting value. The feasibility of the positioning method based on the advanced satellite-based AIS signals is investigated through the experiment. The location error distributions of the TDOA/FDOA joint location model and the proposed localization model are analyzed in this study.


## 2. TDOA/FDOA Localization Based on Least Squares Estimation

In satellite-based AIS, satellites are located at a low orbit from 600 km to 1000 km above the ground. The downlink AIS signals will include the Doppler frequency shift because of the relative satellite-ship velocities, and the frequency shift is up to a maximum of ±4 kHz. In order to achieve localization with the limited number of AIS signals and improve the positioning accuracy as much as possible, the work of ship positioning is carried out on the basis of TDOA/FDOA in this paper.


### 2.1. Principles of TDOA/FDOA

Assuming thatti is the time cost by the ith AIS signal transmitted from satellite to ship, the TDOA between the adjacent signals received by the ship can be expressed as(1)Δti=ti+1-ti=Lsi+1-Lb-Lsi-Lbc+Δnt=1cxsi+1-x2+ysi+1-y2+zsi+1-z21/2-xsi-x2+ysi-y2+zsi-z21/2+Δnt=rtix,y,z+Δnt.The FDOA between the adjacent AIS signals can be expressed as(2)Δfri=fecVsi+1T·Lsi+1-LbLsi+1-Lb-VsiT·Lsi-LbLsi-Lb+Δnf=fecVsxi+1xsi+1-x+Vsyi+1ysi+1-y+Vsxi+1zsi+1-zxsi+1-x2+ysi+1-y2+zsi+1-z21/2-Vsxixsi-x+Vsyiysi-y+Vsxizsi-zxsi-x2+ysi-y2+zsi-z21/2+Δnf=rfix,y,z+Δnf,where fe is the carrier frequency of the AIS signal and c is the signal propagation velocity. Lb=x,y,zT is the vessel position vector in the ECEF reference and Vsi=Vsxi,Vsyi,VsziT and Lsi=xsi,ysi,zsiT are velocity vector and position vector of the satellite when transmitting the ith AIS signal, respectively. Δnt is the difference of noise between the two time measurements and Δnf is the difference of noise between the two frequency measurements.It is assumed that the number of signals received by the ship in the visual time of the satellite isN+1; the localization equation matrix based on (1) and (2) can be written as (3)ΔTΔF=rtx,y,zrfx,y,z+nwith(4)rtx,y,z=rt1x,y,z,rt2x,y,z,…,rtNx,y,zTrfx,y,z=rf1x,y,z,rf2x,y,z,…,rfNx,y,zT,where ΔT=Δt1,Δt2,…,ΔtNT is the TDOA measurement vector obtained by synchronization technique and ΔF=Δfr1,Δfr2,…,ΔfrNT is the FDOA measurement vector. n is the measurement noise matrix.


### 2.2. Calculation Based on Least Squares Estimation

On the basis of least squares criterion, the estimator associated with (3) needs to minimize the differences between the measurements and predictions; the equation to be minimized can be written as(5)Cx,y,z=Υ-qx,y,zTN-1Υ-qx,y,z,where N is the noise covariance matrix, Υ=ΔTΔFT is the measurement vector of TDOA/FDOA, and q=rtrfT.The estimated position valuex^,y^,z^=argminx,y,zCx,y,z can be obtained by Gauss-Newton iteration(6)xj+1,yj+1,zj+1T=xj,yj,zjT+JTxj,yj,zjN-1Jxj,yj,zj-1×JTxj,yj,zjN-1Υ-qxj,yj,zj,where matrix J is(7)J=JtJfwith (8)Jt=rt1x,y,z∂xrt1x,y,z∂yrt1x,y,z∂z⋮⋮⋮rtNx,y,z∂xrtNx,y,z∂yrtNx,y,z∂z,Jf=rf1x,y,z∂xrf1x,y,z∂yrf1x,y,z∂z⋮⋮⋮rfNx,y,z∂xrfNx,y,z∂yrfNx,y,z∂z.The initial position in (6) needs to be defined in advance; the method for determining the initial position (shown in Figure 1) is as follows: (1) setting up a grid with the units of 1°. The grid is centered on the midpoint gλo,gφo of the satellite ground trajectory (the track generated during AIS signals transmission) and the range of grid geodetic coordinates gλ,gφ is gλo-β/2≤gλ≤gλo+β/2,gφo-ε/2≤gφ≤gφo+ε/2, where β and ε are the maximum visible longitude and latitude of satellite, respectively. (2) Connecting the start and end points of the satellite ground trajectory and dividing the grid into S1 and S2 (two parts) by extending the connecting line. (3) Searching within each part of grid and selecting two points with mλi,mφi=argmaxSiE=1/fr-f^2 (i=1,2). In the cost function E, fr is the measured frequency of the AIS signal and f^=fc1-VsiTLsi-Lg/cLsi-Lg is the estimated frequency of received signal at the grid point, where Lg=xg,yg,zgT is the position vector of grid point in ECEF coordinate. In this paper, we select “nearest point” to eliminate the false image which may occur in grid searching, that is, taking the point with the shortest distance from the origin of the ship as the optimal position λ0,φ0. The transformation of vessel location from the geodetic coordinate to ECEF coordinate is defined as follows:(9)x=RNcos⁡λcos⁡φy=RNsin⁡λcos⁡φz=RN1-e2sin⁡φ,where λ and φ are the longitude and latitude coordinates of vessel and x,y,z are the ECEF coordinates of vessel. RN=α/1-e2sin2φ is radius of curvature in prime vertical, where e2=0.00669437999013 and α=6378.137km are the square of the first eccentricity and equatorial radius of the earth defined by WGS-84, respectively.Figure 1
Method of determining the initialization vector.


## 2.1. Principles of TDOA/FDOA

Assuming thatti is the time cost by the ith AIS signal transmitted from satellite to ship, the TDOA between the adjacent signals received by the ship can be expressed as(1)Δti=ti+1-ti=Lsi+1-Lb-Lsi-Lbc+Δnt=1cxsi+1-x2+ysi+1-y2+zsi+1-z21/2-xsi-x2+ysi-y2+zsi-z21/2+Δnt=rtix,y,z+Δnt.The FDOA between the adjacent AIS signals can be expressed as(2)Δfri=fecVsi+1T·Lsi+1-LbLsi+1-Lb-VsiT·Lsi-LbLsi-Lb+Δnf=fecVsxi+1xsi+1-x+Vsyi+1ysi+1-y+Vsxi+1zsi+1-zxsi+1-x2+ysi+1-y2+zsi+1-z21/2-Vsxixsi-x+Vsyiysi-y+Vsxizsi-zxsi-x2+ysi-y2+zsi-z21/2+Δnf=rfix,y,z+Δnf,where fe is the carrier frequency of the AIS signal and c is the signal propagation velocity. Lb=x,y,zT is the vessel position vector in the ECEF reference and Vsi=Vsxi,Vsyi,VsziT and Lsi=xsi,ysi,zsiT are velocity vector and position vector of the satellite when transmitting the ith AIS signal, respectively. Δnt is the difference of noise between the two time measurements and Δnf is the difference of noise between the two frequency measurements.It is assumed that the number of signals received by the ship in the visual time of the satellite isN+1; the localization equation matrix based on (1) and (2) can be written as (3)ΔTΔF=rtx,y,zrfx,y,z+nwith(4)rtx,y,z=rt1x,y,z,rt2x,y,z,…,rtNx,y,zTrfx,y,z=rf1x,y,z,rf2x,y,z,…,rfNx,y,zT,where ΔT=Δt1,Δt2,…,ΔtNT is the TDOA measurement vector obtained by synchronization technique and ΔF=Δfr1,Δfr2,…,ΔfrNT is the FDOA measurement vector. n is the measurement noise matrix.


## 2.2. Calculation Based on Least Squares Estimation

On the basis of least squares criterion, the estimator associated with (3) needs to minimize the differences between the measurements and predictions; the equation to be minimized can be written as(5)Cx,y,z=Υ-qx,y,zTN-1Υ-qx,y,z,where N is the noise covariance matrix, Υ=ΔTΔFT is the measurement vector of TDOA/FDOA, and q=rtrfT.The estimated position valuex^,y^,z^=argminx,y,zCx,y,z can be obtained by Gauss-Newton iteration(6)xj+1,yj+1,zj+1T=xj,yj,zjT+JTxj,yj,zjN-1Jxj,yj,zj-1×JTxj,yj,zjN-1Υ-qxj,yj,zj,where matrix J is(7)J=JtJfwith (8)Jt=rt1x,y,z∂xrt1x,y,z∂yrt1x,y,z∂z⋮⋮⋮rtNx,y,z∂xrtNx,y,z∂yrtNx,y,z∂z,Jf=rf1x,y,z∂xrf1x,y,z∂yrf1x,y,z∂z⋮⋮⋮rfNx,y,z∂xrfNx,y,z∂yrfNx,y,z∂z.The initial position in (6) needs to be defined in advance; the method for determining the initial position (shown in Figure 1) is as follows: (1) setting up a grid with the units of 1°. The grid is centered on the midpoint gλo,gφo of the satellite ground trajectory (the track generated during AIS signals transmission) and the range of grid geodetic coordinates gλ,gφ is gλo-β/2≤gλ≤gλo+β/2,gφo-ε/2≤gφ≤gφo+ε/2, where β and ε are the maximum visible longitude and latitude of satellite, respectively. (2) Connecting the start and end points of the satellite ground trajectory and dividing the grid into S1 and S2 (two parts) by extending the connecting line. (3) Searching within each part of grid and selecting two points with mλi,mφi=argmaxSiE=1/fr-f^2 (i=1,2). In the cost function E, fr is the measured frequency of the AIS signal and f^=fc1-VsiTLsi-Lg/cLsi-Lg is the estimated frequency of received signal at the grid point, where Lg=xg,yg,zgT is the position vector of grid point in ECEF coordinate. In this paper, we select “nearest point” to eliminate the false image which may occur in grid searching, that is, taking the point with the shortest distance from the origin of the ship as the optimal position λ0,φ0. The transformation of vessel location from the geodetic coordinate to ECEF coordinate is defined as follows:(9)x=RNcos⁡λcos⁡φy=RNsin⁡λcos⁡φz=RN1-e2sin⁡φ,where λ and φ are the longitude and latitude coordinates of vessel and x,y,z are the ECEF coordinates of vessel. RN=α/1-e2sin2φ is radius of curvature in prime vertical, where e2=0.00669437999013 and α=6378.137km are the square of the first eccentricity and equatorial radius of the earth defined by WGS-84, respectively.Figure 1
Method of determining the initialization vector.


## 3. Localization Method Combining TDOA/FDOA with the Track Forecast

AIS equipped on vessel can obtain a series of dynamic information, such as speed, heading, and turning rate by connecting external sensors. Assuming that the motion state of ship remains constant during the two adjacent positioning points, it is possible with IMM algorithm to estimate the current ship position by using current dynamic information. Considering that the position of the ship at a moment is related to those of the same ship in the previous moment, in the proposed positioning method, the prediction with IMM is taken as a new measurement and added to the TDOA/FDOA joint positioning measurement.


### 3.1. The Motion Model of Vessel

The ship sailing at sea is a slow maneuvering target, with the consideration of the fact that balance between the model accuracy and computational cost, the constant velocity (CV) model, and the constant turn (CT) model are adopted in this paper.The ship state at timek is defined as(10)Xk=xk,vxk,axk,yk,vyk,ayk,zk,vzk,azkT,where xk, yk, zk are vessel position in ECEF coordinate and vxk, vyk, vzk and axk, ayk, azk are vessel velocity and acceleration in ECEF coordinate, respectively.The CV model equation is shown as follows:(11)Xk=ΦCVXk-1+ΓCVWk-1,where the state transition matrix ΦCV=ΦV000ΦV000ΦV with ΦV=1TS0010000 (TS is time interval between adjacent positioning points). The control input matrix ΓCV=ΓV000ΓV000ΓV with ΓV=TS2/2TS0T and Wk-1 is process noise.The CT model equation is shown as follows:(12)Xk=ΦCTXk-1+ΓCTWk-1,where the control input matrix ΓCT=ΓT000ΓT000ΓT with ΓT=TS3/6TS2/2TST. The state transition matrix ΦCT=ΦT000ΦT000ΦT with ΦT=1sin⁡ωTS/ω1-cos⁡ωTS/ω20cos⁡ωTSsin⁡ωTS/ω0-ωsin⁡ωTScos⁡ωTS and ω is the steering rate of vessel.


### 3.2. Vessel Position Predicting Based on IMM Algorithm

A complete cycle of the IMM consists of four operations, namely, input mixing, model filtering, model probability update, and combination. Taking the recorded position at timek-1 and the dynamic information at time k (motion state is supposed to be unchanged during times k-1 and k) as the initial vessel state of each model, the forecasting process of ship position with IMM algorithm is as follows.Step 1 (input mixing).
p
i
j
i
,
j
=
C
V
,
C
T is defined as the Markov transition probability from model i to model j. The mixing probability is computed as follows:(13)μijk-1∣k-1=pij·μik-1mjkwith(14)mjk=∑i=CV,CTpij·μik-1,where μik-1 is probability of mode i at time k-1
The mixed state estimate for modelj is given by(15)X^0jk-1∣k-1=∑i=CV,CTX^ik-1∣k-1μijk-1∣k-1.The predicted covariance corresponding to the above mixed state estimate is given by(16)P0jk-1∣k-1=∑i=CV,CTμijk-1∣k-1Pik-1∣k-1+X^ik-1∣k-1-X^0jk-1∣k-1×X^ik-1∣k-1-X^0jk-1∣k-1T.Step 2 (model filtering).
The extended Kalman filter (EKF) algorithm is adopted in this stage. The mixed state estimate and the predicted covariance corresponding to CV and CT in (15) and (16) are updated by the following two stages.
Time Update
(17)
X
^
k
∣
k
-
1
=
Φ
X
^
k
-
1
P
k
∣
k
-
1
=
Φ
P
k
-
1
Φ
T
+
Γ
Λ
Γ
T
,where P is the covariance of the vessel location prediction, Φ is the state transition matrix, Γ is the control input matrix, and Λ is the system noise variance matrix.
Measurement Update
(18)
K
k
=
P
k
∣
k
-
1
J
T
X
^
k
∣
k
-
1
×
J
X
^
k
∣
k
-
1
P
k
∣
k
-
1
J
T
X
^
k
∣
k
-
1
+
R
-
1
X
^
k
∣
k
=
X
^
k
∣
k
-
1
+
K
k
Z
k
-
q
X
^
k
∣
k
-
1
P
k
∣
k
=
I
-
K
k
J
X
^
k
∣
k
-
1
P
k
∣
k
-
1
,where K corresponds to the gain matrix, R is the observation noise variance matrix, I is the unit matrix, and J is the Jacobian matrix of the measure function q·, as shown in (7).Step 3 (model probability update).
In the IMM algorithm, the updating of the model probability is carried out by calculating the likelihood function matching the model. However, this method is greatly influenced by the presetting model transition probability. In this work, an optimization method is developed to improve the positioning accuracy and can be described as follows.
(1) The positions estimated using EKF based on the CV model and CT model are denoted byPosCV and PosCT, respectively. Let Pos=μCVPosCV+μCTPosCT and calculate the measurement vector Υp at position Pos according to (1) and (2), where μCV and μCTμCV+μCT=1 correspond to the model probabilities of CV and CT models, respectively.
(2) AssumingμCV=η (in this paper, η=0.05), according to Υp determined by the aforementioned method, the corresponding error of measurement vector denoted by σ is calculated by σ=Υ-ΥpT·Υ-Υp.
(3) DefineμCV=μCV+η and repeat the calculation of σ recursively until μCV=1-η. Select μCVμCT=1-μCV corresponding to the minimum error as output probability.Step 4 (combination).
According to the updated model probability in Step3 and predicted state and covariance in Step 2, the combined state and covariance are represented through the following two equations:(19)X^k∣k=∑i=CV,CTμiX^ik∣kP^k∣k=∑i=CV,CTμiPik∣k+X^ik∣k-X^k∣kX^ik∣k-X^k∣kT.In the process of position predicting with IMM algorithm, it is desirable to have a record of ship trajectory used for initialization. Once there are no records available, the method for initializing vessel state is changed as follows: On the basis of the grid search previously described, set up a new grid centered on λ0,φ0 (the optimal position obtained in grid search); the range of the grid is 2°×2° with the unit of 0.5°. According to the cost function E introduced before, search within the new grid and select a point λ0′,φ0′ corresponding to the maximum of E. The direction of λ0,φ0 pointing to λ0′,φ0′ is taken as the vessel heading and λ0,φ0 is regarded as the initial position.


### 3.3. The Proposed Localization Model Design

The predicted ship position with IMM is defined asL^=x^,y^,z^T. The localization model combining this prediction with TDOA/FDOA measurements can be expressed as(20)ΔTΔFL^=rtx,y,zrfx,y,zx,y,zT+nc.Let Υc=ΔTΔFL^T and qc=rtx,y,zrfx,y,zx,y,zTT; on the basis of least squares criterion, the estimated position of the proposed model is achieved by applying the Gauss-Newton algorithm as follows:(21)xj+1,yj+1,zj+1T=xj,yj,zjT+JcTxj,yj,zjNc-1Jcxj,yj,zj-1×JcTxj,yj,zjNc-1Υc-qcxj,yj,zj,where Jc=JI3×3, Nc=N00NM, and NM is the covariance matrix of ship predicted position obtained in IMM algorithm.


## 3.1. The Motion Model of Vessel

The ship sailing at sea is a slow maneuvering target, with the consideration of the fact that balance between the model accuracy and computational cost, the constant velocity (CV) model, and the constant turn (CT) model are adopted in this paper.The ship state at timek is defined as(10)Xk=xk,vxk,axk,yk,vyk,ayk,zk,vzk,azkT,where xk, yk, zk are vessel position in ECEF coordinate and vxk, vyk, vzk and axk, ayk, azk are vessel velocity and acceleration in ECEF coordinate, respectively.The CV model equation is shown as follows:(11)Xk=ΦCVXk-1+ΓCVWk-1,where the state transition matrix ΦCV=ΦV000ΦV000ΦV with ΦV=1TS0010000 (TS is time interval between adjacent positioning points). The control input matrix ΓCV=ΓV000ΓV000ΓV with ΓV=TS2/2TS0T and Wk-1 is process noise.The CT model equation is shown as follows:(12)Xk=ΦCTXk-1+ΓCTWk-1,where the control input matrix ΓCT=ΓT000ΓT000ΓT with ΓT=TS3/6TS2/2TST. The state transition matrix ΦCT=ΦT000ΦT000ΦT with ΦT=1sin⁡ωTS/ω1-cos⁡ωTS/ω20cos⁡ωTSsin⁡ωTS/ω0-ωsin⁡ωTScos⁡ωTS and ω is the steering rate of vessel.


## 3.2. Vessel Position Predicting Based on IMM Algorithm

A complete cycle of the IMM consists of four operations, namely, input mixing, model filtering, model probability update, and combination. Taking the recorded position at timek-1 and the dynamic information at time k (motion state is supposed to be unchanged during times k-1 and k) as the initial vessel state of each model, the forecasting process of ship position with IMM algorithm is as follows.Step 1 (input mixing).
p
i
j
i
,
j
=
C
V
,
C
T is defined as the Markov transition probability from model i to model j. The mixing probability is computed as follows:(13)μijk-1∣k-1=pij·μik-1mjkwith(14)mjk=∑i=CV,CTpij·μik-1,where μik-1 is probability of mode i at time k-1
The mixed state estimate for modelj is given by(15)X^0jk-1∣k-1=∑i=CV,CTX^ik-1∣k-1μijk-1∣k-1.The predicted covariance corresponding to the above mixed state estimate is given by(16)P0jk-1∣k-1=∑i=CV,CTμijk-1∣k-1Pik-1∣k-1+X^ik-1∣k-1-X^0jk-1∣k-1×X^ik-1∣k-1-X^0jk-1∣k-1T.Step 2 (model filtering).
The extended Kalman filter (EKF) algorithm is adopted in this stage. The mixed state estimate and the predicted covariance corresponding to CV and CT in (15) and (16) are updated by the following two stages.
Time Update
(17)
X
^
k
∣
k
-
1
=
Φ
X
^
k
-
1
P
k
∣
k
-
1
=
Φ
P
k
-
1
Φ
T
+
Γ
Λ
Γ
T
,where P is the covariance of the vessel location prediction, Φ is the state transition matrix, Γ is the control input matrix, and Λ is the system noise variance matrix.
Measurement Update
(18)
K
k
=
P
k
∣
k
-
1
J
T
X
^
k
∣
k
-
1
×
J
X
^
k
∣
k
-
1
P
k
∣
k
-
1
J
T
X
^
k
∣
k
-
1
+
R
-
1
X
^
k
∣
k
=
X
^
k
∣
k
-
1
+
K
k
Z
k
-
q
X
^
k
∣
k
-
1
P
k
∣
k
=
I
-
K
k
J
X
^
k
∣
k
-
1
P
k
∣
k
-
1
,where K corresponds to the gain matrix, R is the observation noise variance matrix, I is the unit matrix, and J is the Jacobian matrix of the measure function q·, as shown in (7).Step 3 (model probability update).
In the IMM algorithm, the updating of the model probability is carried out by calculating the likelihood function matching the model. However, this method is greatly influenced by the presetting model transition probability. In this work, an optimization method is developed to improve the positioning accuracy and can be described as follows.
(1) The positions estimated using EKF based on the CV model and CT model are denoted byPosCV and PosCT, respectively. Let Pos=μCVPosCV+μCTPosCT and calculate the measurement vector Υp at position Pos according to (1) and (2), where μCV and μCTμCV+μCT=1 correspond to the model probabilities of CV and CT models, respectively.
(2) AssumingμCV=η (in this paper, η=0.05), according to Υp determined by the aforementioned method, the corresponding error of measurement vector denoted by σ is calculated by σ=Υ-ΥpT·Υ-Υp.
(3) DefineμCV=μCV+η and repeat the calculation of σ recursively until μCV=1-η. Select μCVμCT=1-μCV corresponding to the minimum error as output probability.Step 4 (combination).
According to the updated model probability in Step3 and predicted state and covariance in Step 2, the combined state and covariance are represented through the following two equations:(19)X^k∣k=∑i=CV,CTμiX^ik∣kP^k∣k=∑i=CV,CTμiPik∣k+X^ik∣k-X^k∣kX^ik∣k-X^k∣kT.In the process of position predicting with IMM algorithm, it is desirable to have a record of ship trajectory used for initialization. Once there are no records available, the method for initializing vessel state is changed as follows: On the basis of the grid search previously described, set up a new grid centered on λ0,φ0 (the optimal position obtained in grid search); the range of the grid is 2°×2° with the unit of 0.5°. According to the cost function E introduced before, search within the new grid and select a point λ0′,φ0′ corresponding to the maximum of E. The direction of λ0,φ0 pointing to λ0′,φ0′ is taken as the vessel heading and λ0,φ0 is regarded as the initial position.


## 3.3. The Proposed Localization Model Design

The predicted ship position with IMM is defined asL^=x^,y^,z^T. The localization model combining this prediction with TDOA/FDOA measurements can be expressed as(20)ΔTΔFL^=rtx,y,zrfx,y,zx,y,zT+nc.Let Υc=ΔTΔFL^T and qc=rtx,y,zrfx,y,zx,y,zTT; on the basis of least squares criterion, the estimated position of the proposed model is achieved by applying the Gauss-Newton algorithm as follows:(21)xj+1,yj+1,zj+1T=xj,yj,zjT+JcTxj,yj,zjNc-1Jcxj,yj,zj-1×JcTxj,yj,zjNc-1Υc-qcxj,yj,zj,where Jc=JI3×3, Nc=N00NM, and NM is the covariance matrix of ship predicted position obtained in IMM algorithm.


## 4. Simulation Research

Due to the limitation of the experimental conditions, the satellite ephemeris information is generated by Satellite Tool Kit. The orbital height of the satellite is set to 1000 km and the inclination angle is50°. Simulations have been conducted to evaluate the estimation performance by using 3 signals for positioning (the change of ship position during the reception of the signals is negligible) and the time intervals separating signals are 20 s and 60 s. The measurement vectors of TDOA and FDOA are calculated by (1) and (2). It is supposed that, in the simulation, the time measurements are affected by an additive Gaussian noise constituting independent samples with zero mean and variance vart=30μs2 and the frequency measurements noise is subjected to the Gaussian distribution of zero mean and variance varf=400Hz2.Figure2 shows the estimated vessel trajectories with TDOA/FDOA positioning model and the proposed positioning model. The reference trajectory of vessel is carried out by AIS message reported by an ocean-going ship from Xiamen towards Long Beach in 24 hours and the average speed of ship is 19 knots. Figure 3 shows examples of reference position and estimated position for records availability scenario. The number of positioning results presented is 21 and the average time separating them is approximately 1 hour. From Figures 2 and 3, it can be seen that, compared with TDOA/FDOA positioning model, the estimated trajectory with the proposed method is more in line with the reference trajectory in the case where a record of previous trajectory is used for initialization.Figure 2
Reference and estimated vessel trajectories.Figure 3
Examples of estimated vessel positions for records availability scenario.In Figure3, it is assumed that the reference positions from low latitude to high latitude correspond to numbers 1 to 21 in order. Figure 4 shows the positioning error of these 21 positions in X, Y, Z direction for records availability scenario and Figure 5 shows the distance between the reference and estimated positions in the same scenario. In Figure 5, the average distance error of TDOA/FDOA positioning model is 30.4947 km, the maximum error occurs at the number 1 position, with the distance of 132.9219 km, and the minimum error is 2.8647 km at the number 20 position. Besides, with the condition of records availability, the average distance error of the proposed method is 22.9933 km and the maximum and the minimum error are 41.6187 km corresponding the number 9 position and 8.0532 km corresponding the number 19 position, respectively. As can be observed from Figures 3 and 5, the stability and accuracy of the proposed model (for records availability scenario) are better than those of the TDOA/FDOA positioning model on the whole. Although there is a case where the performance of TDOA/FDOA positioning is superior to the proposed method in individual positions, the large error values of the TDOA/FDOA method estimated at some positions cannot be ignored; because the TDOA/FDOA method is affected by the relative satellite-ship position, it is unstable in overall positioning accuracy.Figure 4
Positioning error inX, Y, Z direction for records availability scenario.Figure 5
Distance between reference and estimated positions for records availability scenario.The examples of reference position and estimated position for records unavailability scenario are shown in Figure6. In Figures 3 and 6, the reference positions are the same, but it is obvious that the proposed method has a better performance in the case of records availability. Figure 7 shows the positioning error of estimated positions in X, Y, Z direction for records availability scenario.Figure 6
Examples of estimated vessel positions for records unavailability scenario.Figure 7
Positioning error inX, Y, Z direction for records unavailability scenario.Figure8 shows the distance between the estimated position and the reference position in different cases. For the scenario of records unavailability, the average distance error of proposed method is 26.2660 km, the maximum error occurs at the number 18 position, with the distance of 168.1150 km, and the minimum error is 2.9289 km at the number 8 position. According to Figures 5 and 8, it can be seen that, for the same reference positions, the average error with different positioning methods in descending order is as follows: 30.4947 km corresponding to TDOA/FDOA positioning method, 26.2660 km corresponding to the proposed method for records unavailability scenario, and 22.9933 km corresponding to the proposed method for records availability scenario. Overall, the performance of the proposed method in this paper is better than TDOA/FDOA positioning method, and, compared with the scenario of records unavailability, the accuracy of the proposed method is higher in the case where records are available.Figure 8
Distance between reference positions and estimated positions.


## 5. Conclusion

A ship positioning method using AIS signals transmitted from satellite is presented in this paper. In the proposed positioning model, an additional measurement obtained by IMM algorithm is added to the TDOA/FDOA measurements. Besides, a probability update method applied in IMM algorithm is designed in this work. The feasibility of the proposed positioning method is verified by simulations. Regardless of algorithm complexity, the performance of the proposed method is better than TDOA/FDOA positioning model, especially in the case where the record of previous ship track is used for initialization.



---

*Source: 8097187-2017-12-03.xml*