王志崗，周文韜

（1.四川科技職工大學(xué) 應(yīng)急管理系，四川 成都 610000；2.西南科技大學(xué) 環(huán)境與資源學(xué)院，四川 綿陽(yáng) 621010；3.國(guó)家遙感中心綿陽(yáng)科技城分部，四川 綿陽(yáng) 621010）

礦物被開(kāi)采后，巖層內(nèi)部應(yīng)力的平衡狀態(tài)遭到破壞，導(dǎo)致地表產(chǎn)生形變并出現(xiàn)塌陷坑[1]。地表塌陷對(duì)礦區(qū)地下安全開(kāi)采構(gòu)成威脅，也給礦區(qū)周圍環(huán)境、建筑等造成了破壞[2]。為保護(hù)地面建筑并保證地下開(kāi)采工作正常進(jìn)行，監(jiān)測(cè)地表形變十分必要。目前，全球?qū)Ш蕉ㄎ幌到y(tǒng)（GPS）技術(shù)、合成孔徑雷達(dá)干涉測(cè)量（InSAR）技術(shù)已被廣泛應(yīng)用于礦區(qū)地表形變監(jiān)測(cè)中，GPS技術(shù)具有高時(shí)間分辨率和高精度[3]，InSAR技術(shù)具有監(jiān)測(cè)范圍廣、不受天氣影響等優(yōu)點(diǎn)[4]；但這兩種技術(shù)在實(shí)際應(yīng)用中均有不足，即GPS技術(shù)僅能反映離散點(diǎn)的監(jiān)測(cè)結(jié)果，InSAR技術(shù)僅能反映視線向的相對(duì)形變量，若僅將視線向形變投影到垂直方向，該形變往往不能客觀反映實(shí)際形變。王霞迎[5]、曹海坤[6]等基于升降軌InSAR與GPS數(shù)據(jù)，建立了融合函數(shù)模型反演高精度三維形變場(chǎng)；馬艷鴿[7]、丁寧[8]和JI P F[9]等基于先驗(yàn)算法改進(jìn)融合GPS與InSAR的函數(shù)模型，以提高形變監(jiān)測(cè)精度。眾多學(xué)者均是通過(guò)建立融合函數(shù)模型，利用GPS數(shù)據(jù)校正InSAR數(shù)據(jù)來(lái)提高InSAR形變精度[10-12]。然而，當(dāng)?shù)匦吻闆r復(fù)雜時(shí)（如礦區(qū)、冰川等），GPS點(diǎn)的形變不能反映其周圍面的形變，具有較高的偶然性；InSAR技術(shù)的地面分辨率較低，卻能很好地反映面的形變。因此，如何通過(guò)加權(quán)融合，將降低GPS點(diǎn)偶然性和提高InSAR形變監(jiān)測(cè)精度進(jìn)行統(tǒng)一成為亟待解決的問(wèn)題。

鑒于此，本文重點(diǎn)研究地表垂直方向的形變，介紹了不同內(nèi)插GPS點(diǎn)的形變結(jié)果與優(yōu)勢(shì)，再利用變異系數(shù)法（CV）得到組合內(nèi)插GPS形變結(jié)果，最后運(yùn)用算術(shù)平均法將組合內(nèi)插GPS形變結(jié)果與SBAS-InSAR垂直向累積形變結(jié)果進(jìn)行加權(quán)融合，充分發(fā)揮二者優(yōu)勢(shì)，得到GPS-InSAR形變結(jié)果。

1 理論基礎(chǔ)

1.1 空間插值方法

經(jīng)驗(yàn)貝葉斯克里金插值法（EBK）是一種地統(tǒng)計(jì)插值方法，可通過(guò)構(gòu)造子集和模擬的方式自動(dòng)計(jì)算構(gòu)建有效克里金模型過(guò)程中的參數(shù)。與普通克里金插值[13]算法不同，EBK法考慮了半變異函數(shù)估計(jì)的不確定性，可通過(guò)估計(jì)半變異函數(shù)來(lái)說(shuō)明引入的誤差，因此EBK法降低了預(yù)測(cè)的標(biāo)準(zhǔn)誤差[14]。

