周文韜 張文君 繆駿懿 申 銳 訾應(yīng)昆

1 西南科技大學(xué)環(huán)境與資源學(xué)院，四川省綿陽市青龍大道中段59號，621010 2 國家遙感中心綿陽科技城分部，四川省綿陽市青龍大道中段59號，621010 3 四川航遙智測科技有限公司，四川省綿陽市涪金路389號，621010

隨著測繪技術(shù)的不斷更新，礦區(qū)地表形變監(jiān)測方法取得了很大進展[1-2]。傳統(tǒng)GNSS測量方法可得到高精度的單點三維形變數(shù)據(jù)，但僅依靠單點測量無法客觀反映地表真實形變情況[3]。合成孔徑雷達干涉測量(InSAR)技術(shù)具有高空間分辨率、高自動化和廣域監(jiān)測等優(yōu)點，為形變監(jiān)測領(lǐng)域提供了前所未有的機遇[4-5]。然而，InSAR僅能監(jiān)測沿衛(wèi)星視距方向的形變，且由于合成孔徑雷達(SAR)衛(wèi)星近南北向飛行，導(dǎo)致其對南北向形變不敏感[6]。將GNSS與InSAR數(shù)據(jù)融合解算三維形變場是目前國內(nèi)外眾多學(xué)者關(guān)注和研究的熱點[7-9]。目前的融合方法主要分為2種：1)將GNSS南北向形變插值后代入解算InSAR在東西向和垂直向的形變場[10]，但該方法過于依賴GNSS在南北向的形變精度，忽略了InSAR在南北向形變的貢獻；2)通過建立合適的融合模型解算三維形變場[11]，但會遇到先驗方差不準(zhǔn)確及迭代負方差的情況，導(dǎo)致不能完整地表達真實地表三維變化。

本文以金昌市金川西二采礦區(qū)為研究對象，提出一種融合赫爾默特方差分量估計(Helmert variance component estimation, HVCE)和徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)(radial basis function neural network, RBFNN)的三維形變計算法HVCE-RBFNN。利用該方法對同期的GNSS和InSAR數(shù)據(jù)進行迭代定權(quán)，解算出地表三維形變場，并與地下采空區(qū)進行對比分析，以期為礦區(qū)地表三維形變監(jiān)測提供新方法。

1 HVCE-RBFNN模型建立

1.1 InSAR三維幾何分解原理

傳統(tǒng)InSAR技術(shù)僅能獲取沿衛(wèi)星視線(LOS)向的一維形變，然而真實地表形變通常是三維的，因此需要將InSAR單一方向的形變分解為三維方向的形變。根據(jù)衛(wèi)星側(cè)視觀測幾何，沿衛(wèi)星LOS向的形變可分解為地表某一點在東西向、南北向和垂直向的位移[12]，如圖1所示。InSAR的累積形變DLOS可以表示為：

圖1 InSAR三維幾何分解原理Fig.1 3D geometric decomposition schematic of InSAR

DLOS=-DEsinθcosα+DNsinθsinα+DUcosθ

(1)

式中，DE、DN和DU分別為地面點P沿東西向、南北向和垂直向的位移；θ為衛(wèi)星入射角；α為衛(wèi)星方位角，以順時針方向為正。令S=[SESNSU]T、SE=-sinθcosα、SN=sinθsinα、SU=cosθ分別為衛(wèi)星LOS向形變在三維方向的投影。

1.2 HVCE-RBFNN模型

假設(shè)有n組觀測數(shù)據(jù)，且各組數(shù)據(jù)之間互不相關(guān)，根據(jù)間接平差的數(shù)學(xué)模型有：

Vi=Bi-li,i=1,2,…,n

(2)

式中，Vi為改正數(shù)矢量；Bi為系數(shù)矩陣；li為觀測數(shù)據(jù)。法方程表達式為：

=N-1W

(3)

式中，Pi為權(quán)重矩陣，在第1次最小二乘平差中可定義為任意矩陣[13]。

通常情況下，由于很難準(zhǔn)確計算觀測數(shù)據(jù)的先驗方差，因此在第1次平差中給定的權(quán)陣Pi也是不合理的。為此引入HVCE法，設(shè)各組觀測數(shù)據(jù)的單位權(quán)方差為，則有：

=A-1Wθ

(6)

式中，ai為各組觀測數(shù)據(jù)li的觀測量個數(shù)。根據(jù)式(7)計算新的權(quán)重矩陣為：

式中，c為任意常數(shù)。將求得的新權(quán)陣i代入式(2)中計算，直至滿足：

