# 1 2

**Ξ** is the covariance matrix across the rows, and 2 *n n*

, ~ ,, *m n* **X M** *N* **Ξ Ξ** where # *m n*

**X** into the random vector #

**Xs** as a matrix Gaussian distribution, i.e.

**Σ s** is the covariance matrix

**X** is given as

**M M** *vec* is the expectation matrix and

**M** is the

**Ξ** is the

(14)

**X X** *vec* ,

(15)

*T*

characterization is similar to the MGMRF discussed in Section 2.3.

covariance matrix across the columns. Hence, the pdf of #

1 1

**X** be a random matrix distributed as

**3.2.2 Matrix gaussian distribution**

2 2 2 1 2

(Arnold, 1981). Also, if we stack the matrix #

*n m mn p tr*

**Ξ Ξ**

then ~ , *mn* **X M** *N* **Ξ** where

We model the associated noise process #

across the bands, and 2 2 2 *N N*

**3.2.3 Separable covariance structure**

for 1 *m M* where 1 1 1 1 *N N m m*

Fuentes, 2006) has the form

 2 2 2 2 *N N m m kl*

*mn p*

exp 2 2

expectation matrix, 1 *m m*

2 1 *mn mn*

Let #

matrix #

$$\text{cov}\left(X\_{\mathbf{s},(k,l)}, X\_{\mathbf{s},(k,l)} \middle| \mathbf{L}; \Theta\right) = \sigma\_{lk}^{(1)}\left(L\_{\mathbf{s}}\right)\sigma\_{ll}^{(2)}\left(L\_{\mathbf{s}}\right) \tag{19}$$

$$\text{cov}\left(Y\_{\mathbf{s},(k,l)}, Y\_{\mathbf{s},(k,l)} \Big| \mathbf{L}; \Theta\right) = \sigma\_{kk}^{(1)}\left(L\_{\mathbf{s}}\right)\sigma\_{ll}^{(2)}\left(L\_{\mathbf{s}}\right) \,. \tag{20}$$

This corresponds to the product of the variance associated with the reflectance at the *th k* spectral band 1 *kk L***s** and the variance associated with the reflectance at the *th l* temporal slot 2 *ll L***s** . Likewise, the cross-covariance is given as (Arnold, 1981):

$$\text{cov}\left(X\_{\mathbf{s},(i,j)}, X\_{\mathbf{s},(k,l)} \middle| \mathbf{L}; \Theta\right) = \sigma\_{ik}^{(1)}\left(L\_{\mathbf{s}}\right)\sigma\_{jl}^{(2)}\left(L\_{\mathbf{s}}\right) \tag{21}$$

$$\text{cov}\left(Y\_{\mathbf{s},(i,j)}, Y\_{\mathbf{s},(k,l)} \middle| \mathbf{L}; \Theta\right) = \sigma\_{ik}^{(1)}\left(L\_{\mathbf{s}}\right)\sigma\_{jl}^{(2)}\left(L\_{\mathbf{s}}\right) \,. \tag{22}$$

This corresponds to the product of the covariance associated with the reflectance at the *th i* and the *th k* spectral band 1 *ik L***s** and the covariance associated with the reflectance at the *th j* and the *th l* temporal slot 2 *jl L***s** .

The number of parameters in the unpatterned covariance matrix is *N N NN NN* 1 2 12 12 1 2 . On the other hand, the number of parameters for a separable covariance matrix is *NN NN* 11 22 1 12 , which has fewer parameters compared to its non-separable counterpart.

#### **3.2.4 Separable of interaction matrix structure**

We can also model the interaction matrix coefficients with a separable structure for all **r** and 1 *m M* of the form

$$\boldsymbol{\Theta}\_{\mathbf{r}}(m) = \boldsymbol{\Theta}\_{\mathbf{r}}^{(2)}(m) \otimes \boldsymbol{\Theta}\_{\mathbf{r}}^{(1)}(m) \tag{23}$$

where 1 1 1 *N N m* **θr** is the interaction matrix across the bands and 2 2 2 *N N m* **θr** is the interaction matrix across time. In the next section, the interaction matrix coefficient **θr** *m* can be made separable for **r** and1 *m M* provided that **Σ***m* is separable. Furthermore, if **Σ***m* is separable, then the following is the resulting statistical characterization of **Xs** :

$$E\left[\mathbf{X}\_{\sf s} \, \middle| \, \mathbf{L} ; \Theta \right] = \mathbf{0}\_{\sf N \times 1} \tag{24}$$

$$\operatorname{cov}\left(\mathbf{X}\_{\mathsf{s}},\mathbf{X}\_{\mathsf{s}-\mathsf{r}}\,\middle|\,\mathsf{L}\star\Theta\right) = \begin{vmatrix} \mathbf{Z}^{(2)}(\mathrm{L}\_{\mathsf{s}})\otimes\mathbf{Z}^{(1)}(\mathrm{L}\_{\mathsf{s}}) & \mathbf{r} = \mathbf{0}\_{p\times 1} \\ -\left(-\mathsf{G}\_{\mathsf{s}}^{(2)}(\mathrm{L}\_{\mathsf{s}})\boldsymbol{\Sigma}^{(2)}(\mathrm{L}\_{\mathsf{s}})\cdot\mathbf{1}\_{\{\iota\_{\mathsf{s}}=\iota\_{\mathsf{s}\cdot\mathsf{s}}\}}\right)\otimes\left(-\mathsf{G}\_{\mathsf{s}}^{(1)}(\mathrm{L}\_{\mathsf{s}})\boldsymbol{\Sigma}^{(1)}(\mathrm{L}\_{\mathsf{s}})\cdot\mathbf{1}\_{\{\iota\_{\mathsf{s}}=\iota\_{\mathsf{s}\cdot\mathsf{s}}\}}\right) & \mathbf{r}\in\mathcal{N}\;\mathsf{d}\;\mathsf{d}\;\mathsf{d}\;\mathsf{w}\,\mathsf{s} \end{vmatrix} \tag{25}$$

Estimation of the Separable MGMRF Parameters for Thematic Classification 107

This section provides important characterizations enable us to derive the estimators of the GMRF parameters in the next section. We present a simple, yet powerful, method to derive the MPL estimators of the mean and the interaction matrix. Finally, new problems arise in estimating the multivariate observation GMRFs, which were not encountered in the

In this section, a method of deriving the MPL estimators for the mean and the vectorized interaction coefficients are presented regardless of separability. The MPL estimator of the interaction matrix coefficients can be derived by taking the matrix derivative of the log of the pseudo-likelihood function with respect to the interaction matrix coefficient or with

> *f f vec vec*

evaluate. The following proposition provides a simple way of deriving the MPL estimators, where the estimator is either the mean or the vectorized interaction matrix coefficient

mean or the vectorized interaction matrix coefficient. Suppose that **Xs** can be expressed in

covariance matrix **Σ***m* , 1 *m M* is known, then the MPL estimator for **Φ** *m* ,

**Q Σ X 0**

*Proof* From (7) and (9), the log pseudo-likelihood of the image random field conditional

**Θ L Σ X Σ X**

**s s ss**

1 log log 2 log . 2

Taking the gradient of the log pseudo-likelihood function in (33) with respect to **Φ** *m* for

1 . *T*

*q*

1

**s s**

*m*

*m*

**s**

1

*m m*

**s**

*M*

*PL N*

**P P**

**X X** (30)

**X PQ ss s s Φ***L* (31)

(32)

(33)

1

*T*

*L L*

**Q Q s s Θ L** is independent of *L* **Φ s** , and the

**X P** . The latter expression is preferred, since it is easier to

**Φ** , 1 *m M* be a vector of parameters which is either the

**5.1 MPL-based method technique of deriving mean and the interaction matrix**

respect to its vectorized version from the equivalence relation (Neudecker, 1969)

**5. Important MGMRF specifications**

univariate case, are discussed.

where *f* **X** , and , *m n*

*Proposition 1* Let *q* 1

*m*

**P P s s Θ L** , *N q*

1 *m M* is obtained by solving the equation

(Navarro, et. al., 2009).

where *N* 1

to the thematic map is given as

1 *m M* , and equating to *q*1 **0** yields

the form

**estimators**

$$\operatorname{cov}\left(\mathbf{X}\_{\mathbf{s}} \, \mathbf{Y}\_{\mathbf{s}-\mathbf{r}} \, \middle| \, \mathbf{L} \, \middle| \, \Theta\right) = \boldsymbol{\Sigma}^{(2)}\left(L\_{\mathbf{s}}\right) \cdot \mathbf{1}\_{\left\{L\_{\mathbf{s}} = l\_{\mathbf{s}-\mathbf{r}}\right\}} \otimes \boldsymbol{\Sigma}^{(1)}\left(L\_{\mathbf{s}}\right) \cdot \mathbf{1}\_{\left\{l\_{\mathbf{s}} = l\_{\mathbf{s}-\mathbf{r}}\right\}}.\tag{26}$$

The covariance matrix, from the above equation, cov , ; **XX L s sr Θ** has a separable structure between the spectral domain and temporal dimensions. It has a form analogous to that of what is shown in (4) through (6), which is intuitively appealing.

The number of parameters in the unpatterned interaction matrix coefficient is 2 22 1 2 *N NN* . On the other hand, the number of parameters for the separable interaction matrix coefficient is 2 2 *N N* 1 2 , which has fewer parameters compared to its non-separable counterpart.

#### **3.2.5 Separable mean structure**

Likewise, we can also model the mean with a separable structure of the form

$$
\mathfrak{u}(m) = \mathfrak{u}^{(2)}(m) \otimes \mathfrak{u}^{(1)}(m) \tag{27}
$$

for 1 *m M* where 1 1 *N* 1 *m* **μ** is the mean across the bands and 2 2 *N* 1 *m* **μ** is the mean across time. The number of parameters in the unpatterned mean vector is *N NN* 1 2 . On the other hand, the number of parameters for the separable mean vector is *N N* 1 2 which has fewer number of parameters compared to its non-separable counterpart.

#### **3.2.6 Hybrid separable structure**

Finally, we can model the GMRF parameters as having a hybrid separability structure, that is, some of its parameters are separable while the rest are not. Hence, there are eight combinations to consider. As shown in Section 5.2, it is impossible to model a separable interaction matrix with a non-separabable matrix. This leave us six cases to consider in this study.

#### **4. Estimation of thematic map parameters**

The MPLE of **φ** is obtained by taking the derivative of log *PL***φ** with respect to *m* 1 *m M a* and *b* **r r** , then equating to zero (Li, 1995). Accordingly, the estimators are obtained numerically by solving the following set of simultaneous nonlinear equations:

$$\sum\_{\mathbf{s}\in\mathcal{S}} \frac{\exp\left(a\_{m} + \sum\_{\mathbf{r}\in\mathcal{N}} b\_{\mathbf{r}} \cdot V\left(L\_{\mathbf{s}} = m, L\_{\mathbf{s}-\mathbf{r}}\right)\right)}{\sum\_{l=1}^{M} \exp\left(a\_{l} + \sum\_{\mathbf{r}\in\mathcal{N}} b\_{\mathbf{r}} \cdot V\left(L\_{\mathbf{s}} = l, L\_{\mathbf{s}-\mathbf{r}}\right)\right)} = \sum\_{\mathbf{s}\in\mathcal{S}} \mathbf{1}\_{\left\{L\_{\mathbf{s}} = m\right\}} \,\forall a\_{m}, \; 1 \le m \le M\tag{28}$$

$$\sum\_{\mathbf{s}\in\mathcal{S}} \frac{\sum\_{l=1}^{M} \exp\left(a\_{l} + \sum\_{\mathbf{t}\in\mathcal{N}} b\_{\mathbf{t}} \cdot V\left(L\_{\mathbf{s}} = l, L\_{\mathbf{s}-\mathbf{t}}\right)\right) \cdot V\left(L\_{\mathbf{s}} = l, L\_{\mathbf{s}-\mathbf{t}}\right)}{\sum\_{l=1}^{M} \exp\left(a\_{l} + \sum\_{\mathbf{t}\in\mathcal{N}} b\_{\mathbf{t}} \cdot V\left(L\_{\mathbf{s}} = l, L\_{\mathbf{s}-\mathbf{t}}\right)\right)} = \sum\_{\mathbf{s}\in\mathcal{S}} V\left(L\_{\mathbf{s}}, L\_{\mathbf{s}-\mathbf{t}}\right) \quad \forall b\_{\mathbf{r}}, \ \mathbf{r} \in \mathcal{N} \cdot \tag{29}$$

## **5. Important MGMRF specifications**

106 Remote Sensing – Advanced Techniques and Platforms

The covariance matrix, from the above equation, cov , ; **XX L s sr Θ** has a separable structure between the spectral domain and temporal dimensions. It has a form analogous to

The number of parameters in the unpatterned interaction matrix coefficient is 2 22

is 2 2 *N N* 1 2 , which has fewer parameters compared to its non-separable counterpart.

Likewise, we can also model the mean with a separable structure of the form

On the other hand, the number of parameters for the separable interaction matrix coefficient

 2 1

1 1 *N* 1 *m* **μ** is the mean across the bands and

the mean across time. The number of parameters in the unpatterned mean vector is *N NN* 1 2 . On the other hand, the number of parameters for the separable mean vector is *N N* 1 2 which has fewer number of parameters compared to its non-separable counterpart.

Finally, we can model the GMRF parameters as having a hybrid separability structure, that is, some of its parameters are separable while the rest are not. Hence, there are eight combinations to consider. As shown in Section 5.2, it is impossible to model a separable interaction matrix with a non-separabable matrix. This leave us six cases to consider in this

The MPLE of **φ** is obtained by taking the derivative of log *PL***φ** with respect to

obtained numerically by solving the following set of simultaneous nonlinear equations:

**r s sr**

*M L m m*

**r s sr**

**t s st s sr**

**t s st**

 1

, then equating to zero (Li, 1995). Accordingly, the estimators are

**1**

, 1

*a mM*

*VL L b*

**s sr r**

, , .

**r**

(29)

**XY L s sr Θ Σ s s 1 Σ 1** (26)

2 1 cov , ; . *L L L L L L* **s sr s sr**

**μμ μ** *mmm* (27)

1 2 *N NN* .

2 2 *N* 1 *m* **μ** is

(28)

**3.2.5 Separable mean structure**

**3.2.6 Hybrid separable structure**

**4. Estimation of thematic map parameters**

exp ,

**r**

**r**

**t**

*a b V L mL*

*a b V L lL*

*a b V L lL V L lL*

*a b V L lL*

**s s**

**s s**

 **s**

exp ,

exp , ,

exp ,

*m*

*l*

*l*

1

*l*

*l*

**t**

*M*

1

*l*

for 1 *m M* where

study.

*m* 1 *m M a* and *b* **r r**

*M*

*l*

that of what is shown in (4) through (6), which is intuitively appealing.

This section provides important characterizations enable us to derive the estimators of the GMRF parameters in the next section. We present a simple, yet powerful, method to derive the MPL estimators of the mean and the interaction matrix. Finally, new problems arise in estimating the multivariate observation GMRFs, which were not encountered in the univariate case, are discussed.

#### **5.1 MPL-based method technique of deriving mean and the interaction matrix estimators**

In this section, a method of deriving the MPL estimators for the mean and the vectorized interaction coefficients are presented regardless of separability. The MPL estimator of the interaction matrix coefficients can be derived by taking the matrix derivative of the log of the pseudo-likelihood function with respect to the interaction matrix coefficient or with respect to its vectorized version from the equivalence relation (Neudecker, 1969)

$$\frac{\partial f}{\partial \mathbf{X}} = \mathbf{P} \Leftrightarrow \frac{\partial f}{\partial \text{vec}(\mathbf{X})} = \text{vec}(\mathbf{P}) \tag{30}$$

where *f* **X** , and , *m n* **X P** . The latter expression is preferred, since it is easier to evaluate. The following proposition provides a simple way of deriving the MPL estimators, where the estimator is either the mean or the vectorized interaction matrix coefficient (Navarro, et. al., 2009).

*Proposition 1* Let *q* 1 *m* **Φ** , 1 *m M* be a vector of parameters which is either the mean or the vectorized interaction matrix coefficient. Suppose that **Xs** can be expressed in the form

$$\mathbf{X}\_{\mathbf{s}} = \mathbf{P}\_{\mathbf{s}} - \mathbf{Q}\_{\mathbf{s}} \boldsymbol{\Phi} \left( L\_{\mathbf{s}} \right) \tag{31}$$

where *N* 1 **P P s s Θ L** , *N q* **Q Q s s Θ L** is independent of *L* **Φ s** , and the covariance matrix **Σ***m* , 1 *m M* is known, then the MPL estimator for **Φ** *m* , 1 *m M* is obtained by solving the equation

$$\sum\_{\mathbf{s}\in\mathfrak{S}(m)} \mathbf{Q}\_{\mathbf{s}}^{T} \Sigma^{-1}(m) \mathbf{X}\_{\mathbf{s}} = \mathbf{0}\_{q \times 1} \,. \tag{32}$$

