基于提升小波變換的礦區地表形變監測數據處理研究

楊俊山， 劉素潔， 杜春苗

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

隨著人口的增長，人們對礦產資源的需求急劇增加，而礦產資源開采常導致地質沉陷、滑坡等災害，嚴重影響了礦區人民正常的生產生活. 因此，對礦區地表，特別是重點區域地表進行實時形變監測，研究其變化規律，對于及時發現地表形變及災害預報預警具有重要意義.

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

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

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

1 提升小波變換及Kalman濾波動態解算的基本原理

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

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

1.1.1 分解過程

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

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

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

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

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

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

1.1.2 重構過程

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

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

2）反預測. 獲取偶信號序列之后，即可通過預測數據恢復奇信號序列.

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

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

1.2 Kalman濾波動態解算的基本原理

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

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

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

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

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

2 GNSS監測數據的處理分析

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

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

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

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

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

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

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

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

3 結論

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

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

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

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