馬 飛，隋立春，姚頑強(qiáng)，湯伏全

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

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

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

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

1 InSAR技術(shù)獲取礦區(qū)地表沉降量

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

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

Δφ=Δφtopo+Δφdef+

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

(1)

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

圖1 InSAR原理圖

由圖1中幾何關(guān)系推導(dǎo)可知，兩幅影像的相位差為

(2)

其中λ為雷達(dá)波長，ΔRd為地形變化在衛(wèi)星視線向上的投影值；

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

(3)

得到形變相位

(4)

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

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

(5)

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

(6)

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

2 基于GS-SVR的預(yù)計(jì)模型基本思想

SVR的基本思想是通過一個(gè)非線性函數(shù)映射將數(shù)據(jù)轉(zhuǎn)換到高維特征空間，然后對其進(jìn)行線性回歸處理，轉(zhuǎn)化為求解高維特征空間的最優(yōu)決策函數(shù)：

(7)

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

上述問題通過經(jīng)驗(yàn)風(fēng)險(xiǎn)和VC維理論分析轉(zhuǎn)化為

(8)

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

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

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

3 實(shí)驗(yàn)研究

3.1 實(shí)驗(yàn)準(zhǔn)備

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

本文以彬長礦區(qū)工作面為例，選取5景PALSAR數(shù)據(jù)，組成干涉對進(jìn)行方法驗(yàn)證，為減少時(shí)間去相關(guān)的影響，選取的干涉對的時(shí)間間隔盡量小。干涉對組成情況如表1所示。

表1 ALOS PALSAR影像對基本參數(shù)

3.2 礦區(qū)沉降預(yù)測

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

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

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

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

圖3 時(shí)間序列沉降圖

獲取時(shí)間時(shí)間間隔/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 結(jié) 論

將InSAR用于礦區(qū)開采沉陷可以獲得礦區(qū)的整個(gè)開采沉降影響范圍和發(fā)展趨勢，且其監(jiān)測技術(shù)獲得的沉降值與常規(guī)水準(zhǔn)測量方法的精度相當(dāng)，可以為沉陷預(yù)計(jì)模型提供良好的預(yù)測數(shù)據(jù)。

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

表3 預(yù)測結(jié)果與實(shí)測結(jié)果對比

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