局部多項(xiàng)式插值法（LPI）是一種局部加權(quán)最小二乘擬合法，根據(jù)有限的監(jiān)測(cè)數(shù)據(jù)，采用多個(gè)多項(xiàng)式來(lái)擬合表面。它引入了“距離權(quán)”的概念，對(duì)于未知點(diǎn)的計(jì)算，考慮在局部范圍內(nèi)所有已知點(diǎn)對(duì)其的貢獻(xiàn)，距未知點(diǎn)近的點(diǎn)權(quán)重大，反之則權(quán)重小。每個(gè)未知點(diǎn)的預(yù)測(cè)值都對(duì)應(yīng)一個(gè)多項(xiàng)式，每個(gè)多項(xiàng)式都處于特定重疊的鄰近區(qū)域內(nèi)，通過(guò)最小二乘法求解鄰域內(nèi)多項(xiàng)式組成的方程組，從而得到擬合表面[15]。該方法不僅具有趨勢(shì)面法考慮全部數(shù)據(jù)點(diǎn)反映趨勢(shì)性變化的優(yōu)點(diǎn)，而且具有距離法反映局部特征的優(yōu)點(diǎn)。

徑向基函數(shù)插值法（RBF）是一種人工神經(jīng)網(wǎng)絡(luò)方法，根據(jù)有限的監(jiān)測(cè)數(shù)據(jù)，選擇合適的徑向基函數(shù)生成一個(gè)具有最小曲率、且到各樣點(diǎn)的Z值距離最小的曲面。該方法擬合的表面經(jīng)過(guò)所有點(diǎn)數(shù)據(jù)，并可計(jì)算得到高于或低于點(diǎn)Z值的預(yù)測(cè)值，適用于監(jiān)測(cè)點(diǎn)數(shù)據(jù)集大、表面變化平緩的情況[16]。

1.2 CV法

CV法是一種客觀計(jì)算權(quán)重的方法，可直接利用各項(xiàng)評(píng)價(jià)因子包含的信息，通過(guò)計(jì)算得到各項(xiàng)指標(biāo)的權(quán)重。為消除各指標(biāo)量綱不同帶來(lái)的影響，需利用各指標(biāo)的變異系數(shù)來(lái)衡量其取值的差異程度[17]。其計(jì)算公式為：

式中，Vi為第i項(xiàng)指標(biāo)的變異系數(shù)，亦稱標(biāo)準(zhǔn)差系數(shù)；Si為第i項(xiàng)指標(biāo)的標(biāo)準(zhǔn)差；為第i項(xiàng)指標(biāo)的平均值。

根據(jù)式（1）可以得到各項(xiàng)指標(biāo)的權(quán)重，即

1.3 精度評(píng)定指標(biāo)

本文通過(guò)比較GPS點(diǎn)預(yù)測(cè)值與實(shí)測(cè)值的均方根誤差（RMSE）、平均絕對(duì)誤差（MAE）和擬合優(yōu)度（R2）3個(gè)指標(biāo)來(lái)評(píng)定插值精度。其計(jì)算公式為：

式中，xi為第i點(diǎn)的實(shí)測(cè)值；為第i點(diǎn)的預(yù)測(cè)值；為實(shí)測(cè)值x的平均值。RMSE和MAE值越小、R2越趨近1，表示預(yù)測(cè)精度越高。

2 礦區(qū)地表形變監(jiān)測(cè)實(shí)驗(yàn)分析

本文以金川某礦西二采區(qū)為研究區(qū)（以下簡(jiǎn)稱西二采），位于我國(guó)甘肅省金昌市中部。礦區(qū)地勢(shì)較平坦，平均海拔約為1 500～1 800 m；地表水系不發(fā)育，常年干枯，屬溫帶大陸性氣候[18]。近年來(lái)，由于地下不斷開(kāi)采，導(dǎo)致地表塌陷、裂縫明顯，為銅鎳礦資源開(kāi)采留下了極大的安全隱患。

2.1 數(shù)據(jù)基礎(chǔ)

