Model-Based Adaptive Iterative Hard Thresholding Compressive Sensing in Sensor Network for Volcanic Earthquake Detection

Abstract

Recent years have witnessed pilot deployments of inexpensive wireless sensor networks (WSNs) for volcanic eruption detection, where the volcano-seismic signals were collected and processed by sensor nodes. However, it is faced with the limitation of energy resources and the transmission bottleneck of sensors in WSN. In this paper, a Model-Based Adaptive Iterative Hard Thresholding (MAIHT) compressive sensing scheme is developed, where a large number of inexpensive sensors are used to collect fine-grained, real-time volcano-seismic signals while a small number of powerful coordinator nodes process and pick arrival times of primary waves (i.e., P-phases). The paper contribution is two-fold. Firstly, a sparse measurement matrix with theoretical analysis of its restricted isometry property (RIP) is designed to simplify the acquisition process, thereby reducing required storage space and computational demands in sensors. Secondly, a compressive sensing reconstruction algorithm with theoretical analysis of its error bound is presented. Experimental results based on real volcano-seismic data collected from a volcano show that our method can recover the original seismic signal and achieve accurate P-phase picking based on the reconstructed seismic signal.

1. Introduction

Recent years have witnessed pilot deployments of inexpensive wireless sensor networks (WSNs) for active volcano monitoring [1–3]. In such applications, the sensor nodes with limited computation resources were used to collect seismic signals, while a few coordinators process and detect the arrival time of primary waves (i.e., P-phase). These researches mostly focus on improving system robustness, time synchronization, network efficiency, and communication performance issue. However, the energy of each node (generally depending on battery power) and transmission bandwidth are limited in a large-scale wireless sensor network as the sensor signals are sampled at high frequency; sometimes, it is up to 200 Hz. Compressive sensing (CS) theory demonstrates that it is possible to reduce the amount of sampling and reconstruct the original signal accurately with high probability [4–9].

There are two challenges in integrating CS into volcano application for active volcano monitoring. Firstly, in sensor node, it is impossible to generate a random Gaussian matrix due to its limited storage capacity. Meanwhile, the multiplication operation of measurement matrix and the original signal is too complex to be done in the low-end sensor node, especially when the measurement matrix is nonsparse. Mamaghanian et al. proposed a sparse binary matrix with alternative restricted isometry property (RIP) which reduces the operation time for electrocardiogram processing [10, 11]. However, the size of the measurement matrix is over the storage capacity limit of a sensor, such as the popular TelosB platform, which was developed and published to the research community by UC Berkeley. TelosB platform integrates TI MSP430 as the main processing unit and Chipcon CC2420 as radio transceiver. It has a 10 KB internal RAM and a 48 KB program flash memory. Zhang et al. proposed a sparse random matrix for real-time compressive tracking [12]. It can be used for real-time object tracking. However, the parameter s for generating sparse random matrix can only be set to 2 or 3 when considering the RIP of measurement matrix. That is to say, this type of matrix with $s=2$ or 3 satisfies restricted isometry property. Secondly, although a large amount of CS reconstruction algorithms—such as Gradient Projection for Sparse Reconstruction (GRSR) method [13], Iterative Hard Thresholding (IHT) method [14], Matching Pursuit (MP) algorithm [15], Orthogonal Matching Pursuit (OMP) algorithm [16], and Compression Sampling Matching Pursuit (CoSaMP) method [17]—has been proposed, fast convergence CS algorithm with lower reconstruction error bound is still an open issue. The wavelet tree model in [18] is fast. However its accuracy of reconstruction is not ideal. Previous studies have shown that condensing sort and select algorithm (CSSA) [19, 20] can reduce time complexity by modeling the sparse signal into a tree model. Compared with IHT, CoSaMP has more strict requirements for RIP and less accuracy. Unfortunately, the fixed step of CoSaMP bounds its range of application and speed of convergence.

In this paper, we summarize the deficiencies and propose solutions specific to the problems mentioned above. In our preliminary research [21], the measurement matrix and CS reconstruction were simply integrated into volcanic earthquake detection without details. It demonstrated the probability of applying CS in sensors. However, due to the length limitation, further discussion is not presented. Based on our former study [21], here, we give the theoretical analysis of the RIP of the sparse measurement matrix. The MAIHT scheme is developed, which has the following advantages: (1) combining IHT and wavelet tree model algorithms to obtain higher reconstruction accuracy and (2) using an adaptive parameter μ to speed up the convergence. The error bound of the CS reconstruction algorithm was deduced by theoretical analysis. Extensive simulations based on real data collected from a volcano were employed to examine the effectiveness of the proposed scheme.

This paper is organized as follows. Section 2 describes the proposed matrix and presents the analysis of its RIP. In Section 3 the recovery algorithm and the analysis of its error bound are described. Experiments are conducted in Section 4 and conclusions in Section 5.

2. Sparse Measurement Matrix

(A) Sparse Sensing. Element ${\theta }_i$ of the M-dimension measurement matrix Ф is designed as follows:

