面向礦區地表形變的GPS-InSAR融合方法研究

陳月娟 栗向蕾 黃平平 蔡秋歡 譚維賢 尹 博

(1.內蒙古自治區雷達技術與應用重點實驗室,內蒙古 呼和浩特 010051;2.內蒙古工業大學信息工程學院,內蒙古 呼和浩特 010051;3.內蒙古工業大學資源與環境工程學院,內蒙古 呼和浩特 010051)

礦產資源對于國家的經濟發展至關重要[1],不僅是社會發展的基礎,還是衡量一個國家實力的重要指標。隨著工業化和城市化進程的加快,人類對煤炭資源的開采量也隨之增大。對地下煤炭過度開采容易導致礦區地面發生形變[2],從而破壞礦區的地表結構和周圍環境。為保護礦區地質環境并保證地下開采工作正常進行,對礦區開采沉陷引發的地面變形進行監測并及時防治尤為重要。

目前,全球定位系統(Global Positioning System,GPS)與合成孔徑雷達干涉(Synthetic Aperture Radar Interferometry,InSAR)測量技術已被廣泛應用于礦區地表形變監測中[3]。GPS技術具有高時間分辨率和高精度,可以實現全天候的連續監測[4]。但其空間分辨率低,且難以實現大范圍的形變情況監測。In-SAR技術具有監測范圍廣、不受天氣影響等特點,但該技術的測量主要依靠衛星環繞地球的周期,因此其時間分辨率較低。GPS技術與InSAR技術具有良好的互補性[5],將GPS與InSAR兩種技術得到的數據進行融合不僅可以獲取監測區完整的形變情況,同時可以實現對特殊點的有效監測,進一步提高對地監測的精度[6]。馬艷鴿等[7]提出利用函數模型與隨機模型雙約束的GPS-InSAR融合技術,該模型主要利用了GPS監測技術高分辨率、高精度的優勢,將GPS監測結果作為先驗信息,通過添加隨機模型約束,建立新的三維變形模型。然而對于地形復雜的區域,GPS點不能很好地反映其周圍面的形變,具有較高的偶然性。隨后,高晨等[8]利用克里金(Kriging)插值算法的線性變異函數模型對GPS數據進行插值,實現了對GPS速度場與InSAR技術的雷達視線向形變速率的綜合求解,并解決了信息缺失的問題。但克里金插值算法在確定模型參數時具有較強的偶然性,從而影響監測精度。因此,如何對GPS與InSAR數據進行融合并提高監測的精度對礦山地表形變引發的災害預防至關重要[9]。本文對GPS與InSAR數據進行融合,探究融合對地表形變的監測精度的影響。

1 融合GPS-InSAR數據的流程

GPS技術對地監測易于實現且操作簡單,具有非常高的時間分辨率。但針對地表變化復雜、高度較大的區域,衛星信號反饋的數據存在一定誤差,其監測精度也很難滿足地面沉降變形研究工作的需求。In-SAR技術能提供整個區域連續的形變信息,但所提供圖像具有一定的周期性,因此其時間分辨率較低[10]。另外InSAR技術在進行地表形變監測時,由于大氣時延與空間失相干的影響,在形變梯度較大的區域容易出現信息缺失的情況[11]。針對以上問題,融合GPS與InSAR數據并進行克里金插值,可進一步掌握研究區完整的形變信息。如圖1所示,GPS與InSAR數據融合的具體步驟如下。

(1)數據準備,融合GPS與InSAR監測數據集。由于GPS技術所獲取的數據為研究區的高程值,In-SAR技術可以直接監測研究區年平均沉降速率,因此在融合數據之前需要將相同(相近)時間段內的GPS高程觀測值轉為年平均沉降速率。

(2)數據處理,對融合后的數據集進行數據處理與統一。為方便后續模型的分析,需要對數據集進行離散值的剔除。隨后統一GPS數據集與InSAR數據集的坐標基準將其融合成一個總的數據集。