當(dāng)滿足式(11)時，可輸出計算結(jié)果。各類觀測數(shù)據(jù)差異較大時，往往會導(dǎo)致式(11)出現(xiàn)負方差的情況，這是因為觀測方程系數(shù)矩陣秩虧，導(dǎo)致矩陣N病態(tài)或完全秩虧，使得N-1→∞，N-1與法矩陣Ni的乘積的跡tr(N-1Ni)→∞，從而導(dǎo)致迭代出現(xiàn)負方差?；诖?，采用RBFNN算法對負方差的點進行計算。

RBFNN具有自主學(xué)習(xí)、自主組合和自主適應(yīng)等特點，可對差異較大的數(shù)據(jù)進行訓(xùn)練，從而達到數(shù)據(jù)融合的目的，不僅解決了計算效率的問題，還可完整地表達各組數(shù)據(jù)在融合中的貢獻[14]。RBFNN由3層前向網(wǎng)絡(luò)構(gòu)成，第1層為輸入層，第2層為隱含層，第3層為輸出層，其數(shù)學(xué)模型表示為：

(12)

式中，wij為權(quán)值；X為訓(xùn)練樣本；xi為基函數(shù)中心；φ(‖X-xi‖2)為基函數(shù)；k為輸出單元數(shù)。

RBF函數(shù)中心確定的方法不同，RBFNN的學(xué)習(xí)策略也不同。根據(jù)各組觀測數(shù)據(jù)的特點，采用隨機選取固定中心的學(xué)習(xí)策略，使基函數(shù)中心和標(biāo)準(zhǔn)差恒定不變。當(dāng)各組數(shù)據(jù)比較典型、具有代表性時，這種策略的學(xué)習(xí)效率會大幅提升。傳遞函數(shù)選擇高斯分布函數(shù)：

式中，r為模糊半徑；σ為基函數(shù)的標(biāo)準(zhǔn)差。為防止RBF函數(shù)過尖或過平，對σ進行選取：

式中，dmax為選取中心之間的距離；n為隱含節(jié)點個數(shù)。得到基函數(shù)為：

i=1,2,…,n

(15)

最后采用偽逆法計算權(quán)值：

ωij=φ(‖X-xi‖2)dkj

(16)

式中，dkj為期望輸出值。

2 實驗分析

2.1 研究區(qū)概況及數(shù)據(jù)介紹

金川西二采礦區(qū)位于甘肅省金昌市(圖2)，地處河西走廊中段、祁連山北麓，平均海拔1 700 ～2 700 m。該區(qū)屬大陸性溫帶干旱氣候，光照充足，氣候干燥，降水少且蒸發(fā)強，地下水系不發(fā)育，總體生態(tài)環(huán)境較弱[15]。礦區(qū)地下開采范圍南北長420 m，東西寬230 m，總面積約43 512 m2。由于長期受地下開采和斷層活動影響，礦區(qū)地表塌陷明顯，安全隱患極大。

(a)礦區(qū)地理位置；(b)實驗SAR影像覆蓋情況；(c)監(jiān)測點位分布；(d)地面監(jiān)測樁圖2 研究區(qū)概況Fig.2 Overview of the study area

本研究在地表均勻布設(shè)139個監(jiān)測點，得到2019-04～2020-06期間GNSS三維形變數(shù)據(jù)。通過歐洲航天局官網(wǎng)下載同期67景C波段升降軌Sentinel-1A斜距單視復(fù)數(shù)(SLC)影像，衛(wèi)星重訪周期為12 d，分辨率為5 m×20 m，為提高數(shù)據(jù)處理效率，裁剪影像至合適區(qū)域(圖2(b))，實驗參數(shù)詳見表1。

表1 Sentinel-1A實驗參數(shù)Tab.1 The experimental parameters of Sentinel-1A

根據(jù)前文可知，受SAR衛(wèi)星飛行軌道影響，InSAR在各方向的敏感程度不同。提取表1中衛(wèi)星升降軌入射角和方位角，可由式(1)表示為：

由式(17)可知，InSAR數(shù)據(jù)對垂直向形變最敏感，東西向次之，南北向不敏感。

2.2 實驗過程與分析

針對GNSS監(jiān)測數(shù)據(jù)，采用克里金(Kriging)插值法將點數(shù)據(jù)內(nèi)插至與InSAR相同像元的面，得到GNSS三維形變場(圖3)。采用短基線集(SBAS)方法處理InSAR數(shù)據(jù)，并引入AUX_POEORB精密定軌星歷和30 m分辨率的SRTM DEM數(shù)據(jù)去除地形相位，得到升降軌LOS向累積形變場(圖4，其中東、北和上為正方向)。

圖3 Kriging插值的GNSS三維形變場Fig.3 GNSS 3D deformation fields with Kriging interpolation

圖4 升降軌InSAR LOS向形變場Fig.4 InSAR LOS direction deformation fields of ascending and descending

利用GNSS三維形變與InSAR升降軌LOS向形變結(jié)果提取融合后的三維形變場，根據(jù)式(2)構(gòu)建誤差方程：