*Proof* From (7) and (9), the log pseudo-likelihood of the image random field conditional to the thematic map is given as

$$\log \text{PL}\left(\Theta \middle| \mathbf{L}\right) = -\frac{1}{2} \sum\_{m=1}^{M} \sum\_{\mathbf{s} \in \mathcal{S}(m)} \left[ N \log 2\pi + \log \left| \Sigma(L\_{\mathbf{s}}) \right| + \mathbf{X}\_{\mathbf{s}}^{T} \Sigma^{-1}(L\_{\mathbf{s}}) \mathbf{X}\_{\mathbf{s}} \right]. \tag{33}$$

Taking the gradient of the log pseudo-likelihood function in (33) with respect to **Φ** *m* for 1 *m M* , and equating to *q*1 **0** yields

Estimation of the Separable MGMRF Parameters for Thematic Classification 109

1 1 2 12 1 2 1

*mmm m m m*

*T*

**θ Σθ Σ r r** (43)

**θ Σθ Σ r r** , (44)

,1 *m M* and thus these cases are not possible.

(42)

,

(45)

1

*mmm mmm*

*T T*

 1 1 11 1 *T m mmm*

 1 2 22 2 *T m mmm*

By considering the hybrid separability cases which involve a separable interaction matrix and a non-separable covariance matrix, the expression *T* 1 *mm m* **Σθ Σ r** is not separable, in general. This implies that **θr** *m* cannot be expressed in the form

This section proposes an estimation procedure for the GMRF parameters for both separable

**Proposition 2**Assume that the interaction matrix coefficients **θr** *m* for **r**

*m m mm*

**s r r**

**s r r μ I θ 1 Σ I θ 1**

1 *m M* and the covariance matrices **Σ***m* for 1 *m M* are known.

*T T*

**Σ Σθ θ Σ Σ**

**θ θ Σ Σθ θ Σ Σ**

The identification of **θr** *m* is completely specified from (39) if we take

**r r r r**

which is analogous to the relation in (39).

2 1 **θθ θ rr r** *mmm* for *S* **r**

**6. GMRF parameter estimation**

**6.1 Mean parameter estimation**

parameters are estimated as follows:

a. Non-Separable Case:

for 1 *m M* .

b. Separable Case:

and non-separable cases based on the MPL.

*m*

*m*

In addition, if we assume the following for 1 *m M* :

1 ˆ

1 2 1 2 12 1 2 1

*m m m mm m m m*

$$\mathbf{0}\_{q\times 1} = \frac{\partial}{\partial \mathbf{\Theta}(m)} \log P L\{\boldsymbol{\Theta} | \mathbf{L}\} = -\frac{1}{2} \sum\_{l=1}^{M} \sum\_{\mathbf{s} \in \mathbf{\mathcal{S}}(l)} \frac{\partial}{\partial \mathbf{\Theta}(m)} \mathbf{X}\_{\mathbf{s}}^{T} \boldsymbol{\Sigma}^{-1}(L\_{\mathbf{s}}) \mathbf{X}\_{\mathbf{s}}.\tag{34}$$

Since

$$\begin{split} \mathbf{X}\_{\sf s}^{T} \boldsymbol{\Sigma}^{-1} \left( \mathbf{L}\_{\sf s} \right) \mathbf{X}\_{\sf s} &= \left( \mathbf{P}\_{\sf s} - \mathbf{Q}\_{\sf s} \boldsymbol{\Phi} \left( \mathbf{L}\_{\sf s} \right) \right)^{\sf T} \boldsymbol{\Sigma}^{-1} \left( \mathbf{L}\_{\sf s} \right) \left( \mathbf{P}\_{\sf s} - \mathbf{Q}\_{\sf s} \boldsymbol{\Phi} \left( \mathbf{L}\_{\sf s} \right) \right) \\ &= \mathbf{P}\_{\sf s}^{T} \boldsymbol{\Sigma}^{-1} \left( \mathbf{L}\_{\sf s} \right) \mathbf{P}\_{\sf s} - 2 \mathbf{P}\_{\sf s}^{T} \boldsymbol{\Sigma}^{-1} \left( \mathbf{L}\_{\sf s} \right) \mathbf{Q}\_{\sf s} \boldsymbol{\Phi} \left( \mathbf{L}\_{\sf s} \right) + \boldsymbol{\Phi}^{T} \left( \mathbf{L}\_{\sf s} \right) \mathbf{Q}\_{\sf s}^{T} \boldsymbol{\Sigma}^{-1} \left( \mathbf{L}\_{\sf s} \right) \mathbf{Q}\_{\sf s} \boldsymbol{\Phi} \left( \mathbf{L}\_{\sf s} \right) \end{split} \tag{35}$$

then taking the gradient in (34) with respect to **Φ** yields

$$\begin{split} \frac{\partial}{\partial \mathbf{Q}} \mathbf{X}\_{\mathbf{s}}^{T} \boldsymbol{\Sigma}^{-1} (\boldsymbol{L}\_{\mathbf{s}}) \mathbf{X}\_{\mathbf{s}} &= 2 \mathbf{Q}\_{\mathbf{s}}^{T} \boldsymbol{\Sigma}^{-1} (\boldsymbol{L}\_{\mathbf{s}}) \mathbf{P}\_{\mathbf{s}} \mathbf{1}\_{\{\boldsymbol{L}\_{\mathbf{s}} = \boldsymbol{m}\}} - 2 \mathbf{Q}\_{\mathbf{s}}^{T} \boldsymbol{\Sigma}^{-1} (\boldsymbol{L}\_{\mathbf{s}}) \mathbf{Q}\_{\mathbf{s}} \boldsymbol{\Theta} (\boldsymbol{L}\_{\mathbf{s}}) \mathbf{1}\_{\{\boldsymbol{L}\_{\mathbf{s}} = \boldsymbol{m}\}} \\ &= 2 \mathbf{Q}\_{\mathbf{s}}^{T} \boldsymbol{\Sigma}^{-1} (\boldsymbol{L}\_{\mathbf{s}}) \left( \mathbf{P}\_{\mathbf{s}} - \mathbf{Q}\_{\mathbf{s}} \boldsymbol{\Theta} (\boldsymbol{L}\_{\mathbf{s}}) \right) \mathbf{1}\_{\{\boldsymbol{L}\_{\mathbf{s}} = \boldsymbol{m}\}} \\ &= 2 \mathbf{Q}\_{\mathbf{s}}^{T} \boldsymbol{\Sigma}^{-1} (\boldsymbol{L}\_{\mathbf{s}}) \mathbf{X}\_{\mathbf{s}} \mathbf{1}\_{\{\boldsymbol{L}\_{\mathbf{s}} = \boldsymbol{m}\}} \end{split} \tag{36}$$

Finally, substituting the result of (36) into (34) gives us the identity

$$\mathbf{O}\_{q\times 1} = \sum\_{l=1}^{M} \sum\_{\mathbf{s}\in\mathcal{S}(l)} \mathbf{Q}\_{\mathbf{s}}^{T} \boldsymbol{\Sigma}^{-1} \left(L\_{\mathbf{s}}\right) \mathbf{X}\_{\mathbf{s}} \mathbf{1}\_{\{l\_{\mathbf{s}} = m\}} = \sum\_{\mathbf{s}\in\mathcal{S}(m)} \mathbf{Q}\_{\mathbf{s}}^{T} \boldsymbol{\Sigma}^{-1} \left(m\right) \mathbf{X}\_{\mathbf{s}}.\tag{37}$$

#### **5.2 Interaction matrix identities**

From the covariance identity

$$\operatorname{cov}\left(\mathbf{X}\_{\text{s}},\mathbf{X}\_{\text{s}-\text{r}} \middle| \mathbf{L};\Theta\right) = \operatorname{cov}^{T}\left(\mathbf{X}\_{\text{s}-\text{r}},\mathbf{X}\_{\text{s}} \middle| \mathbf{L};\Theta\right) \tag{38}$$

(Ravishanker and Dey, 2002), from (5), we obtain the following relationship:

$$
\boldsymbol{\Theta}\_{-\mathbf{r}}\left(\boldsymbol{L}\_{\mathbf{s}}\right) = \boldsymbol{\Sigma}\left(\boldsymbol{L}\_{\mathbf{s}}\right)\boldsymbol{\Theta}\_{\mathbf{r}}^{T}\left(\boldsymbol{L}\_{\mathbf{s}}\right)\boldsymbol{\Sigma}^{-1}\left(\boldsymbol{L}\_{\mathbf{s}}\right).\tag{39}
$$

One consequence of this result is that **Xs** can be written as follows:

$$\mathbf{X}\_{\mathbf{s}} = \left(\mathbf{Y}\_{\mathbf{s}} - \boldsymbol{\mathfrak{m}}(L\_{\mathbf{s}})\right) - \sum\_{\mathbf{r}\mathbf{s}, \boldsymbol{\mathcal{N}}\_{\mathbf{s}}} \left[\boldsymbol{\mathfrak{e}}\_{\mathbf{r}}(L\_{\mathbf{s}})\mathbf{1}\_{\left(L\_{\mathbf{s}} = l\_{\mathbf{s} + \mathbf{r}}\right)}\left(\mathbf{Y}\_{\mathbf{s} - \mathbf{r}} - \boldsymbol{\mathfrak{m}}(L\_{\mathbf{s}})\right) + \boldsymbol{\mathfrak{s}}\,\Sigma(L\_{\mathbf{s}})\boldsymbol{\mathfrak{e}}\_{\mathbf{r}}^{T}(L\_{\mathbf{s}})\Sigma^{-1}(L\_{\mathbf{s}})\mathbf{1}\_{\left(L\_{\mathbf{s}} = l\_{\mathbf{s} + \mathbf{r}}\right)}\left(\mathbf{Y}\_{\mathbf{s} + \mathbf{r}} - \boldsymbol{\mathfrak{m}}(L\_{\mathbf{s}})\right)\right] \tag{40}$$

where *S* , a subset of which represents the symmetric neighborhood set (Kashyap and Chellappa, 1983), is defined as follows: *S S* **r r** and . *S S* � **r r**

Another consequence of (39) are the specifications of the interaction matrices in the separable case. If the interaction matrices are modeled as separable, then by (39), we obtain

$$\boldsymbol{\Theta}\_{-\mathbf{r}}\left(m\right) = \boldsymbol{\Theta}\_{-\mathbf{r}}^{(2)}\left(m\right) \otimes \boldsymbol{\Theta}\_{-\mathbf{r}}^{(1)}\left(m\right) = \boldsymbol{\Sigma}\left(m\right)\left(\boldsymbol{\Theta}\_{\mathbf{r}}^{(2)}\left(m\right) \otimes \boldsymbol{\Theta}\_{\mathbf{r}}^{(1)}\left(m\right)\right)^{T}\boldsymbol{\Sigma}^{-1}\left(m\right) = \boldsymbol{\Sigma}\left(m\right)\boldsymbol{\Theta}\_{\mathbf{r}}^{T}\left(m\right)\boldsymbol{\Sigma}^{-1}\left(m\right) \tag{41}$$

for 1 *m M* . The RHS of (40) can be made separable if **Σ** *m* is also separable. Hence,

$$\begin{split} \boldsymbol{\Theta}\_{\mathbf{r}}^{(2)}(m) \otimes \boldsymbol{\Phi}\_{\mathbf{r}}^{(1)}(m) &= \left(\boldsymbol{\Sigma}^{(2)}(m) \otimes \boldsymbol{\Sigma}^{(1)}(m)\right) \left(\boldsymbol{\Theta}\_{\mathbf{r}}^{(2)}(m) \otimes \boldsymbol{\Phi}\_{\mathbf{r}}^{(1)}(m)\right)^{\mathrm{T}} \left(\boldsymbol{\Sigma}^{(2)}(m) \otimes \boldsymbol{\Sigma}^{(1)}(m)\right)^{-1} \\ &= \left(\boldsymbol{\Sigma}^{(2)}(m) \otimes \boldsymbol{\Sigma}^{(1)}(m)\right) \left(\boldsymbol{\Phi}\_{\mathbf{r}}^{(2)T}(m) \otimes \boldsymbol{\Phi}\_{\mathbf{r}}^{(1)T}(m)\right) \left(\boldsymbol{\Sigma}^{(2)}(m)\right)^{-1} \otimes \left(\boldsymbol{\Sigma}^{(1)}(m)\right)^{-1} \\ &= \boldsymbol{\Sigma}^{(2)}(m) \boldsymbol{\Phi}\_{\mathbf{r}}^{(2)T}(m) \left(\boldsymbol{\Sigma}^{(2)}(m)\right)^{-1} \otimes \boldsymbol{\Sigma}^{(1)}(m) \boldsymbol{\Phi}\_{\mathbf{r}}^{(1)T}(m) \left(\boldsymbol{\Sigma}^{(1)}(m)\right)^{-1}. \end{split} \tag{42}$$

The identification of **θr** *m* is completely specified from (39) if we take

$$\boldsymbol{\Theta}\_{-\mathbf{r}}^{(1)}(m) = \boldsymbol{\Sigma}^{(1)}(m)\boldsymbol{\Theta}\_{\mathbf{r}}^{(1)T}(m)\left(\boldsymbol{\Sigma}^{(1)}(m)\right)^{-1} \tag{43}$$

$$
\Theta\_{-\mathbf{r}}^{(2)}(m) = \boldsymbol{\Sigma}^{(2)}(m)\boldsymbol{\Theta}\_{\mathbf{r}}^{(2)\text{Tr}}\left(m\right)\left(\boldsymbol{\Sigma}^{(2)}(m)\right)^{-1}\tag{44}
$$

which is analogous to the relation in (39).

108 Remote Sensing – Advanced Techniques and Platforms

*m m*

**Φ Φ**

11 1

2 2

2 .

*l l m*

**s s**

*L*

**s ss**

**Q Σ X 1**

1 1

*TT T*

*L LL L*

**s ss s s s s s s s**

**X Σ X PQ Φ Σ P Q Φ**

2

*T*

*T*

Finally, substituting the result of (36) into (34) gives us the identity

*q L m*

**s ss s**

(Ravishanker and Dey, 2002), from (5), we obtain the following relationship:

One consequence of this result is that **Xs** can be written as follows:

1

*M*

**0 Θ L X Σ X**

1

1 1

then taking the gradient in (34) with respect to **Φ** yields

*T T*

1

**5.2 Interaction matrix identities**

*S*

Chellappa, 1983), is defined as follows: *S S* **r r**

**r**

*S* , a subset of

where

From the covariance identity

*q*

**Φ**

Since

*M*

1 1 log . 2

1 1 1 2

**P Σ P P Σ Q Φ Φ Q Σ Q Φ**

1 1

*L m*

*L m*

*T T*

**s**

*L L L L*

*L L*

*T T T T*

1

**s sr s sr s s s r s sr s s r s s sr s**

Another consequence of (39) are the specifications of the interaction matrices in the separable case.  The RHS of (40) can be made separable if **Σ** *m* is also separable.  *S S*

**r r**

*L m*

**s**

cov , ; cov , ; *T* **XX L s sr Θ X XL sr s Θ** (38)

*T*

