基于提升小波變換的礦區(qū)地表形變監(jiān)測數(shù)據(jù)處理研究

楊俊山， 劉素潔， 杜春苗

（1.河南省地球物理空間信息研究院，鄭州 450000； 2.河南省地質(zhì)物探工程技術(shù)研究中心，鄭州 450000）

隨著人口的增長，人們對礦產(chǎn)資源的需求急劇增加，而礦產(chǎn)資源開采常導(dǎo)致地質(zhì)沉陷、滑坡等災(zāi)害，嚴重影響了礦區(qū)人民正常的生產(chǎn)生活. 因此，對礦區(qū)地表，特別是重點區(qū)域地表進行實時形變監(jiān)測，研究其變化規(guī)律，對于及時發(fā)現(xiàn)地表形變及災(zāi)害預(yù)報預(yù)警具有重要意義.

在礦區(qū)地表形變監(jiān)測及數(shù)據(jù)處理方面，研究人員開展了富有成效的工作. 欒元重和韓李濤［1］通過對礦區(qū)GPS變形監(jiān)測網(wǎng)的建立、實時監(jiān)測、基線平差、變形分析及分形特征等問題進行研究指出，地表點的移動具有較強的分形增長規(guī)律，GPS 變形監(jiān)測技術(shù)能夠揭示地表移動的非線性特征. 范洪冬［2］利用DInSAR、PS-DInSAR、SBAS等技術(shù)對天津市的失水沉降及某礦的開采沉陷監(jiān)測進行了試驗研究，闡述了SAR、InSAR、DInSAR 等技術(shù)的基本原理，分析了DInSAR 技術(shù)在沉降監(jiān)測應(yīng)用中的局限性. 然而，開采沉陷是一種大變形的地質(zhì)災(zāi)害，僅僅采用DInSAR技術(shù)無法獲得開采沉陷的主值以及完整的沉陷信息. 成曉倩等［3］為了彌補DInSAR提取大形變的不足，以DInSAR和PIM結(jié)果組成的混合數(shù)據(jù)為數(shù)據(jù)源，采用Gauss模型重建開采沉陷主斷面，實現(xiàn)了下沉特征曲線數(shù)學模型的重建. 王偉等［4］利用衛(wèi)星導(dǎo)航定位基準站（CORS）網(wǎng)全球衛(wèi)星導(dǎo)航系統(tǒng)（GNSS）對區(qū)域地面大地高、地面重力和地傾斜變化進行了連續(xù)監(jiān)測，同時通過與部分地質(zhì)災(zāi)害信息進行時空分析，確定了用于地面穩(wěn)定性變化分析的要素和權(quán)重，進而定量分析了地面穩(wěn)定性變化. 雖然GPS技術(shù)的精度較高［5］，但是其觀測點的空間密度低、觀測周期長，并且GPS 測站及觀測點布設(shè)與維護較為困難，布設(shè)的觀測點極易遭到破壞. DInSAR 技術(shù)主要用于短時期的形變研究，而在進行長時間、緩慢地表形變監(jiān)測時，會受到時/空基線失相干以及大氣效應(yīng)等因素的影響，存在精度低、結(jié)果可靠性較差的問題［6-7］，同時DInSAR 技術(shù)常需要有針對性地購買區(qū)域高精度影像數(shù)據(jù)，監(jiān)測成本較高. GNSS是利用接收上空的衛(wèi)星信號實現(xiàn)定位的，靜態(tài)測量精度能達到厘米級甚至毫米級，在復(fù)雜山區(qū)具有獨特的優(yōu)勢. GNSS形變監(jiān)測具有測站間無須通視和全天候自動觀測等優(yōu)點，但由于GNSS數(shù)據(jù)采集過程中易受天氣、樹木遮擋、多路徑效應(yīng)等多種因素影響，因此GNSS形變序列中常包含測量噪聲，這些噪聲會嚴重影響有效形變信息的提取和測量精度，導(dǎo)致監(jiān)測物的形變信息不準確［8-9］. 因此，如何從嚴重干擾的信號序列中提取出有效的形變信息、提高GNSS形變監(jiān)測解算的精度和穩(wěn)定性至關(guān)重要.

