Blind separation of complex-valued satellite-AIS data for marine surveillance: a spatial quadratic time-frequency domain approach

In this paper, the problem of the blind separation of complex-valued Satellite-AIS data for marine surveillance is addressed. Due to the speciﬁc properties of the sources under consideration: they are cyclo-stationary signals with two close cyclic frequencies, we opt for spatial quadratic time-frequency domain methods. The use of an additional diversity, the time delay, is aimed at making it possible to undo the mixing of signals at the multi-sensor receiver. The suggested method involves three main stages. First, the spatial generalized mean Ambiguity function of the observations across the array is constructed. Second, in the Ambiguity plane, Delay-Doppler regions of high magnitude are determined and Delay-Doppler points of peaky values are selected. Third, the mixing matrix is estimated from these Delay-Doppler regions using our proposed non-unitary joint zero-(block) diagonalization algorithms as to perform separation.


INTRODUCTION
This paper concerns the spatial automatic identification system (S-AIS) dedicated to marine surveillance by satellite.It covers a larger area than the terrestrial automatic identification system [1], [2].The idea of switching to satellite monitoring was introduced because of the fast and dynamic development of the marine traffic [3][4][5].It was an emergency to adopt a method that operates a global monitoring with reliability, efficiency and robustness.However, this generalization to space involves several phenomena.Among these phenomena, we found: 1.The speed of the satellite movement generates the Doppler effect which creates frequency offsets at the S-AIS signals [6], 2. The propagation delay of the signals and their attenuation due to the satellite altitude [7], 3. When a wide area is covered by the satellite, it certainly includes several traditional AIS cells.In fact, the time propagation of signals transmitted from vessels to the satellite vary according to the position of the Journal Homepage: http://www.iaescore.com/journals/index.php/IJECE IJECE ISSN: 2088-8708 1733 ships and the maximum coverage area of the satellite antenna.Due to these two problems, it mainly affects the organizational mechanism of S-AIS signals.It results a collision data, as illustrated in the Figure 1, issued by vessels located in different AIS cells but received at the antenna of the same satellite [8], [9].
For this reason, we present new approaches to address this problem where the Doppler effect and the propagation delay are also taken into consideration.In fact, to solve the collision problem, few works have focused on blind separation of sources (BSS) methods [10], [11].In [11], Zhou et al. present a multi-user receiver equipped with an array of antennas embedded in Low Orbit Earth (LEO) satellite.The principle of this receiver is to exploit spatial multiplexing in a nonstationary asynchronous context.Indeed, the authors consider the equation below: where is the Schur-Hadamard operator, X = [x 1 , x 2 , . . ., is the matrix of sources and , where ϕ k = e j2π∆f k Ts contains the Doppler frequencies of the sources.
The principle of this method is based on joint diagonalization (JD) of matrices in order to reconstruct the S-AIS sources from separation matrix estimation [12].However, because of the very specific properties of the S-AIS signals (complex and cyclo-stationary with two close cyclic frequencies), we opt for spatial quadratic time-frequency domain methods.Our aim is reshaping the collision problem into BSS problem more simpler than (1).We will show how another type of decomposition matrix named joint zero-diagonalization (JZD) of matrices set resulting from spatial quadratic time-frequency distributions allows the restitution of S-AIS sources.

TRANSMISSION SCHEME 2.1. AIS Frame
The AIS frame is a length of 256 bits and occupies one minute.It is divided into 2250 time slots where one slot equals 26.67 ms [13].Its structure as illustrated in Figure 2  Check Sequence (FCS) (or 16 bits Cyclic Redundancy Code (CRC)) is added to the data information (168 bits) in which a zero is inserted after every five continuous one.The binary sequence {a k } 0≤k≤K of the AIS frame takes the values {−1, +1} since the NRZI encoding is used.Moreover, the modulation specified by S-AIS standard is Gaussian Minimum Shift Keying (GMSK) [14].The encoded message is modulated and transmitted at 9600 bps on 161.975MHz and 162.025MHz frequencies carrier.

GMSK modulation
The resulting sequence after the bit stuffing and NRZI coding procedure is modulated with GMSK which is a frequency-shift keying modulation producing constant-envelope and continuous-phase.Hence, the signal can be written as s g (t) = +∞ k=0 a k g(t − kT s ), where a k are the transmitted symboles, T s is the symbol period and g(t) = 2π log 2 B exp − 2 log 2 (πBt) 2 represents the shaping Gaussian filter where B is the bandwidth of the Gaussian filter.The GMSK modulation is described by the bandwidth-time (BT) product where S-AIS uses BT= 0.4 and T s = 1 9600 s).Making the signal on one of the frequencies carrier f c , produces a signal of spectral characteristic which is adapted to the band-pass channel transmission.The GMSK signal is, thus, expressed as : s(t) = {e −j(2πfct+φ(t)) } = I(t) cos(2πf c t) − Q(t) sin(2πf c t), where {•} is the real part of a complex number, φ(t) = 2πh +∞ k=0 a k g(t − kT s ) is the instantaneous phase of s g (t) where, in the AIS system, the modulation index is theoretically equal to h = 0.5 [15], I(t) (resp.Q(t)) modulates the frequency carrier in phase (resp. in phase quadrature).All steps of the GMSK modulation can be presented in the Figure 3.