$$\begin{array}{cc}& {\theta }_i=\sqrt{\frac{\alpha }{M}}\times \left\{\begin{array}{cc}\mathrm{1}\hfill & \text{with\hspace{0.17em}\hspace{0.17em}probability\hspace{0.17em}\hspace{0.17em}}\frac{\mathrm{1}}{\mathrm{2}\alpha }\hfill \\ \mathrm{0}\hfill & \text{with\hspace{0.17em}\hspace{0.17em}probability\hspace{0.17em}\hspace{0.17em}}\mathrm{1}-\frac{\mathrm{1}}{\alpha }\hfill \\ -\mathrm{1}\hfill & \text{with\hspace{0.17em}\hspace{0.17em}probability\hspace{0.17em}\hspace{0.17em}}\frac{\mathrm{1}}{\mathrm{2}\alpha }.\hfill \end{array}\right.\end{array}$$

(1)

It has zero-mean value, $E\left[{\theta }_i^{\mathrm{2}}\right]=\mathrm{1}/M$ , and $\left|{\theta }_i\right|\le \sqrt{\alpha /M}$ .

Obviously, when the value of α increases, the number of “0” elements increases, and the matrix Ф becomes sparser. The key steps of generating sparse random matrix can be summarized as follows:

• (1)

Firstly, generating two “coordinate” matrices randomly. Each random element is the positive integer in interval $[1,M\times N]$ and not duplicate, which means the probability of “1” and “ $-1$ ” in whole elements of Ф is $\mathrm{1}/\mathrm{2}\alpha$ , respectively.

(2)

Multiplying the measurement matrix and the sparse vector and finding out the elements’ positions in sparse vector x, which is corresponding to the line j of “coordinate” matrix. The “1” elements will conduct addition operation and the “ $-1$ ” will conduct subtraction. Finally, multiply the matrix with a constant $\sqrt{\alpha /M}$ .

This matrix is very easy to compute by a uniform random generator.

(B) RIP Analysis of the Proposed Matrix. As described in [22], currently there is no known way to test in polynomial time whether a given matrix satisfies RIP or not. Certain probabilistic constructions of matrices can be shown to satisfy RIP with high probability. Particularly, the matrix Ф is considered to satisfy the RIP ( $K,\epsilon$ ) if for order K and parameter ε there is

$$\begin{array}{cc}& \left(\mathrm{1}-\epsilon \right){‖x‖}_{\mathrm{2}}^{\mathrm{2}}\le {‖\Phi x‖}_{\mathrm{2}}^{\mathrm{2}}\le \left(\mathrm{1}+\epsilon \right){‖x‖}_{\mathrm{2}}^{\mathrm{2}},\end{array}$$

(2)

where x is a K-sparse vector and $\delta >\mathrm{0}$ . In other words, the RIP ( $K,\epsilon$ ) guideline, for a matrix Ф, is defined as the singular values of all submatrices satisfying in the range $\left(\sqrt{\mathrm{1}-\epsilon },\sqrt{\mathrm{1}+\epsilon }\right)$ [22]. We follow this idea and show that the following theorem stands.

Theorem 1.

Let parameter $\epsilon >0$ be a universal constant. Given a sparsity level K and the length of signal N, for any $c<{\epsilon }^{\mathrm{2}}/\mathrm{8}$ , if the number of measurements satisfies

$$\begin{array}{cc}& M\ge \frac{\mathrm{24}{\alpha }^{\mathrm{2}}}{{\epsilon }^{\mathrm{2}}-\mathrm{8}c}{K}^{\mathrm{2}}\log N,\end{array}$$

(3)

where $\alpha >\mathrm{0}$ is defined in (1), the RIP in (2) holds with probability at least $\mathrm{1}-\mathrm{e}\mathrm{x}\mathrm{p}\left(-cM/{\alpha }^{\mathrm{2}}{K}^{\mathrm{2}}\right)$ .

Proof.

Before the proof, three lemmas will be useful.

Lemma 2 (Gershgorin disc theorem [23]).

Let $A\in C^{n\times n}$ , $A=\left(a_{ij}\right)$ , so the n eigenvalues will fall into the union $\left(\cup D_i\left(A\right)\right)$ of n disks

$$\begin{array}{cc}& D_i\left(A\right)=\left\{z:\left|z-a_{ij}\right|\le R_i\right\},i=\mathrm{1,2},\dots ,n,\end{array}$$

(4)

and $R_i={\sum }_{j=\mathrm{1},j\ne i}^n\left|a_{ij}\right|$ .

Lemma 3 (see [22]).

Let $x_i,\hspace{0.17em}\hspace{0.17em}i=\mathrm{1,2},\dots ,k$ , be a sequence of i.i.d. zero-mean bounded random variables such that $\left|x_i\right|\le a$ , and with variance $E\left[x_i^{\mathrm{2}}\right]={\sigma }^{\mathrm{2}}$ , then,

$$\begin{array}{cc}& P\left(\left|{\displaystyle \sum }_{i=\mathrm{1}}^k{x}_i^{\mathrm{2}}-k{\sigma }^{\mathrm{2}}\right|\ge t\right)\le \mathrm{2}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\mathrm{2}{t}^{\mathrm{2}}}{k{a}^{\mathrm{4}}}\right).\end{array}$$