傅里葉變換能夠把信號映射到頻域內(nèi)，提取信號的頻譜，然后利用信號的頻譜特性分析時域內(nèi)難以看清的問題，但其無法反映某一局部時間內(nèi)信號的頻譜特性，因此可能會忽略短時間內(nèi)的信號變化［10-13］. 小波變換繼承并發(fā)展了短時傅里葉變換局部化的思想，能夠克服窗口大小不隨頻率變化的缺點，是進行信號時頻分析處理的理想工具，當前已應(yīng)用于信號去噪及信息提取等方面［14-16］. 小波變換能夠充分突出問題某些方面的特征，實現(xiàn)對時空頻率的局部化分析，并通過伸縮和平移運算對信號進行多尺度分解與細化，獲取信號的細節(jié). 小波去噪其實就是抑制信號中的無用部分、恢復(fù)信號中的有用部分的過程［17-18］. 在小波變換的基礎(chǔ)上，提升小波變換能夠通過提升模式構(gòu)造出小波函數(shù)，減小常規(guī)小波變換計算的復(fù)雜度，提高運算速度［19-21］，因而被廣泛應(yīng)用于大地測量中.

為了削弱GNSS監(jiān)測數(shù)據(jù)中的干擾噪聲、提取有效的監(jiān)測形變信息、提高數(shù)據(jù)解算精度，本研究以河南某礦區(qū)的GNSS監(jiān)測數(shù)據(jù)為研究對象，將提升小波變換方法引入到GNSS監(jiān)測數(shù)據(jù)處理中，對該礦區(qū)多個監(jiān)測點的GNSS監(jiān)測數(shù)據(jù)進行了提升小波降噪及分解解算處理，并對比分析了常規(guī)Kalman濾波動態(tài)解算方法和提升小波降噪解算方法對GNSS監(jiān)測數(shù)據(jù)解算的精度和穩(wěn)定性.

1 提升小波變換及Kalman濾波動態(tài)解算的基本原理

1.1 提升小波變換的基本原理

提升小波變換算法是對傳統(tǒng)小波變換的一種改進，小波函數(shù)也不再由函數(shù)的平移和伸縮而產(chǎn)生，所有的運算都在時域上進行，能夠在一定程度上減少計算量，提高監(jiān)測數(shù)據(jù)的解算效率. 提升小波變換應(yīng)用于信號去噪時，首先是對分解后得到的低頻近似信號和高頻細節(jié)信號進行一些處理，然后再進行完全的反向重構(gòu)，即提升小波變換包含了分解和重構(gòu)兩個過程.

1.1.1 分解過程

提升小波變換的分解過程包括3個步驟，即分裂、預(yù)測、更新［20］. 具體過程如下：

1）分裂. 將獲取的監(jiān)測信號序列Sn分解成奇信號序列S2n+1和偶信號序列S2n，n=1，2，3，…. 除按奇偶信號序列劃分外，還可按其他方法進行分裂.

2）預(yù)測. 在分裂完成之后即可進行預(yù)測. 預(yù)測是在初始監(jiān)測信號序列相關(guān)性的基礎(chǔ)上，通過預(yù)測算子P利用偶信號序列預(yù)測奇信號序列，用預(yù)測奇信號序列值與實際奇信號序列值之差（即預(yù)測誤差δn）代替奇信號序列，以此來表示信號序列細節(jié)信息.

預(yù)測誤差δn可視為原信號序列的高頻部分.

3）更新. 更新過程利用更新算子U來處理預(yù)測誤差δn，然后將處理結(jié)果疊加至偶信號序列以對其進行更新，得到近似信號（又稱尺度系數(shù)T）.

通過重復(fù)上述3個步驟即可實現(xiàn)多尺度、多分辨率的提升小波變換.

1.1.2 重構(gòu)過程

提升小波變換重構(gòu)過程包含反更新、反預(yù)測及合并3個步驟，過程如下：

1）反更新. 給定信號Tn和預(yù)測誤差δn，則可恢復(fù)偶信號序列.

2）反預(yù)測. 獲取偶信號序列之后，即可通過預(yù)測數(shù)據(jù)恢復(fù)奇信號序列.

3）合并. 利用加法運算即可得到原信號序列.

在提升小波變換中，預(yù)測與更新是算法的核心，利用預(yù)測過程能夠獲取監(jiān)測信號中的高頻信息，而利用更新過程則可獲取監(jiān)測信號中的低頻信息. 因此，在變換時首先可將原始監(jiān)測信號分解為尺度系數(shù)及新的小波系數(shù)，根據(jù)有效形變信息及噪聲信息對應(yīng)的系數(shù)確定合適的閾值，然后利用提升小波算法對系數(shù)進行反變換，即可完成對監(jiān)測信號的降噪. 降噪之后，對監(jiān)測信號進行小波多尺度分解，得到分解后的信號在各個頻段的分布狀況，根據(jù)前期統(tǒng)計分析結(jié)果，可提高分解后的處理效率，削弱干擾噪聲的影響，提取有效的監(jiān)測信息.