which represents the symmetric neighborhood set (Kashyap and

1 . *T L LL L* **θ Σθ Σ rs sr s s** (39)

�

**s s**

By considering the hybrid separability cases which involve a separable interaction matrix and a non-separable covariance matrix, the expression *T* 1 *mm m* **Σθ Σ r** is not separable, in general. This implies that **θr** *m* cannot be expressed in the form 2 1 **θθ θ rr r** *mmm* for *S* **r** ,1 *m M* and thus these cases are not possible.

#### **6. GMRF parameter estimation**

This section proposes an estimation procedure for the GMRF parameters for both separable and non-separable cases based on the MPL.

#### **6.1 Mean parameter estimation**

**Proposition 2**Assume that the interaction matrix coefficients **θr** *m* for **r** , 1 *m M* and the covariance matrices **Σ***m* for 1 *m M* are known. Then the mean parameters are estimated as follows:

a. Non-Separable Case:

$$
\hat{\boldsymbol{\mu}}\left(\boldsymbol{m}\right) = \left[\sum\_{\boldsymbol{\kappa}\in\mathcal{S}(m)} \left(\mathbf{I}\_N - \sum\_{\boldsymbol{\kappa}\prec\boldsymbol{\kappa}'} \boldsymbol{\Phi}\_{\boldsymbol{\kappa}}\left(m\right) \mathbf{1}\_{\{\boldsymbol{\ell}\_{\boldsymbol{\kappa}\_{\boldsymbol{\kappa}}-\boldsymbol{m}}\}}\right)^T \boldsymbol{\Sigma}^{-1}\left(m\right) \left(\mathbf{I}\_N - \sum\_{\boldsymbol{\kappa}\prec\boldsymbol{\kappa}'} \boldsymbol{\Phi}\_{\boldsymbol{\kappa}}\left(m\right) \mathbf{1}\_{\{\boldsymbol{\ell}\_{\boldsymbol{\kappa}\_{\boldsymbol{\kappa}}}-\boldsymbol{m}\}}\right)\right]^{-1}.\tag{45}
$$

$$
\left[\sum\_{\boldsymbol{\kappa}\in\mathcal{S}(m)} \left(\mathbf{I}\_N - \sum\_{\boldsymbol{\kappa}\prec\boldsymbol{\kappa}'} \boldsymbol{\Phi}\_{\boldsymbol{\kappa}}\left(m\right) \mathbf{1}\_{\{\boldsymbol{\ell}\_{\boldsymbol{\kappa}\_{\boldsymbol{\kappa}}}-\boldsymbol{m}\}}\right)^T \boldsymbol{\Sigma}^{-1}\left(m\right) \left(\mathbf{Y}\_{\boldsymbol{\kappa}} - \sum\_{\boldsymbol{\kappa}\prec\boldsymbol{\kappa}'} \boldsymbol{\Phi}\_{\boldsymbol{\kappa}}\left(m\right) \mathbf{1}\_{\{\boldsymbol{\ell}\_{\boldsymbol{\kappa}\_{\boldsymbol{\kappa}}}-\boldsymbol{m}\}}\mathbf{Y}\_{\boldsymbol{\kappa}-\boldsymbol{\kappa}}\right)\right]$$

for 1 *m M* .

b.  *L L L L*

 2 1 *N N* . *L L L L*

**H***mm m* **Ψ Γ** (53)

*S S*

  (54)

(55)

(56)

(51)

(52)

1

*m m row col mm m*

*m*

**s**

matrix estimators is estimated as follows:

**s**

**r**

by applying Preposition 1 and rearranging terms, we obtain (47).

by solving the simultaneous linear equations given as follows:

**r**

By applying Preposition 1 and rearranging terms, we obtain (46).

 2

For this case, we recognize the following from (50):

For this case from (50), we recognize

a. Non-Separable Case:

where

and

**6.2 Interaction matrix parameter estimation**

 1 **μ** *m* is estimated, given that 2 **μ** *m* is known 2 **μ** *m* is estimated, given that 1 **μ** *m* is known. Thus

$$\begin{split} & \hat{\mathbf{p}}^{(1)}(m) \\ = & \left[ \sum\_{\boldsymbol{\omega} \in \mathcal{S}(m)} \left( \boldsymbol{\mu}^{(1)}(m) \otimes \mathbf{I}\_{\boldsymbol{\omega}\_{1}} \right)^{\intercal} \left( \mathbf{I}\_{\boldsymbol{\aleph}} - \sum\_{\boldsymbol{\kappa} \in \boldsymbol{\boldsymbol{\nu}}^{\boldsymbol{\prime}}} \boldsymbol{\Phi}\_{\boldsymbol{\epsilon}}(m) \mathbf{1}\_{\{\boldsymbol{\ell}\_{\boldsymbol{\zeta}\_{\boldsymbol{\omega}\_{\boldsymbol{\mathsf{m}}}} \mathbf{m} \}} \right)^{\intercal} \boldsymbol{\Sigma}^{-1}(m) \left( \mathbf{I}\_{\boldsymbol{\aleph}} - \sum\_{\boldsymbol{\kappa} \in \boldsymbol{\boldsymbol{\nu}}^{\boldsymbol{\prime}}} \boldsymbol{\Phi}\_{\boldsymbol{\epsilon}}(m) \mathbf{1}\_{\{\boldsymbol{\ell}\_{\boldsymbol{\sigma}\_{\boldsymbol{\mathsf{m}}} \mathbf{m} \}} \right) \left( \mathbf{p}^{(2)}(m) \otimes \mathbf{I}\_{\boldsymbol{\aleph}\_{\boldsymbol{\mathsf{m}}}} \right)^{\intercal} \\ = & \left[ \sum\_{\boldsymbol{\omega} \in \mathcal{S}(m)} \left( \boldsymbol{\mu}^{(2)}(m) \otimes \mathbf{I}\_{\boldsymbol{\aleph}\_{\boldsymbol{\mathsf{m}}}} \right)^{\intercal} \left( \mathbf{I}\_{\boldsymbol{\aleph}} - \sum\_{\boldsymbol{\kappa} \in \boldsymbol{\boldsymbol{\nu}}^{\boldsymbol{\prime}}} \boldsymbol{\Phi}\_{\boldsymbol{\epsilon}}(m) \mathbf{1}\_{\{\boldsymbol{\ell}\_{\boldsymbol{\sigma}\_{\boldsymbol{\mathsf{m}}} \mathbf{m}} \}} \right)^{\intercal} \boldsymbol{\Sigma}^{-1}(m) \left( \mathbf{I}\_{\boldsymbol{\aleph}\_{\boldsymbol{\mathsf{m}}}} - \sum\_{\boldsymbol{\kappa} \in \boldsymbol{\boldsymbol{\$$

for 1 *m M* .

*Proof*

$$\mathbf{X}\_{\mathbf{s}} = \left(\mathbf{Y}\_{\mathbf{s}} - \sum\_{\mathbf{r} \in \mathcal{N}} \boldsymbol{\Theta}\_{\mathbf{r}} \left(\boldsymbol{L}\_{\mathbf{s}}\right) \mathbf{1}\_{\{\boldsymbol{L}\_{\mathbf{s}\leftarrow\mathbf{r}} = \boldsymbol{L}\_{\mathbf{s}}\}} \mathbf{Y}\_{\mathbf{s}-\mathbf{r}}\right) - \left(\mathbf{I}\_{N} - \sum\_{\mathbf{r} \in \mathcal{N}} \boldsymbol{\Theta}\_{\mathbf{r}} \left(\boldsymbol{L}\_{\mathbf{s}}\right) \mathbf{1}\_{\{\boldsymbol{L}\_{\mathbf{s}\leftarrow\mathbf{r}} = \boldsymbol{L}\_{\mathbf{s}}\}}\right) \boldsymbol{\upmu}\left(\boldsymbol{L}\_{\mathbf{s}}\right). \tag{48}$$

For the separable case, the mean can be written as follows:

$$\begin{split} \mathfrak{p}\left(m\right) &= \mathfrak{p}^{\left(2\right)}\left(m\right) \otimes \mathfrak{p}^{\left(1\right)}\left(m\right) \\ &= \left(\mathfrak{p}^{\left(2\right)}\left(m\right) \otimes \mathbf{I}\_{N\_1}\right)\left(\mathbf{1} \otimes \mathfrak{p}^{\left(1\right)}\left(m\right)\right) = \left(\mathfrak{p}^{\left(2\right)}\left(m\right) \otimes \mathbf{I}\_{N\_1}\right)\mathfrak{p}^{\left(1\right)}\left(m\right) \\ &= \left(\mathbf{I}\_{N\_2} \otimes \mathfrak{p}^{\left(1\right)}\left(m\right)\right)\left(\mathfrak{p}^{\left(2\right)}\left(m\right) \otimes \mathbf{1}\right) = \left(\mathbf{I}\_{N\_2} \otimes \mathfrak{p}^{\left(1\right)}\left(m\right)\right)\mathfrak{p}^{\left(2\right)}\left(m\right). \end{split} \tag{49}$$

Plugging the results of (49) into (48) yields

$$\begin{split} \mathbf{X\_{s}} &= \left(\mathbf{Y\_{s}} - \sum\_{\mathbf{re},\boldsymbol{\mathcal{V}}} \boldsymbol{\Theta}\_{\mathbf{r}}(L\_{\mathbf{s}}) \mathbf{1}\_{\{l\_{\star,\boldsymbol{\bullet}}=l\_{\star}\}} \mathbf{Y\_{s-\mathbf{r}}}\right) - \left(\mathbf{I\_{N}} - \sum\_{\mathbf{re},\boldsymbol{\mathcal{V}}} \boldsymbol{\Theta}\_{\mathbf{r}}(L\_{\mathbf{s}}) \mathbf{1}\_{\{l\_{\star,\boldsymbol{\bullet}}=l\_{\star}\}}\right) \left(\mathbf{\mu}^{(2)}(L\_{\mathbf{s}}) \otimes \mathbf{I}\_{N\_{\mathbf{r}}}\right) \mathbf{\mu}^{(1)}(L\_{\mathbf{s}}) \\ &= \left(\mathbf{Y\_{s}} - \sum\_{\mathbf{re},\boldsymbol{\mathcal{V}}} \boldsymbol{\Theta}\_{\mathbf{r}}(L\_{\mathbf{s}}) \mathbf{1}\_{\{l\_{\star,\boldsymbol{\bullet}}=l\_{\star}\}} \mathbf{Y\_{s-\mathbf{r}}}\right) - \left(\mathbf{I\_{N}} - \sum\_{\mathbf{re},\boldsymbol{\mathcal{V}}} \boldsymbol{\Theta}\_{\mathbf{r}}(L\_{\mathbf{s}}) \mathbf{1}\_{\{l\_{\star,\boldsymbol{\bullet}}=l\_{\star}\}}\right) \left(\mathbf{I\_{N\_{\mathbf{s}}}} \otimes \mathbf{\mu}^{(1)}(L\_{\mathbf{s}})\right) \mathbf{\mu}^{(2)}(L\_{\mathbf{s}}). \end{split} \tag{50}$$

$$\left(\mathbf{1}^{\cdot}\right)\_{\cdot} \qquad \Phi\left(m\right) = \mu^{\left(\cdot\right)}\left(m\right), \ 1 \le m \le M$$

For this case, we recognize the following from (50):

$$\mathbf{Q\_s} = \left(\mathbf{I\_N} - \sum\_{\mathbf{r} \neq \boldsymbol{\mathcal{N}}} \mathbf{\theta\_r}(L\_\mathbf{s}) \mathbf{1}\_{\{L\_\mathbf{s} = \boldsymbol{I}\_\mathbf{s}\}}\right) \left(\mathbf{p}^{(2)}(L\_\mathbf{s}) \otimes \mathbf{I}\_{N\_1}\right). \tag{51}$$

By applying Preposition 1 and rearranging terms, we obtain (46).

$$\left(\begin{array}{c} \begin{pmatrix} \circ \\ \end{pmatrix} \right) \qquad \Phi\left(\begin{pmatrix} m \\ \end{pmatrix}\right) = \mathfrak{h}^{\left(\begin{smallmatrix} \\ \end{smallmatrix}\right)}\left(\begin{smallmatrix} m \\ \end{smallmatrix}\right) , \; 1 \le m \le M \end{smallmatrix}$$

For this case from (50), we recognize

$$\mathbf{Q\_s} = \left(\mathbf{I\_N} - \sum\_{\mathbf{r} \neq \angle \forall} \mathbf{\hat{e}\_r}(L\_\mathbf{s}) \mathbf{1}\_{\{L\_\mathbf{s-r} = L\_\mathbf{s}\}}\right) \left(\mathbf{I\_{N\_2}} \otimes \mathbf{\mu^{(1)}}(L\_\mathbf{s})\right). \tag{52}$$

by applying Preposition 1 and rearranging terms, we obtain (47).

#### **6.2 Interaction matrix parameter estimation**

**Proposition 3** Assume that the mean vectors **μ** *m* for 1 *m M* and the covariance matrices **Σ** *m* for 1 *m M* are known, then interaction matrix parameters are estimated by solving the simultaneous linear equations given as follows:

a.    From (3), **Xs** can be written as follows:

Plugging the results of (49) into (48) yields

For the separable case, the mean can be written as follows:

*mmm*

**r r**

**r r**

**μμ μ**

2 1

**sr r**

*T T*

*T T*

**sr r**

1

$$\mathbf{H}(m) = row \left\{ col \left( \sum\_{\mathbf{s} \in \mathcal{S}(m)} \mathbf{A}\_{\mathbf{s},\mathbf{t}}(m) \boldsymbol{\Sigma}^{-1}(m) \mathbf{A}\_{\mathbf{s},\mathbf{t}}^{T}(m), \mathbf{r} \in \mathcal{N}\_{S} \right), \mathbf{t} \in \mathcal{N}\_{S} \right\} \tag{54}$$

$$\Gamma(m) = row \left( \sum\_{s \in \mathcal{S}(m)} \mathbf{A}\_{\mathbf{s},\mathbf{t}}(m) \boldsymbol{\Sigma}^{-1}(m) \big( \mathbf{Y}\_{\mathbf{s}} - \boldsymbol{\mathfrak{u}}(m) \big), \ t \in \mathcal{N}\_S \right) \tag{55}$$

$$\mathbf{WP}\left(m\right) = row\left(vec\left(\hat{\boldsymbol{\theta}}\_{\mathbf{r}}\left(m\right)\right), \ \mathbf{r} \in \mathcal{N}\_S\right) \tag{56}$$

and

$$\mathbf{A}\_{\mathbf{s},\mathbf{r}}(m) = \left( \left( \mathbf{Y}\_{\mathbf{s}\cdot\mathbf{r}} - \boldsymbol{\mathfrak{p}}\left(m\right) \right) \mathbf{1}\_{\{l\_{\mathbf{s}\cdot\mathbf{r}}=m\}} \otimes \mathbf{I}\_{N} \right) + \mathbf{K}\_{N,N} \left( \boldsymbol{\Sigma}^{-1}(m) \otimes \boldsymbol{\mathfrak{s}}(m) \right) \left( \left( \mathbf{Y}\_{\mathbf{s}+\mathbf{r}} - \boldsymbol{\mathfrak{p}}\left(m\right) \right) \mathbf{1}\_{\{l\_{\mathbf{s}\cdot\mathbf{r}}=m\}} \otimes \mathbf{I}\_{N} \right). \tag{57}$$

From the invariance property of the MPL, the complete set of non-separable interaction matrix estimators is estimated as follows:

Estimation of the Separable MGMRF Parameters for Thematic Classification 113

a. The proof for the non-separable case is derived by applying Proposition 1 (Navarro, et.

*T vec L L vec L* **s s s s sr s r s r**

The above expression can also be written using the following matrix identities (Magnus and

, *T*

where *m n* **A** .

(69)

(74)

(75)

(76)

for *S* **r**

al., 2009).

Neudecker, 1999)

where *m n*

Furthermore, since

then,

*Proof*

and 1 *m M* .

b. From (3), **Xs** can be written as

where *m n* **A** , *np* **B** , and . *p q* **C**

**A** and *p q*

$$\hat{\boldsymbol{\Theta}}\_{\mathbf{r}}(m) = \text{reshape}\left(\text{vec}\left(\hat{\boldsymbol{\Theta}}\_{\mathbf{r}}(m)\right), \boldsymbol{N}, \boldsymbol{N}\right) \tag{58}$$

$$
\hat{\boldsymbol{\Theta}}\_{-\mathbf{r}}(m) = \boldsymbol{\Sigma}(m)\hat{\boldsymbol{\Theta}}\_{\mathbf{r}}^T(m)\left(\boldsymbol{\Sigma}(m)\right)^{-1} \tag{59}
$$

for *S* **r** , 1 *m M* .

#### b. Separable Case:

In addition, if we assume the following for *S* **r** and 1 *m M* :

then

$$\mathbf{H}^{(k)}(m)\mathbf{W}^{(k)}(m) = \mathbf{I}^{(k)}(m) \tag{60}$$

where

$$\mathbf{H}^{(k)}(m) = row \left\{ col \left( \sum\_{\mathbf{s} \in \mathcal{S}(m)} \mathbf{A}\_{\mathbf{s},\mathbf{t}}^{(k)}(m) \Sigma^{-1}(m) \mathbf{A}\_{\mathbf{s},\mathbf{t}}^{(k)T}(m), \ \mathbf{r} \in \mathcal{N}\_S \right), \ \mathbf{t} \in \mathcal{N}\_S \right\} \tag{61}$$

$$\Gamma^{(k)}(m) = row \left( \sum\_{\mathbf{s} \in \mathcal{S}(m)} \mathbf{A}\_{\mathbf{s},\mathbf{t}}^{(k)}(m) \Sigma^{-1}(m) (\mathbf{Y}\_{\mathbf{s}} - \mathfrak{u}(m))\_{\mathbf{t}} \ \mathbf{t} \in \mathcal{N}\_S \right) \tag{62}$$

$$\mathbf{W}^{(k)}(m) = \text{row}\left(\text{vec}\left(\hat{\mathsf{H}}\_{\mathbf{r}}^{(k)}(m)\right), \ \mathbf{r} \in \mathcal{N}\_S\right) \tag{63}$$

for 1 2 *k* and

$$\mathbf{A}\_{\mathbf{s},\mathbf{r}}^{(1)}(m) = \left(\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(2)}(m)\right) \otimes \mathbf{I}\_{N\_1}\right)^{T} \left(\mathbf{I}\_{N\_2} \otimes \mathbf{K}\_{N\_1,N\_2} \otimes \mathbf{I}\_{N\_1}\right)^{T} \mathbf{A}\_{\mathbf{s},\mathbf{r}}(m) \tag{64}$$

$$\mathbf{A}\_{\mathbf{s},\mathbf{r}}^{(2)}(m) = \left(\mathbf{I}\_{N\_2} \otimes vec\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(1)}(m)\right)\right)^{T} \left(\mathbf{I}\_{N\_2} \otimes \mathbf{K}\_{N\_1,N\_2} \otimes \mathbf{I}\_{N\_1}\right)^{T} \mathbf{A}\_{\mathbf{s},\mathbf{r}}(m). \tag{65}$$

From the invariance property of the MPL, the complete set of separable interaction matrix estimators is estimated as follows for *S* **r** , 1 *m M* , 1 2 *k* :

$$\hat{\boldsymbol{\Theta}}\_{\mathbf{r}}^{(k)}(m) = reshape\left(vec\left(\hat{\boldsymbol{\Theta}}\_{\mathbf{r}}^{(k)}\left(m\right)\right), \boldsymbol{N}\_{k'} \boldsymbol{N}\_{k}\right) \tag{66}$$

$$
\hat{\boldsymbol{\Theta}}\_{\mathbf{r}}^{(k)}(m) = \boldsymbol{\Sigma}^{(k)}(m)\hat{\boldsymbol{\Theta}}\_{\mathbf{r}}^{(k)T}\left(m\right)\left(\boldsymbol{\Sigma}^{(k)}(m)\right)^{-1} \tag{67}
$$

and also

$$
\hat{\boldsymbol{\Theta}}\_{\mathbf{r}}(m) = \hat{\boldsymbol{\Theta}}\_{\mathbf{r}}^{(2)}(m) \otimes \hat{\boldsymbol{\Theta}}\_{\mathbf{r}}^{(1)}(m) \tag{68}
$$

for *S* **r** and 1 *m M* .

*Proof*

112 Remote Sensing – Advanced Techniques and Platforms

1 , , , , *k k k T*

**s t s r**

**H A Σ Art**

, , *k k*

**s t s**

ˆ , *k k* **Ψ θ** *m row vec m* **r** *S* **r**

1 2 12 1

2 2 1 2 1

From the invariance property of the MPL, the complete set of separable interaction matrix

ˆ ˆ , , *k k*

 1 ˆ ˆ *k kk T k m mmm*

, , , . *T T*

, , ,

*m row mm m*

**Γ A Σ Y μ t**

1

 2 **θr** *m* is known

 1 **θr** *m* is known

*m m row col mm m*

*m*

**s**

1 2

estimators is estimated as follows for *S* **r**

2 1

**s**

for *S* **r**

then

where

b. Separable Case:

1

2

for 1 2 *k* and

and also

, 1 *m M* .

 In addition, from the identity (Magnus and Neudecker, 1999)

$$vec{ec}\left(\mathbf{A}\otimes\mathbf{B}\right) = \left(\mathbf{I}\_n\otimes\mathbf{K}\_{q,m}\otimes\mathbf{I}\_p\right)\cdot\left(vec\left(\mathbf{A}\right)\otimes vec\left(\mathbf{B}\right)\right)\tag{73}$$

where *m n* **A** and *p q* **B** , it follows that

$$\begin{split} \text{vec}\left(\boldsymbol{\Phi}\_{\mathbf{r}}\left(m\right)\right) &= \text{vec}\left(\boldsymbol{\Phi}\_{\mathbf{r}}^{(2)}\left(m\right) \otimes \boldsymbol{\Phi}\_{\mathbf{r}}^{(1)}\left(m\right)\right) \\ &= \left(\mathbf{I}\_{N\_{2}} \otimes \mathbf{K}\_{N\_{1}, N\_{2}} \otimes \mathbf{I}\_{N\_{1}}\right) \cdot \left(\text{vec}\left(\boldsymbol{\Phi}\_{\mathbf{r}}^{(2)}\left(m\right)\right) \otimes \text{vec}\left(\boldsymbol{\Phi}\_{\mathbf{r}}^{(1)}\left(m\right)\right)\right). \end{split} \tag{74}$$

Furthermore, since

$$\begin{aligned} &\left(\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(2)}(m)\right)\otimes\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(1)}(m)\right)\right) \\ &= \left(\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(2)}(m)\right)\otimes\mathbf{I}\_{N\_{1}}\right)\left(\mathbf{1}\otimes\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(1)}(m)\right)\right) = \left(\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(2)}(m)\right)\otimes\mathbf{I}\_{N\_{1}}\right)\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(1)}(m)\right) \\ &= \left(\mathbf{I}\_{N\_{2}}\otimes\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(1)}(m)\right)\right)\left(\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(2)}(m)\right)\otimes\mathbf{1}\right) = \left(\mathbf{I}\_{N\_{2}}\otimes\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(1)}(m)\right)\right)\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(2)}(m)\right) \end{aligned} \tag{75}$$