(5)

Lemma 4 (see [22]).

Let $x_i,\hspace{0.17em}\hspace{0.17em}i=\mathrm{1,2},\dots ,k$ , and $y_i,\hspace{0.17em}\hspace{0.17em}i=\mathrm{1,2},\dots ,k$ , be sequences of i.i.d. zero-mean bounded random variables satisfying $\left|x_i\right|\le a$ , $\left|x_i{y}_i\right|\le a^{\mathrm{2}}$ ; then,

$$\begin{array}{cc}& P\left(\left|{\displaystyle \sum }_{i=\mathrm{1}}^k{x}_i{y}_i\right|\ge t\right)\le \mathrm{2}\exp \left(-\frac{t^{\mathrm{2}}}{\mathrm{2}k{a}^{\mathrm{4}}}\right).\end{array}$$

(6)

According to the results in [22], there is a Gramian matrix $G={\Phi }^T\Phi$ in which the nonnegative square roots of its eigenvalues are the singular values of Ф, so that all eigenvalues of the matrix satisfied RIP will be in the range $\left(\mathrm{1}-\epsilon ,\mathrm{1}+\epsilon \right)$ . Bounding up with Lemma 2, let $\epsilon ={\epsilon }_{\mathrm{1}}+{\epsilon }_{\mathrm{2}}$ , ${\epsilon }_{\mathrm{1}}>\mathrm{0}$ , ${\epsilon }_{\mathrm{2}}>\mathrm{0}$ . The distance between center of Gershgorin disk and 1 is less than ${\epsilon }_{\mathrm{1}}$ , and radius of all circles is less than ${\epsilon }_{\mathrm{2}}$ . Based on the above descriptions, the following result can be obtained that if matrix Ф satisfies RIP, the diagonal element $g_{ii}$ of its Gramian matrix $G={\Phi }^T\Phi$ satisfies $\left|g_{ii}-1\right|<{\epsilon }_1$ , and off-diagonal elements satisfy $\left|g_{ij}\right|<{\epsilon }_2/K$ .

In the following proof, we will show that the proposed matrix will satisfy RIP with high probability. Applying Lemma 3, each diagonal element $g_{ii}={\sum }_{j=\mathrm{1}}^M{\theta }_j^{\mathrm{2}}$ of Gramian matrix $G={\Phi }^T\Phi$ satisfies

$$\begin{array}{cc}& P\left(\left|g_{ii}-\mathrm{1}\right|\ge {\epsilon }_{\mathrm{1}}\right)\le \mathrm{2}\exp \left(-\frac{\mathrm{2}M{\epsilon }_{\mathrm{1}}^{\mathrm{2}}}{{\alpha }^{\mathrm{2}}}\right),\\ P\left({\displaystyle \bigcup }_{i=\mathrm{1}}^N\left\{\left|g_{ii}-\mathrm{1}\right|\ge {\epsilon }_{\mathrm{1}}\right\}\right)\le \mathrm{2}N\exp \left(-\frac{\mathrm{2}M{\epsilon }_{\mathrm{1}}^{\mathrm{2}}}{{\alpha }^{\mathrm{2}}}\right).\end{array}$$

(7)

For the off-diagonal element given by $g_{ij}$ , such as $g_{12}={\theta }_1{\theta }_2+{\theta }_{n+1}{\theta }_{n+2}+{\theta }_{2n+1}{\theta }_{2n+2}+\cdots$ , we recognize that the entries of the sum are independent. From Lemma 4, we can get

$$\begin{array}{cc}& P\left(\left|g_{ij}\right|\ge \frac{{\epsilon }_{\mathrm{2}}}{K}\right)\le \mathrm{2}\exp \left(-\frac{M{\epsilon }_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{2}{\alpha }^{\mathrm{2}}{K}^{\mathrm{2}}}\right).\end{array}$$

(8)

Due to the fact that $g_{ij}=g_{ji}$ , the number of off-diagonal elements $g_{ij}$ is $N^{\mathrm{2}}-N<N^{\mathrm{2}}$ . We have

$$\begin{array}{cc}& P\left({\displaystyle \bigcup }_{i=1}^N{\displaystyle \bigcup }_{\begin{array}{c}j=1\\ j\ne i\end{array}}^N\left|g_{ij}\right|\ge \frac{{\epsilon }_{\mathrm{2}}}{K}\right)\le \mathrm{2}{N}^{\mathrm{2}}\exp \left(-\frac{M{\epsilon }_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{2}{\alpha }^{\mathrm{2}}{K}^{\mathrm{2}}}\right).\end{array}$$

(9)

Let $T=\{{\Phi }_{M\times N}$ , which satisfies RIP $(K,\epsilon )\}$ , and $\stackrel{-}{T}$ is the complement of T; we have

