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

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

(1.內(nèi)蒙古自治區(qū)雷達(dá)技術(shù)與應(yīng)用重點實驗室,內(nèi)蒙古 呼和浩特 010051;2.內(nèi)蒙古工業(yè)大學(xué)信息工程學(xué)院,內(nèi)蒙古 呼和浩特 010051;3.內(nèi)蒙古工業(yè)大學(xué)資源與環(huán)境工程學(xué)院,內(nèi)蒙古 呼和浩特 010051)

礦產(chǎn)資源對于國家的經(jīng)濟發(fā)展至關(guān)重要[1],不僅是社會發(fā)展的基礎(chǔ),還是衡量一個國家實力的重要指標(biāo)。隨著工業(yè)化和城市化進(jìn)程的加快,人類對煤炭資源的開采量也隨之增大。對地下煤炭過度開采容易導(dǎo)致礦區(qū)地面發(fā)生形變[2],從而破壞礦區(qū)的地表結(jié)構(gòu)和周圍環(huán)境。為保護(hù)礦區(qū)地質(zhì)環(huán)境并保證地下開采工作正常進(jìn)行,對礦區(qū)開采沉陷引發(fā)的地面變形進(jìn)行監(jiān)測并及時防治尤為重要。

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

1 融合GPS-InSAR數(shù)據(jù)的流程

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

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

(2)數(shù)據(jù)處理,對融合后的數(shù)據(jù)集進(jìn)行數(shù)據(jù)處理與統(tǒng)一。為方便后續(xù)模型的分析,需要對數(shù)據(jù)集進(jìn)行離散值的剔除。隨后統(tǒng)一GPS數(shù)據(jù)集與InSAR數(shù)據(jù)集的坐標(biāo)基準(zhǔn)將其融合成一個總的數(shù)據(jù)集。

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

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

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

2.1 克里金插值方法介紹

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

式中,為點(x0,y0)處的估計值;λi為第i個監(jiān)測點的權(quán)重系數(shù)。

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

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

球狀模型:

指數(shù)模型:

高斯模型:

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

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

圖2 球狀模型Fig.2 Spherical model

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

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

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

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

圖3 指數(shù)模型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°的空間相關(guān)性最好,其走勢基本與高斯模型一致,因此選取高斯模型的135°方位角進(jìn)行年平均形變速率的克里金插值,這樣所組合的樣本點相較于隨機選取的樣本點更符合高斯模型,其估計值也更符合實際測量值。

2.3 遺傳算法優(yōu)化克里金插值

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

在變異函數(shù)模型確定的前提下,克里金插值是否準(zhǔn)確的本質(zhì)在于函數(shù)模型參數(shù)的確定[20-21]?？死锝鸩逯底儺惡瘮?shù)模型的相關(guān)參數(shù)是通過函數(shù)擬合獲得的,因此具有較大的隨意性。采用遺傳算法確定克里金插值中的變異函數(shù)模型的相關(guān)參數(shù),即c0、c1與a的值,可以減小函數(shù)模型構(gòu)建的主觀性與偶然性。

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

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

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

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

(2)依據(jù)選取的變異函數(shù)確定相關(guān)參數(shù)c0、c1、a及約束條件,確定適應(yīng)度函數(shù)。

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

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

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

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

3.1 遺傳-克里金插值的實驗結(jié)果

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

圖10 遺傳優(yōu)化算法適應(yīng)度曲線Fig.10 Genetic optimization algorithm fitness curve

圖11 遺傳優(yōu)化算法擬合結(jié)果Fig.11 Genetic optimization algorithm fitting results

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

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

表2 模型參數(shù)對比Table 2 Comparison of model parameters

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

圖12 監(jiān)測結(jié)果對比Fig.12 Comparison of monitoring results

圖13 誤差結(jié)果對比Fig.13 Comparison of error results

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

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

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

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

4 結(jié) 論

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

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

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

(3)克里金插值模型能夠?qū)Χ嘧宰兞坑成渲g的高度非線性關(guān)系進(jìn)行很好的擬合,其擬合精度主要依賴于參數(shù)c0、c1與a的確定。通過遺傳算法優(yōu)化克里金插值半變異函數(shù)模型,最終得到高斯模型的最優(yōu)參數(shù)為c0=8.617 902、c1=1 498.095 4、a=493.357 2,利用得到的參數(shù)建立優(yōu)化模型并進(jìn)行綜合分析。通過對比InSAR、融合GPS數(shù)據(jù)的克里金插值、融合GPS數(shù)據(jù)的遺傳-克里金插值3種模型的監(jiān)測結(jié)果與誤差結(jié)果,發(fā)現(xiàn)融合GPS數(shù)據(jù)的遺傳-克里金插值模型得到的年平均形變速率與GPS實測數(shù)據(jù)最為接近,且相對誤差要明顯小于其他2種模型,說明該模型的模擬效果與精度得到了顯著的改善。