1.2 Kalman濾波動態(tài)解算的基本原理

在GNSS監(jiān)測數(shù)據(jù)處理中，常采用Kalman 濾波算法進行動態(tài)解算. 在Kalman 濾波中，如果系統(tǒng)狀態(tài)噪聲和觀測噪聲均符合Gaussian分布，則離散線性狀態(tài)模型如公式（6）所示.

式中：Xk表示狀態(tài)向量；Φk/k-1表示狀態(tài)轉(zhuǎn)移矩陣；Zk表示觀測向量；Hk為觀測矩陣；Wk和Vk分別表示狀態(tài)噪聲和測量噪聲. Kalman濾波方差準則如式（7）所示.

式中：vk和vk/k-1分別表示觀測向量的殘差和預(yù)測狀態(tài)向量的殘差；Pk/k-1表示預(yù)測狀態(tài)向量的協(xié)方差矩陣；Rk表示觀測噪聲的協(xié)方差矩陣；min表示最小.

在Kalman濾波方差準則下可得到Kalman濾波解，即：

式中：Xk/k-1和Pk/k-1分別表示預(yù)測狀態(tài)向量及其協(xié)方差矩陣；Rk表示觀測噪聲的協(xié)方差矩陣；Kk表示增益矩陣；Xk表示狀態(tài)向量；Hk為觀測矩陣；Pk為Xk的協(xié)方差矩陣；Zk為觀測向量.

2 GNSS監(jiān)測數(shù)據(jù)的處理分析

河南某礦區(qū)位于河南西部地區(qū)，東起沙河與汝河交匯帶的洛崗斷層，西抵紅石山附近的郟縣斷層，南至湛河北岸的煤層露頭，北達汝河附近的襄郟斷層，東西長約40 km，南北寬約20 km，含煤面積約650 km2. 受地勢及環(huán)境影響，該礦區(qū)地形復(fù)雜，難以采用傳統(tǒng)水準測量方法對其地表進行形變監(jiān)測. 另外，由于DInSAR技術(shù)存在成本較高的問題，難以適用大范圍區(qū)域長期的地表形變監(jiān)測工作，因此近年來一直是在該礦區(qū)建設(shè)的CORS站的基礎(chǔ)上利用GNSS技術(shù)對其進行地表形變監(jiān)測. 但在長期的監(jiān)測工作中發(fā)現(xiàn)，監(jiān)測解算出的形變量與實際形變量常存在不符合的情況，有的監(jiān)測點監(jiān)測解算出的形變量與實際形變量的差距還相對較大，無法忽略不計. 通過分析得知，這是因為得到的GNSS監(jiān)測數(shù)據(jù)受噪聲干擾較為嚴重，易淹沒有效形變信息. 因此本研究以該礦區(qū)多個監(jiān)測點的GNSS監(jiān)測數(shù)據(jù)為研究對象，先采用提升小波變換方法對該礦區(qū)的GNSS監(jiān)測數(shù)據(jù)進行降噪處理后再對其進行解算，以期獲取更為精確的地表形變信息. 同時，為了驗證提升小波變換在GNSS形變監(jiān)測數(shù)據(jù)處理中的有效性，分別對該礦區(qū)的GNSS監(jiān)測數(shù)據(jù)進行了常規(guī)Kalman濾波動態(tài)解算（算法一）和提升小波降噪解算（算法二）. 以監(jiān)測點BD25為例，圖1給出了通過不同算法解算得到的該點的位置誤差.

圖1 通過不同算法解算得到的BD25監(jiān)測點的位置誤差Fig.1 Position errors of the monitoring point BD25 obtained by different algorithms

由圖1可知，通過提升小波降噪解算方法得到的位置誤差整體比通過常規(guī)Kalman濾波動態(tài)解算方法得到的位置誤差小，這說明經(jīng)過提升小波變換處理之后，干擾噪聲被明顯削弱，GNSS 監(jiān)測信號得到了明顯改善，所以最終的定位解算精度明顯得到提升. 從圖1還可以看出：通過常規(guī)Kalman 濾波動態(tài)解算方法得到的位置誤差在1.5 cm以內(nèi)，絕大部分歷元的位置誤差在1.0 cm以內(nèi)；通過提升小波降噪解算方法得到的位置誤差在1.0 cm以內(nèi)，精度相對較高. 因此，在實際的礦區(qū)地表形變監(jiān)測中，建議使用提升小波降噪解算方法對GNSS監(jiān)測數(shù)據(jù)進行解算，因為通過該方法能夠有效削弱干擾噪聲的影響，獲取有效的監(jiān)測位移信息，且監(jiān)測精度能夠滿足基本需求.