$$\begin{array}{cc}& P\left(\stackrel{-}{T}\right)\le \mathrm{2}N\exp \left(-\frac{\mathrm{2}M{\epsilon }_{\mathrm{1}}^{\mathrm{2}}}{{\alpha }^{\mathrm{2}}}\right)+\mathrm{2}{N}^{\mathrm{2}}\exp \left(-\frac{M{\epsilon }_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{2}{\alpha }^{\mathrm{2}}{K}^{\mathrm{2}}}\right).\end{array}$$

(10)

Assuming $N\ge \mathrm{3}$ and ${\epsilon }_{\mathrm{1}}={\epsilon }_{\mathrm{2}}=\epsilon /\mathrm{2}$ , we obtain

$$\begin{array}{l}P\left(\stackrel{-}{T}\right)\le \mathrm{2}N\exp \left(-\frac{M{\epsilon }^{\mathrm{2}}}{\mathrm{2}{\alpha }^{\mathrm{2}}}\right)+\mathrm{2}{N}^{\mathrm{2}}\exp \left(-\frac{M{\epsilon }^{\mathrm{2}}}{\mathrm{8}{\alpha }^{\mathrm{2}}{K}^{\mathrm{2}}}\right)\hfill \\ \le \mathrm{3}{N}^{\mathrm{2}}\exp \left(-\frac{M{\epsilon }^{\mathrm{2}}}{\mathrm{8}{\alpha }^{\mathrm{2}}{K}^{\mathrm{2}}}\right).\hfill \end{array}$$

(11)

Then we have $c<{\epsilon }^{\mathrm{2}}/\mathrm{8}$ , $M\ge \left(\mathrm{24}{\alpha }^{\mathrm{2}}/\left({\epsilon }^{\mathrm{2}}-\mathrm{8}c\right)\right)K^{\mathrm{2}}\log N$ , and

$$\begin{array}{cc}& P\left(\stackrel{-}{T}\right)\le \exp \left(-\frac{cM}{{\alpha }^{\mathrm{2}}{K}^{\mathrm{2}}}\right).\end{array}$$

(12)

So for any $c<{\epsilon }^{\mathrm{2}}/\mathrm{8}$ and $M\ge \left(\mathrm{24}{\alpha }^{\mathrm{2}}/\left({\epsilon }^{\mathrm{2}}-\mathrm{8}c\right)\right)K^{\mathrm{2}}\log N$ , the measurement matrix ${\Phi }_{M\times N}$ satisfies RIP ( $K,\epsilon$ ) with probability at least $\mathrm{1}-\mathrm{e}\mathrm{x}\mathrm{p}\left(-cM/{\alpha }^{\mathrm{2}}{K}^{\mathrm{2}}\right)$ .

So far, we recognize that the RIP of measurement matrix Ф is verified by the above proof.

3. Model-Based Adaptive IHT

(A) Overview of the Algorithm. The proposed algorithm is summarized as follows:

• (1)

Initialization: $x_{\mathrm{1}}=\mathrm{0}$ , residual $\mathit{r}=\mathit{y}$ , iterations $n=0$ ;

 

$x_n$ is the result of original signal which converted by condensing sort and select algorithm (CSSA) with parameter m;

 

Supports $\Lambda =Ø$ , $J=Ø$ ;

(2)

while halting criterion false do

 

$n=n+1$ ;

 

$x_{n+\mathrm{1}}=x_n+\mu {\Phi }^{T}r$ , $\mu ={‖x_k-x_{k-1}‖}_2^2/{‖\Phi \left(x_k-x_{k-1}\right)‖}_2^2$ ;

 

$J\leftarrow \mathfrak{M}\hspace{0.17em}\hspace{0.17em}\left(x_{n+\mathrm{1}},\mathrm{2}m\right)$ (getting support according to signal model $\mathfrak{M}$ , 2m maximum indexes);

 

$T=\sup p\left(x_n\right)\cup J$ (merging supports);

 

${\left.b\right|}_T={\Phi }_T^{†}y$ (signal estimating);

 

$x_{n+\mathrm{1}}\leftarrow \mathfrak{M}\hspace{0.17em}\hspace{0.17em}\left(b,m\right)$ (pruning signal estimate b with m maximum indexes);

 

${\mathbf{r}}_{\mathrm{n}\mathrm{e}\mathrm{w}}=\mathbf{y}-\mathbf{\Phi }{x}_{n+\mathrm{1}}$ (updating measurement residual);

(3)

End while

  

return $x=x_n$ .

(B) Error Bound Analysis of the Algorithm. Starting with the basic principle of IHT, the iteration is

$$\begin{array}{cc}& x_{n+\mathrm{1}}=H_S\left(x_n+{\Phi }^T\left(y-\Phi x_n\right)\right),\end{array}$$

(13)

with initial conditions $x_{\mathrm{0}}=\mathrm{0}$ , where $H_S\left(\alpha \right)$ is a nonlinear operator. Based on a predefined order, it remains the maximum amplitude part s in α and sets the rest to be zero. IHT has good convergence performance [14]. Literature [24] shows that the iterative formula of IHT method can be rewritten as