在研究區(qū)范圍內(nèi)，布設(shè)GPS點(diǎn)進(jìn)行靜態(tài)觀測(cè)，監(jiān)測(cè)時(shí)間為2019年4月—2020年6月，監(jiān)測(cè)周期為60 d，共進(jìn)行8期監(jiān)測(cè)。為與GPS監(jiān)測(cè)時(shí)間統(tǒng)一，實(shí)驗(yàn)選取覆蓋研究區(qū)的38景C波段Sentinel-1A升軌影像，時(shí)間跨度為2019年3月22日—2020年6月8日，設(shè)置臨界基線閾值為45%，時(shí)間基線閾值為120 d，生成325個(gè)干涉對(duì)，利用小基線（SBAS）方法提取垂直向累積形變結(jié)果。Sentinel-1A升軌影像覆蓋范圍和GPS點(diǎn)位分布情況如圖1所示。

圖1 西二采影像覆蓋范圍和GPS點(diǎn)位分布圖

2.2 地表形變監(jiān)測(cè)實(shí)驗(yàn)與分析

2.2.1 基于SBAS-InSAR技術(shù)的地表形變監(jiān)測(cè)

本文利用SBAS-InSAR技術(shù)獲取了研究區(qū)2019年3月22日—2020年6月8日的地表垂直向累積形變量，如圖2所示，可以看出，礦區(qū)地表沉降區(qū)域主要分布在礦體西南部，該區(qū)域受地下開(kāi)采影響較大；在15個(gè)月的時(shí)間里，礦區(qū)內(nèi)最大沉降量達(dá)到-168.41 mm，地表整體呈下降趨勢(shì)，沉降面積約為0.267 km2。

圖2 SBAS-InSAR垂直向累積形變量

2.2.2 基于CV的三插GPS形變實(shí)驗(yàn)與分析

為融合GPS數(shù)據(jù)與InSAR數(shù)據(jù)，需先將GPS點(diǎn)元數(shù)據(jù)內(nèi)插為面元數(shù)據(jù)。EBK法克服了傳統(tǒng)克里金法需人機(jī)交互式建模來(lái)尋找準(zhǔn)確結(jié)果的缺點(diǎn)，提高了預(yù)測(cè)精度，并可準(zhǔn)確預(yù)測(cè)不夠穩(wěn)定的數(shù)據(jù)；但處理時(shí)間隨監(jiān)測(cè)點(diǎn)數(shù)、子集大小或重疊系數(shù)的增加而快速增加，降低了計(jì)算效率，且對(duì)于含有異常值的數(shù)據(jù)，可能會(huì)得到大于或小于輸入點(diǎn)值若干個(gè)數(shù)量級(jí)的預(yù)測(cè)結(jié)果，從而影響插值結(jié)果。因此，本文結(jié)合兩種精確插值方法來(lái)提高運(yùn)算效率，降低異常數(shù)據(jù)對(duì)插值結(jié)果的影響。LPI法突出點(diǎn)數(shù)據(jù)反映趨勢(shì)性變化，同時(shí)通過(guò)距離法強(qiáng)調(diào)局部特征。RBF法可預(yù)測(cè)大于最大測(cè)量值和小于最小測(cè)量值的值，適用于表面變化平緩的情況。上述3種插值方法的結(jié)果如圖3～5所示，其中EBK法插值后的最大累積沉降量為-320.88 mm，最大抬升量為123.47 mm；LPI法插值后的最大累積沉降量為-428.49 mm，最大抬升量為565.24 mm；RBF法插值后的最大累積沉降量為-368.34 mm，最大抬升量為264.89 mm；說(shuō)明3種插值方法內(nèi)插后的面數(shù)據(jù)整體形變趨勢(shì)較一致，3～8行的沉降量大，9～12行出現(xiàn)輕微抬升。本文綜合3種插值方法的優(yōu)點(diǎn)，根據(jù)式（2）得到EBK-LPI-RBF法的插值結(jié)果，如圖6所示。3種插值方法的權(quán)重分配如表1所示。

表1 CV法權(quán)重分配

圖3 EBK法插值結(jié)果

圖4 LPI法插值結(jié)果

圖5 RBF法插值結(jié)果

圖6 EBK-LPI-RBF法插值結(jié)果