為驗證HVCE-RBFNN法的有效性，分別用3種方法進行解算：1)結(jié)合GNSS南北向形變與InSAR升降軌LOS向形變，解算東西向和垂直向形變場；2)結(jié)合GNSS三維形變與InSAR升降軌LOS向形變，采用等權(quán)估計法定權(quán)，利用最小二乘法解算融合后的東西向、南北向和垂直向形變場；3)采用HVCE-RBFNN法融合解算方法2中的5組形變數(shù)據(jù)的東西向、南北向和垂直向形變場。由方法1可以解算東西向和垂直向形變場，南北向形變場由GNSS插值代替，方法2和方法3可解算東西向、南北向和垂直向形變場。將139個監(jiān)測點的觀測量作為真值，分析3種方法三維形變的精度，結(jié)果如表2所示。

表2 不同方法的三維形變誤差Tab.2 The 3D deformation errors of different methods

由表2可以看出，方法1南北向RMSE為0 mm，但東西向RMSE達到50.22 mm，垂直向RMSE達到75.63 mm，監(jiān)測精度較低。對比方法1和方法2發(fā)現(xiàn)，方法2在綜合考慮了GNSS三維形變和InSAR升降軌LOS向形變的情況下，東西向和垂直向的形變精度顯著提升，南北向形變精度雖不及方法1，但顧及了InSAR監(jiān)測在南北向的貢獻，精度可達mm級。對比方法2和方法3發(fā)現(xiàn)，后者在三維方向的精度略優(yōu)于前者，盡管二者的精度差別不大，但在一般情況下，方法2利用等權(quán)法決定權(quán)重難以將各組數(shù)據(jù)進行合理融合。綜上，通過方法3迭代定權(quán)解算的三維形變結(jié)果精度優(yōu)于方法1和方法2。

2.3 地表三維形變結(jié)果及分析

圖5為139個監(jiān)測點GNSS三維形變結(jié)果與方法3得到的結(jié)果的對比。不難發(fā)現(xiàn)，由于InSAR對垂直向形變最為敏感，對南北向形變最不敏感，因此圖5(c)中融合形變較GNSS形變有一定差異，但曲線走勢基本一致；圖5(b)中2種形變曲線高度一致，驗證了InSAR對南北向形變貢獻小的問題；圖5(a)中曲線走勢介于圖5(b)和圖5(c)之間。整體來看，利用HVCE-RBFNN法得到的三維形變結(jié)果與GNSS三維形變結(jié)果較為一致。

圖5 GNSS與HVCE-RBFNN法的三維形變結(jié)果對比Fig.5 Comparison of 3D deformation results between GNSS and HVCE-RBFNN method

提取西二采礦區(qū)融合后的三維形變場，并與地下采空區(qū)位置進行疊加，以分析地下開采對地表的影響。由圖6(a)可知，采空區(qū)以西區(qū)域向東偏移，最大偏移量為228 mm；采空區(qū)以東區(qū)域向西偏移，最大偏移量為62 mm。由圖6(b)可知，采空區(qū)以南區(qū)域向北偏移，最大偏移量為300 mm；采空區(qū)以北區(qū)域向南偏移，最大偏移量為99 mm。由圖6(c)可知，礦區(qū)最大沉降量為193 mm，最大抬升量64 mm，沉降中心與采空區(qū)中心重疊，在地表形成沉降盆地。

圖6 基于HVCE-RBFNN法的三維形變場Fig.6 3D deformation fields based on HVCE-RBFNN method

采空區(qū)地表東西向和南北向的形變量在采空區(qū)中心表現(xiàn)平穩(wěn)，并由中心向兩側(cè)先增大后減小，垂直向形變量由采空區(qū)中心向四周逐漸減小。整體來看，研究區(qū)地表形變受地下開采和斷層影響，但其三維形變結(jié)果與采空區(qū)基本一致，符合開采沉陷規(guī)律。

3 結(jié) 語

本文以金昌市金川西二采礦區(qū)為研究對象，分別采用GNSS和SBAS-InSAR方法對礦區(qū)地表進行監(jiān)測，得到GNSS在三維方向和升降軌InSAR在LOS向的形變數(shù)據(jù)；然后根據(jù)GNSS和InSAR的優(yōu)勢互補性提出HVCE-RBFNN方法解算礦區(qū)地表三維形變，并通過3種不同的融合方法驗證其有效性。計算表明，本文方法解算的三維形變結(jié)果在東西向、南北向和垂直向的RMSE分別為20.85 mm、7.41 mm和34.47 mm，優(yōu)于其他2種方法。同時，利用本文方法獲取的三維形變場與采空區(qū)基本一致，符合開采沉陷的基本規(guī)律。由此可知，本文提出的方法可用于礦區(qū)地表的三維形變監(jiān)測。