$$\begin{array}{cc}& x_{n+\mathrm{1}}=H_S\left(x_n+\mu {\Phi }^T\left(y-\Phi x_n\right)\right),\end{array}$$

(14)

where the step $\mu =g_{{\Gamma }^n}^T{g}_{{\Gamma }^n}/g_{{\Gamma }^n}^T{\Phi }_{{\Gamma }^n}^T{\Phi }_{{\Gamma }^n}{g}_{{\Gamma }^n}$ and $\mathit{g}={\Phi }^T(\mathit{y}-\Phi \mathit{x})$ . This parameter μ can improve the rate of convergence. However, it increases the complexity of calculation. Furthermore, the step μ is derived from steepest descent method, which is sensitive to ill-conditioning. In order to reduce the complexity and maintain the good reliability, we use two-point step size gradient method [25] to increase the rate of convergence and give the error bound. In the following, we will analytically demonstrate the CS reconstruction error bound.

Theorem 5.

Given a noisy measurement $y=\Phi x^s+e$ and $x^s$ is a s-sparse vector, measurement matrix Ф satisfies RIP principle with $\delta <\mathrm{0.15}$ ; then at iteration k, Model-Based Adaptive IHT will recover an approximation $x^k$ satisfying

$$\begin{array}{cc}& {‖x^s-x_k‖}_2\le 2^{-k}{‖x^s‖}_2+5.04{‖e‖}_2.\end{array}$$

(15)

Proof.

From the following equations

$$\begin{array}{cc}& x_{n+1}=x_n-\mu g,\\ \mu =\arg \min {‖\Delta x-\mu \Delta g‖}^{\mathrm{2}},\end{array}$$

(16)

we have

$$\begin{array}{cc}& \mu =\frac{\Delta x^T\Delta x}{\Delta x^T\Delta g},\end{array}$$

(17)

where

$$\begin{array}{cc}& \Delta x=x_{n+1}-x_n,\\ \Delta g=g_{n+1}-g_n.\end{array}$$

(18)

In the proposed method, (17) was applied to (14). It provides faster convergence rate and simple calculation of μ. Given a noisy observation $y=\Phi x^s+e$ and $x^s$ is a s-sparse vector. The matrix Ф satisfies the RIP of order s, if there exists ${\delta }_s<\mathrm{1}$ satisfying

$$\begin{array}{cc}& \left(\mathrm{1}-{\delta }_s\right){‖x^s‖}_{\mathrm{2}}^{\mathrm{2}}\le {‖\Phi x^s‖}_{\mathrm{2}}^{\mathrm{2}}\le \left(\mathrm{1}+{\delta }_s\right){‖x^s‖}_{\mathrm{2}}^{\mathrm{2}}.\end{array}$$

(19)

At the iteration $k-1$ and k, the recovery approximations in $x_{k-1}$ and $x_k$ have the same properties with $x^s$ . They are s-sparse in the same index. We can infer that $x_k-x_{k-1}$ is s-sparse and applicable to the RIP of order s. Therefore we have

$$\begin{array}{l}\left(1-{\delta }_s\right){‖x_k-x_{k-1}‖}_2^2\le {‖\Phi \left(x_k-x_{k-1}\right)‖}_2^2\hfill \\ \le \left(1+{\delta }_s\right){‖x_k-x_{k-1}‖}_2^2,\hfill \\ 1-{\delta }_s\le \frac{{‖\Phi \left(x_k-x_{k-1}\right)‖}_2^2}{{‖x_k-x_{k-1}‖}_2^2}\le 1+{\delta }_s.\hfill \end{array}$$

(20)

Based on the equation $x_{n+1}=x_n-\mu g$ , we recognize that the formula at iteration $n+1$ is $x_{n+1}=x_n+\mu {\Phi }^T(y-\Phi x_n)$ . Let $g={\Phi }^T\left(\Phi x-y\right)$ be the negative gradient of ${‖y-\Phi x‖}_{\mathrm{2}}^{\mathrm{2}}$ . Substituting g into the term of μ,

$$\begin{array}{l}\mu =\frac{\Delta x^T\Delta x}{\Delta x^T\Delta g}\\ =\frac{{\left(x_k-x_{k-\mathrm{1}}\right)}^T\left(x_k-x_{k-\mathrm{1}}\right)}{{\left(x_k-x_{k-\mathrm{1}}\right)}^T\left({\Phi }^T\Phi x_k-{\Phi }^{T}y-{\Phi }^T\Phi x_{k-\mathrm{1}}+{\Phi }^{T}y\right)}\\ =\frac{{\left(x_k-x_{k-\mathrm{1}}\right)}^T\left(x_k-x_{k-\mathrm{1}}\right)}{{\left(x_k-x_{k-\mathrm{1}}\right)}^T{\Phi }^T\Phi \left(x_k-x_{k-\mathrm{1}}\right)}\\ =\frac{{‖x_k-x_{k-\mathrm{1}}‖}_{\mathrm{2}}^{\mathrm{2}}}{{‖\Phi \left(x_k-x_{k-\mathrm{1}}\right)‖}_{\mathrm{2}}^{\mathrm{2}}},\end{array}$$