PROBLEM STATEMENT: COLLISION & BSS IN INSTANTANEOUS CONTEXT 3.1. Mathematical model of collision problem
The collision problem can be simply expressed as follows: IJECE Vol. 9, No. where x(t) is the received signal by the satellite, s j is the transmitted signal by the j − th vessel, h j , τ j and ∆f j are respectively the coefficients of the channel, the delay and the Doppler shift corresponding to the j − th vessel with J is the number of vessels and n(t) is an additive stationary white Gaussian noise, mutually uncorrelated, independent from the s j , with the variance σ 2 n .

Reshape the collision problem into BSS problem
We show, here, that (2) can be written in BSS nomenclature in which the delay and the Doppler shift caused by the satellite speed are considered.However, before any reformulation, we notice that the mixing matrix considered for S-AIS application is an instantaneous mixture due to the absence of obstacles in the ocean.Thus, we set J = n, the collision problem can be easily modeled in a BSS problem as follows: where T are respectively the (m × 1) observations and noises vectors.The superscript (.) T denotes the transpose operator.
Our developments are based on the following assumptions: Assumption A: The noises n j (t) for all j = 1, . . ., m are stationary, white, zero-mean, mutually uncorrelated random signals and independent from the sources with variance σ 2 n .Assumption B: For each s j of the n sources, there Delay-Doppler points of only one source is present in the Ambiguity plane.Assumption C: The number of sensors m and the number of sources n are both known and m ≥ n to deal with an over-determined model (the under-determined case is outside of the scope in this paper).

PRINCIPLE OF THE PROPOSED METHODS BASED ON THE SPATIAL GENERALIZED MEAN AMBIGUITY FUNCTION
We show, here, how the algorithms proposed in [16], [17] adresses the problem of the separation of instantaneous mixtures of S-AIS data.The principle of the proposed methods are based on three main steps: first, the SGMAF of the observations across the array is constructed.Second, in the Ambiguity plane, Delay-Doppler regions of high magnitude are determined and Delay-Doppler points of peaky values are selected.Third, the mixing matrix is estimated from these Delay-Doppler regions so as to perform separation and to undo the mixing of signals at the multi-sensor receiver.

The Spatial Generalized Mean Ambiguity Function
With regard to BSS, it has been shown that spatial time-frequency distributions are an effective tool when signature of the sources differ in certain points of the time-frequency plan [18].However, in the cyclo-stationary sources case, the Delay-Doppler frequency domain seems to be a more natural field for the re-estimation of sources than the time-frequency domain.As mentioned in [19], the approaches based on information derived from spatial Ambiguity function (SAF) or on SGMAF should be used.In fact, for any vectorial complex signal z(t), the SGMAF is expressed as [20][21][22]: where (s τ,ν z) is the operator of elementary Delay-Doppler translations of z defined by , where R z (t, τ ) stands for the correlation matrix of z(t), E {.} stands for the mathematical expectation operator and superscript (.) H denotes the conjugate transpose operator.Āz (ν, τ ) characterizes the average correlation of all pairs separated by τ in time and by ν in frequency [21], [22].Notice that the diagonal terms of the matrix Āz (ν, τ ) are called auto-terms, while the other ones are called cross-terms.

