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

Authors: Guojin Liu; Qian Zhang; Yuyuan Yang; Zhenzhi Yin; Bin Zhu
Journal: International Journal of Distributed Sensor Networks
 (2015)

---

## 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.


---

## Body

## 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. Section2 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:(1)θi=αM×1with  probability12α0with  probability1-1α-1with  probability12α.It has zero-mean value, Eθi2=1/M, and θi≤α/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×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 vectorx, 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 α/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(2)1-εx22≤Φx22≤1+εx22,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 1-ε,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(3)M≥24α2ε2-8cK2log⁡N,where α>0 is defined in (1), the RIP in (2) holds with probability at least 1-exp-cM/α2K2.Proof.
Before the proof, three lemmas will be useful.Lemma 2 (Gershgorin disc theorem [23]).
LetA∈Cn×n, A=aij, so the n eigenvalues will fall into the union ∪DiA of n disks(4)DiA=z:z-aij≤Ri,i=1,2,…,n,and Ri=∑j=1,j≠inaij.Lemma 3 (see [22]).
Letxi,i=1,2,…,k, be a sequence of i.i.d. zero-mean bounded random variables such that xi≤a, and with variance Exi2=σ2, then,(5)P∑i=1kxi2-kσ2≥t≤2exp-2t2ka4.Lemma 4 (see [22]).
Letxi,i=1,2,…,k, and yi,i=1,2,…,k, be sequences of i.i.d. zero-mean bounded random variables satisfying xi≤a, xiyi≤a2; then,(6)P∑i=1kxiyi≥t≤2exp⁡-t22ka4.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 gii of its Gramian matrix G=ΦTΦ satisfies gii-1<ε1, and off-diagonal elements satisfy gij<ε2/K.In the following proof, we will show that the proposed matrix will satisfy RIP with high probability. Applying Lemma3, each diagonal element gii=∑j=1Mθj2 of Gramian matrix G=ΦTΦ satisfies(7)Pgii-1≥ε1≤2exp⁡-2Mε12α2,P⋃i=1Ngii-1≥ε1≤2Nexp⁡-2Mε12α2.For the off-diagonal element given by gij, such as g12=θ1θ2+θn+1θn+2+θ2n+1θ2n+2+⋯, we recognize that the entries of the sum are independent. From Lemma 4, we can get(8)Pgij≥ε2K≤2exp⁡-Mε222α2K2.Due to the fact thatgij=gji, the number of off-diagonal elements gij is N2-N<N2. We have(9)P⋃i=1N⋃j=1j≠iNgij≥ε2K≤2N2exp⁡-Mε222α2K2.LetT={ΦM×N, which satisfies RIP (K,ε)}, and T- is the complement of T; we have (10)PT-≤2Nexp⁡-2Mε12α2+2N2exp⁡-Mε222α2K2.AssumingN≥3 and ε1=ε2=ε/2, we obtain(11)PT-≤2Nexp⁡-Mε22α2+2N2exp⁡-Mε28α2K2≤3N2exp⁡-Mε28α2K2.Then we havec<ε2/8, M≥24α2/ε2-8cK2log⁡N, and(12)PT-≤exp⁡-cMα2K2.So for anyc<ε2/8 and M≥24α2/ε2-8cK2log⁡N, the measurement matrix ΦM×N satisfies RIP (K,ε) with probability at least 1-exp-cM/α2K2.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:x1=0, residual r=y, iterations n=0;xn 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 don=n+1;xn+1=xn+μΦTr, μ=xk-xk-122/Φxk-xk-122;J←Mxn+1,2m (getting support according to signal model M, 2m maximum indexes);T=sup⁡pxn∪J (merging supports);bT=ΦT†y (signal estimating);xn+1←Mb,m (pruning signal estimate b with m maximum indexes);rnew=y-Φxn+1 (updating measurement residual);(3)
End whilereturnx=xn.(B) Error Bound Analysis of the Algorithm. Starting with the basic principle of IHT, the iteration is(13)xn+1=HSxn+ΦTy-Φxn,with initial conditions x0=0, where HSα 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(14)xn+1=HSxn+μΦTy-Φxn,where the step μ=gΓnTgΓn/gΓnTΦΓnTΦΓngΓ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 measurementy=Φxs+e and xs 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 xk satisfying(15)xs-xk2≤2-kxs2+5.04e2.Proof.
From the following equations(16)xn+1=xn-μg,μ=arg⁡min⁡Δx-μΔg2,we have(17)μ=ΔxTΔxΔxTΔg,where(18)Δx=xn+1-xn,Δg=gn+1-gn.
In the proposed method, (17) was applied to (14). It provides faster convergence rate and simple calculation of μ. Given a noisy observation y=Φxs+e and xs is a s-sparse vector. The matrix Ф satisfies the RIP of order s, if there exists δs<1 satisfying(19)1-δsxs22≤Φxs22≤1+δsxs22.
At the iterationk-1 and k, the recovery approximations in xk-1 and xk have the same properties with xs. They are s-sparse in the same index. We can infer that xk-xk-1 is s-sparse and applicable to the RIP of order s. Therefore we have (20)1-δsxk-xk-122≤Φxk-xk-122≤1+δsxk-xk-122,1-δs≤Φxk-xk-122xk-xk-122≤1+δs.Based on the equation xn+1=xn-μg, we recognize that the formula at iteration n+1 is xn+1=xn+μΦT(y-Φxn). Let g=ΦTΦx-y be the negative gradient of y-Φx22. Substituting g into the term of μ,(21)μ=ΔxTΔxΔxTΔg=xk-xk-1Txk-xk-1xk-xk-1TΦTΦxk-ΦTy-ΦTΦxk-1+ΦTy=xk-xk-1Txk-xk-1xk-xk-1TΦTΦxk-xk-1=xk-xk-122Φxk-xk-122,we obtain(22)μ=xk-xk-122Φxk-xk-122.Replacing μ in (20), we have(23)11+δ≤μ≤11-δ.
Firstly, we use the triangle inequality to bound(24)xs-xn+12≤xs-a2+xn+1-a2,and we know that xn+1 is a better s-term approximation to a than xs, so that we obtain(25)xs-xn+12≤2xs-a2.
Then expandinga from the iteration, we have xn+1=xn+μΦT(y-Φxn). We recognize that the following equation stands:(26)xn+1=xn+μΦTΦxn+1+e-Φxn.Assuming r is residual, we have(27)a=xn+μΦTΦr+e=xn+μΦTΦr+μΦTe.With these positions, we obtain(28)xs-xn+12≤2xs-xn-μΦTΦrn-μΦTe2≤2rBn+1n-μΦBn+1TΦrn2+2μΦBn+1Te2=2I-μΦBn+1TΦBn+1rBn+1n-μΦBn+1TΦBn∖Bn+1rBn∖Bn+1n2+2μΦBn+1Te2≤2I-μΦBn+1TΦBn+1rBn+1n2+2μΦBn+1TΦBn∖Bn+1rBn∖Bn+1n2+2μΦBn+1Te2,where the sets Bn∖Bn+1 and Bn+1 are incoherent and Bn∪Bn+1≤3s (3s refers to the sparsity).
From (19), we have(29)ΦΓTy2≤1+δy2,(30)1-δΓxΓ2≤ΦΓTΦΓxΓ2≤1+δΓxΓ2.Combining the above two equations with (23), we can get (31)I-μΦΓTΦΓxΓ2≤δΓ1-δΓxΓ2.For two disjoint sets Γ and Λ, all Φ for which the RIP holds with s=Γ∪Λ,(32)μΦΓTΦΛxΛ2≤δΛ1-δΛxΛ2.Combining (29), (31), and (32), we can get (33)ΦΓTy2≤1+δy2,I-μΦΓTΦΓxΓ2≤δΓ1-δΓxΓ2,μΦΓTΦΛxΛ2≤δΛ1-δΛxΛ2.Substituting (33) into (28), we have(34)xs-xn+12≤2δ2s1-δ2srBn+1n2+2δ3s1-δ3srBn∖Bn+1n2+21+δ1-δe2.
Sinceδ2s and δ3s refer to the parameters from RIP of orders 2s and 3s, respectively, we have the bound δ2s≤δ3s.
FurthermorerBn+1n and rBn∖Bn+1n are orthogonal, so that rBn+1n2+rBn∖Bn+1n2≤2r2. According to (34), we have(35)rn+12≤22δ3s1-δ3sr2+21+δ1-δe2.If δ<0.15, the following inequality stands:(36)rn+12≤0.5rn2+2.52e2.Iterating the relationship (36), and realizing that 1+1/2+⋯+1/2k×2.52≤5.04 with the initial condition x0=0, finally, we obtain the error bound described by Theorem 5(37)rk2≤2-kxs2+5.04e2.


## 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 matrixA. 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 α=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 ofD (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 ofD.Figure2 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 Table2. 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
-
4Figure 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).Figure4 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 Figure5. 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.



---

*Source: 769246-2015-07-13.xml*