(21)

we obtain

$$\begin{array}{cc}& \mu =\frac{{‖x_k-x_{k-1}‖}_2^2}{{‖\Phi \left(x_k-x_{k-1}\right)‖}_2^2}.\end{array}$$

(22)

Replacing μ in (20), we have

$$\begin{array}{cc}& \frac{\mathrm{1}}{\mathrm{1}+\delta }\le \mu \le \frac{\mathrm{1}}{\mathrm{1}-\delta }.\end{array}$$

(23)

Firstly, we use the triangle inequality to bound

$$\begin{array}{cc}& {‖x^s-x_{n+1}‖}_2\le {‖x^s-a‖}_2+{‖x_{n+1}-a‖}_2,\end{array}$$

(24)

and we know that $x^{n+1}$ is a better s -term approximation to a than $x^s$ , so that we obtain

$$\begin{array}{cc}& {‖x^s-x_{n+1}‖}_2\le 2{‖x^s-a‖}_2.\end{array}$$

(25)

Then expanding a from the iteration, we have $x_{n+1}=x_n+\mu {\Phi }^T(y-\Phi x_n)$ . We recognize that the following equation stands:

$$\begin{array}{cc}& x_{n+\mathrm{1}}=x_n+\mu {\Phi }^T\left(\Phi \left(x_{n+\mathrm{1}}+e\right)-\Phi x_n\right).\end{array}$$

(26)

Assuming r is residual, we have

$$\begin{array}{cc}& a=x_n+\mu {\Phi }^T\left(\Phi r+e\right)=x_n+\mu {\Phi }^T\Phi r+\mu {\Phi }^{T}e.\end{array}$$

(27)

With these positions, we obtain

$$\begin{array}{l}{\Vert x^s-x_{n+1}\Vert }_2\le 2{\Vert x^s-x_n-u{\Phi }^T\Phi r_n-u{\Phi }^{T}e\Vert }_2\hfill \\ \le 2{\Vert r_{B^{N+1}}^n-u{\Phi }_{B^{n+1}}^T\Phi r_n\Vert }_2+2{\Vert u{\Phi }_{B^{n+1}}^{T}e\Vert }_2\hfill \\ =\Vert \left(I-u{\Phi }_{B^{m+1}}^T{\Phi }_{B^{n+1}}\right)r_{B^{n+1}}^n\hfill \\ -u{\Phi }_{B^{n+1}}^T{\Phi }_{B^n\backslash B^{n+1}}{r}_{B^n\backslash B^{n+1}}^n|{|}_2+2{\Vert u{\Phi }_{B^{n+1}}^{T}e\Vert }_2\hfill \\ \le 2{\Vert \left(I-u{\Phi }_{B^{n+1}}^T{\Phi }_{B^{n+1}}\right)r_{B^{n+1}}^n\Vert }_2\hfill \\ +2{\Vert u{\Phi }_{B^{n+1}}^T{\Phi }_{B^n\backslash B^{n+1}}{r}_{B^n\backslash B^{n+1}}^n\Vert }_2+2{\Vert u{\Phi }_{B^{n+1}}^{T}e\Vert }_2,\hfill \end{array}$$

(28)

where the sets $B^n\setminus B^{n+1}$ and $B^{n+1}$ are incoherent and $\left|B^n\cup B^{n+1}\right|\le 3s$ ( $3s$ refers to the sparsity).

From (19), we have

$$\begin{array}{cc}& {‖{\Phi }_{\Gamma }^{T}y‖}_2\le \sqrt{1+\delta }{‖y‖}_2,\end{array}$$

(29)

$$\begin{array}{cc}& \left(1-{\delta }_{\Gamma }\right){‖x_{\Gamma }‖}_2\le {‖{\Phi }_{\Gamma }^T{\Phi }_{\Gamma }{x}_{\Gamma }‖}_2\le \left(1+{\delta }_{\Gamma }\right){‖x_{\Gamma }‖}_2.\end{array}$$

(30)

Combining the above two equations with (23), we can get

$$\begin{array}{cc}& {‖\left(I-\mu {\Phi }_{\Gamma }^T{\Phi }_{\Gamma }\right)x_{\Gamma }‖}_2\le \frac{{\delta }_{\Gamma }}{1-{\delta }_{\Gamma }}{‖x_{\Gamma }‖}_2.\end{array}$$

(31)

For two disjoint sets Γ and Λ, all Φ for which the RIP holds with $s=\left|\Gamma \cup \Lambda \right|$ ,

$$\begin{array}{cc}& {‖\mu {\Phi }_{\Gamma }^T{\Phi }_{\Lambda }{x}_{\Lambda }‖}_2\le \frac{{\delta }_{\Lambda }}{1-{\delta }_{\Lambda }}{‖x_{\Lambda }‖}_2.\end{array}$$

(32)

Combining (29), (31), and (32), we can get