Selection of peaky Delay-Doppler points
Under the linear data model in (3), the SGMAF of observations across the array at a given Delay-Doppler point is a (m × m) matrix admits the following decomposition: where Ās (ν, τ ) represents the (n × n) SGMAF of sources defined similarly to Āz (ν, τ ) in ( 4) and dt and I m is the m × m identity matrix.It is known that the matrix Ās (ν, τ ) for any τ and ν has no special structure.However, there are some Delay-Doppler points where this matrix has a specific algebraic structure : • Diagonal, for points where each of them corresponds to a single auto-source term for all source signals, • Zero-diagonal for points where each of them correspond to all two by two cross-source term (this structure is exploited because the signature of the sources differ in certain points of the Delay-Doppler plan on the zero-diagonal part (as shown in section 5.).
Our aim is to take advantage of these properties of the Āx (t, ν) in ( 5) since the element of this is no more (zero) diagonal matrices due to the mixing effect in order to estimate the separation matrix B (the pseudo-inverse of matrix H) and restore the unknown sources.We use the detector suggested in [23] (denoted C Ins ) for the instantaneous mixture considered without pre-whitening of the observations.The idea is to find "useful" Delay-Doppler points which consists in keeping Delay-Doppler points with a sufficient energy, then using the rank-one property to detect single cross-source terms (we don't make any assumptions on the knowledge of cyclic frequencies) in the following way: where 1 , 2 are (sufficiently) small positif values and λ max {.} is the largest eigenvalue of a matrix.

Non-unitary joint zero-(block) diagonalization algorithms (NU − JZ(B)D)
The matrices belonging to the set M (whose size is denoted by N m (N m ∈ N * )) all admit a particular structure since they can be decomposed into H Ās (ν, τ )H H with Ās (ν, τ ) a zero-diagonal matrix with only one non null term on its zero-diagonal.One possible way to recover the mixing matrix B is to directly joint zerodiagonalize the matrix set M. It has to be noticed that the recovered sources (after multiplying the observations vector by the estimated matrix B) are obtained up to a permutation (among the classical indetermination of the BSS).Hence, two BSS methods can be derived.The first called JZD CG DD algorithm based on conjugate gradient approach [16].The second JZD LM DD algorithm based on Levenbreg-Marquardt scheme [17].
To tackle that problem, we propose here, to consider the following cost function [16], [17], , where the matrix operator Bdiag (n) {.} is defined as follows: where M is a N × N (N = n(L + L ) where L is the order of the FIR filter and L is the number of delays considered when the convolutif mixture is considered) square matrix whose block components M ij for all i, j = 1, . . ., r are n i × n j matrices (and n 1 + . . .+ n r = N ) denoting by n = (n 1 , n 2 , . . ., n r ).Note that when L = 0, L = 1 we find the instantaneous model since Āx are no more matrices but scalars.Thus, it leads to the minimization of the following cost function: where M i = Āx i is the i − th of the N m matrices belonging to M. We suggest to use conjugate gradient and Levenberg-Marquardt algorithms [16], [17] to minimize the cost function given by Equation .(7) in order to estimate the matrix B ∈ C n×m .It means that B is re-estimated at each iteration m (denoted B (m) or b (m) when the vector b (m) = vec B (m) is considered).The matrix B (or the vector b) is updated according to the following adaptation rule for all m = 1, 2, . . .

Conjugate gradient approach
where µ is a positive small factor called the step-size, d B is the direction of search, β is an exact line search and [16] how the optimal step-size µ opt , ∇ a C ZD (B) and β are calculated at each iteration).
where [.] −1 denotes the inverse of a matrix, λ is positive a small damping factor, is the Hessian matrix of C ZD (B) composed of four complex matrices with: where the operator ⊗ denotes the Kronecker product, K N,M is a square commutation matrix of size N M × N M and T Boff = I N 2 − T Diag , is the N 2 × N 2 "transformation" matrix, with T Diag = diag{vec(BDiag{1 N })}, 1 N is the N × N matrix whose components are all ones, diag{a} is a square diagonal matrix whose diagonal elements are the elements of the vector a, I N 2 = Diag{1 N 2 } is the N 2 × N 2 identity matrix, and Diag{A} is the square diagonal matrix with the same diagonal elements as A.

COMPUTER SIMULATIONS
Computer simulations are performed to illustrate the good behavior of the suggested methods and to compare them with the same kind of existing approach denoted by JZD Chabriel DD proposed in [24] with the Delay-Doppler point C Ins detector.We consider m = 3 mixtures of n = 2 frames of 256 bits correspond to two vessels with different characteristics.The frames are generated according to the S-AIS recommendation as mentioned in the Figure 2 (see also [11], [10]).These frames are encoded with NRZI and modulated in GMSK with a The real part and the imaginary part of their SGMAF is given on the left and on the right of the Figure 4 respectively.Then, the SGMAF of the observations x is then calculated by ( 5) and finally the 100 resulting SGMAF are averaged.We have chosen 1 = 0.07 and 2 = 0.08 for the detector C Ins in order to construct the set M to be joint zero-diagonalized.The signal-to-noise ratio SNR is defined by SNR = 10 log( 1 ) of mean 0 and variance σ 2 n .The selected Delay-Doppler points using the proposed detector are represented in the Figure 5 for SNR = 10 dB and 100 dB.To measure the quality of the estimation, the ensuing error index is used [25] : The SGMAF of the source 1 where (T ) i,j for all i, j ∈ 1, . . ., n is the (i, j)-th element of T = BH.The separation is perfect when the error index I(•) is close to 0 in a linear scale (−∞ in a logarithmic scale).All the displayed results have been averaged over 30 Monte-Carlo trials.We plot, in the Figure 6, the evolution of the error index versus the SNR in order to emphasize the influence of this in the quality of the estimation.All algorithms are initialized using the same initialization suggested in [24].First, we can deduce from the Figure 4 that the diversity in the Delay-Doppler regions is obtained on the zero-diagonal part which supports the use of zero diagonalization algorithms.Then, our analysis are examined on the Figure 6 according to noiseless and noisy contexts .For the noiseless context (when SNR=100 dB), the JZD CG DD and JZD LM DD reach approximately -64 dB and -60 comparing with JZD Chabriel DD method which reaches -20 dB.From this comparison, we have checked the validity of the good behavior of JZD CG DD and JZD LM DD compared to the JZD Chabriel DD approach.Moreover, we observe that the JZD LM DD based on the computation of exact Hessian matrices is more efficient than the JZD CG DD approach.Even in a difficult (noisy) context (for example SNR=15 dB), we note that the best results are generally obtained using the JZD LM DD (-36 dB) then JZD CG DD (-33 dB) especially the JZD LM DD algorithm based on the computation of exact Hessian matrices.It may be concluded that the approaches exploiting the Delay-Doppler diversity of S-AIS signals seem rather promising.Due to its robustness to the noise, it seems to be able to solve the problem of BSS (i.e the collision problem) in a marine surveillance context.

Figure 1 .
Figure 1.Collision problem: The AIS signals from two different SO-TDMA cells received to the satellite antenna at the same time.
contains a training sequence (TS) consisting zero and one which takes 24 bits.The start flag (SF) and the end flag (EF) for information takes 8 bits.A Frame Blind separation of complex-valued satellite-AIS data ... (Omar Cherrak) ISSN: 2088-8708

4. 3 .
Construction of M (set of Delay-Doppler matrices of the observations across the array at the chosen Delay-Doppler points)

Figure 4 .
Figure 4. Left : The SGMAF real part of the S-AIS sources.Right : The SGMAF imaginary part of the S-AIS sources.

Figure 6 .
Figure 6.Comparison of the different methods: evolution of the error index I(T) in dB versus SNR.