由于礦區(qū)的環(huán)境復(fù)雜，因此很多情況下，監(jiān)測信號受到的干擾是未知的，無法有效判斷干擾源，而通過分析監(jiān)測點的位置誤差振幅和頻率可在一定程度上為未知干擾的判斷提供依據(jù). 以BD25監(jiān)測點為例，通過繪制該監(jiān)測點提升小波降噪解算方法的頻譜圖，得到了該監(jiān)測點的位置誤差振幅和頻率，如圖2所示. 從圖2可以看出，大部分歷元解算的位置誤差在5 mm 以內(nèi)，而這些信號的頻率多集中在0～0.5 Hz 之間，這為后續(xù)有效形變監(jiān)測位移信息的獲取提供了幫助.

圖2 BD25監(jiān)測點提升小波降噪解算方法的頻譜圖Fig.2 Spectrum diagram of lifting wavelet denoising method for the monitoring point BD25

為進一步對比分析不同算法的性能，隨機選取該礦區(qū)6 個監(jiān)測點（AD6、AD8、AD25、BD25、BD28、CJ15）的監(jiān)測數(shù)據(jù)為研究對象，分別用兩種算法對這6個監(jiān)測點的GNSS監(jiān)測數(shù)據(jù)進行解算，得到各監(jiān)測點的位置誤差如表1所示.

表1 通過不同算法解算得到的各個監(jiān)測點的位置誤差Tab.1 Position errors of each monitoring point calculated by different algorithms

從表1可以看出，通過提升小波降噪解算方法解算獲取的位置誤差明顯小于通過常規(guī)Kalman濾波動態(tài)解算方法解算獲取的位置誤差，這是由于提升小波變換解算方法能夠有效削弱干擾噪聲的影響. 通過設(shè)定合理的閾值，能夠排除相關(guān)頻段的噪聲信號，從而有助于獲取有效的監(jiān)測信息，這不僅能夠為復(fù)雜環(huán)境下礦區(qū)地表GNSS靜態(tài)形變監(jiān)測技術(shù)提供幫助，還能為動態(tài)定位解算及其他領(lǐng)域的提升小波變換應(yīng)用提供借鑒.從表1還可以看出：通過常規(guī)Kalman 濾波動態(tài)解算方法解算出的各監(jiān)測點的位置誤差相差較大，這與各個監(jiān)測點的位置及周邊監(jiān)測環(huán)境有關(guān)；通過提升小波降噪解算方法解算出的各監(jiān)測點的位置誤差整體相差較小，但仍存在個別位置誤差較大的點，這可能是因為在預(yù)處理解算中并未完全剔除觀測粗差，導(dǎo)致解算結(jié)果出現(xiàn)較大偏差. 因此，在前期預(yù)處理解算過程中，可根據(jù)實際解算效果合理調(diào)整粗差探測的閾值，有效削弱觀測粗差的影響，以便在后續(xù)處理中進一步提高解算的精度.

3 結(jié)論

針對礦區(qū)地表形變監(jiān)測實踐中遇到的復(fù)雜噪聲干擾問題，利用提升小波變換對GNSS監(jiān)測數(shù)據(jù)進行了處理，并通過對比試驗驗證了提升小波降噪解算算法的有效性，得出結(jié)論如下：

1）提升小波變換能夠有效降低GNSS監(jiān)測信號的噪聲水平，削弱干擾噪聲的影響.

2）通過頻譜分析能夠獲取有效信號的頻率范圍，從而有助于利用提升小波變換在噪聲干擾環(huán)境下快速獲取有效的形變信息.

3）通過提升小波降噪解算算法對GNSS監(jiān)測數(shù)據(jù)進行解算得到的點位誤差可達到毫米級，但解算結(jié)果會受到異常觀測值的影響，因此，解算前應(yīng)加強GNSS監(jiān)測數(shù)據(jù)的預(yù)處理工作，有效剔除異常觀測數(shù)據(jù)，以便在后續(xù)處理中進一步提高解算的精度.