$$\begin{array}{cc}& {‖{\Phi }_{\Gamma }^{T}y‖}_2\le \sqrt{1+\delta }{‖y‖}_2,\\ {‖\left(I-\mu {\Phi }_{\Gamma }^T{\Phi }_{\Gamma }\right)x_{\Gamma }‖}_2\le \frac{{\delta }_{\Gamma }}{1-{\delta }_{\Gamma }}{‖x_{\Gamma }‖}_2,\\ {‖\mu {\Phi }_{\Gamma }^T{\Phi }_{\Lambda }{x}_{\Lambda }‖}_2\le \frac{{\delta }_{\Lambda }}{1-{\delta }_{\Lambda }}{‖x_{\Lambda }‖}_2.\end{array}$$

(33)

Substituting (33) into (28), we have

$$\begin{array}{l}{‖x^s-x_{n+1}‖}_2\le \frac{2{\delta }_{2s}}{1-{\delta }_{2s}}{‖r_{B^{n+1}}^n‖}_2+\frac{2{\delta }_{3s}}{1-{\delta }_{3s}}{‖r_{B^n\setminus B^{n+1}}^n‖}_2\\ +\frac{2\sqrt{1+\delta }}{1-\delta }{‖e‖}_2.\end{array}$$

(34)

Since ${\delta }_{\mathrm{2}s}$ and ${\delta }_{3s}$ refer to the parameters from RIP of orders $2s$ and $3s$ , respectively, we have the bound ${\delta }_{\mathrm{2}s}\le {\delta }_{\mathrm{3}s}$ .

Furthermore $r_{B^{n+1}}^n$ and $r_{B^n\setminus B^{n+1}}^n$ are orthogonal, so that ${‖r_{B^{n+1}}^n‖}_2+{‖r_{B^n\setminus B^{n+1}}^n‖}_2\le \sqrt{2}{‖r‖}_2$ . According to (34), we have

$$\begin{array}{cc}& {‖r^{n+1}‖}_2\le \frac{2\sqrt{2}{\delta }_{3s}}{1-{\delta }_{3s}}{‖r‖}_2+\frac{2\sqrt{1+\delta }}{1-\delta }{‖e‖}_2.\end{array}$$

(35)

If $\delta <\mathrm{0.15}$ , the following inequality stands:

$$\begin{array}{cc}& {‖r^{n+1}‖}_2\le 0.5{‖r^n‖}_2+2.52{‖e‖}_2.\end{array}$$

(36)

Iterating the relationship (36), and realizing that $\left[\mathrm{1}+\mathrm{1}/\mathrm{2}+\cdots +{\left(\mathrm{1}/\mathrm{2}\right)}^k\right]\times \mathrm{2.52}\le \mathrm{5.04}$ with the initial condition $x_0=\mathrm{0}$ , finally, we obtain the error bound described by Theorem 5

$$\begin{array}{cc}& {‖r_k‖}_2\le 2^{-k}{‖x^s‖}_2+5.04{‖e‖}_2.\end{array}$$

(37)

4. Experimental Results and Discussion

(A) Experiments on Sensor Nodes. In this part, we conduct testbed experiments and extensive simulations based on real data traces collected by 12 nodes on Mount St. Helens in the OASIS project [3]. Our system implementation is based on TelosB motes. Our current implementation of CS uses predefined binary random matrices for sensors. We use a laptop computer to simulate the coordinator.

We evaluate the computation and storage overhead of the algorithms running on low-end sensors in a testbed of 12 TelosB motes. The min execution time and max execution time of CS are 0.20 s and 0.88 s, respectively, while the RAM and ROM requirement of CS in TelosB are 3 and 22 kb, respectively. We can see that the accumulated execution time does not exceed 1 second, which introduces tiny workload and energy consumption to the sensors. The variation of execution time is mainly caused by sparsity of the seismic signal and CS. As the seismic signals at sensors have different sparsity, sensors have different number of rows in the project matrix A. This leads to the variable execution time of CS.

Due to the high computation and storage overhead requirement, random Gaussian matrix and sparse binary matrix [16] cannot be run in TelosB. We simulate random Gaussian matrix and sparse binary matrix [16] in a laptop with an Intel i3-M480 CPU and 4 GB RAM for comparison. We test the simulation with the data length $N=1024$ , measurements $M=4K$ (K represents the sparse of signal), $D=4$ , and $\alpha =M/\mathrm{4}$ for 1000 times. The results are shown in Table 1, where the space complexity refers to the storage space. The processing time is the interval from the generation of measurement matrix to the generation of measurement vector. The averages results over the 1000 times were used for further comparisons. From Table 1, it can be seen that the storage requirement of random Gaussian matrix and sparse binary matrix outstrips the capacity of TelosB.

Table 1

Comparison of complexity of random Gaussian matrix and sparse binary matrix.

	Space complexity (bytes)	Processing time (s)
Random Gaussian matrix	2990080	0.0142
Sparse binary matrix	2990080	0.0089