不同插值方法的精度比較如表2所示，可以看出，LPI法和RBF法的精度略高于EBK法，基于CV法的EBK-LPI-RBF法 的RMSE為8.68 mm，MAE為5.52 mm，R2為0.99，滿足精度要求。

表2 不同插值方法的精度比較

2.2.3 加權(quán)融合GPS-InSAR的地表形變實(shí)驗(yàn)與分析

加權(quán)融合的方法有很多，如變權(quán)組合法、算術(shù)平均法（等權(quán)組合法）、非線性組合法、最優(yōu)加權(quán)法等，本文為實(shí)現(xiàn)降低GPS點(diǎn)偶然性和提高InSAR形變監(jiān)測(cè)精度的統(tǒng)一，采用算術(shù)平均法融合GPS與InSAR數(shù)據(jù)獲取礦區(qū)地表形變場(chǎng)。根據(jù)EBK-LPI-RBF法和InSAR形變結(jié)果圖，統(tǒng)一像元分辨率（16.9×16.9）和坐標(biāo)基準(zhǔn)，得到加權(quán)融合的GPS-InSAR形變結(jié)果，如圖7所示，可以看出，沉降區(qū)域分布在礦體西南部，與InSAR和GPS監(jiān)測(cè)結(jié)果保持一致。

圖7 GPS-InSAR累積形變結(jié)果

為驗(yàn)證本文提出方法的有效性，分別采用兩種方案進(jìn)行精度驗(yàn)證：①方案1，將GPS監(jiān)測(cè)值作為真值，分析InSAR、EBK-LPI-RBF和GPS-InSAR三種方法在垂直向的形變精度；②方案2，將InSAR監(jiān)測(cè)值作為真值，分析GPS、EBK-LPI-RBF和GPS-InSAR三種方法在垂直向的形變精度，結(jié)果如表3所示。

通過(guò)分析表3和諸圖可以看出：

表3 不同研究方法的精度比較

1）方案1將GPS監(jiān)測(cè)值作為真值，EBK-LPI-RBF法的內(nèi)插結(jié)果精度達(dá)到毫米級(jí)，RMSE為8.68 mm，MAE為5.52 mm，R2高達(dá)0.99；但I(xiàn)nSAR形變結(jié)果的精度較低，其RMSE為73.28 mm，R2僅為0.45。方案2將InSAR監(jiān)測(cè)值作為真值，EBK-LPI-RBF法的內(nèi)插結(jié)果精度略優(yōu)于GPS監(jiān)測(cè)結(jié)果。

2）對(duì)比兩個(gè)方案發(fā)現(xiàn)，雖然GPS-InSAR法精度低于EBK-LPI-RBF法，但卻很好地降低了離散點(diǎn)的偶然性，提升了InSAR的形變監(jiān)測(cè)精度。兩種方案中GPS-InSAR法的R2分別為0.88和0.77，很好地說(shuō)明了該方法的可靠性。

3 結(jié)語(yǔ)

傳統(tǒng)GPS與InSAR數(shù)據(jù)融合僅依據(jù)融合函數(shù)模型來(lái)提升InSAR形變精度，未考慮兩種數(shù)據(jù)融合時(shí)本身性質(zhì)不同所造成的影響。鑒于此，本文首先基于SBAS-InSAR技術(shù)獲取研究區(qū)地表2019年4月—2020年6月的形變場(chǎng)，提取垂直向累積形變量；再針對(duì)GPS離散點(diǎn)形變反映面形變具有偶然性和InSAR監(jiān)測(cè)不能客觀反映實(shí)際形變的問(wèn)題，綜合了3種插值方法的優(yōu)點(diǎn)，利用基于CV法的EBK-LPI-RBF法內(nèi)插GPS點(diǎn)數(shù)據(jù)，并計(jì)算得到加權(quán)融合的GPS-InSAR地表形變結(jié)果。由精度評(píng)定指標(biāo)可知，該方法有效降低了GPS點(diǎn)數(shù)據(jù)的偶然性，提高了InSAR監(jiān)測(cè)精度，監(jiān)測(cè)結(jié)果滿足精度要求，為GPS與InSAR數(shù)據(jù)融合提供了新方法。