(3)插值處理,其核心為模型搭建,計算總數據集中所有觀測數據之間的分離距離與半方差。依據距離-半方差圖確定克里金插值算法的最佳變異函數模型,并依據任意空間距離計算出相應的半方差。計算所有已知點的半方差并得到系數矩陣,求解方程組得權重系數,通過對已知點生成的權重進行加權求和預測未知點的估計值。

(4)效果評價,將未參與融合的GPS點作為檢驗點,分別與該點的InSAR監測結果、融合GPS-InSAR數據的插值結果進行對比,對融合模型的精度進行評價。

2 融合GPS-InSAR插值模型的建立

2.1 克里金插值方法介紹

對于克里金插值是在考慮研究變量的空間相關性與空間屬性變異分布的基礎上,通過已知數據對未知數據進行最優無偏估計的空間自協方差內插。假設研究區共有n個監測點,待插點為z0=z(x0,y0),待估點的值可表示為空間周圍觀測點的線性組合?？死锝鸩逯捣?

式中,為點(x0,y0)處的估計值;λi為第i個監測點的權重系數。

2.2 克里金插值半變異函數的分析與選取

克里金插值的關鍵在于變異函數的選取[12-14],因此在進行克里金插值時需要有依據地確定變異函數模型及其相關參數。本研究主要選用球狀模型、高斯模型和指數模型3種半變異函數模型作為理論變異函數對插值效果進行分析,從而選擇擬合度較好的半變異函數模型。變異函數的半方差函數如下。

球狀模型:

指數模型:

高斯模型:

式中,γ(h)為半方差值;c0為塊金值;c1為偏基臺值;a為變程。

采用3種模型對處理后的GPS-InSAR數據進行研究區域半變異函數各向同性曲線的繪制,并通過擬合優度R2與殘差平方和RSS選擇擬合度最好的變異函數模型[15]。以分離距離為橫坐標,以半方差為縱坐標繪制年均形變速率各向同性曲線。圖2～圖4為3種模型的繪制結果。

圖2 球狀模型Fig.2 Spherical model

半變異函數在度量空間自相關性時,反映的是全局自相關程度。為進一步確定半變異函數模型,以R2與RSS作為評價指標,對3種模型的擬合精度進行分析,對比結果如表1所示。

表1 不同模型的精度對比Table 1 Accuracy comparison of different models

R2在回歸分析中可用于衡量回歸方程對樣本觀測點的整體擬合程度。R2的值越接近1,說明回歸直線對樣本觀測點的擬合程度越好。RSS是一種用于評價擬合數據質量的指標,用于衡量擬合曲線的泛化能力。因此,RSS的值越小,說明擬合曲線與實際數據之間的偏差越小,即擬合數據的質量越高。表1結果表明:在擬合優度方面,高斯模型的擬合優度最好,球狀模型次之,指數模型的擬合優度最差。在殘差平方和方面,高斯模型的殘差平方和明顯小于其他2種模型,即擬合數據的質量最高。綜上所述,在進行融合GPS與InSAR數據的克里金插值算法中確定高斯模型為變異函數模型。

通過圖2、圖3、圖4所繪制的研究區域年均形變速率各向同性曲線圖可確定變異函數模型,之后需要通過年均形變速率的各向異性曲線選擇擬合的方位角進行研究區域的切分[16]。高斯模型的各向異性曲線如圖5～圖8所示。

圖3 指數模型Fig.3 Exponential model

圖4 高斯模型Fig.4 Gaussian model

圖5 0°方位角Fig.5 0° azimuth angle

圖6 45°方位角Fig.6 45° azimuth angle

圖8 135°方位角Fig.8 135° azimuth angle

由圖5～圖8可以看出,年均形變速率的各向異性是不同的。從0°、45°、90°、135°4個方位來看,135°的空間相關性最好,其走勢基本與高斯模型一致,因此選取高斯模型的135°方位角進行年平均形變速率的克里金插值,這樣所組合的樣本點相較于隨機選取的樣本點更符合高斯模型,其估計值也更符合實際測量值。

2.3 遺傳算法優化克里金插值

