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 $θ_{i}$ of the M-dimension measurement matrix Ф is designed as follows:

\begin{matrix} θ_{i} = \sqrt{\frac{α}{M}} \times \{\begin{cases} 1 & with  probability \frac{1}{2 α} \\ 0 & with  probability 1 - \frac{1}{α} \\ - 1 & with  probability \frac{1}{2 α} . \end{cases} \end{matrix}

(1)

It has zero-mean value,

E [θ_{i}^{2}] = 1 / M

, and

|θ_{i}| \leq \sqrt{α / 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 $1 / 2 α$ , 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{α / 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, ε$ ) if for order K and parameter ε there is

\begin{matrix} (1 - ε) {‖x‖}_{2}^{2} \leq {‖Φ x‖}_{2}^{2} \leq (1 + ε) {‖x‖}_{2}^{2}, \end{matrix}

(2)

where x is a K-sparse vector and

δ > 0

. In other words, the RIP (

K, ε

) guideline, for a matrix Ф, is defined as the singular values of all submatrices satisfying in the range

(\sqrt{1 - ε}, \sqrt{1 + ε})

[22]. We follow this idea and show that the following theorem stands.

Theorem 1.

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

\begin{matrix} M \geq \frac{24 α^{2}}{ε^{2} - 8 c} K^{2} \log N, \end{matrix}

(3)

where

α > 0

is defined in (1), the RIP in (2) holds with probability at least

1 - e x p (- c M / α^{2} K^{2})

Proof.

Before the proof, three lemmas will be useful.

Lemma 2 (Gershgorin disc theorem [23]).

Let $A \in C^{n \times n}$ , $A = (a_{i j})$ , so the n eigenvalues will fall into the union $(\cup D_{i} (A))$ of n disks

\begin{matrix} D_{i} (A) = \{z : |z - a_{i j}| \leq R_{i}\}, i = 1,2, \dots, n, \end{matrix}

(4)

and

R_{i} = \sum_{j = 1, j \neq i}^{n} |a_{i j}|

Lemma 3 (see [22]).

Let $x_{i}, i = 1,2, \dots, k$ , be a sequence of i.i.d. zero-mean bounded random variables such that $|x_{i}| \leq a$ , and with variance $E [x_{i}^{2}] = σ^{2}$ , then,

\begin{matrix} P (|\sum_{i = 1}^{k} x_{i}^{2} - k σ^{2}| \geq t) \leq 2 e x p (- \frac{2 t^{2}}{k a^{4}}) . \end{matrix}

(5)

Lemma 4 (see [22]).

Let $x_{i}, i = 1,2, \dots, k$ , and $y_{i}, i = 1,2, \dots, k$ , be sequences of i.i.d. zero-mean bounded random variables satisfying $|x_{i}| \leq a$ , $|x_{i} y_{i}| \leq a^{2}$ ; then,

\begin{matrix} P (|\sum_{i = 1}^{k} x_{i} y_{i}| \geq t) \leq 2 \exp (- \frac{t^{2}}{2 k a^{4}}) . \end{matrix}

(6)

According to the results in [22], there is a Gramian matrix $G = Φ^{T} Φ$ 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 $(1 - ε, 1 + ε)$ . Bounding up with Lemma 2, let $ε = ε_{1} + ε_{2}$ , $ε_{1} > 0$ , $ε_{2} > 0$ . The distance between center of Gershgorin disk and 1 is less than $ε_{1}$ , and radius of all circles is less than $ε_{2}$ . Based on the above descriptions, the following result can be obtained that if matrix Ф satisfies RIP, the diagonal element $g_{i i}$ of its Gramian matrix $G = Φ^{T} Φ$ satisfies $|g_{i i} - 1| < ε_{1}$ , and off-diagonal elements satisfy $|g_{i j}| < ε_{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_{i i} = \sum_{j = 1}^{M} θ_{j}^{2}$ of Gramian matrix $G = Φ^{T} Φ$ satisfies

\begin{matrix} P (|g_{i i} - 1| \geq ε_{1}) \leq 2 \exp (- \frac{2 M ε_{1}^{2}}{α^{2}}), \\ P (⋃_{i = 1}^{N} \{|g_{i i} - 1| \geq ε_{1}\}) \leq 2 N \exp (- \frac{2 M ε_{1}^{2}}{α^{2}}) . \end{matrix}

(7)

For the off-diagonal element given by

g_{i j}

, such as

g_{12} = θ_{1} θ_{2} + θ_{n + 1} θ_{n + 2} + θ_{2 n + 1} θ_{2 n + 2} + \dots

, we recognize that the entries of the sum are independent. From Lemma 4, we can get

\begin{matrix} P (|g_{i j}| \geq \frac{ε_{2}}{K}) \leq 2 \exp (- \frac{M ε_{2}^{2}}{2 α^{2} K^{2}}) . \end{matrix}

(8)

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

\begin{matrix} P (⋃_{i = 1}^{N} ⋃_{\begin{matrix} j = 1 \\ j \neq i \end{matrix}}^{N} |g_{i j}| \geq \frac{ε_{2}}{K}) \leq 2 N^{2} \exp (- \frac{M ε_{2}^{2}}{2 α^{2} K^{2}}) . \end{matrix}

(9)

Let $T = {Φ_{M \times N}$ , which satisfies RIP $(K, ε)}$ , and $\bar{T}$ is the complement of T; we have

\begin{matrix} P (\bar{T}) \leq 2 N \exp (- \frac{2 M ε_{1}^{2}}{α^{2}}) + 2 N^{2} \exp (- \frac{M ε_{2}^{2}}{2 α^{2} K^{2}}) . \end{matrix}

(10)

Assuming $N \geq 3$ and $ε_{1} = ε_{2} = ε / 2$ , we obtain

\begin{array}{l} P (\bar{T}) \leq 2 N \exp (- \frac{M ε^{2}}{2 α^{2}}) + 2 N^{2} \exp (- \frac{M ε^{2}}{8 α^{2} K^{2}}) \\ \leq 3 N^{2} \exp (- \frac{M ε^{2}}{8 α^{2} K^{2}}) . \end{array}

(11)

Then we have $c < ε^{2} / 8$ , $M \geq (24 α^{2} / (ε^{2} - 8 c)) K^{2} \log N$ , and

\begin{matrix} P (\bar{T}) \leq \exp (- \frac{c M}{α^{2} K^{2}}) . \end{matrix}

(12)

So for any $c < ε^{2} / 8$ and $M \geq (24 α^{2} / (ε^{2} - 8 c)) K^{2} \log N$ , the measurement matrix $Φ_{M \times N}$ satisfies RIP ( $K, ε$ ) with probability at least $1 - e x p (- c M / α^{2} K^{2})$ .

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_{1} = 0$ , residual $r = 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 $Λ = Ø$ , $J = Ø$ ;

(2)

while halting criterion false do

$n = n + 1$ ;

$x_{n + 1} = x_{n} + μ Φ^{T} r$ , $μ = {‖x_{k} - x_{k - 1}‖}_{2}^{2} / {‖Φ (x_{k} - x_{k - 1})‖}_{2}^{2}$ ;

$J \leftarrow M (x_{n + 1}, 2 m)$ (getting support according to signal model $M$ , 2m maximum indexes);

$T = \sup p (x_{n}) \cup J$ (merging supports);

${b|}_{T} = Φ_{T}^{†} y$ (signal estimating);

$x_{n + 1} \leftarrow M (b, m)$ (pruning signal estimate b with m maximum indexes);

$r_{n e w} = y - Φ x_{n + 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{matrix} x_{n + 1} = H_{S} (x_{n} + Φ^{T} (y - Φ x_{n})), \end{matrix}

(13)

with initial conditions

x_{0} = 0

, where

H_{S} (α)

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{matrix} x_{n + 1} = H_{S} (x_{n} + μ Φ^{T} (y - Φ x_{n})), \end{matrix}

(14)

where the step

μ = g_{Γ^{n}}^{T} g_{Γ^{n}} / g_{Γ^{n}}^{T} Φ_{Γ^{n}}^{T} Φ_{Γ^{n}} g_{Γ^{n}}

and

g = Φ^{T} (y - Φ 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 = Φ x^{s} + e$ and $x^{s}$ is a s-sparse vector, measurement matrix Ф satisfies RIP principle with $δ < 0.15$ ; then at iteration k, Model-Based Adaptive IHT will recover an approximation $x^{k}$ satisfying

\begin{matrix} {‖x^{s} - x_{k}‖}_{2} \leq 2^{- k} {‖x^{s}‖}_{2} + 5.04 {‖e‖}_{2} . \end{matrix}

(15)

Proof.

From the following equations

\begin{matrix} x_{n + 1} = x_{n} - μ g, \\ μ = \arg \min {‖Δ x - μ Δ g‖}^{2}, \end{matrix}

(16)

we have

\begin{matrix} μ = \frac{Δ x^{T} Δ x}{Δ x^{T} Δ g}, \end{matrix}

(17)

where

\begin{matrix} Δ x = x_{n + 1} - x_{n}, \\ Δ g = g_{n + 1} - g_{n} . \end{matrix}

(18)

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

\begin{matrix} (1 - δ_{s}) {‖x^{s}‖}_{2}^{2} \leq {‖Φ x^{s}‖}_{2}^{2} \leq (1 + δ_{s}) {‖x^{s}‖}_{2}^{2} . \end{matrix}

(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} (1 - δ_{s}) {‖x_{k} - x_{k - 1}‖}_{2}^{2} \leq {‖Φ (x_{k} - x_{k - 1})‖}_{2}^{2} \\ \leq (1 + δ_{s}) {‖x_{k} - x_{k - 1}‖}_{2}^{2}, \\ 1 - δ_{s} \leq \frac{{‖Φ (x_{k} - x_{k - 1})‖}_{2}^{2}}{{‖x_{k} - x_{k - 1}‖}_{2}^{2}} \leq 1 + δ_{s} . \end{array}

(20)

Based on the equation

x_{n + 1} = x_{n} - μ g

, we recognize that the formula at iteration

n + 1

x_{n + 1} = x_{n} + μ Φ^{T} (y - Φ x_{n})

. Let

g = Φ^{T} (Φ x - y)

be the negative gradient of

{‖y - Φ x‖}_{2}^{2}

. Substituting g into the term of μ,

\begin{array}{l} μ = \frac{Δ x^{T} Δ x}{Δ x^{T} Δ g} \\ = \frac{{(x_{k} - x_{k - 1})}^{T} (x_{k} - x_{k - 1})}{{(x_{k} - x_{k - 1})}^{T} (Φ^{T} Φ x_{k} - Φ^{T} y - Φ^{T} Φ x_{k - 1} + Φ^{T} y)} \\ = \frac{{(x_{k} - x_{k - 1})}^{T} (x_{k} - x_{k - 1})}{{(x_{k} - x_{k - 1})}^{T} Φ^{T} Φ (x_{k} - x_{k - 1})} \\ = \frac{{‖x_{k} - x_{k - 1}‖}_{2}^{2}}{{‖Φ (x_{k} - x_{k - 1})‖}_{2}^{2}}, \end{array}

(21)

we obtain

\begin{matrix} μ = \frac{{‖x_{k} - x_{k - 1}‖}_{2}^{2}}{{‖Φ (x_{k} - x_{k - 1})‖}_{2}^{2}} . \end{matrix}

(22)

Replacing μ in (20), we have

\begin{matrix} \frac{1}{1 + δ} \leq μ \leq \frac{1}{1 - δ} . \end{matrix}

(23)

Firstly, we use the triangle inequality to bound

\begin{matrix} {‖x^{s} - x_{n + 1}‖}_{2} \leq {‖x^{s} - a‖}_{2} + {‖x_{n + 1} - a‖}_{2}, \end{matrix}

(24)

and we know that

x^{n + 1}

is a better s -term approximation to a than

x^{s}

, so that we obtain

\begin{matrix} {‖x^{s} - x_{n + 1}‖}_{2} \leq 2 {‖x^{s} - a‖}_{2} . \end{matrix}

(25)

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

\begin{matrix} x_{n + 1} = x_{n} + μ Φ^{T} (Φ (x_{n + 1} + e) - Φ x_{n}) . \end{matrix}

(26)

Assuming r is residual, we have

\begin{matrix} a = x_{n} + μ Φ^{T} (Φ r + e) = x_{n} + μ Φ^{T} Φ r + μ Φ^{T} e . \end{matrix}

(27)

With these positions, we obtain

\begin{array}{l} {‖ x^{s} - x_{n + 1} ‖}_{2} \leq 2 {‖ x^{s} - x_{n} - u Φ^{T} Φ r_{n} - u Φ^{T} e ‖}_{2} \\ \leq 2 {‖ r_{B^{N + 1}}^{n} - u Φ_{B^{n + 1}}^{T} Φ r_{n} ‖}_{2} + 2 {‖ u Φ_{B^{n + 1}}^{T} e ‖}_{2} \\ = ‖ (I - u Φ_{B^{m + 1}}^{T} Φ_{B^{n + 1}}) r_{B^{n + 1}}^{n} \\ - u Φ_{B^{n + 1}}^{T} Φ_{B^{n} \ B^{n + 1}} r_{B^{n} \ B^{n + 1}}^{n} | |_{2} + 2 {‖ u Φ_{B^{n + 1}}^{T} e ‖}_{2} \\ \leq 2 {‖ (I - u Φ_{B^{n + 1}}^{T} Φ_{B^{n + 1}}) r_{B^{n + 1}}^{n} ‖}_{2} \\ + 2 {‖ u Φ_{B^{n + 1}}^{T} Φ_{B^{n} \ B^{n + 1}} r_{B^{n} \ B^{n + 1}}^{n} ‖}_{2} + 2 {‖ u Φ_{B^{n + 1}}^{T} e ‖}_{2}, \end{array}

(28)

where the sets

B^{n} ∖ B^{n + 1}

and

B^{n + 1}

are incoherent and

|B^{n} \cup B^{n + 1}| \leq 3 s

(

3 s

refers to the sparsity).

From (19), we have

\begin{matrix} {‖Φ_{Γ}^{T} y‖}_{2} \leq \sqrt{1 + δ} {‖y‖}_{2}, \end{matrix}

(29)

\begin{matrix} (1 - δ_{Γ}) {‖x_{Γ}‖}_{2} \leq {‖Φ_{Γ}^{T} Φ_{Γ} x_{Γ}‖}_{2} \leq (1 + δ_{Γ}) {‖x_{Γ}‖}_{2} . \end{matrix}

(30)

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

\begin{matrix} {‖(I - μ Φ_{Γ}^{T} Φ_{Γ}) x_{Γ}‖}_{2} \leq \frac{δ_{Γ}}{1 - δ_{Γ}} {‖x_{Γ}‖}_{2} . \end{matrix}

(31)

For two disjoint sets Γ and Λ, all Φ for which the RIP holds with

s = |Γ \cup Λ|

\begin{matrix} {‖μ Φ_{Γ}^{T} Φ_{Λ} x_{Λ}‖}_{2} \leq \frac{δ_{Λ}}{1 - δ_{Λ}} {‖x_{Λ}‖}_{2} . \end{matrix}

(32)

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

\begin{matrix} {‖Φ_{Γ}^{T} y‖}_{2} \leq \sqrt{1 + δ} {‖y‖}_{2}, \\ {‖(I - μ Φ_{Γ}^{T} Φ_{Γ}) x_{Γ}‖}_{2} \leq \frac{δ_{Γ}}{1 - δ_{Γ}} {‖x_{Γ}‖}_{2}, \\ {‖μ Φ_{Γ}^{T} Φ_{Λ} x_{Λ}‖}_{2} \leq \frac{δ_{Λ}}{1 - δ_{Λ}} {‖x_{Λ}‖}_{2} . \end{matrix}

(33)

Substituting (33) into (28), we have

\begin{array}{l} {‖x^{s} - x_{n + 1}‖}_{2} \leq \frac{2 δ_{2 s}}{1 - δ_{2 s}} {‖r_{B^{n + 1}}^{n}‖}_{2} + \frac{2 δ_{3 s}}{1 - δ_{3 s}} {‖r_{B^{n} ∖ B^{n + 1}}^{n}‖}_{2} \\ + \frac{2 \sqrt{1 + δ}}{1 - δ} {‖e‖}_{2} . \end{array}

(34)

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

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

\begin{matrix} {‖r^{n + 1}‖}_{2} \leq \frac{2 \sqrt{2} δ_{3 s}}{1 - δ_{3 s}} {‖r‖}_{2} + \frac{2 \sqrt{1 + δ}}{1 - δ} {‖e‖}_{2} . \end{matrix}

(35)

δ < 0.15

, the following inequality stands:

\begin{matrix} {‖r^{n + 1}‖}_{2} \leq 0.5 {‖r^{n}‖}_{2} + 2.52 {‖e‖}_{2} . \end{matrix}

(36)

Iterating the relationship (36), and realizing that

[1 + 1 / 2 + \dots + {(1 / 2)}^{k}] \times 2.52 \leq 5.04

with the initial condition

x_{0} = 0

, finally, we obtain the error bound described by Theorem 5

\begin{matrix} {‖r_{k}‖}_{2} \leq 2^{- k} {‖x^{s}‖}_{2} + 5.04 {‖e‖}_{2} . \end{matrix}

(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 = 4 K$ (K represents the sparse of signal), $D = 4$ , and $α = M / 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	$1.92 * e - 3$	$2.86 * e - 4$	$1.65 * e - 4$
Total time	1.984	0.014	$6.82 * e - 4$	$2.58 * e - 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.

Footnotes

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Thanks are due to Michigan State University for providing the real monitoring data of St. Helens volcano while the first author was with Michigan State University as visiting scholar. This work was supported in part by National Natural Science Foundation of China under Grant 61171089 and in part by the Fundamental Research Funds for the Central Universities under Grant CDJZR10160005. This paper is financially supported by the 2013 Innovative Team Construction Project of Chongqing Universities.

References

Werner-Allen

Johnson

Ruiz

Lees

Welsh

Monitoring volcanic eruptions with a wireless sensor network

Proceedings of the 2nd European Workshop onWireless Sensor Networks (EWSN ′05)

February 2005

Istanbul, Turkey

108 120

10.1109/ewsn.2005.1462003

2-s2.0-27744578667

Werner-Allen

Lorincz

Johnson

Lees

Welsh

Fidelity and yield in a volcano monitoring sensor network

Proceedings of the 7th USENIX Symposium on Operating Systems Design and Implementation (OSDI ′06)

November 2006

Seattle, Wash, USA

USENIX Association

381 396

Song

W.-Z.

Huang

Shirazi

LaHusen

Air-dropped sensor network for real-time high-fidelity volcano monitoring

Proceedings of the 7th ACM International Conference on Mobile Systems, Applications, and Services (MobiSys ′09)

June 2009

ACM

305 318

10.1145/1555816.1555847

2-s2.0-70450241180

Baraniuk

R. G.

Compressive sensing

IEEE Signal Processing Magazine 2007 24 4 118 124

10.1109/MSP.2007.4286571

2-s2.0-34548253373

Donoho

D. L.

Tsaig

Extensions of compressed sensing

Signal Processing 2006 86 3 533 548

Candès

E. J.

Eldar

Y. C.

Needell

Randall

Compressed sensing with coherent and redundant dictionaries

Applied and Computational Harmonic Analysis 2011 31 1 59 73

10.1016/j.acha.2010.10.002

MR2795875

2-s2.0-79955778301

Candes

E. J.

Romberg

Tao

Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information

IEEE Transactions on Information Theory 2006 52 2 489 509

10.1109/tit.2005.862083

MR2236170

2-s2.0-31744440684

Bajwa

Haupt

Sayeed

Nowak

Compressive wireless sensing

Proceedings of the 5th ACM/IEEE International Conference on Information Processing in Sensor Networks (IPSN ′06)

April 2006

Nashville, Tenn, USA

ACM

134 142

10.1145/1127777.1127801

2-s2.0-34247357879

Cevher

Boufounos

Baraniuk

R. G.

Gilbert

A. C.

Strauss

M. J.

Near-optimal Bayesian localization via incoherence and sparsity

Proceedings of the International Conference on Information Processing in Sensor Networks (IPSN ′09)

April 2009

205 216

2-s2.0-71049194736

10.

Mamaghanian

Khaled

Atienza

Vandergheynst

Compressed sensing for real-time energy-efficient ECG compression on wireless body sensor nodes

IEEE Transactions on Biomedical Engineering 2011 58 9 2456 2466

10.1109/TBME.2011.2156795

2-s2.0-80052049053

11.

Berinde

Gilbert

A. C.

Indyk

Karloff

Strauss

M. J.

Combining geometry and combinatorics: a unified approach to sparse signal recovery

Proceedings of the 46th Annual Allerton Conference on Communication, Control, and Computing

September 2008

798 805

10.1109/allerton.2008.4797639

2-s2.0-64549162832

12.

Zhang

Yang

M.-H.

Real-time compressive tracking

Computer Vision—ECCV 2012 2012 7574

Berlin, Germany

Springer

864 877 Lecture Notes in Computer Science

10.1007/978-3-642-33712-3_62

13.

Figueiredo

M. A. T.

Nowak

R. D.

Wright

S. J.

Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems

IEEE Journal of Selected Topics in Signal Processing 2007 1 4 586 597

10.1109/jstsp.2007.910281

2-s2.0-39449126969

14.

Blumensath

Davies

M. E.

Iterative hard thresholding for compressed sensing

Applied and Computational Harmonic Analysis 2009 27 3 265 274

10.1016/j.acha.2009.04.002

15.

Mallat

S. G.

Zhang

Matching pursuits with time-frequency dictionaries

IEEE Transactions on Signal Processing 1993 41 12 3397 3415

10.1109/78.258082

ZBL0842.94004

2-s2.0-0027842081

16.

Tropp

J. A.

Gilbert

A. C.

Signal recovery from random measurements via orthogonal matching pursuit

IEEE Transactions on Information Theory 2007 53 12 4655 4666

10.1109/tit.2007.909108

MR2446929

2-s2.0-64649083745

17.

Needell

Tropp

J. A.

CoSaMP: iterative signal recovery from incomplete and inaccurate samples

Applied and Computational Harmonic Analysis 2009 26 3 301 321

10.1016/j.acha.2008.07.002

MR2502366

2-s2.0-62749175137

18.

Baraniuk

R. G.

Optimal tree approximation with wavelets

3813

Wavelet Applications in Signal and Image Processing VII

July 1999

Denver, Colo, USA

196 207 Proceedings of SPIE

10.1117/12.366780

19.

Baraniuk

R. G.

Jones

D. L.

A signal-dependent time-frequency representation: fast algorithm for optimal kernel design

IEEE Transactions on Signal Processing 1994 42 1 134 146

10.1109/78.258128

2-s2.0-0028312573

20.

Baraniuk

R. G.

Cevher

Duarte

M. F.

Hegde

Model-based compressive sensing

IEEE Transactions on Information Theory 2010 56 4 1982 2001

10.1109/tit.2010.2040894

MR2654489

2-s2.0-77950244328

21.

Liu

Tan

Zhou

Xing

Song

W.-Z.

Lees

J. M.

Volcanic earthquake timing using wireless sensor networks

Proceedings of the 12th ACM/IEEE International Conference on Information Processing in Sensor Networks (IPSN ′13)

April 2013

Philadelphia, Pa, USA

ACM

91 102

10.1145/2461381.2461396

2-s2.0-84876753952

22.

Haupt

Bajwa

W. U.

Raz

Nowak

Toeplitz compressed sensing matrices with applications to sparse channel estimation

IEEE Transactions on Information Theory 2010 56 11 5862 5875

10.1109/tit.2010.2070191

MR2808938

2-s2.0-77956563404

23.

Varga

R. S.

Geršgorin and his Circles 2004 36

Berlin, Germany

Springer

10.1007/978-3-642-17798-9

MR2093409

24.

Blumensath

Davies

M. E.

Normalized iterative hard thresholding: guaranteed stability and performance

IEEE Journal on Selected Topics in Signal Processing 2010 4 2 298 309

10.1109/jstsp.2010.2042411

2-s2.0-77949703326

25.

Barzilai

Borwein

J. M.

Two-point step size gradient methods

IMA Journal of Numerical Analysis 1988 8 1 141 148

10.1093/imanum/8.1.141

MR967848

2-s2.0-0001531895

26.

Sleeman

van Eck

Robust automatic P-phase picking: an on-line implementation in the analysis of broadband seismogram recordings

Physics of the Earth and Planetary Interiors 1999 113 1–4 265 275

10.1016/s0031-9201(99)00007-2

2-s2.0-0033393307