結合SBAS-InSAR與支持向量回歸的開采沉陷監測與預測

師蕓，李杰，呂杰，馬東暉

(1.西安科技大學 測繪科學技術學院，西安 710054；2.自然資源部煤炭資源勘查與綜合利用重點實驗室，西安 710021)

0 引言

煤炭是我國的主要消耗能源，開采后采空區周圍巖體失去原來的平衡狀態而發生移動，誘發滑坡、塌陷、地裂縫等地質災害[1]。煤炭開采引起地面沉降、塌陷等礦區地質災害問題亟待當地政府和企業解決。分析預測開采沉陷的可能性對礦區安全開采具有積極的意義[2]。

合成孔徑雷達差分干涉測量技術(differential interferometric synthetic aperture radar，D-InSAR)對礦區進行沉降監測已被眾多國內外學者證實具有可行性，其監測精度可達到cm甚至mm級[3-4]。然而，傳統D-InSAR技術觀測單一，易受時空基線的限制導致失相干現象[5]，從而無法獲得礦區內時序沉降量。為了克服上述問題，時序InSAR技術[6-7]逐漸發展了起來，以小基線集InSAR(small baseline subset InSAR，SBAS-InSAR)為代表的多主影像時序InSAR技術能夠進行大尺度長時間序列的地表形變監測，廣泛應用于城市地表形變[8-9]、地質災害監測[10]、礦區地表形變等領域[11]。

目前，常用的礦區沉降預測方法主要有概率積分法[12]、神經網絡(neural networks，NN)[13]、灰色模型(grey models，GM)[14]和支持向量機(support vector machine，SVM)[15]等。支持向量機是由Vapnik提出的一種機器學習算法，具有估算精度高、可解決小樣本和非線性問題、泛化能力強等特點。然而，SVM的性能受其參數影響較大，模型參數選取存在一定的主觀性與隨意性使其效率低下。遺傳算法(genetic algorithm，GA)[16]是一種通過模擬自然界遺傳機制和生物進化論來進行隨機優化的算法，它具有較強的全局搜索能力，不依賴特定的求解模型，但基本遺傳算法(sample genetic algorithm，SGA)存在欺騙問題，即進化過程緩慢，易早熟收斂。針對這一問題，Srinivas等[17]提出了一種根據種群適應度自動調整的自適應遺傳算法(adaptive genetic algorithm，AGA)，對SVM模型中的參數進行尋優，可以避免選參的盲目性從而得到全局最優解，同時提高預測效率和準確度。

本文提出將AGA-SVR模型應用于開采沉陷預測，利用17景Sentinel-1A影像，采用SBAS-InSAR技術實現對礦區的長時間序列的地表沉降監測，將獲得的沉降結果作為支持向量回歸(support vector regression，SVR)的學習訓練樣本，采用AGA對SVR模型中的參數進行尋優并得到預測結果。

1 礦區沉降監測與預測原理

1.1 SBAS-InSAR基本原理

SBAS-InSAR技術是在傳統InSAR技術的基礎上發展起來的，其主要思想是采用多主影像的方式，通過設置合理的時間與空間基線閾值將SAR影像組合成若干個小基線集合從而獲取干涉對，其次在每個子集內利用最小二乘法(least squares，LS)求解對應的地面沉降值，對不同子集間采用奇異值分解法(singular value decomposition，SVD)進行聯合求解，可以獲取觀測時間段內礦區的沉降時間序列。

假設在時間t0，t1，…，tn內獲取了N+1幅同一研究區域的SAR影像，選擇一幅作為主影像，將其余影像統一配準至該影像上，設置合適的時空基線后生成M幅干涉圖，M滿足式(1)。

(N+1)/2≤M≤N(N+1)/2

(1)

