基于InSAR技術和GS-SVR算法的礦區地表開采沉陷預計

馬 飛，隋立春，姚頑強，湯伏全

(1.長安大學 地質工程與測繪學院,陜西 西安 710054;2.西安科技大學 測繪科學與技術學院,陜西 西安 710054;3.地理國情監測國家測繪地理信息局 工程技術研究中心,陜西 西安 710054)

煤炭開采引起的礦區地表沉陷災害已經引起了人們越來越多的關注。然而礦區沉降的常規監測方法大多是重復采集被測地區的數據[1-5]，通過計算得出被監測地區在不同時期的沉降速率。但是西部黃土高原地區，土質疏松，塬梁破碎，溝壑縱橫，難以尋找穩定的控制點，使得常規的監測方法耗費人力、物力巨大，效率低下。

近年來，隨著衛星遙感領域的不斷發展，隨著大量的商用SAR數據的獲取，合成孔徑雷達理論技術也在不斷地完善，差分雷達干涉測量技術(D-InSAR)在地表形變監測方面的研究應用得到很大發展。D-InSAR技術所獲取的不是離散點的信息，而是大面積區域連續的地表形變信息，與傳統測量技術相比，其覆蓋范圍大、成本低、空間分辨率高，可以全天候工作，已成為獲取地表形變信息的新手段[6-10]。至今，劉廣、黃寶偉、范洪冬等學者利用該技術在城市地表形變監測、礦區沉陷監測等方面獲得了大量成果[13-15]。

支持向量機回歸(SVR)是支持向量在函數回歸領域的應用，它可以高效處理小樣本、非線性、高維數據問題。相比于最小二乘法、灰色模型、神經網絡等算法具有更高的運算效率和預測精度。因此，本文通過InSAR技術獲得彬長礦區的沉降區域范圍，監測隨著時間推移礦區地表累計沉降變化量。將沉降結果作為SVR算法的訓練學習樣本，建立礦區地表沉降預測模型，利用格網搜索法(GS)優化選取模型參數，對礦區地表沉降進行預測估計；結果表明，InSAR技術能夠準確獲取礦區的沉降范圍，監測得到的沉降值與常規水準數據用于礦區沉陷預計可以獲取相當的預測精度。

1 InSAR技術獲取礦區地表沉降量

InSAR(合成孔徑雷達干涉測量技術)是以波的干涉為基礎，用衛星兩次飛過同一地區(重復軌道方式)上空所獲得的兩幅微波圖像，兩幅圖像滿足相干條件，對其進行相位干涉處理，產生干涉條紋，反映出相位的變化，這種圖像叫做干涉圖。理想狀態，如果地面沒有變形或其它影響，通過解纏處理，解算出每一點的相位值，從而計算得出地面點到雷達的斜距以及地面點的高程。

二軌法首先利用一對跨越形變期的SAR圖像進行干涉處理(見圖1)，獲得的干涉相位可以表示為

Δφ=Δφtopo+Δφdef+

Δφflat+Δφatm+Δφnoise.

(1)

上式中Δφ表示干涉圖的纏繞相位，Δφtopo表示地形相位，Δφdef表示形變相位，Δφflat表示平地相位，Δφatm表示大氣影響，Δφnoise表示噪聲相位(包括系統熱噪聲、時間與空間失相關噪聲等)。為了得到準確的形變相位Δφdef，(1)式右邊其余各項應當逐項消除。利用二軌差分法結合已有的外部DEM模擬地形信息從而實現地形相位的去除[1]。大氣影響Δφatm是最主要的誤差源之一，尤其是電離層延遲誤差對L波段的SAR數據干涉測量影響不容忽視[3]，本文利用GPS數據進行電離層延遲誤差改正；平地相位Δφflat和噪聲相位Δφnoise分別通過基線估計和自適應濾波進行去除。各項誤差消除之后通過相位解纏得到形變相位Δφdef，進而計算出地表的形變量。

圖1 InSAR原理圖

由圖1中幾何關系推導可知，兩幅影像的相位差為

(2)

其中λ為雷達波長，ΔRd為地形變化在衛星視線向上的投影值；

由已知DEM反演得到的地形相位為

(3)

得到形變相位

(4)

式(4)給出了差分相位對地形形變的敏感度，ENVISAT/ASAR C波段，波長為5.6 cm，一般情況下，2.8 cm的斜距向變化即可引起一個2π的相位變化，ALOS/PALSAR L波段，波長為23.5 cm，11.7 cm的斜距向變化即可引起一個2π的相位變化。

在二軌法差分中，地形誤差對差分相位的影響主要取決于外部DEM的精度。

(5)

由上式可知，地形誤差對干涉相位的影響為

(6)

其中B⊥為基線垂直于衛星視線方向的分量，θ為主圖像視角。

2 基于GS-SVR的預計模型基本思想

SVR的基本思想是通過一個非線性函數映射將數據轉換到高維特征空間，然后對其進行線性回歸處理，轉化為求解高維特征空間的最優決策函數：

(7)

其中，X→Rn,w∈f，b為閾值。

上述問題通過經驗風險和VC維理論分析轉化為