遺傳算法是模仿遺傳學自然選擇進化過程建立起來的一種概率搜索算法[17-19],在非線性、多模型、多目標的函數優化方面具有較強的優勢。遺傳算法的基本思想是通過優勝劣汰的方式對基因進行選擇、交叉和變異等遺傳操作,逐步對問題進行優化,最終確定全局最優解。

在變異函數模型確定的前提下,克里金插值是否準確的本質在于函數模型參數的確定[20-21]。克里金插值變異函數模型的相關參數是通過函數擬合獲得的,因此具有較大的隨意性。采用遺傳算法確定克里金插值中的變異函數模型的相關參數,即c0、c1與a的值,可以減小函數模型構建的主觀性與偶然性。

考慮到變異函數模型最優擬合的目的是使實際變異函數與理論變異函數更為接近,因此在遺傳算法尋優過程中,將各種距離條件下實際變異函數與理論變異函數的差值之和作為目標函數[22]。遺傳優化算法中的目標函數:

圖9為遺傳算法優化克里金插值模型的流程圖,遺傳算法優化克里金插值方法的具體步驟如下:

圖9 遺傳-克里金算法流程Fig.9 Genetic-Kriging algorithm process

(1)載入GPS-InSAR數據集。

(2)依據選取的變異函數確定相關參數c0、c1、a及約束條件,確定適應度函數。

(3)初始化種群,確定個體適應度的量化評價方法并計算個體適應度值。種群大小決定了算法的搜索空間。群體規模太小,會出現近親交配產生病態基因;太大會導致結果難以收斂且穩健性下降。種群大小一般設置為50～300,這樣既保證了種群中有足夠的個體可以參與優化也不會浪費資源。

(4)進行個體選擇、交叉、變異等遺傳運算,確定遺傳算法的相關運行參數。變異概率決定了每一代進行變異的個體數,變異概率太大容易導致丟失有效基因,太小會導致過早收斂,變異概率一般設置為0.01～0.1。交叉概率決定了參與交叉的個體數量,交叉概率過大會導致優良的個體在交叉過程中被過度破壞,個體的優秀特性難以傳遞給下一代,過小則不能及時對種群進行更新,其一般設置為0.6～0.9。

(5)判斷是否滿足精度要求或達到最大迭代次數,輸出變異函數模型的最優參數,擬合出最優變異函數模型。

3 遺傳-克里金插值的實驗結果及驗證分析

3.1 遺傳-克里金插值的實驗結果

首先基于GPS-InSAR數據計算采樣點之間的協方差與實際變異函數值,根據數據集的大小對遺傳算法控制參數進行初始化。本研究參與遺傳運算的個體規模較大,為保證較快獲取最優解,設置遺傳算法的種群規模為300,變異概率為0.1,交叉概率為0.9,最大迭代次數為100次,迭代過程中目標函數最優值隨迭代次數的變化如圖10所示?？梢钥闯鲭S著迭代次數的增加,目標函數下降的速率逐漸趨于平緩,迭代次數達到15次時就已經得到最優值?；谶z傳算法擬合得到的理論變異函數與實際變異函數如圖11所示。可以看出理論變異函數與實際變異函數值的曲線走向一致且誤差較小,表明改進模型的變異函數擬合效果很好。

圖10 遺傳優化算法適應度曲線Fig.10 Genetic optimization algorithm fitness curve

圖11 遺傳優化算法擬合結果Fig.11 Genetic optimization algorithm fitting results

3.2 遺傳-克里金插值的驗證分析

由上文確定選用的變異函數模型為高斯模型,依據融合后的GPS-InSAR數據集,通過遺傳算法的全局尋優實現對克里金插值模型相關參數的解算,表2為克里金插值模型與優化模型解算出的c0、c1與a的值。

表2 模型參數對比Table 2 Comparison of model parameters

為驗證融合GPS與InSAR數據集插值模型在地表監測方面相對于InSAR監測技術的優越性,將高斯模型及表2中求解的2種模型的相關參數作為理論變異函數模型進行仿真實驗。未參與融合的12個GPS站點的年平均形變速率值作為核驗數據,分別與InSAR數據、融合后的克里金插值結果、融合后的遺傳-克里金插值結果進行對比。圖12為年均形變速率值對比曲線,圖13為誤差對比情況。