對于由ta和tb(ta

(2)

式中：λ代表雷達中心波長；dtb、dta分別代表以t0為初始時刻，對應于tb與ta時刻在雷達視線方向(line of sight，LOS)的累積形變量；Δφtop，j、Δφatm，j、Δφnoise，j分別代表殘余的地形相位差、大氣相位差、噪聲相位差。去除形變量以外的相位后，干涉相位簡化由式(3)至式(4)表示。

(3)

dtb-dta=νi(tb-ta)

(4)

式中：νi代表ta至tb時間段雷達視線方向平均形變速率。由此解纏后的差分干涉圖相位可由矩陣表示為式(5)。

Aν=δφ

(5)

式中：A是一個M×N的秩虧矩陣，需利用SVD法進行求解，從而對得到的各時間段內的平均速率進行積分最終獲取時間序列形變量。考慮到礦區地表以下沉為主，且水平移動對雷達視線方向的貢獻遠小于下沉值，所以可將LOS向形變dLOS直接轉換為垂直沉降值W(式(6))。

W=dLOS/cosθ

(6)

式中：θ為雷達入射角[18]。

1.2 SVR算法原理

SVR是建立在SVM基礎上的回歸算法，在預測回歸方面具有良好的泛化能力。其基本思想是通過核函數定義的非線性特征，將低維樣本數據映射至高維空間中，并在高維空間中對樣本進行線性回歸，獲取最優決策函數f(x)(式(7))。其中，輸入樣本為(x1，y1)，(x2，y2)，…，(xn，yn)。

f(x)=ωφ(x)+b

(7)

式中：ω為權值矢量；x為輸入變量；b為閾值。其優化問題用式(8)表示。

(8)

滿足：

(9)

引入拉格朗日函數，并將式(9)轉換為對偶形式，同時選取徑向基函數作為SVR的核函數，滿足KKT(Karush-Kuhn-Tucker，KKT)約束條件后，可求解出SVR回歸函數模型(式(10))。

(10)

(11)

1.3 自適應遺傳算法參數尋優

在SVR模型中，C、ε、g3個參數的選取對預測精度有很大的影響，對SVR進行參數尋優能有效避免算法陷入局部最優解。首先，建立SVR目標函數與AGA適應度函數的關系，在全局范圍內求解目標參數的最優解；然后，隨機創建初始種群并開始最優解迭代搜索，計算種群個體的適應度；最后，以最大適應度為標準自動調整交叉和變異的概率并生成下一代種群。循環迭代這一過程直到滿足優化終止條件，得到最優個體。算法流程如圖1所示。

圖1 自適應遺傳算法優化流程圖

適應度是描述種群個體的主要指標，根據適應度的大小對個體進行優勝劣汰。適應度函數是根據目標函數確定的判別種群中個體優劣性的評價標準。選擇SVR目標函數的倒數作為適應度函數，適應度函數表示為式(12)。

(12)

式中：e為避免分母為零所加的一個常數。

遺傳算法的主要操作包括選擇、交叉、變異等。選擇操作對種群中的個體按適應度大小進行排序，根據個體的適應度從上一代種群中選擇優良個體復制到下一代。交叉是遺傳算法的核心操作，即模擬遺傳學的雜交原理，隨機結合種群中的2個個體，交換部分染色體從而形成新的個體。變異操作以變異概率改變個體的基因組，從而維持種群的多樣性，避免算法過早收斂。交叉概率Pc和變異概率Pm分別如式(13)、式(14)所示。

(13)

(14)

式中：fmax為每代種群的最大適應度；favg為每代種群的平均適應度；f′為要交叉的2個個體中較大的適應度；f為變異個體的適應度。

自適應遺傳算法的核心思想是在算法前期保持較大的交叉和變異概率，以豐富種群的多樣性；在算法后期降低交叉和變異的概率，保證優良的個體不受破壞。同時，對于適應度高于種群平均適應度的個體，減小Pc和Pm，使得優良個體能夠進入下一代；對于低于平均適應度的個體，增大Pc和Pm，淘汰劣質個體。

2 實驗結果與分析

2.1 研究區概況和實驗數據

李家壕煤礦位于內蒙古自治區鄂爾多斯中南部，地處鄂爾多斯高原腹地。井田東西長約8.0 km，南北約9.0 km，于2011年5月開始試生產。井田構造為向西南傾斜的單斜構造，傾向220°～260°，煤層傾角小于5°，采用綜合機械化采煤方法進行綜采，其中12113工作面寬度300 m，采高2.3 m，埋深170 m，2017年可推進長度630 m，采煤量可達52萬噸。本文數據來源為哨兵1號衛星。哨兵衛星是歐洲航天局哥白尼計劃中的地球觀測衛星，能夠獲取C波段的雷達影像，其中Sentinel-1A重訪周期為12 d。選用2017年10月2日至2018年4月24日共17景Sentinel-1A影像，影像模式為干涉寬幅模式IW(interferometric wide)，分辨率為5 m×20 m，條帶寬度240 km，極化方式為HH。同時，使用精密軌道文件去除平地效應和大氣相位的影響。

2.2 時序InSAR監測及結果分析

時序InSAR通過對多時相的SAR影像進行參數估計，可以有效減少平地效應和大氣延遲的影響。SBAS技術是一種時序技術，按照空間基線的大小將所有數據組成若干個集合，每個集合的時間序列形變量可以通過最小二乘法求得，使用奇異值分解將多個小基線集進行聯合求解就可以獲得整個觀測周期形變信息，最終利用入射角將視線向形變量轉換為垂直沉降信息。

圖2表示2017年10月至2018年4月時間段該礦區年平均沉降速率圖，該工作面是自東向西開采。由圖2可以看出，整個工作面自東向西呈現逐漸下沉的趨勢，這是因為東側要比西側先開采，同時比西側更早進行回填。最大沉降速率為360 mm/a。

圖3表示以2017年10月2日為時間參考基準的累計沉降量時間序列圖。由圖3可知，隨著時間變化，該工作面的沉降量不斷增大，同時沉降范圍沿著工作面開采方向由東向西不斷擴大，形成明顯的沉降漏斗，截至2018年4月24日，最大的累計沉降量為163 mm。

圖2 年平均沉降速率圖

圖3 研究區時序累計沉降量

為了量化分析該工作面的時序累計沉降值，在該工作面中心做剖面提取若干像點進行分析，圖4為走向時序累計沉降值折線圖。由圖4可知，該工作面自東向西進行開采并隨著時間的累計沉降值不斷變大，出現一個快速下沉盆地，開采中心沉降值從開始的20 mm以360 mm/a平均速率增加到163 mm，說明該工作面還未穩定，隨著時間推移沉降值還會不斷增大且沉降中心持續向西移動。結合礦區資料可知，自2018年1月起對該工作面開采結束的部分進行回填，導致沉降的區域有所抬升。故圖4中 35至55像元處出現明顯的抬升區域。

圖4 走向時序累計沉降

為了驗證SBAS-InSAR技術在礦區形變監測的可靠性，將GPS獲得的三維形變投影至雷達視線方向并于與SBAS-InSAR監測結果進行對照。由于GPS的觀測數據表示一個點的形變，而InSAR的監測結果表示一個像元分辨率大小的形變，當地表變形不劇烈時，可認為GPS和InSAR獲得的沉降值近似相等[19]。圖5、圖6為SBAS-InSAR技術和GPS技術在觀測點的形變序列結果，G1點GPS與InSAR獲得的沉降量最大絕對誤差、均方誤差分別為4.7 mm和3.4 mm；G2點最大絕對誤差、均方誤差分別為4.8 mm和3.6 mm，符合礦區沉降監測的要求[20]。

圖5 G1點雷達視線向SBAS-InSAR與GPS形變結果對比

圖6 G2點雷達視線向SBAS-InSAR與GPS形變結果對比

2.3 礦區沉降預測及結果分析

本文取17景哨兵1號影像，共獲取16期沉降量，將前12期數據作為訓練樣本，后4期數據作為測試樣本進行預測。采用AGA對SVR模型中C、g、ε3個參數進行優化，初始種群個數為30，最大進化次數為100，C的尋優范圍為[0，100]，g的尋優范圍為[0，0.01]，ε的尋優范圍為[0，0.1]。迭代100次后，AGA優化后的參數C=80，g=0.001，ε=0.05。為了驗證本文預測模型的可靠性，與GM、SVR、GA-SVR 3種算法進行對比，預測結果如圖7所示。

圖7 各模型預測結果

由圖7可知，在前6期的觀測中觀測點位沉降較為緩慢，此時該工作面剛開始開采，從第7期開始，該點位隨著工作面的掘進沉降速率逐漸增大。4種預測模型對于前12期觀測數據擬合較為理想，但從13期開始，GM和SVR模型在預測時出現了較大的偏差，沉降速率變化趨勢與原始值不相符，而AGA-SVR模型與原始值反映的沉降速率基本一致，能夠更好地預測礦區地表沉降。各模型預測精度如表1所示。由表1可知，GM預測結果較差，13期和15期2期預測結果的絕對誤差均大于10 mm，而SVR模型的15期預測結果絕對誤差超過15 mm，與SBAS-InSAR結果偏差較大。與未經參數優化的SVR算法相比，GA-SVR模型預測效果相對較好，但仍存在過早收斂等問題，采用自適應遺傳算法進行參數尋優的SVR算法取得了良好的預測結果，4期數據的絕對誤差均在5 mm以內，與GA-SVR算法相比，絕對誤差分別減小了37.56%、26.31%、24.73%、31.40%。

表1 各預測模型精度對比 mm

分別計算各預測模型的平均絕對誤差(mean absolute error，MAE)、均方根誤差(root mean square error，RMSE)和運行時間，如表2所示。GM和SVR模型雖然運算時間短但精度較差，MAE和RMSE均大于11 mm，GA-SVR模型的精度有所提升，但運算時間也大幅增加，可能與陷入局部最優解有關。相比其他3種模型，本文提出的AGA-SVR模型表現出更優的精度，MAE和RMSE均小于5 mm，同時運算時間相比 GA-SVR 算法有所下降，避免了遺傳算法中存在的欺騙問題，提高了參數尋優效率。

表2 各個預測模型精度對比

為了進一步驗證AGA-SVR算法的有效性，在圖4剖面自沉降邊緣至沉降中心等間距選取3個點A1、A2、A3進行預測。剖面上3點預測結果如表3所示。由表3可知，不同沉降區域的3個點位均取得了較好的預測結果，絕對誤差均在5 mm以內。

表3 剖面上所選點位預測結果 mm

從時間距離和預測精度的關系考慮，3個點的預測精度均與時間距離成反比，即隨著觀測期數的增加，各點預測精度逐漸降低。本文對第13～16期數據采用了滾動預測的方法，即先用第1～12期的觀測數據預測第13期的數據，然后將第13期的預測結果編入訓練樣本，用第2～13期的數據對下一期數據進行預測，以此類推完成所有的預測。由于第1～12期的數據均為獨立觀測獲得，而預測第14～16期數據所用的訓練樣本中均含有預測值，且越往后各個訓練樣本間的相關性越強，導致預測模型訓練結果較差，因此時間距離越短，預測的精度相對越高。

3個點的預測精度如表4所示。由于A3點的累計沉降量最大，因此A3點的相對誤差最小。各個點位的MAE和RMSE均在5 mm以內，取得了較好的預測結果，進一步驗證了AGA-SVR算法是一種較為可靠的預測方法。

表4 各個點位預測精度

3 結束語

針對傳統D-InSAR技術存在的失相干等問題，利用SBAS-InSAR技術對內蒙古李家壕煤礦某工作面進行開采沉降監測，獲得了礦區連續地表形變特征，在此基礎上使用AGA-SVR模型對礦區沉降值進行預測。監測結果表明，自2017年10月該工作面開始開采，該地區存在連續下沉現象，截至2018年4月，該工作面最大累計垂直沉降量為163 mm，最大沉降速率為360 mm/a。SBAS-InSAR技術監測礦區沉降的最大絕對誤差為4.8 mm。與其他預測模型相比，AGA-SVR模型能夠避免算法陷入局部最優解，顯著提高了預測精度，最大絕對誤差為5.0 mm，最大相對誤差為4.08%，MAE和RMSE均優于4.5 mm，監測精度和預測精度均滿足工程實踐需要。

在礦區地表形變監測中，利用InSAR技術獲取開采沉陷的影響范圍與發展趨勢，對于形變劇烈的區域，應考慮在實地布設GPS監測點并定期進行水準觀測。同時，建立預測模型，將InSAR技術提取的形變量作為訓練樣本進行開采沉陷動態預測，當預測結果的形變速率顯著增大時，應結合其他監測手段進行分析并對礦區開采活動發出預警，以保障礦區生產活動的平穩進行。下一步的工作將考慮增加數據量并采用不同極化方式的數據以提取礦區的三維形變信息。