then,

$$\begin{split} \text{vec}\left(\boldsymbol{\Phi\_{\mathbf{r}}}\left(m\right)\right) &= \left(\mathbf{I}\_{N\_{2}} \otimes \mathbf{K}\_{N\_{1},N\_{2}} \otimes \mathbf{I}\_{N\_{1}}\right) \cdot \left(\text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(2)}\left(m\right)\right) \otimes \mathbf{I}\_{N\_{1}}\right) \text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(1)}\left(m\right)\right) \\ &= \left(\mathbf{I}\_{N\_{2}} \otimes \mathbf{K}\_{N\_{1},N\_{2}} \otimes \mathbf{I}\_{N\_{1}}\right) \cdot \left(\mathbf{I}\_{N\_{2}} \otimes \text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(1)}\left(m\right)\right)\right) \text{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(2)}\left(m\right)\right). \end{split} \tag{76}$$

Estimation of the Separable MGMRF Parameters for Thematic Classification 115

**Σ** *m* is known

**Σ** *m* is known

 1 1 2 # #

 1 2 1 # #

The above estimators are not in their closed form. The estimators can be solved iteratively

The multispectral and multitemporal satellite image under consideration is the 'Butuan' image acquired from the LANDSAT TM. The image shows the scenery of Butuan City and its surroundings in Northeastern Mindanao, Philippines. It consists of six spectral bands and four temporal slots with a dynamic range of 8 bits. The images were captured chronologically on the following dates: August 1, 1992, August 7, 2000, May 22, 2001, and December 3, 2002. The images were radiometrically corrected, geometrically co-registered with each other, and have been resized to 600 x 800 pixels.  Separable Case:

**7. Data preparation**

1

2

then

In addition, if we assume the following for 1 *m M* :

**Σ** *m* is estimated, given that 2

**Σ** *m* is estimated, given that 1

using the flip-flop algorithm (Dutilleul, 1999).

scaled RGB realization captured on May 22, 2001.

Fig. 1. RGB image of 'Butuan' captured on May 22, 2001.

Plugging the results of (76) into (69) yields

 1 1 , 2 2 , *S S T T L L vec L L L vec L* **ss s sr s r s r s s sr s r s r X Y μ A θ Y μ A θ** (77) 1 1 *m vec m* **Φ θt** , , 1 **t** *S m M*

For this case, we recognize from (77),

$$\mathbf{Q}\_{\mathbf{s}} = \mathbf{A}\_{\mathbf{s},\mathbf{t}}^{(1)T}(m) \,. \tag{78}$$

By applying Preposition 1 and rearranging terms, we obtain the following expression

$$\sum\_{\mathbf{s}\in\mathcal{N}\_{\mathbf{s}}}\sum\_{\mathbf{s}\in\mathcal{S}(m)}\mathbf{A}\_{\mathbf{s},\mathbf{t}}^{(1)}(m)\Sigma^{-1}(m)\mathbf{A}\_{\mathbf{s},\mathbf{t}}^{(1)T}(m)\operatorname{vec}\Big(\boldsymbol{\theta}\_{\mathbf{r}}^{(1)}(m)\Big) = \sum\_{\mathbf{s}\in\mathcal{S}(m)}\mathbf{A}\_{\mathbf{s},\mathbf{t}}^{(1)}(m)\Sigma^{-1}(m)\left(\mathbf{Y}\_{\mathbf{s}}-\boldsymbol{\mathfrak{y}}(m)\right). \tag{79}$$

By aggregating the equation in (79) for *S* **t** , the interaction matrix coefficients are estimated by solving the simultaneous linear equations in (60) for 1. *k*

$$\mathfrak{a}\left(2^{\cdot}\right) \qquad \mathfrak{sp}\left(m\right) = \mathsf{vec}\left(\mathfrak{G}\_{\ast}^{(\cdot)}\left(m\right)\right), \ \mathsf{t} \in \mathscr{N}\_{\mathbb{S}}, \ \mathsf{l} \le m \le M$$

For this case, we recognize from (77)

$$\mathbf{Q}\_{\mathbf{s}} = \mathbf{A}\_{\mathbf{s},\mathbf{t}}^{\{2\}^{\rm T}}(m). \tag{80}$$

By applying Preposition 1 and rearranging terms, we obtain the following expression

$$\sum\_{\mathbf{r}\in\mathcal{N}\_{\mathbf{s}}}\sum\_{\mathbf{s}\in\mathcal{S}(m)}\mathbf{A}\_{\mathbf{s},\mathbf{t}}^{(2)}(m)\Sigma^{-1}(m)\mathbf{A}\_{\mathbf{s},\mathbf{t}}^{(2)T}(m)\mathrm{vec}\left(\boldsymbol{\theta}\_{\mathbf{r}}^{(2)}(m)\right) = \sum\_{\mathbf{s}\in\mathcal{S}(m)}\mathbf{A}\_{\mathbf{s},\mathbf{t}}^{(2)}(m)\Sigma^{-1}(m)(\mathbf{Y}\_{\mathbf{s}}-\boldsymbol{\mathfrak{y}}(m)). \tag{81}$$

By aggregating the equations in (79) for *S* **t** , the interaction matrix coefficients are estimated by solving the simultaneous linear equations in (60) for 2. *k*

#### **6.3 Covariance matrix parameter estimation**

Since **Xs** is dependent on a covariance matrix in finding the MPL estimator of **Σ** *m* , for all 1 *m M* is cumbersome to derive. As an alternative, we estimate the covariance matrix as the sample covariance matrix given that the mean vectors **μ***m* for 1 *m M* and the interaction matrix coefficients *m* **r θ** , for **r** , 1 *m M* are known, then the covariance matrix parameters are estimated as follows:

a. Non-Separable Case:

$$
\hat{\Sigma}(m) = \frac{1}{r(m)} \sum\_{\mathbf{s} \in \mathcal{S}(m)} \mathbf{X}\_{\mathbf{s}} \mathbf{X}\_{\mathbf{s}}^T \tag{82}
$$

#### b. Separable Case:

In addition, if we assume the following for 1 *m M* :

then

114 Remote Sensing – Advanced Techniques and Platforms

*S*

**ss s sr s r s r**

**X Y μ A θ**

**Y μ A θ**

**r s s**

By aggregating the equation in (79) for *S* **t**

 2

For this case, we recognize from (77)

*T m m*

estimated by solving the simultaneous linear equations in (60) for 1. *k*

**r s s**

By aggregating the equations in (79) for *S* **t**

**6.3 Covariance matrix parameter estimation**

interaction matrix coefficients *m* **r θ** , for **r**

matrix parameters are estimated as follows:

*T m m*

estimated by solving the simultaneous linear equations in (60) for 2. *k*

**r**

By applying Preposition 1 and rearranging terms, we obtain the following expression

1 11 1 1 1

**s t s r r s t s**

By applying Preposition 1 and rearranging terms, we obtain the following expression

2 22 2 1 1

**s t s r r s t s**

*S*

 1 *m vec m* **Φ θt** , , 1 **t**

,

*T*

*L L vec L*

*T*

*L L vec L*

1 1

2 2

*S m M*

*T* **Q A s st** *m* . (78)

*S m M*

*T* **Q A s st** *m* .   Non-Separable Case:

*S*

*S*