(B) Simulation of Seismic Signal Reconstruction Algorithm. To quantitatively evaluate the reconstruction performance, we compare the proposed method with common algorithms in CS area, that is, Orthogonal Matching Pursuit (OMP), Compressive Sampling Matching Pursuit (CoSaMP), Model-Based Recovery, Iterative Hard Thresholding (IHT), and Normalized Iterative Hard Thresholding (μ-IHT). We quantified an algorithm's performance in terms of two quality indexes: root mean square error (RMSE) and average execution time. The simulation data is the real monitoring data of the St. Helens volcano in USA, which has an obvious vibration in case of the seismic wave arrival. In the practical implementation, a large number of iterations were set (e.g., 300 in this paper). The RMSE and time index are temporarily stored. Final result is the mean value of them.

First of all, we analyze the impact of D (the number of nonzeros in each column) on RMSE, where the length of test random data is $N=1024$ , M is the number of measurements, K is sparsity of signal, and D is the number of nonzeros in each column. The result is shown in Figure 1. It is observed that the greater the value of M, the smaller the error. From the three lines of 8, 6, and 4 of D, it can be seen that efficient reconstruction is guaranteed when M is small. On the contrary, $D=2$ is too low to achieve an accurate reconstruction. As we know the larger the value of D, the sparser the measurement matrix. $D=4$ is used in this paper.

Figure 1

Validation of D.

Figure 2 reflects the effectiveness of the proposed method. From the result, it can be seen that although the original volcano signal has a lot of vibration, the proposed method can reconstruct the original signal. Figure 3 displays the reconstruction results for four different traces in volcano dataset, respectively. It is obvious to see that the proposed method has a better reconstruction result.

Figure 2

Reconstructed results of real monitoring data with our method between original signal (red bold solid) and reconstructed result (blue solid). The length of test signal is 1600 ( $N=1600$ ).

Figure 3

Reconstructed results of four different monitoring data with our method between the original signal (red bold solid) and reconstructed result (blue solid). The length of test signal is 640 ( $N=640$ ).

Next, the analysis of computational complexity and performance comparison in terms of RMSE of the following four methods, CoSaMp, IHT, Model-Based Recovery, and our method, is shown in Table 2. It shows that our method is faster than the others. Our method is about two times faster than Model-Based method. When compared with IHT method, our method is about 50 times faster than it. The accuracy of each method is shown in Figures 4 and 5. It indicates our method is better.

Table 2

Comparison of computational complexity performance of recovery.

Time (seconds)	CoSaMP	IHT	Model-Based	Our method
Convergence time	0.043	$$\mathit{1.92}\ast e-\mathit{3}$$	$$\mathit{2.86}\ast e-\mathit{4}$$	$$\mathit{1.65}\ast e-\mathit{4}$$
Total time	1.984	0.014	$$\mathit{6.82}\ast e-\mathit{4}$$	$$\mathit{2.58}\ast e-\mathit{4}$$

Figure 4

Reconstructed results between the original signal (red bold solid) and reconstructed result (blue solid) for five different recovery algorithms using the same data. The length of test signal is 640 ( $N=640$ ).

Figure 5

Comparison of different algorithms. Normalised IHT (solid), CoSaMP (dashed), Model-Based Recovery (dash-dotted), Model-Based μ-IHT (solid-ring), and proposed method (solid-asterisk).

Figure 4 shows the recovery results for CoSaMP, IHT, μ-IHT, Model-Based method, and our method, respectively. It can be seen that CoSaMP cannot reconstruct the small variation in the signal, while IHT and μ-IHT can recover most of the variation. Model-Based method and the proposed method have higher accuracy in recovery than the other three algorithms.

The effectiveness of the proposed method is confirmed by the plots of the RMSE shown in Figure 5. Simulation results show that the proposed algorithm has a nice performance, and the RMSE are improved in the same precondition. With the increase of measurement numbers, our method can reach the minimal RMSE firstly.

(C) Accuracy Impact on the P-Phase Picking. In this section, we use the reconstruction trace to pick the arrival time of the seismic signal to verify the earthquake detection. We evaluate the accuracy of the P-phase picked by the AR-AIC picker [26] on the original seismic signal and reconstructed seismic signal, respectively. Figure 6 shows the original and reconstructed signals received by a sensor in one volcano earthquake. We can see from the figure that the signals can be reconstructed and the P-phase can be well preserved.

Figure 6

Original and reconstructed seismic signals with picks. Vertical dashed/solid lines represent the picks by AR-AIC algorithm on the original/reconstructed seismic signals.

5. Conclusion and Perspectives

In this paper, a CS-based system consists of a sparse measurement matrix and a new Model-Based Adaptive Iterative Hard Thresholding recovery algorithm was proposed. We integrate the CS scheme into the volcano detection application. Theoretical proof was presented to specify the performance of the method. Simulation results show that the measurement matrix can satisfy RIP and the reconstruction algorithm has a good convergence. Testbed experiments and extensive simulations based on real data collected on a volcano show that our method can achieve P-phase picking based on the reconstructed seismic signal.

In the future, questions concerning more effective reconstruction algorithms and sparse binary matrices will be studied.