(8)

式中，C為懲罰因子，ζ為誤差，ε為損失函數參數。

損失函數參數ε用于控制支持向量的個數和泛化能力，取值越小，精度越低，則支持向量越少，為了達到最優的擬合效果一般取值為(0.000 1～0.01)；懲罰因子C用于控制模型的復雜度，一般取值為(1～1 000)。為了選擇合適的預計模型參數，前人提出了多種參數優化選取方法，本文采用網格優化算法(GS)對模型參數進行尋優。首先確定C和ε的初始值，然后基于網格法全局搜索，獲取最優的C和ε值，確定預測方程。

根據于廣明[23]、于學義[24]等人的研究成果，通過大量實測資料分析礦區地表沉陷，發現沉陷非線性機理導致地表點的下沉過程在時間上沉陷呈不規則性，地表點下沉量在相等時間內大小不等呈非線性。本文通過InSAR技術獲取一系列與時間相關的沉降值。這些沉降值表現為在時間上非線性關系：{xi}={x1,x2,…,xn}，取前n-p個數據作為GS-SVR算法的訓練數據，構造如式(7)所示的預計函數，通過對式(8)進行計算分析，求解得到預計函數；對p組數據進行預測，將其與實測值比較分析，從而確定預測模型的精度。

3 實驗研究

3.1 實驗準備

彬長礦區位于陜西省關中西北部長武和彬縣境內，是國家規劃的黃隴基地的主力礦區之一。本區地處渭北黃土高原塬梁溝壑區，地勢從黃土塬梁向中間涇河谷地傾斜。礦區水土流失較為嚴重，土壤主要是黑壚土、黃綿土、紅土、淤土、潮土等，植被類型以闊葉落葉灌叢和草本植被為主。礦區東西長46 km，南北寬36.5 km，規劃面積978 km2，地質儲量為67.29億t，整個礦區生產能力達5 000萬t。其中，大佛寺煤礦2006年建成，大佛寺40301工作面為首采區，開采煤層厚度平均為11 m。

本文以彬長礦區工作面為例，選取5景PALSAR數據，組成干涉對進行方法驗證，為減少時間去相關的影響，選取的干涉對的時間間隔盡量小。干涉對組成情況如表1所示。

表1 ALOS PALSAR影像對基本參數

3.2 礦區沉降預測

本文數據處理采用二軌法差分干涉的方法，包括影像的預處理、主輔影像配準、生成干涉圖、干涉圖濾波、去平地效應、相位解纏、基線參數計算、去地形相位、生成差分干涉圖、地理編碼等環節，最終得到相干性分析圖、相位干涉圖、差分干涉相位圖、形變圖，具體處理流程見圖2。按照表1中的干涉對可獲取4組沉降圖，以第一組為參考將后續得到的沉降圖依次疊加即可得到在監測時間段內的累積沉降量，具體如圖3所示。

由圖3可知，隨著時間的推移，礦區的累積沉降量越來越大，逐漸形成為一個沉陷盆地， 本文選取沉降中心的一個沉降點(O)以及沿走向線(T1、T2)和傾向線(Q1、Q2)方向各兩個點為研究對象，對第3節提出的方法進行驗證。取前4組數據作為測試樣本，對第5組沉降值進行預測，選取的沉降信息如表2所示。按照式(7)和式(8)建立預測函數，預測結果如表3所示：表3中MAE表示絕對誤差，MRE表示相對誤差。

由表3可知，預測結果與GPS技術測得的結果較為一致，最大絕對誤差為3 mm，最大相對誤差為5.9%，其預測精度滿足工程實踐的需求，將InSAR技術與SVR算法相結合可以用于礦區沉降預測的實際應用。分析表3可知在沉降量大的區域其預測精度相對較高，所以在絕對誤差相等時，O點處的相對誤差最小。

圖2 二軌法差分干涉處理流程圖

圖3 時間序列沉降圖

獲取時間時間間隔/dO/mmT1/mmT2/mmQ1/mmQ2/mm2007年7月0000002007年8月46-15-10-12-12-102007年10月46-26-17-20-18-142007年11月46-57-37-39-37-322008年1月46-81-52-51-48-46

4 結 論

將InSAR用于礦區開采沉陷可以獲得礦區的整個開采沉降影響范圍和發展趨勢，且其監測技術獲得的沉降值與常規水準測量方法的精度相當，可以為沉陷預計模型提供良好的預測數據。

本文提出用InSAR技術與GS-SVR算法相結合預測礦區沉降值，在樣本數據量較少的情況下也可以準確的預測礦區的沉降值。本實驗采用5組數據進行實驗，結果表明：本文提出的預測模型獲得的沉降值與GPS技術監測得到的沉降值最大誤差為3 mm，最大相對誤差為5.9%，其預測精度符合工程的應用需求。

表3 預測結果與實測結果對比

分析表3中的絕對誤差和相對誤差可知，二者在絕對誤差相同的情況下，沉降量較大值的相對誤差反而小。而表3中的S0點為實驗區的沉降中心點。所以該方法用于沉降中心的沉降預測比非沉降中心可獲得更高的預測精度。