$$
\hat{\boldsymbol{\Sigma}}^{(1)}(m) = \frac{1}{r(m)N\_{2}} \sum\_{\mathbf{s} \in \mathcal{S}(m)} \mathbf{X}\_{\mathbf{s}}^{\#} \left(\hat{\boldsymbol{\Sigma}}^{(2)}(m)\right)^{-1} \mathbf{X}\_{\mathbf{s}}^{\#T} \tag{83}
$$

$$
\hat{\boldsymbol{\Sigma}}^{(2)}(m) = \frac{1}{r\left(m\right)N\_{1}} \sum\_{\mathbf{s} \in \mathcal{S}\left(m\right)} \mathbf{X}\_{\mathbf{s}}^{\#\mathcal{T}} \left(\boldsymbol{\Sigma}^{(1)}(m)\right)^{-1} \mathbf{X}\_{\mathbf{s}}^{\#}.\tag{84}
$$

The above estimators are not in their closed form. The estimators can be solved iteratively using the flip-flop algorithm (Dutilleul, 1999).

## **7. Data preparation**

The multispectral and multitemporal satellite image under consideration is the 'Butuan' image acquired from the LANDSAT TM. The image shows the scenery of Butuan City and its surroundings in Northeastern Mindanao, Philippines. It consists of six spectral bands and four temporal slots with a dynamic range of 8 bits. The images were captured chronologically on the following dates: August 1, 1992, August 7, 2000, May 22, 2001, and December 3, 2002. The images were radiometrically corrected, geometrically co-registered with each other, and have been resized to 600 x 800 pixels. The image in Fig. 1 is a grayscaled RGB realization captured on May 22, 2001.

Fig. 1. RGB image of 'Butuan' captured on May 22, 2001.

Estimation of the Separable MGMRF Parameters for Thematic Classification 117

(Hazel, 2000) which in general, does not hold the multivariate case. This relation, however, holds in the univariate case (Kashyap and Chellappa, 1983) as well as the Rellier's GMRF. On the other hand, anisotropic models, such as Rellier's GMRF, and our model exhibited a substantially better classification performance as compared to the GSC. Since the covariance matrix estimators used a sub-optimal alternative, some slight performance degradation has

Denote *S***μ** , *S***θ** , and *S***Σ** to be the separable indicators for the mean, interaction matrix, and

Since the GSC model is a degenerate form of our MGMRF with zero interaction matrices, the separability structure of the mean and covariance matrices are examined. The results are presented in Table 3 showed that no improvement in the classification performance,

> *S***Σ** *S***μ** Accuracy 0 0 55.3% 0 1 54.2% 0 0 54.3% 1 1 54.1%

The hybrid separable anisotropic MGMRF shows the separability of the covariance matrix has a slight improvement in performance over a non-separable spectro-temporal observation. As discussed in Section 5.2, the hybrid separable model with separable interaction matrix, together with a non-separable matrix, were excluded in the model performance as these modes are not possible. The classification accuracy is presented in

> *S***Σ** *S***θ** *S***μ** Accuracy 0 0 0 84.3% 0 0 1 84.6%

resulted.

Table 4.

**8.2 Hybrid separable case**

covariance matrix, respectively.

**8.2.1 Hybrid separable GSC model**

regardless of separability of the parameters.

Table 3. Classification Accuracy of Hybrid Separable GSC models.

0 1 0 0 1 1

**8.2.2 Hybrid separable anisotropic GMRF model**

**θ θ rs rs** *L L* (85)

The thematic classes were established by employing the k-means algorithm (Richards and Jia, 2006). The thematic classes were identified and their mean reflectance vector form the training data are shown in Table 1.

Table 1. Average reflectances from the training data.

Training and verification sites were obtained from a random sample of 1200 sites. The firstorder neighborhood system in the MRF modeling of the thematic map and the image were used.

## **8. Discussion**

## **8.1 Non-separable case**

The classification performance of our model with non-separable MGMRF parameters, as compared to the GSC, Hazel's, and Rellier's models are presented in Table 2.

Table 2. Classification Accuracy of Different MGMRF models.

The GSC model has a low accuracy compared to the remaining MGMRF models. It substantiates that Markov dependence would yield a better accuracy to the thematic map classification than to the site independence model.

It is noticeable that Hazel's GMRF presents a relatively poor classification accuracy which is attributed to the bilateral symmetry imposed into the interaction matrices, that is,

$$\boldsymbol{\Theta}\_{\mathbf{r}}\left(\boldsymbol{L}\_{\mathbf{s}}\right) = \boldsymbol{\Theta}\_{-\mathbf{r}}\left(\boldsymbol{L}\_{\mathbf{s}}\right) \tag{85}$$

(Hazel, 2000) which in general, does not hold the multivariate case. This relation, however, holds in the univariate case (Kashyap and Chellappa, 1983) as well as the Rellier's GMRF.

On the other hand, anisotropic models, such as Rellier's GMRF, and our model exhibited a substantially better classification performance as compared to the GSC. Since the covariance matrix estimators used a sub-optimal alternative, some slight performance degradation has resulted.

## **8.2 Hybrid separable case**

116 Remote Sensing – Advanced Techniques and Platforms

The thematic classes were established by employing the k-means algorithm (Richards and Jia, 2006). The thematic classes were identified and their mean reflectance vector form the

1 Thick Vegetation 62 48 33 91 69 29

2 Sparse Vegetation 70 58 43 99 83 37

3 Built Up Areas 77 63 54 75 78 41

4 Body of Water 72 41 29 12 13 11

5 Thin Clouds 104 84 76 88 85 53

6 Thick Clouds 197 190 190 144 167 115

Training and verification sites were obtained from a random sample of 1200 sites. The firstorder neighborhood system in the MRF modeling of the thematic map and the image were

The classification performance of our model with non-separable MGMRF parameters, as

Model Accuracy

GSC 55.3%

Hazel's GMRF 45.6%

Rellier's GMRF 83.1%

Our Model 84.3%

The GSC model has a low accuracy compared to the remaining MGMRF models. It substantiates that Markov dependence would yield a better accuracy to the thematic map

It is noticeable that Hazel's GMRF presents a relatively poor classification accuracy which is

attributed to the bilateral symmetry imposed into the interaction matrices, that is,

compared to the GSC, Hazel's, and Rellier's models are presented in Table 2.

Table 2. Classification Accuracy of Different MGMRF models.

classification than to the site independence model.

*1 2 3 4 5 7*

M Thematic Class Landsat TM Band Number

training data are shown in Table 1.

Table 1. Average reflectances from the training data.

used.

**8. Discussion**

**8.1 Non-separable case**

Denote *S***μ** , *S***θ** , and *S***Σ** to be the separable indicators for the mean, interaction matrix, and covariance matrix, respectively.

## **8.2.1 Hybrid separable GSC model**

Since the GSC model is a degenerate form of our MGMRF with zero interaction matrices, the separability structure of the mean and covariance matrices are examined. The results are presented in Table 3 showed that no improvement in the classification performance, regardless of separability of the parameters.

Table 3. Classification Accuracy of Hybrid Separable GSC models.

## **8.2.2 Hybrid separable anisotropic GMRF model**

The hybrid separable anisotropic MGMRF shows the separability of the covariance matrix has a slight improvement in performance over a non-separable spectro-temporal observation. As discussed in Section 5.2, the hybrid separable model with separable interaction matrix, together with a non-separable matrix, were excluded in the model performance as these modes are not possible. The classification accuracy is presented in Table 4.

Estimation of the Separable MGMRF Parameters for Thematic Classification 119

Fig. 3. Thematic Map – GSC with separable mean and covariance matrix

Fig. 2. Thematic Map – Hazel's MGMRF

Table 4. Classification Accuracy of Hybrid Separable Anisotropic MGMRF models.

## **8.3 Thematic maps**

Some of the thematic map labels are presented in Figs. 2 to 4, based on the May 22, 2001 satellite image. For clarity of visual presentation, thematic map labels were based on the gray-scaled average RGB reflectance of the training data.

#### Fig. 2. Thematic Map – Hazel's MGMRF

118 Remote Sensing – Advanced Techniques and Platforms

*S***Σ** *S***θ** *S***μ** Accuracy 1 0 0 84.5% 1 0 1 86.6% 1 1 0 83.8% 1 1 1 86.2%

Table 4. Classification Accuracy of Hybrid Separable Anisotropic MGMRF models.

gray-scaled average RGB reflectance of the training data.

Some of the thematic map labels are presented in Figs. 2 to 4, based on the May 22, 2001 satellite image. For clarity of visual presentation, thematic map labels were based on the

**8.3 Thematic maps**

Fig. 3. Thematic Map – GSC with separable mean and covariance matrix

Estimation of the Separable MGMRF Parameters for Thematic Classification 121

Casella, G. & Berger, R. L. (2002). *Statistical Inference 2nd ed.*, Wadswoth Group, ISBN 978-053-

Dutilleul, P. (1999). The MLE Algorithm for the Matrix Normal Distribution. *Journal of Statistical Computation and Simulation*, Vol. 64, No. 2 , ISSN 0094-9655 Fuentes, M. (2006). Testing for Separability of Spatial-Temporal Covariance Functions. *Journal of Statistical Planning and Inference*, Vol. 136, pp. 447-466, ISSN 0378-3758 Geman, S. & Graffigne, C. (1987). Markov Random Field Models and Their Applications to

Hazel, G. G. (2000). Multivariate Gaussian MRF for Multispectral Scene Segmentation and

Huizenga, H., Munck, J., Waldorp, R., Grasman, R. (2002). Spatiotemporal EEG/MEG

Kashyap, R. and Chellappa, R. (1983). Estimation and Choice of Neighbors in Spatial-

Kreyzig, E. (2005). *Advanced Engineering Mathematics, 8th. ed.*, Wiley, ISBN<PHONE_NUMBER>-85-9,

Kyriakidis, P. C. & Journel, A. G. (1999). Geostatistical Space-Time Models: A Review. *Mathematical Geology*, Vol. 31, No. 6, (August 1999), pp. 651-684, ISSN 0882-8121 Li, S. Z. (1995) *Markov Random Field Modeling in Computer Vision*, Springer-Verlag, ISBN 978-

Lu, N. & Zimmerman, D. (2005). The Likelihood Ratio Test for a Separable Covariance

Magnus, J. R. & H. Neudecker (1999). *Matrix Differential Calculus with Applications in Statistics*

Moura, J. M. F. & Balram N. (1993). Chapter 15: Statistical Algorithms for Noncausal Markov

Naik, D. N. & Rao, S. S. (2001). Analysis of Multivariate Repeated Measures Data with a

Navarro, R. D. Jr., Magadia, J. C., & Paringit, E. C. (2009). Estimating the Gauss-Markov

Neudecker, H. (1969). Some Theorems on Matrix Differentiation with Special Reference to

*and Econometrics 2nd ed.*, Wiley, ISBN<PHONE_NUMBER>-33-1, Chichester

623-691, North Holland, ISBN<PHONE_NUMBER>-05-8, Amsterdam

28, No. 1, (January 2001), pp. 91-105, ISSN 0013-1644;

No. 327, (September 1969), pp. 953-963, ISSN 0162-1459

Matrix. *Statistics and Probability Letters*, Vol. 73, No. 4, (July 2005), pp. 449-457, ISSN

Random Fields, In: *Handbook of Statistics Volume 10*, Bose, N. K. & Rao, C. R., pp.

Kronecker Product Structured Covariance Matrix. *Journal of Applied Statistics*, Vol.

Random Field Parameters for Remote Sensing Image Textures, Proceedings of TENCON 2009 - 2009 IEEE Region 10 Conference, ISBN<PHONE_NUMBER>-46-2,

Kronecker Matrix Products. Journal of American Statistical Association, Vol. 64,

ISBN<PHONE_NUMBER>-10-9,Berkeley, CA, August, 1986

3, (May 2000), pp. 1199–1211, ISSN 0196-2892

1, (January 1983), pp. 60-72, ISSN 0018-9448

Computer Vision, Proceedings of the International Congress of Mathematicians,

Anomaly Detection. *IEEE Transactions on Geoscience and Remote Sensing*, Vol. 38, No.

Source Analysis Based on a Parametric Noise Covariance Model. *IEEE Transactions on Biomedical Engineering*, Vol. 49, No. 6, (June 2002), pp. 533-539, ISSN 0018-9294 Jeng, F. and Woods, J. (1991). Compound Gauss-Markov Random Fields for Image

Estimation. *IEEE Transactions on Signal Processing*, Vol. 39, No. 3, (March 1991), pp.

Interaction Models of Images. *IEEE Transactions on Information Theory*, Vol. 29, No.

4243-12-8, Pacific Grove, CA

683–697, ISSN 1053-587X

4431701453, New York

Singapore, November, 2009

New York

0167-7152

Fig. 4. Thematic Map – Anisotropic MGMRF separable mean, interaction matrices, and covaraince matrices

## **9. Summary, conclusions, and recommendations**

This study presents a parameter estimation procedure based on the MPL for an anisotropic MGMRF with hybrid-separable parameters. Although the MGMRF is a natural extension of its univariate counterpart, the interaction matrix relationship is, in general, dependent on the covariance matrix. In an effort to make the estimation and classification procedure more tractable to compute, some sub-optimal approximations were incorporated. This resulted in a slight degradation in the classification performance. The classification performance based on our model performed well when compared to the GSC model and Hazel's MGMRF. Nonetheless, its performance is comparable to the Rellier's MGMRF. Moreover, for spectrotemporal observations, the separability of the interaction matrix as well as the covariance matrix improved the classification performance. Computational capabilities are foreseen to further advance in the near future following the improvement of numerical estimation and classification procedures.

This study presents a parameter estimation procedure based on the MPL for anisotropic MGMRF with hybrid-separable parameters. Although the MGMRF is a natural extension of its univariate counterpart, the interaction matrix relationship is, in general, dependent on the covariance matrix. In an effort to make the estimation and classification procedure more tractable to compute, some sub-optimal approximations were incorporated in the process. This resulted in a slight degradation in the classification performance. The classification performance, based on our model, has performed well, as compared to the GSC model and Hazel's MGMRF. Furthermore, its performance is comparable to Rellier's MGMRF. In terms of spectro-temporal observations, the separability of the covariance matrix has improved the classification performance. This study can be improved even more with numerical estimation and classification procedure as computational capabilities. This is foreseen to further advance in the near future.

## **10. Acknowledgment**

We acknowledge the invaluable support of extended by the Statistical Training and Research Center of the Philippine Statistical System.

## **11. References**

This study presents a parameter estimation procedure based on the MPL for an anisotropic MGMRF with hybrid-separable parameters. Although the MGMRF is a natural extension of its univariate counterpart, the interaction matrix relationship is, in general, dependent on the covariance matrix. In an effort to make the estimation and classification procedure more tractable to compute, some sub-optimal approximations were incorporated. This resulted in a slight degradation in the classification performance. The classification performance based on our model performed well when compared to the GSC model and Hazel's MGMRF. Nonetheless, its performance is comparable to the Rellier's MGMRF. Moreover, for spectrotemporal observations, the separability of the interaction matrix as well as the covariance matrix improved the classification performance. Computational capabilities are foreseen to further advance in the near future following the improvement of numerical estimation and

This study presents a parameter estimation procedure based on the MPL for anisotropic MGMRF with hybrid-separable parameters. Although the MGMRF is a natural extension of its univariate counterpart, the interaction matrix relationship is, in general, dependent on the covariance matrix. In an effort to make the estimation and classification procedure more tractable to compute, some sub-optimal approximations were incorporated in the process. This resulted in a slight degradation in the classification performance. The classification performance, based on our model, has performed well, as compared to the GSC model and Hazel's MGMRF. Furthermore, its performance is comparable to Rellier's MGMRF. In terms of spectro-temporal observations, the separability of the covariance matrix has improved the classification performance. This study can be improved even more with numerical estimation and classification procedure as computational capabilities. This is foreseen to

We acknowledge the invaluable support of extended by the Statistical Training and

Aarts, E. and Korts, J. (1987). *Simulated Annealing and Boltzmann Machines*, Wiley, ISBN 978-

Arnold, S. F. (1981). *The Theory of Linear Models and Multivariate Analysis*, Wiley, ISBN 978-

Besag, J. (1986). On the Statistical Analysis of Dirty Pictures (with discussions). *Journal of Royal Statistical Society B*. Vol. 48, No. 3, pp. 259-302, ISSN 0035-9246 Campbell, N. A. and Kiiveri, H. T. (1988). Spectral-Temporal Indices for Discrimination.

*Applied Statistics*, Vol. 37, No. 1, pp. 51-62, ISSN 0035-9254

Fig. 4. Thematic Map – Anisotropic MGMRF separable mean, interaction matrices, and

**9. Summary, conclusions, and recommendations**

covaraince matrices

classification procedures.

further advance in the near future.

Research Center of the Philippine Statistical System.

047-1921-46-2, New York

047-1050-65-0, New York

**10. Acknowledgment**

**11. References**

**0**

**6**

*Spain*

**Low Rate High Frequency Data Transmission**

This chapter deals with the difficulties of transmitting data gathered from sensors placed in very remote areas where energy supplies are scarce. The data link is established by means of the ionosphere, a layer of the upper atmosphere that is ionized by solar radiation. Communications through the ionosphere have persisted, although the use of artificial repeaters, such as satellites, has provided more reliable communication. In spite of being random, noisy and susceptible to interference, ionospheric transmission still has favorable characteristics (e.g. low cost equipment, worldwide coverage, invulnerability, etc.) that

The Research Group in Electromagnetism and Communications (GRECO) from La Salle - Universitat Ramon Llull (Spain) is investigating techniques for the improvement of remote sensing and skywave digital communications. The GRECO has focused its attention on the link between Antarctica and Spain. The main objectives of this study are: to implement a long-haul oblique ionospheric sounder and to transmit data from sensors located at the

The SAS is located on Livingston Island (62.7 ◦S, 299.6 ◦E; geomagnetic latitude 52.6 ◦S) in the South Shetlands archipelago. Spanish research is focused on the study of the biological and geological environment, and also the physical geography. Many of the research activities undertaken at the SAS collect data on temperature, position, magnetic field, height, etc. which is temporarily stored in data loggers on-site. Part of this data is then transmitted to research laboratories in Spain. Even though the SAS is only manned during the austral summer, data collection never stops. While the station is left unmanned, the sets of data are stored in memory devices, and are not downloaded until the next Antarctic season. The information that has to be analyzed in almost real-time is transmitted to Spain through a satellite link. The skywave digital communication system, presented here, is intended to transmit the information from the Antarctic sensors as a backup, or even as an alternative

Antarctica is a continent of great scientific interest in terms of remote sensing experiments related to physics and geology. Due to the peculiarities of Antarctica, some of these experiments cannot be conducted anywhere else on the Earth and this fact might oblige the

to the satellite, without depending on other entities for support or funding.

**1. Introduction**

appeal to current communications engineering.

Spanish Antarctic Station (SAS) Juan Carlos I to Spain.

**from Very Remote Sensors**

and Joan Ramon Regué *La Salle - Universitat Ramon Llull*

Pau Bergada, Rosa Ma Alsina-Pages, Carles Vilella

## **Low Rate High Frequency Data Transmission from Very Remote Sensors**

Pau Bergada, Rosa Ma Alsina-Pages, Carles Vilella and Joan Ramon Regué *La Salle - Universitat Ramon Llull Spain*

#### **1. Introduction**

122 Remote Sensing – Advanced Techniques and Platforms

Ravishanker, N. & Dey, D. K. (2002). *A First Course in Linear Model Theory*, CRC Press LLC,

Richards, J. A. & Jia, X. (1999). *Remote Sensing Image Analysis: An Introduction, 4th ed*.,

Rellier, G., Descombes, X., Falzon, F., & Zerubia, J. (2004). Texture Feature Analysis Using a

Winkler, G. (2006). *Image Analysis, Random Fields and Dynamic Monte Carlo Methods: A Mathematical Introduction 2nd ed.*, Springer-Verlag. ISBN<PHONE_NUMBER>-13-4, Berlin

Gauss–Markov Model in Hyperspectral Image Segmentation. *IEEE Transactions on Geoscience and Remote Sensing*, Vol. 42, No. 7, (July 2004), pp. 1543-1551, ISSN 0196-

ISBN<PHONE_NUMBER>-47-3, Boca Raton, FL

2892

Springer-Verlag. ISBN<PHONE_NUMBER>-28-6, Berlin

This chapter deals with the difficulties of transmitting data gathered from sensors placed in very remote areas where energy supplies are scarce. The data link is established by means of the ionosphere, a layer of the upper atmosphere that is ionized by solar radiation. Communications through the ionosphere have persisted, although the use of artificial repeaters, such as satellites, has provided more reliable communication. In spite of being random, noisy and susceptible to interference, ionospheric transmission still has favorable characteristics (e.g. low cost equipment, worldwide coverage, invulnerability, etc.) that appeal to current communications engineering.

The Research Group in Electromagnetism and Communications (GRECO) from La Salle - Universitat Ramon Llull (Spain) is investigating techniques for the improvement of remote sensing and skywave digital communications. The GRECO has focused its attention on the link between Antarctica and Spain. The main objectives of this study are: to implement a long-haul oblique ionospheric sounder and to transmit data from sensors located at the Spanish Antarctic Station (SAS) Juan Carlos I to Spain.

The SAS is located on Livingston Island (62.7 ◦S, 299.6 ◦E; geomagnetic latitude 52.6 ◦S) in the South Shetlands archipelago. Spanish research is focused on the study of the biological and geological environment, and also the physical geography. Many of the research activities undertaken at the SAS collect data on temperature, position, magnetic field, height, etc. which is temporarily stored in data loggers on-site. Part of this data is then transmitted to research laboratories in Spain. Even though the SAS is only manned during the austral summer, data collection never stops. While the station is left unmanned, the sets of data are stored in memory devices, and are not downloaded until the next Antarctic season. The information that has to be analyzed in almost real-time is transmitted to Spain through a satellite link. The skywave digital communication system, presented here, is intended to transmit the information from the Antarctic sensors as a backup, or even as an alternative to the satellite, without depending on other entities for support or funding.

Antarctica is a continent of great scientific interest in terms of remote sensing experiments related to physics and geology. Due to the peculiarities of Antarctica, some of these experiments cannot be conducted anywhere else on the Earth and this fact might oblige the

for data transmission from the SAS to Spain. Data provided by the VIS is used to conduct ionospheric research, mainly to characterize the climatology of the ionospheric characteristics and to investigate the ionospheric effects caused during geomagnetically disturbed periods

Low Rate High Frequency Data Transmission from Very Remote Sensors 125

The oblique ionosonde monitors various parameters to model the HF radiolink between the SAS and Spain (Vilella et al., 2008). These parameters include link availability, power delay profile and frequency dispersion of the channel. The sounder includes a transmitter, placed on the premises of the SAS and a receiver deployed in Spain. The main drawback of the oblique sounder is the difficulty in establishing the ionospheric link. Firstly, the long distance of the link (12700 km) requires four hops to reach the receiver. And secondly, the transmitted signal

The ionosphere study can be approximated from several points of view. Vertical incidence soundings provide accurate information about electron density profiles below the peak electron density. However, when using this technique the electron profile must be extrapolated from the peak point to the upper limit of the ionosphere. Moreover, the low density of vertical ionosondes, especially in oceans and remote areas, is a serious impairment. GNSS receivers constitute a high temporal and spatial resolution sounding network which despite gaps over oceans and remote regions, can be used to study fast perturbations affecting local regions, such as Travelling Ionospheric Disturbances and scintillations, or wider regions such as solar flares. Data gathered from GNSS receivers can provide information about the Total Electron Content (TEC) between receivers and satellites by means of proper tomographic modeling approaches. Spatial and temporal variations of the main ionospheric events can be monitored by means of GNSS receivers, especially those placed in the Antarctic Region, which is considered the entrance point of many ionospheric disturbances coming from Solar events.

This chapter will study, analyze and experimentally verify a possible candidate for the physical layer of a long-haul ionospheric data link, focusing on the case SAS-Spain. Preliminary studies of data transmission feasibility over this link were already performed

The first application of this link is the transmission of data generated by a geomagnetic sensor installed at the SAS. Future applications may include sending information of another nature

The minimum requirements regarding the geomagnetic sensor data transmissions from the

Furthermore, TEC reaches its highest variability peaks in the Antarctica area.

in (Deumal et al., 2006) and (Bergada et al., 2009), with encouraging results.

• The system should support a data throughput of 5120 bits per hour. • The maximum delivery delay of the data should not exceed 24 hours.

such as temperature, glacier movements, seismic activity, etc.

(see (Solé et al., 2006) and (Vilella et al., 2009)).

has to cross the equator and four different time zones.

**1.1.4 Global Navigation Satellite Systems**

**1.1.3 Oblique ionosonde**

**1.2 Data transmission**

SAS to Spain are:

researchers to transmit gathered data to laboratories placed on other continents for intensive study. Because of the remoteness of the transmitter placed at the SAS, the system suffers from power restrictions mainly during austral winter. Therefore, maintaining the radio link, even at a reduced throughput, is a challenge. One possible solution to increase data rate, with minimal power, is to improve the spectral efficiency of the physical layer of the radio link while maintaining acceptable performance. The outcomes and conclusions of this research work may be extrapolated to other environments where communication is scarcely possible due to economic or coverage problems. Therefore, the solutions presented in this study may be adopted in other situations, such as communications in developing countries or in any other remote area.

#### **1.1 Remote sensors at the SAS**

In this section we describe the main sensors located at the SAS, including a geomagnetic sensor, a vertical incidence ionosonde, an oblique incidence ionosonde and a Global Navigation Satellite System (GNSS) receiver. They have all been deployed in the premises of the SAS by engineers of the GRECO and scientist of the Observatori de l'Ebre. The geomagnetic sensor, the vertical incidence ionosonde and the GNSS receiver are commercial solutions. The oblique incidence ionosonde, used to sound the ionospheric channel between Antarctica and Spain, was developed by the GRECO in the framework of this research work.

#### **1.1.1 Geomagnetic sensor**

Ground-based geomagnetic observatories provide a time series of accurate measurements of the natural magnetic field vector in a particular location on the Earth's surface. This data is used for several scientific and practical purposes, including the synthesis and updates of global magnetic field models, the study of the solar-terrestrial relationships and the Earth's space environment, and support for other types of geophysical studies.

Once the raw observatory data is processed, it is sent to the World Data Centers, where the worldwide scientific community can access them. International Real-time Magnetic Observatory Network (INTERMAGNET) provides means to access the data by an almost real-time satellite link. The data is packed, sent to the geostationary satellites, and collected by Geomagnetic Information Nodes (GINs), where the information can be accessed freely. However, experience has shown that the satellite link is not 100% reliable, and it is preferable to have alternative means to retrieve the geomagnetic data.

There are three main reasons for designing a transmission backup system by skywave. Firstly, visibility problems appear when trying to reach geostationary satellites from polar latitudes. Secondly, end-to-end reliability can be increased by transmitting each frame repeatedly throughout the day. And finally, the ionospheric channel is freely accessed anywhere, whereas satellite communications have operational costs.

#### **1.1.2 Ionosonde: vertical incidence soundings of the ionosphere**

A vertical incidence ionospheric sounder (VIS) (Zuccheretti et al., 2003) was installed in order to have a sensor providing ionospheric monitoring in this remote region. This ionosonde is also being used to provide information for the High Frequency (HF) radio link employed for data transmission from the SAS to Spain. Data provided by the VIS is used to conduct ionospheric research, mainly to characterize the climatology of the ionospheric characteristics and to investigate the ionospheric effects caused during geomagnetically disturbed periods (see (Solé et al., 2006) and (Vilella et al., 2009)).

## **1.1.3 Oblique ionosonde**

2 Will-be-set-by-IN-TECH

researchers to transmit gathered data to laboratories placed on other continents for intensive study. Because of the remoteness of the transmitter placed at the SAS, the system suffers from power restrictions mainly during austral winter. Therefore, maintaining the radio link, even at a reduced throughput, is a challenge. One possible solution to increase data rate, with minimal power, is to improve the spectral efficiency of the physical layer of the radio link while maintaining acceptable performance. The outcomes and conclusions of this research work may be extrapolated to other environments where communication is scarcely possible due to economic or coverage problems. Therefore, the solutions presented in this study may be adopted in other situations, such as communications in developing countries or in any

In this section we describe the main sensors located at the SAS, including a geomagnetic sensor, a vertical incidence ionosonde, an oblique incidence ionosonde and a Global Navigation Satellite System (GNSS) receiver. They have all been deployed in the premises of the SAS by engineers of the GRECO and scientist of the Observatori de l'Ebre. The geomagnetic sensor, the vertical incidence ionosonde and the GNSS receiver are commercial solutions. The oblique incidence ionosonde, used to sound the ionospheric channel between Antarctica and Spain, was developed by the GRECO in the framework of this research work.

Ground-based geomagnetic observatories provide a time series of accurate measurements of the natural magnetic field vector in a particular location on the Earth's surface. This data is used for several scientific and practical purposes, including the synthesis and updates of global magnetic field models, the study of the solar-terrestrial relationships and the Earth's

Once the raw observatory data is processed, it is sent to the World Data Centers, where the worldwide scientific community can access them. International Real-time Magnetic Observatory Network (INTERMAGNET) provides means to access the data by an almost real-time satellite link. The data is packed, sent to the geostationary satellites, and collected by Geomagnetic Information Nodes (GINs), where the information can be accessed freely. However, experience has shown that the satellite link is not 100% reliable, and it is preferable

There are three main reasons for designing a transmission backup system by skywave. Firstly, visibility problems appear when trying to reach geostationary satellites from polar latitudes. Secondly, end-to-end reliability can be increased by transmitting each frame repeatedly throughout the day. And finally, the ionospheric channel is freely accessed anywhere, whereas

A vertical incidence ionospheric sounder (VIS) (Zuccheretti et al., 2003) was installed in order to have a sensor providing ionospheric monitoring in this remote region. This ionosonde is also being used to provide information for the High Frequency (HF) radio link employed

space environment, and support for other types of geophysical studies.

to have alternative means to retrieve the geomagnetic data.

**1.1.2 Ionosonde: vertical incidence soundings of the ionosphere**

satellite communications have operational costs.

other remote area.

**1.1 Remote sensors at the SAS**

**1.1.1 Geomagnetic sensor**

The oblique ionosonde monitors various parameters to model the HF radiolink between the SAS and Spain (Vilella et al., 2008). These parameters include link availability, power delay profile and frequency dispersion of the channel. The sounder includes a transmitter, placed on the premises of the SAS and a receiver deployed in Spain. The main drawback of the oblique sounder is the difficulty in establishing the ionospheric link. Firstly, the long distance of the link (12700 km) requires four hops to reach the receiver. And secondly, the transmitted signal has to cross the equator and four different time zones.

## **1.1.4 Global Navigation Satellite Systems**

The ionosphere study can be approximated from several points of view. Vertical incidence soundings provide accurate information about electron density profiles below the peak electron density. However, when using this technique the electron profile must be extrapolated from the peak point to the upper limit of the ionosphere. Moreover, the low density of vertical ionosondes, especially in oceans and remote areas, is a serious impairment.

GNSS receivers constitute a high temporal and spatial resolution sounding network which despite gaps over oceans and remote regions, can be used to study fast perturbations affecting local regions, such as Travelling Ionospheric Disturbances and scintillations, or wider regions such as solar flares. Data gathered from GNSS receivers can provide information about the Total Electron Content (TEC) between receivers and satellites by means of proper tomographic modeling approaches. Spatial and temporal variations of the main ionospheric events can be monitored by means of GNSS receivers, especially those placed in the Antarctic Region, which is considered the entrance point of many ionospheric disturbances coming from Solar events. Furthermore, TEC reaches its highest variability peaks in the Antarctica area.

## **1.2 Data transmission**

This chapter will study, analyze and experimentally verify a possible candidate for the physical layer of a long-haul ionospheric data link, focusing on the case SAS-Spain. Preliminary studies of data transmission feasibility over this link were already performed in (Deumal et al., 2006) and (Bergada et al., 2009), with encouraging results.

The first application of this link is the transmission of data generated by a geomagnetic sensor installed at the SAS. Future applications may include sending information of another nature such as temperature, glacier movements, seismic activity, etc.

The minimum requirements regarding the geomagnetic sensor data transmissions from the SAS to Spain are:

**1.2.3 HF communication standards**

2000).

its non-suitability for the purposes of this project.

Regarding the interests of this work it is noted that:

robust to interference modulations.

power and antennas (Vilella et al., 2008).

transmitting in free bands (IEEE802.11, 2007).

impairments of the channel, the environment and other services.

synchronized code in the receiver to despread and retrieve data.

systems is not considered.

with other users.

In this section we briefly review the current communication standards for HF and we justify

Low Rate High Frequency Data Transmission from Very Remote Sensors 127

Due to the proliferation of modems in the field of HF communications, interoperability between equipment from different manufacturers became a problem (NTIA, 1998). Hence the need to standardize communication protocols. Worldwide, there are three organizations proposing standards regarding HF communications: (*i*) the U.S. Department of Defense proposes the Military Standards (MIL-STD-188-110A, 1991; MIL-STD-188-110B, 2000; MIL-STD-188-141A, 1991), (*ii*) the Institute for Telecommunications Science (ITS), which depends on the U.S. Department of Commerce, writes the Federal Standard (FED-STD) and (*iii*) NATO proposes the Standardization Agreements (STANAG-4406, 1999; STANAG-5066,

• The standard modes are designed for primary or secondary services. Therefore:

**–** The bandwidth of the channels is standardized (3 kHz or multiples). Interference reduction, i.e. minimize the output power spectral density, with other transmitting

**–** No modes are considered based on short sporadic burst transfers to reduce interference

**–** There are anti-jamming techniques (see MIL-STD-188-148) for additional application on a appropriate communication standard, but the proposals are not based on intrinsically

• Robust configurations require a minimum signal to noise ratio (SNR) of 0 dB at 3 kHz bandwidth, which is not common in this link under the specified conditions of transmitted

We conclude that the configurations proposed by current standards do not meet the desirable characteristics for the type of communication that is required in this work, and consequently, a new proposal should be suggested. In this chapter, we study a number of alternatives based on the use of Direct Sequence Spread Spectrum (DS-SS) techniques in order to cope with the

Spread Spectrum (SS) techniques are described by (Pickholtz et al., 1982) as a kind of transmission in which signal occupies a greater bandwidth than the necessary bandwidth to send the information; bandwidth spreading is achieved by an independent data source, and a

SS began to be developed especially for military purposes in the mid twentieth century and has continued in the forefront of research to present, which is, nowadays, a key point for the 3G mobile cellular systems (Third Generation Partnership Project, 1999) and wireless systems

**2. Data transmission with Direct Sequence Spread Spectrum techniques**

#### **1.2.1 Constraints**

The extreme conditions prevailing at the SAS impose a number of restrictions that affect the transmission system. We highlight the following ones:

#### **1.2.2 Approach**

This section justifies the need for a new data communication system adapted to the requirements of the project and presents the main ideas of this proposal. Firstly, we review the mechanisms that exist worldwide regarding the regulation of occupation of the radio spectrum. Then we review the features of current standards of HF data communications and discuss the non-suitability of these to the requirements of the project.

The International Telecommunication Union (ITU) is responsible for regulating the use of radio spectrum. From the point of view of frequency allocation, it has divided the world into three regions. Broadly speaking, region 1 comprises Europe and Africa, Asia and Oceania constitute region 2 and North and South America region 3.

In each region, the ITU recommends the allocation of each frequency band to one or several services. When multiple services are attributed to the same frequency band in the same region, these fall into two categories: primary or secondary. The ones that are classified as secondary services can not cause interference with the primary services and can not claim protection from interference from the primary services; however, they can demand protection from interference from other secondary services attributed afterwards.

Given these considerations, we propose a system transmission with the following guidelines:

To meet these requirements, we propose a system with the following characteristics:

Moreover, given the ionospheric channel measures described in (Vilella et al., 2008) the following additional features are required:

## **1.2.3 HF communication standards**

4 Will-be-set-by-IN-TECH

The extreme conditions prevailing at the SAS impose a number of restrictions that affect the

• The transmission power should be minimal. It is noted that the SAS is inhabited only during the austral summer, approximately from November to March. During this period there is no limitation regarding the maximum power consumption. However, the transmission system is designed to continue operating during the austral winter, when energy is obtained entirely from batteries powered by wind generators and solar panels.

• Environmental regulations applicable at the site advise against the installation of large

This section justifies the need for a new data communication system adapted to the requirements of the project and presents the main ideas of this proposal. Firstly, we review the mechanisms that exist worldwide regarding the regulation of occupation of the radio spectrum. Then we review the features of current standards of HF data communications and

The International Telecommunication Union (ITU) is responsible for regulating the use of radio spectrum. From the point of view of frequency allocation, it has divided the world into three regions. Broadly speaking, region 1 comprises Europe and Africa, Asia and Oceania

In each region, the ITU recommends the allocation of each frequency band to one or several services. When multiple services are attributed to the same frequency band in the same region, these fall into two categories: primary or secondary. The ones that are classified as secondary services can not cause interference with the primary services and can not claim protection from interference from the primary services; however, they can demand protection

Given these considerations, we propose a system transmission with the following guidelines: • It can not cause harmful interference to any other service stations (primary or secondary).

Moreover, given the ionospheric channel measures described in (Vilella et al., 2008) the

• Robustness against noise (possibility of working with negative signal to noise ratio).

To meet these requirements, we propose a system with the following characteristics: • Reduced transmission power (accordingly with the consumption constrains).

structures that would be needed to install certain types of directive antennas.

**1.2.1 Constraints**

**1.2.2 Approach**

transmission system. We highlight the following ones:

Hence the power amplifier is set to a maximum of only 250 watts.

discuss the non-suitability of these to the requirements of the project.

from interference from other secondary services attributed afterwards.

• It can not claim protection from interference from other services.

• Robustness against time and frequency dispersive channels.

• Low power spectral density. • Robustness to interference.

• Sporadic communications.

• Burst transmissions (few seconds).

following additional features are required:

constitute region 2 and North and South America region 3.

In this section we briefly review the current communication standards for HF and we justify its non-suitability for the purposes of this project.

Due to the proliferation of modems in the field of HF communications, interoperability between equipment from different manufacturers became a problem (NTIA, 1998). Hence the need to standardize communication protocols. Worldwide, there are three organizations proposing standards regarding HF communications: (*i*) the U.S. Department of Defense proposes the Military Standards (MIL-STD-188-110A, 1991; MIL-STD-188-110B, 2000; MIL-STD-188-141A, 1991), (*ii*) the Institute for Telecommunications Science (ITS), which depends on the U.S. Department of Commerce, writes the Federal Standard (FED-STD) and (*iii*) NATO proposes the Standardization Agreements (STANAG-4406, 1999; STANAG-5066, 2000).

Regarding the interests of this work it is noted that:

	- **–** The bandwidth of the channels is standardized (3 kHz or multiples). Interference reduction, i.e. minimize the output power spectral density, with other transmitting systems is not considered.
	- **–** No modes are considered based on short sporadic burst transfers to reduce interference with other users.
	- **–** There are anti-jamming techniques (see MIL-STD-188-148) for additional application on a appropriate communication standard, but the proposals are not based on intrinsically robust to interference modulations.

We conclude that the configurations proposed by current standards do not meet the desirable characteristics for the type of communication that is required in this work, and consequently, a new proposal should be suggested. In this chapter, we study a number of alternatives based on the use of Direct Sequence Spread Spectrum (DS-SS) techniques in order to cope with the impairments of the channel, the environment and other services.

## **2. Data transmission with Direct Sequence Spread Spectrum techniques**

Spread Spectrum (SS) techniques are described by (Pickholtz et al., 1982) as a kind of transmission in which signal occupies a greater bandwidth than the necessary bandwidth to send the information; bandwidth spreading is achieved by an independent data source, and a synchronized code in the receiver to despread and retrieve data.

SS began to be developed especially for military purposes in the mid twentieth century and has continued in the forefront of research to present, which is, nowadays, a key point for the 3G mobile cellular systems (Third Generation Partnership Project, 1999) and wireless systems transmitting in free bands (IEEE802.11, 2007).

**2.1 Robustness against interference**

against interference.

Then (Pickholtz et al., 1982),

 *Eb No*

It is noted that when *Bss* increases *Eb*

of *Eb No i*(*t*)

is achieved.

Ionospheric communications have global coverage range. Consequently, any system operating in a given area might potentially interfere with other remote systems operating at the same frequency band. Hence the transmission system proposed in this work might be interfered with primary or secondary services that are assigned the same frequency band. For these reasons it is appropriate to review the characteristics of DS-SS regarding robustness

Low Rate High Frequency Data Transmission from Very Remote Sensors 129

Let a DS-SS based system that transmits *Rb* bits per second in a bandwidth *Bss* (*Bss* � *Rb*) in the presence of additive white Gaussian noise *z*(*t*) with power spectral density *No* [*W*/*Hz*]

> *Pn Pn* + *Pi*

where *Pn* = *BssNo* is the noise power within the transmission bandwidth and *P* = *EbRb* is the signal power. We can deduce from Equation 6 that we can reduce the effect of interfering signals by increasing *Bss*. In other words, as *Bss* = *L* · *Rb*, the larger the spreading factor the

*Gp* and is a measure of the robustness of a spread spectrum system against interference. In

in the same proportion, whereas an increase of *Bss* implies an equivalent improvement

improvement regarding narrowband interfering signals whereas no improvement over noise

Feasibility studies of DS-SS systems with different types of interference can be found in the literature. See, for instance, (Schilling et al., 1980) when the interference is a narrowband

According to the analysis described in (Vilella et al., 2008), the ionospheric channel established between the SAS and Spain shows a maximum multipath delay spread (*τmax*) that varies, depending on time and frequency, between 0.5 ms and 2.5 ms. Therefore, the coherence bandwidth of the channel, which can be considered as approximately the inverse of the maximum multipath delay spread (Proakis, 1995), can be narrower than 400 Hz. In case of transmitting with a wider bandwidth the channel would be frequency selective and distortion due to multipath would arise. Below, the properties of DS-SS against multipath are discussed.

Let a DS-SS based system with bandwidth *Bss* in a channel with coherence bandwidth *Wc* ∼ 1 *τmax* � *Bss*. If symbol time *Ts* � *τmax* intersymbol interference due to multipath can be

, as *Pi* is unchanged. To summarize, the use of spread spectrum provides

DS-SS systems the processing gain coincides with the spreading sequence length (*L*).

*No z*(*t*) = *P Pn Bss Rb*

*rss*(*t*) = *sss*(*t*) + *i*(*t*) + *z*(*t*). (5)

*N*0 *N*0 + *Pi Bss*

does not change because *Pn* = *NoBss* increases

, (6)

*Rb* is called the process gain

and narrowband interference *i*(*t*) with power *Pi*. At the receiver side:

= *P Pn Bss Rb*

*z*(*t*),*i*(*t*)

lower the degradation due to interfering signals. The quotient *Bss*

signal and (Milstein, 1988) for multiple interfering signals.

**2.2 Robustness against multipath channels**

In the field of HF communications new techniques have always been slowly introduced due to a widespread sense that reliable communications were not feasible in this frequency band, while improvements of its implementation would be irrelevant. However, SS techniques have been suggested several times as suitable for the lower band of frequencies (i.e. LF, MF and, by extension, HF) (see (Enge & Sarwate, 1987)), since the intrinsic characteristics of SS systems to cope with multipath and interference (typical ionospheric channel characteristics).

There are three types of spread spectrum systems (Peterson et al., 1995): Direct Sequence, Frequency Hopping and hybrid systems composed by a mixture of both. In this study we will focus on Direct Sequence schemes.

DS-SS systems spread the spectrum by multiplying the information data by a spreading sequence. Consider the following model (Proakis, 1995):

$$s\_{\rm ss}(t) = \sum\_{i=0}^{N\_s - 1} d\_l c(t - iT\_s)\_\prime \ c(t) = \sum\_{l=0}^{L-1} c\_l p(t - lT\_\varepsilon)\_\prime \tag{1}$$

where *di* denotes the *ith* symbol, of length *Ts*, of a modulated signal:

$$d = \{d\_0, d\_1, \dots, d\_{N\_s - 1}\} \tag{2}$$

and *cl* are the chips 1, of length *Tc*, of a spreading sequence of length *L*:

$$\overline{\mathfrak{C}} = \{c\_{0}, c\_{1}, \dots, c\_{L-1}\} \tag{3}$$

and *p*(*t*) is a pulse shaping defined as

$$p(t) = \begin{cases} 1, t \in [0, T\_{\mathcal{C}}) \\ 0 \Rightarrow \text{otherwise} \end{cases} \tag{4}$$

In addition, it holds that *LTc* = *Ts* and thus if the base band signal is formed by the symbols *di* and occupies a bandwidth of 1 *Ts* the spread spectrum signal *sss*(*t*) occupies a bandwidth of 1 *Tc* = *L* 1 *Ts* .

The spreading sequence *c* should have good properties of autocorrelation and cross-correlation in order to ease the detection and synchronization at the receiver side.

Some of the main advantages of a system based on DS-SS are: (*i*) jamming and interference robustness, (*ii*) privacy, (*iii*) ability to use Code Division Multiple Access (CDMA) and (*iv*) robustness against multipath and time variant channels. On the other hand, the drawbacks of this technique are: (*i*) bandwidth inefficiency and (*ii*) receiver complexity: chip-level synchronization, symbol despreading (DS-SS signaling) and channel estimation and detection (RAKE receiver) (Viterbi, 1995).

Throughout the following sections we will discuss the most important considerations that justify the choice of DS-SS; as well as the technical basis to design the data modem for the ionospheric link between the SAS and Spain.

1 The bits of a spreading sequence are called chips

#### **2.1 Robustness against interference**

6 Will-be-set-by-IN-TECH

In the field of HF communications new techniques have always been slowly introduced due to a widespread sense that reliable communications were not feasible in this frequency band, while improvements of its implementation would be irrelevant. However, SS techniques have been suggested several times as suitable for the lower band of frequencies (i.e. LF, MF and, by extension, HF) (see (Enge & Sarwate, 1987)), since the intrinsic characteristics of SS systems to

There are three types of spread spectrum systems (Peterson et al., 1995): Direct Sequence, Frequency Hopping and hybrid systems composed by a mixture of both. In this study we

DS-SS systems spread the spectrum by multiplying the information data by a spreading

*L*−1 ∑ *l*=0

*d* = {*d*0, *d*1, ..., *dNs*−1} (2)

*c* = {*c*0, *c*1, ..., *cL*−1} (3)

*Ts* the spread spectrum signal *sss*(*t*) occupies a bandwidth of

0 ⇒ otherwise (4)

*cl p*(*t* − *lTc*), (1)

*dic*(*t* − *iTs*), *c*(*t*) =

1, *t* ∈ [0, *Tc*)

In addition, it holds that *LTc* = *Ts* and thus if the base band signal is formed by the symbols

The spreading sequence *c* should have good properties of autocorrelation and cross-correlation in order to ease the detection and synchronization at the receiver side.

Some of the main advantages of a system based on DS-SS are: (*i*) jamming and interference robustness, (*ii*) privacy, (*iii*) ability to use Code Division Multiple Access (CDMA) and (*iv*) robustness against multipath and time variant channels. On the other hand, the drawbacks of this technique are: (*i*) bandwidth inefficiency and (*ii*) receiver complexity: chip-level synchronization, symbol despreading (DS-SS signaling) and channel estimation and detection

Throughout the following sections we will discuss the most important considerations that justify the choice of DS-SS; as well as the technical basis to design the data modem for the

cope with multipath and interference (typical ionospheric channel characteristics).

will focus on Direct Sequence schemes.

and *p*(*t*) is a pulse shaping defined as

*di* and occupies a bandwidth of 1

(RAKE receiver) (Viterbi, 1995).

ionospheric link between the SAS and Spain.

1 The bits of a spreading sequence are called chips

1 *Tc* = *L* 1 *Ts* .

sequence. Consider the following model (Proakis, 1995):

*Ns*−1 ∑ *i*=0

where *di* denotes the *ith* symbol, of length *Ts*, of a modulated signal:

and *cl* are the chips 1, of length *Tc*, of a spreading sequence of length *L*:

*p*(*t*) =

*sss*(*t*) =

Ionospheric communications have global coverage range. Consequently, any system operating in a given area might potentially interfere with other remote systems operating at the same frequency band. Hence the transmission system proposed in this work might be interfered with primary or secondary services that are assigned the same frequency band. For these reasons it is appropriate to review the characteristics of DS-SS regarding robustness against interference.

Let a DS-SS based system that transmits *Rb* bits per second in a bandwidth *Bss* (*Bss* � *Rb*) in the presence of additive white Gaussian noise *z*(*t*) with power spectral density *No* [*W*/*Hz*] and narrowband interference *i*(*t*) with power *Pi*. At the receiver side:

$$r\_{\rm ss}(t) = s\_{\rm ss}(t) + i(t) + z(t). \tag{5}$$

Then (Pickholtz et al., 1982),

$$\left(\frac{E\_b}{N\_0}\right)\_{z(t),i(t)} = \frac{P}{P\_n}\frac{B\_{ss}}{R\_b}\frac{P\_n}{P\_n + P\_i} = \frac{P}{P\_n}\frac{B\_{ss}}{R\_b}\frac{N\_0}{N\_0 + \frac{P\_i}{B\_{ss}}}\tag{6}$$

where *Pn* = *BssNo* is the noise power within the transmission bandwidth and *P* = *EbRb* is the signal power. We can deduce from Equation 6 that we can reduce the effect of interfering signals by increasing *Bss*. In other words, as *Bss* = *L* · *Rb*, the larger the spreading factor the lower the degradation due to interfering signals. The quotient *Bss Rb* is called the process gain *Gp* and is a measure of the robustness of a spread spectrum system against interference. In DS-SS systems the processing gain coincides with the spreading sequence length (*L*).

It is noted that when *Bss* increases *Eb No z*(*t*) does not change because *Pn* = *NoBss* increases in the same proportion, whereas an increase of *Bss* implies an equivalent improvement of *Eb No i*(*t*) , as *Pi* is unchanged. To summarize, the use of spread spectrum provides improvement regarding narrowband interfering signals whereas no improvement over noise is achieved.

Feasibility studies of DS-SS systems with different types of interference can be found in the literature. See, for instance, (Schilling et al., 1980) when the interference is a narrowband signal and (Milstein, 1988) for multiple interfering signals.

#### **2.2 Robustness against multipath channels**

According to the analysis described in (Vilella et al., 2008), the ionospheric channel established between the SAS and Spain shows a maximum multipath delay spread (*τmax*) that varies, depending on time and frequency, between 0.5 ms and 2.5 ms. Therefore, the coherence bandwidth of the channel, which can be considered as approximately the inverse of the maximum multipath delay spread (Proakis, 1995), can be narrower than 400 Hz. In case of transmitting with a wider bandwidth the channel would be frequency selective and distortion due to multipath would arise. Below, the properties of DS-SS against multipath are discussed.

Let a DS-SS based system with bandwidth *Bss* in a channel with coherence bandwidth *Wc* ∼ 1 *τmax* � *Bss*. If symbol time *Ts* � *τmax* intersymbol interference due to multipath can be

**2.4 Flexibility regarding spectral efficiency**

system.

*k* is the corresponding number of bits per symbol).

and we note the spectral efficiency of each of them.

**2.4.1 DS-SS M-ary signaling**

Let a set of *M* spreading sequences *Q* =

*sss*(*t*) =

*Ps* =

in the demodulation of the coded bits in *d*:

*M*−1 ∑ *p*=1

*P*2 = 1 2

(−1)*p*+1

certain sequence can be computed from the following expression (Proakis, 1995):

*Ps* + (1 − *Ps*) *Q*

*Ns*−1 ∑ *i*=0 *di L*−1 ∑ *l*=0 *c* (*v*)

log2(*M*) bits of information. Then

The signal model expressed by Equation 1 is able to transmit *k* = log2 *K* bits (*b*

� *c*(1) , *c*(2)

modulated in *di* during a period *Ts* (*K* is the number of possible modulation symbols in *d* and

Low Rate High Frequency Data Transmission from Very Remote Sensors 131

The spectral efficiency (*Css* = *k*/(*Ts* · *Bss*)), expressed in [*bits*/*s*/*Hz*] and defined as the ratio between bit-rate and transmission bandwidth, is *Gp* times lower than the non spreading

There are several alternatives to increase spectral efficiency without decreasing process gain (and hence, robustness to interference) at the expense of increasing computational cost of the receiver. In the following sections we describe two of them: DS-SS M-ary signaling and quadriphase spreading. We briefly present the signal model, a study of the probability of error

relationship (orthogonal or nearly orthogonal according to Equation 8). Suppose that a certain sequence *v* from the previous set (*v* ∈ [1, *M*]) is transmitted depending on the value of *m* =

This technique is called DS-SS M-ary signaling (see, for example, (Enge & Sarwate, 1987) for orthogonal sequences). On the receiver side, the optimum demodulator correlates the received signal with a replica of each of the *M* possible sequences belonging to the set *Q*. A noncoherent detector will make a decision based on the computation of the maximum likelihood of the *M* envelopes at the output of each correlator. The probability *Ps* of detecting an incorrect sequence in the presence of only additive white noise is given by (Proakis, 1995):

> � *M* − 1 *p*

*P*1 = 2*m*−1 2*m* − 1

The probability *P*1 of making an error in the demodulation of coded bits transmitted in a

Once the sequence is detected we proceed to compute the probability *P*2 of making an error

⎛ ⎝ � 2 *k* � *Eb No* �� ⎞

� 1 *p* + 1 *e* − *p*

*l p*(*t* − *iTs* − *lTC*) =

,..., *c*(*M*)

�

*Ns*−1 ∑ *i*=0

*p*+1 (*m*+*k*) *Eb*

*Ps*. (12)

*dic*(*v*)

(0) 0 ... *b*

that satisfy a certain correlation

(*t*). (10)

*No* . (11)

⎠ , (13)

(0) *k*−1)

neglected. Moreover, if *Ts* � 1 *υmax* (where *υmax* denotes the maximum Doppler spread) the channel is almost invariant during a symbol time. Under these conditions it can be shown that (Proakis, 1995):

$$r\_{\rm ss}^{(k)}(t) = \sum\_{n=1}^{N} h\left(\frac{n}{B\_{\rm ss}}\right) s\_{\rm ss}^{(k)}\left(t - \frac{n}{B\_{\rm ss}}\right) + z(t),\tag{7}$$

where *h n Bss* denotes a coefficient of the equivalent low-pass of the channel impulse response, *sss*(*t*) corresponds to the base band spread signal defined in Equation 1, (*k*) denotes the contribution due to symbol *k*, *N* = *τmaxBss* is the number of non zero channel taps (since *sss*(*t*) has a limited bandwidth of *Bss*) and *z*(*t*) is additive white Gaussian noise. In consequence, the signal reception is formed by delayed replicas of the transmitted signal. Then, we substitute Equation 1 in Equation 7 and apply an array of correlators to correlate the received signal with *N* copies of the spreading sequence *c* (each of them delayed a chip time). Let *c* be a sequence with good properties of circular autocorrelation:

$$\rho(m) = \sum\_{l=1}^{L} c\_l c\_{l+m} \begin{cases} 1 & m = 0 \\ 0 \Rightarrow \text{otherwise} \end{cases} \tag{8}$$

The output *Um* of each correlator can be expressed as:

$$d\mathcal{U}\_{\rm ll} = d\_k \hbar \left(\frac{m}{B\_{\rm ss}}\right) + \int\_0^{T\_s} c(t)z(t + \frac{m}{B\_{\rm ss}}) \, dt, \; m \in [0, N-1]. \tag{9}$$

Therefore, at the output of each correlator we obtain each transmitted symbol (*dk*) multiplied by a channel coefficient *h m Bss* plus a noise term. Hence, the use of DS-SS can take advantage of different replicas of the signal if correctly combined. The most general linear combination is the criterion of Maximal Ratio Combining that chooses the coefficients that maximize instantaneous SNR (Peterson et al., 1995). To properly apply this method it is mandatory to know the coefficients of the channel. Alternatively, the outputs of the correlators can be equally weighed (Equal Gain Combining), thus simplifying the receiver at the expense of worse performance.

#### **2.3 Transmission with low spectral density power**

One of the requirements of the proposed transmission system consists in causing minimal interference with primary and secondary services. For this purpose we propose alternatives to minimize power spectral density. We should note that as process gain increases, DS-SS techniques enable transmission with arbitrarily low power density. Suppose the transmission of a data stream *d* using a bandwidth *Bd* and power *P*. Then the average power spectral density is *D* = *P Bd W Hz* . Under the same conditions of power consider the transmission of the same data stream with DS-SS (*sss*(*t*)) by means of a spreading sequence *c*(*t*) of length *L*. Then, the spectral occupancy of *sss*(*t*) will be at least *L* · *Bd* and the average power spectral density will be *Dss* = *P L*·*Bd W Hz* .

Therefore, the use of DS-SS involves an average reduction of power spectral density by a factor equal to the process gain *Gp* = *L*. Then, the spectral occupancy proportionally increases; however, it is not an inconvenience in this case since there is no limitation in this regard.

#### **2.4 Flexibility regarding spectral efficiency**

8 Will-be-set-by-IN-TECH

channel is almost invariant during a symbol time. Under these conditions it can be shown

response, *sss*(*t*) corresponds to the base band spread signal defined in Equation 1, (*k*) denotes the contribution due to symbol *k*, *N* = *τmaxBss* is the number of non zero channel taps (since *sss*(*t*) has a limited bandwidth of *Bss*) and *z*(*t*) is additive white Gaussian noise. In consequence, the signal reception is formed by delayed replicas of the transmitted signal. Then, we substitute Equation 1 in Equation 7 and apply an array of correlators to correlate the received signal with *N* copies of the spreading sequence *c* (each of them delayed a chip time).

denotes a coefficient of the equivalent low-pass of the channel impulse

1, *m* = 0

*m Bss*

plus a noise term. Hence, the use of DS-SS can take advantage

. Under the same conditions of power consider the transmission of

*υmax* (where *υmax* denotes the maximum Doppler spread) the

+ *z*(*t*), (7)

0 ⇒ otherwise (8)

) *dt*, *m* ∈ [0, *N* − 1]. (9)

neglected. Moreover, if *Ts* � 1

*r* (*k*) *ss* (*t*) =

*N* ∑ *n*=1 *h n Bss s* (*k*) *ss t* − *n Bss*

Let *c* be a sequence with good properties of circular autocorrelation:

*L* ∑ *l*=1

*cl cl*+*m*

*c*(*t*)*z*(*t* +

Therefore, at the output of each correlator we obtain each transmitted symbol (*dk*) multiplied

of different replicas of the signal if correctly combined. The most general linear combination is the criterion of Maximal Ratio Combining that chooses the coefficients that maximize instantaneous SNR (Peterson et al., 1995). To properly apply this method it is mandatory to know the coefficients of the channel. Alternatively, the outputs of the correlators can be equally weighed (Equal Gain Combining), thus simplifying the receiver at the expense of

One of the requirements of the proposed transmission system consists in causing minimal interference with primary and secondary services. For this purpose we propose alternatives to minimize power spectral density. We should note that as process gain increases, DS-SS techniques enable transmission with arbitrarily low power density. Suppose the transmission of a data stream *d* using a bandwidth *Bd* and power *P*. Then the average power spectral

the same data stream with DS-SS (*sss*(*t*)) by means of a spreading sequence *c*(*t*) of length *L*. Then, the spectral occupancy of *sss*(*t*) will be at least *L* · *Bd* and the average power spectral

Therefore, the use of DS-SS involves an average reduction of power spectral density by a factor equal to the process gain *Gp* = *L*. Then, the spectral occupancy proportionally increases; however, it is not an inconvenience in this case since there is no limitation in this regard.

*ρ*(*m*) =

 *m Bss* + *Ts* 0

The output *Um* of each correlator can be expressed as:

 *m Bss*

**2.3 Transmission with low spectral density power**

*L*·*Bd*

 *W Hz* .

*Um* = *dkh*

that (Proakis, 1995):

 *n Bss*

by a channel coefficient *h*

worse performance.

density is *D* = *P*

density will be *Dss* = *P*

*Bd W Hz*

where *h*

The signal model expressed by Equation 1 is able to transmit *k* = log2 *K* bits (*b* (0) 0 ... *b* (0) *k*−1) modulated in *di* during a period *Ts* (*K* is the number of possible modulation symbols in *d* and *k* is the corresponding number of bits per symbol).

The spectral efficiency (*Css* = *k*/(*Ts* · *Bss*)), expressed in [*bits*/*s*/*Hz*] and defined as the ratio between bit-rate and transmission bandwidth, is *Gp* times lower than the non spreading system.

There are several alternatives to increase spectral efficiency without decreasing process gain (and hence, robustness to interference) at the expense of increasing computational cost of the receiver. In the following sections we describe two of them: DS-SS M-ary signaling and quadriphase spreading. We briefly present the signal model, a study of the probability of error and we note the spectral efficiency of each of them.

#### **2.4.1 DS-SS M-ary signaling**

Let a set of *M* spreading sequences *Q* = � *c*(1) , *c*(2) ,..., *c*(*M*) � that satisfy a certain correlation relationship (orthogonal or nearly orthogonal according to Equation 8). Suppose that a certain sequence *v* from the previous set (*v* ∈ [1, *M*]) is transmitted depending on the value of *m* = log2(*M*) bits of information. Then

$$s\_{\rm ss}(t) = \sum\_{i=0}^{N\_\circ - 1} d\_i \sum\_{l=0}^{L-1} c\_l^{(v)} p(t - iT\_\circ - lT\_\circ) = \sum\_{i=0}^{N\_\circ - 1} d\_i c^{(v)}(t). \tag{10}$$

This technique is called DS-SS M-ary signaling (see, for example, (Enge & Sarwate, 1987) for orthogonal sequences). On the receiver side, the optimum demodulator correlates the received signal with a replica of each of the *M* possible sequences belonging to the set *Q*. A noncoherent detector will make a decision based on the computation of the maximum likelihood of the *M* envelopes at the output of each correlator. The probability *Ps* of detecting an incorrect sequence in the presence of only additive white noise is given by (Proakis, 1995):

$$P\_{\mathbf{S}} = \sum\_{p=1}^{M-1} (-1)^{p+1} \binom{M-1}{p} \frac{1}{p+1} e^{-\frac{p}{p+1}(m+k)\frac{\mathbb{E}\_k}{\mathbb{N}\_0}}.\tag{11}$$

The probability *P*1 of making an error in the demodulation of coded bits transmitted in a certain sequence can be computed from the following expression (Proakis, 1995):

$$P\_1 = \frac{2^{m-1}}{2^m - 1} P\_s.\tag{12}$$

Once the sequence is detected we proceed to compute the probability *P*2 of making an error in the demodulation of the coded bits in *d*:

$$P\_2 = \frac{1}{2}P\_s + (1 - P\_s) \bigotimes \left(\sqrt{\frac{2}{k}\left(\frac{E\_b}{N\_o}\right)'}\right) \tag{13}$$

In a symbol time *Ts* we send *k* + *m* bits (*b*

the computational cost of the receiver.

{*c*(1)

*sss*(*t*) =

, ··· , *c*(*M*/2)

*Ns*−1 ∑ *i*=0

**2.4.2 DS-SS M-ary signaling + Quadriphase**

� {*di*} *c*(*v*1)

symbol encoded in the modulation of *di*).

the presence of only additive white Gaussian noise is:

It is worth noting that the factor 1/2 multiplying *Eb*

*M*/2−1 ∑ *p*=1

Then, the probability *P*2 for BPSK and QPSK, respectively, is:

*PsPs* + 2*Ps*(1 − *Ps*)*Q*�

*P*2 = 1 2 (−1)*p*+1

*Ps* =

sequence. Consequently, the spectral efficiency is:

*C* = *k* + *m Gp*

} on one side and *Qi* = {*c*(*M*/2+1)

DS-SS M-ary signaling on both the real and the imaginary part of *d*. Thus,

(*t*) + *j* · � {*di*} *c*(*v*2)

(1) 0 ··· *b*

= *Css* +

the bits modulated with *d* and *m* corresponds to the bits used in the choice of the spreading

Low Rate High Frequency Data Transmission from Very Remote Sensors 133

Therefore, the larger the *M* the lower the BER for a given SNR per bit and the greater the spectral efficiency at the expense of greater computational cost of the receiver. In addition, the larger the number of bits per symbol the better the BER for a given SNR per bit and the lower

Another alternative is to divide the set of *M* sequences into two subsets: *Qr* =

(*t*)

This variant is called quadriphase chip spreading and permits us to send *M* = 2 log2 (*M*/2) bits per symbol by choosing a sequence from each of the two sets (plus *k* additional bits per

At the receiver side, the demodulator correlates the received signal with a replica of each of the *M* possible sequences. The detector will decide on the envelopes computed at the output of correlators corresponding to the sequences of the subset *Qr* and a similar decision on the subset *Qi*. The probability of incorrectly detecting a sequence from both the set *Qr* and *Qi* in

> *M*/2 − 1 *p*

energy is equally distributed between real and imaginary parts (see Equation 16). The probability *P*1 of incorrectly demodulating the bits coded in a sequence that belongs to the

> *P*1 = 2(*m*/2)−1 2(*m*/2) − 1

We then discuss the probability *P*2 of error on the bits modulated in *d* in the case of BPSK and QPSK. For BPSK, the decision is based on the sign of the sum of the two outputs of the correlators for the two sequences detected in the previous step. The case of QPSK is equivalent to two BPSK with half the SNR per each bit, both independently demodulated from the detection of two sequences (the corresponding to subsets *Qr* and *Qi*, respectively).

subset *Qr* or *Qi* can be obtained by applying an equation analogous to Equation 12:

 1 *p* + 1 *e* − *p p*+1 ( *m*+*k* 2 ) *Eb*

*m Gp*

(1) *k*−1*b* (1) 0 ··· *b*

= *Css* + log2 *M*

, ··· , *c*(*M*)

(1)

*Gp*

*m*−1), where *k* corresponds to

. (15)

} on the other side. Then, apply

*No* . (17)

, *v*1 ∈ [1, *M*/2], *v*2 ∈ [*M*/2 + 1, *M*]. (16)

*No* comes from considering that the symbol

*Ps*. (18)

*bpsk* + (1 − *Ps*)(1 − *Ps*)*Qbpsk* , (19)

Fig. 1. Probability of error versus SNR per bit using DS-SS M-ary signaling for various values of *M* (32, 64, 128) and *k* (*k* = 0 : no modulation, *k* = 1: BPSK, *k* = 2: QPSK). Probability is analytically (A) computed and derived from algorithm simulations (S)

where *Eb No* � = *Eb No* (*m* + *k*). Finally, the joint probability *Pb* of bit error considering the contribution of both mechanisms is:

$$P\_b = \frac{m \cdot P\_1 + k \cdot P\_2}{m + k}.\tag{14}$$

Figure 1 shows that the higher the *M*, the lower the SNR per bit required to obtain a certain BER. It can be explained by the fact that *L* increases as *M* (in a DS-SS system) and so does the process gain. It can be shown that the minimum SNR per bit required to obtain an arbitrarily small BER when *M* → ∞ is -1.6 dB. Figure 1 also shows that the larger the *k*, the smaller the SNR per bit required to achieve a given BER. This apparent contradiction can be derived from the following two arguments: (*i*) for a given bit-rate, a high value of *k* enables the reduction of transmission bandwidth (and thus reduction of noise) and hence, improve the probability of finding the transmitted sequence (Equation 13). (*ii*) The probability *Pb* of total error (Equation 14) is a balance between *P*1 and *P*2. The second term in *P*2 (Equation 13) derives from the probability of error in demodulating the bits in *d* once the sequence is successfully detected. So, if this term is lower than both the first term in Equation 13 and *P*1 (Equation 12) the use of any kind of modulation will not result in significant degradation in *Pb*.

In a symbol time *Ts* we send *k* + *m* bits (*b* (1) 0 ··· *b* (1) *k*−1*b* (1) 0 ··· *b* (1) *m*−1), where *k* corresponds to the bits modulated with *d* and *m* corresponds to the bits used in the choice of the spreading sequence. Consequently, the spectral efficiency is:

$$\mathcal{C} = \frac{k+m}{\mathcal{G}\_p} = \mathcal{C}\_{ss} + \frac{m}{\mathcal{G}\_p} = \mathcal{C}\_{ss} + \frac{\log\_2 M}{\mathcal{G}\_p}.\tag{15}$$

Therefore, the larger the *M* the lower the BER for a given SNR per bit and the greater the spectral efficiency at the expense of greater computational cost of the receiver. In addition, the larger the number of bits per symbol the better the BER for a given SNR per bit and the lower the computational cost of the receiver.

#### **2.4.2 DS-SS M-ary signaling + Quadriphase**

10 Will-be-set-by-IN-TECH

DS−SS Signaling

(A) 32−0 (S) 32−0 (A) 64−0 (S) 64−0 (S) 128−0 (A) 32−1 (S) 32−1 (A) 32−2 (S) 32−2

0 2 4 6 8 10

Eb/No [dB]

(*m* + *k*). Finally, the joint probability *Pb* of bit error considering the

*m* + *k* . (14)

Fig. 1. Probability of error versus SNR per bit using DS-SS M-ary signaling for various values of *M* (32, 64, 128) and *k* (*k* = 0 : no modulation, *k* = 1: BPSK, *k* = 2: QPSK). Probability is

*Pb* = *m* · *P*1 + *k* · *P*2

Figure 1 shows that the higher the *M*, the lower the SNR per bit required to obtain a certain BER. It can be explained by the fact that *L* increases as *M* (in a DS-SS system) and so does the process gain. It can be shown that the minimum SNR per bit required to obtain an arbitrarily small BER when *M* → ∞ is -1.6 dB. Figure 1 also shows that the larger the *k*, the smaller the SNR per bit required to achieve a given BER. This apparent contradiction can be derived from the following two arguments: (*i*) for a given bit-rate, a high value of *k* enables the reduction of transmission bandwidth (and thus reduction of noise) and hence, improve the probability of finding the transmitted sequence (Equation 13). (*ii*) The probability *Pb* of total error (Equation 14) is a balance between *P*1 and *P*2. The second term in *P*2 (Equation 13) derives from the probability of error in demodulating the bits in *d* once the sequence is successfully detected. So, if this term is lower than both the first term in Equation 13 and *P*1 (Equation 12) the use of

analytically (A) computed and derived from algorithm simulations (S)

any kind of modulation will not result in significant degradation in *Pb*.

10−4

 *Eb No* � = *Eb No*

contribution of both mechanisms is:

where

10−3

10−2

BER

10−1

100

Another alternative is to divide the set of *M* sequences into two subsets: *Qr* = {*c*(1) , ··· , *c*(*M*/2) } on one side and *Qi* = {*c*(*M*/2+1) , ··· , *c*(*M*) } on the other side. Then, apply DS-SS M-ary signaling on both the real and the imaginary part of *d*. Thus,

$$s\_{\rm ss}(t) = \sum\_{i=0}^{N\_t - 1} \left( \Re \left\{ d\_i \right\} c^{\left( v\_1 \right)}(t) + j \cdot \Im \left\{ d\_i \right\} c^{\left( v\_2 \right)}(t) \right), \ v\_1 \in [1, M/2], \\ v\_2 \in [M/2 + 1, M]. \tag{16}$$

This variant is called quadriphase chip spreading and permits us to send *M* = 2 log2 (*M*/2) bits per symbol by choosing a sequence from each of the two sets (plus *k* additional bits per symbol encoded in the modulation of *di*).