圖12 監測結果對比Fig.12 Comparison of monitoring results

圖13 誤差結果對比Fig.13 Comparison of error results

由圖12可知:InSAR監測數據與GPS實測值的曲線走勢一致,說明InSAR技術能很好地對地表形變進行大范圍的監測,但部分監測點的數據不穩定且與實測值存在較大的差距,因此InSAR技術對點測量方面精度不高。融合GPS-InSAR數據集并進行插值處理,其獲得的監測結果與實測值較為接近,說明對2種數據進行融合可以很好地對2種技術進行互補,從而提高監測的精度。但仍存在部分點偏離實測值的情況,這是因為普通克里金插值算法在進行變異函數模型相關參數的選擇方面存在主觀性。經遺傳算法優化后的克里金插值結果與實測值曲線發展趨勢一致且擬合度最高,表明該模型的效果最好,在礦區地面變形監測領域更加適用。

由圖13可知:將3種模型的監測結果與GPS實際觀測值之間的誤差進行對比,說明InSAR技術監測結果與GPS實測數據之間的誤差較大。融合GPS數據的克里金插值模型精度有所提高但仍有偏差較大的情況,說明內插的精度不夠理想。經遺傳算法優化后的克里金插值模型的結果與GPS實測值之間的誤差最小,其相對誤差曲線基本在零基準線附近擺動,說明其監測結果與真實值較為接近。由此可見,該模型在對地面形變速率監測領域具有很好的實用性。

為了定量地對3種模型的監測精度進行評價,利用平均絕對誤差(MAE)、均方誤差(MSE)、均方根誤差(RMSE)3個評價指標進行對比分析。表3為3種模型的性能比較結果。結果表明:遺傳-克里金插值模型在3個指標上的計算結果明顯小于其他2種模型,說明該模型的整體性能最好。

表3 模型性能比較Table 3 Comparison of model performance

4 結 論

針對GPS技術與InSAR技術在進行地表形變監測方面存在的不足,提出了融合GPS-InSAR數據的遺傳-克里金插值方法,利用未參與插值的GPS數據對插值結果進行驗證,可以得到以下結論:

(1)分析克里金插值方法的原理及GPS與In-SAR數據融合的插值流程,對克里金插值算法中球狀模型、指數模型、高斯模型3種模型的半變異函數進行各向同性分析。通過對比3種函數的擬合優度R2與殘差平方和RSS,選定使用高斯模型進行插值,隨后,對高斯模型在0°、45°、90°、135°4個方位角上的各向異性曲線進行分析,依據曲線擬合程度選用135°方位角劃分研究區域,從而解決InSAR技術對點監測精度不高的問題,完整地反映研究區域的形變情況。

(2)比較InSAR技術、融合GPS數據的克里金插值得到的年平均形變速率,發現融合GPS數據的插值結果與GPS站點的實測結果較為接近。對比2種數據的誤差曲線,并將MAE、MSE、RMSE作為精度評定指標,發現融合GPS數據的克里金插值模型得到的形變速率結果相對誤差要明顯低于InSAR技術,證明經克里金插值后的精度更高,但相較于實際的測量結果某些點位還是存在較大的誤差,因此該模型還需要進一步的完善。

(3)克里金插值模型能夠對多自變量映射之間的高度非線性關系進行很好的擬合,其擬合精度主要依賴于參數c0、c1與a的確定。通過遺傳算法優化克里金插值半變異函數模型,最終得到高斯模型的最優參數為c0=8.617 902、c1=1 498.095 4、a=493.357 2,利用得到的參數建立優化模型并進行綜合分析。通過對比InSAR、融合GPS數據的克里金插值、融合GPS數據的遺傳-克里金插值3種模型的監測結果與誤差結果,發現融合GPS數據的遺傳-克里金插值模型得到的年平均形變速率與GPS實測數據最為接近,且相對誤差要明顯小于其他2種模型,說明該模型的模擬效果與精度得到了顯著的改善。