楊 雋 汪 娟 王 瑞 郭金城 馮 瀚 向淇文

（1.貴州省第一測繪院；2.貴州省北斗導(dǎo)航位置服務(wù)中心；3.贛南科技學(xué)院）

InSAR 是一種新興的雷達(dá)遙感技術(shù)，近20 a 來得到了迅速發(fā)展，能夠高精度地獲取地表形變信息［1］。相較于傳統(tǒng)方法，該方法具有空間分辨率高、覆蓋范圍廣、全天候監(jiān)測和不受天氣影響等優(yōu)勢，可以有效地揭示區(qū)域性地表宏觀趨勢［2］。許多國內(nèi)外學(xué)者已經(jīng)對InSAR在礦區(qū)監(jiān)測方面進(jìn)行了大量研究，比如楊澤發(fā)等［3］利用水準(zhǔn)測量和Levenberg-Marquard 算法模擬礦區(qū)動態(tài)沉降過程，陳磊［4］提出了基于Konthe模型的礦區(qū)大尺度形變信息的預(yù)計(jì)方法。

然而，在礦山開采過程中地表沉陷是三維的，傳統(tǒng)的InSAR只對視線向的形變敏感，對于東西向和南北向變形監(jiān)測精度不高。為獲取地表真實(shí)變形，本文提出了一種單視線向解算方法，結(jié)合開采沉陷梯度關(guān)系來提高對東西向和南北向變形監(jiān)測的精度。該方法利用開采引起的水平移動與傾斜之間的數(shù)學(xué)關(guān)系，對反演三維形變場提供約束條件，并構(gòu)建解算模型，通過該模型可以對開采礦區(qū)地表進(jìn)行三維形變提取，不僅能夠探索礦山開采過程中地表形變的規(guī)律，而且為礦山的安全生產(chǎn)和地表災(zāi)害的預(yù)測預(yù)防提供了可靠的技術(shù)支撐。

1 InSAR視線向模糊問題分析

合成孔徑雷達(dá)自身側(cè)視成像的工作模式存在著視線向模糊的問題，即地表東西（E-W）、南北（S-N）、豎直（U）三維方向上的真實(shí)的地表形變分量為視線向（LOS）上的投影，這使得觀測到的視線向形變并不是地表的真實(shí)形變［5］。將地表點(diǎn)的變形看成在垂直、東西、南北（U，E，N）3個方向分量的組合，并且3個方向的分量對于LOS 形變的貢獻(xiàn)各不相同，圖1 為In-SAR視線向形變量與地表三維形變場之間的關(guān)系，αh表示衛(wèi)星飛行方向與正北方向的夾角，θ表示雷達(dá)的入射方位角。

假設(shè)地面目標(biāo)視線向形變遠(yuǎn)離雷達(dá)為正，靠近雷達(dá)為負(fù)時，則dLOS與東西dE、南北dN和垂直dU向形變量之間的關(guān)系可以有：

由式（1）可知，利用一維形變求解三維形變會使得計(jì)算矩陣產(chǎn)生秩虧，因此，需要有不同視角的雷達(dá)觀測量進(jìn)行聯(lián)合解算，但是在實(shí)際應(yīng)用中存在許多制約和不足，如需要3個不同SAR傳感器的觀測值比較難以滿足；此外該方法還假定觀測目標(biāo)在多平臺觀測的時間間隔內(nèi)是不變或者線性變化的，這對于礦區(qū)的非線性形變來說，受到了一定的限制。

地下開采引起的巖層及地表移動過程是一個復(fù)雜的時間-空間過程，在移動盆地內(nèi)地表各點(diǎn)的移動方向和移動大小都不相同，巖層的移動會逐漸傳遞直至地表，使得地表發(fā)生下沉和水平移動，形成一個三維形變場。這導(dǎo)致了InSAR 在礦區(qū)地表形變監(jiān)測中無法準(zhǔn)確地獲得地面形變特征和規(guī)律。

2 基于開采沉陷知識的三維形變提取方法

2.1 開采引起的水平移動與傾斜的關(guān)系

本方法將結(jié)合采礦區(qū)開采沉陷規(guī)律，利用開采引起的水平位移與垂直位移梯度之間的函數(shù)關(guān)系作為約束條件，更好地將差分干涉圖反演為地表真實(shí)的三維形變場。

隨機(jī)介質(zhì)模型是開采沉陷工程中預(yù)測采動位移最常用的方法之一［6］，在這個模型中，一個地下采掘面可以分成許多無窮小的采掘單元，巖體內(nèi)部（x，z）點(diǎn)受單元開采影響產(chǎn)生的水平移動量Ue（x，z）和下沉量We（x，z）如圖2所示。

根據(jù)Knothe，地表單元下沉函數(shù)服從正態(tài)分布：

式中，z為開采深度，m；x為開采邊界投影到地表的距離，m；rz表示在z處開采影響的球體半徑，m。

根據(jù)隨機(jī)介質(zhì)模型可知，無窮小的采掘單元引起的地表位移等于其位移之和，假設(shè)工作面板上的巖體不可壓縮，則在二維平面上x、z方向上的法向應(yīng)變之和為零，通過數(shù)學(xué)推導(dǎo)可得到

式中，b為水平位移常數(shù)；β水平位移主要影響角，（°）；T(x,H)為沿x方向的傾斜。

式（3）表明，在x-z平面上，地表點(diǎn)的水平位移Ue(x,H)與垂直下沉的水平梯度成比例，同理，該方法在y-z平面也同樣適用。

2.2 SAR三維形變場解算模型

對于礦區(qū)而言，開采所影響的地表點(diǎn)在水平位移與垂直下沉之間存在一定的函數(shù)關(guān)系，可以使用單對干涉差分圖反演地表三維變形，但是由于合成孔徑雷達(dá)失相干導(dǎo)致視線向形變場的缺失，將導(dǎo)致解算模型失去連續(xù)的梯度從而給結(jié)果引入誤差，因此需要通過克里金插值來補(bǔ)充形變場的空值區(qū)域［7］。

假設(shè)地理編碼后的視線向形變場維數(shù)為n×m，則點(diǎn)（i，j）（i=1，…n；j=1，…m）的東西TE（i，j）和南北向TN（i，j）的傾斜可以表示為

式中，W（i，j）為地表點(diǎn)的下沉；ΔE、ΔN分別為差分干涉圖在東西向和南北向的空間分辨率。

由于當(dāng)前SAR 衛(wèi)星大多都是側(cè)視成像，所以SAR影像的像元陣列如圖3所示。

聯(lián)立式（3）、式（4），可以進(jìn)一步推導(dǎo)EW、SN向的水平位移為

式中，H（i，j）為開采深度。

在解算地表三維變形的過程中，首先利用dLOS求解W（n，m），再將其代入到前面的方程直到求解出W（1，1），然后可以求解出東西向和南北向的水平位移。單軌干涉合成孔徑雷達(dá)反演礦區(qū)地表三維位移場的處理流程如圖4所示。

3 應(yīng)用實(shí)例

3.1 研究區(qū)概況

萬年礦位于峰峰礦區(qū)北部，該礦區(qū)橫跨武安市伯延、磁山、午汲3個鄉(xiāng)鎮(zhèn)，北距武安市10 km，東距邯鄲35 km。礦井的設(shè)計(jì)生產(chǎn)能力為180萬t/a。在井田范圍內(nèi)，井田南北長約9 km，東西寬1～4 km，面積約21.167 km2，礦井的開采標(biāo)高為325～-750 m。研究區(qū)在萬年礦范圍內(nèi)，屬于伯延鎮(zhèn)轄區(qū)內(nèi)，如圖5所示。

3.2 試驗(yàn)處理

本次試驗(yàn)使用的是Sentinel-1A 衛(wèi)星于2020 年2—4 月間獲取的峰峰地區(qū)的兩景數(shù)據(jù)，SAR 數(shù)據(jù)的基本信息見表1，由于VV 極化方式對于地表的穿透更強(qiáng)，因此試驗(yàn)采用了該極化方式成像的數(shù)據(jù)。

試驗(yàn)中選擇2019 年11 月的數(shù)據(jù)作為主影像，2020 年2 月的作為從影像，2 次影像獲取的時間基線為48 d，經(jīng)過主輔影像配準(zhǔn)、多視、干涉圖生成、去平地效應(yīng)、噪聲濾波、相位解纏、地理編碼［8-10］等處理步驟，最終得到研究區(qū)地表視線向形變，由于受地表不同地物相干性的影響，導(dǎo)致相干性低的區(qū)域出現(xiàn)形變場的缺失，這會使得圖像像素失去連續(xù)性，從而影響模型的解算精度。因此，在反演研究區(qū)的三維形變前需要對LOS 向形變?nèi)笔^(qū)域進(jìn)行插值，本文選用克里金插值法，利用方差圖作為權(quán)重函數(shù)，基于最小二乘法對數(shù)據(jù)點(diǎn)進(jìn)行空間建模和插值，如圖6 所示。

3.3 三維變形反演

研究區(qū)位于萬年礦東邊，靠近井田邊界，含煤地層主要為太原組、二疊系的山西組和石炭系的本溪組，根據(jù)萬年礦以往的開采實(shí)測資料，可以獲得該地質(zhì)條件下的開采沉陷參數(shù)見表2。

通過使用基于開采沉陷知識的三維形變提取算法，將單視線向LOS 形變代入模型計(jì)算，并對其進(jìn)行迭代求解得到三維形變場（圖7）。

圖7（a）為南北向水平移動，向北方向移動為負(fù)值，最大移動量為31 mm，向南方方向移動為正值，最大移動量為42 mm；圖7（b）為東西向水平移動，向東方向?yàn)樨?fù)值，最大移動量為35 mm，向西方向?yàn)檎?，最大移動量?7 mm；圖7（c）為垂直向下沉，最大下沉為61 mm?？梢钥闯?，利用該模型提取的三維形變場基本符合礦區(qū)的開采沉陷規(guī)律［11］，在沉陷盆地中，東西向和南北向的變形是向盆地中心進(jìn)行移動，而垂直向下沉為一個明顯的沉降漏斗。

在垂直向下沉中，出現(xiàn)了較為清晰的沉降阻隔線，在東西向水平移動變形圖中，該區(qū)域在遠(yuǎn)離盆地的情況下出現(xiàn)大幅向西方向的移動，將在監(jiān)測期間獲取的三維形變場進(jìn)行矢量合成（圖8）。結(jié)合地質(zhì)資料合理推斷是由于工作面在采動過程中巖層的靜態(tài)結(jié)構(gòu)發(fā)生變化，影響到地質(zhì)弱面斷層的穩(wěn)定性所導(dǎo)致的，斷層作為巖層中的弱面，力學(xué)強(qiáng)度遠(yuǎn)遠(yuǎn)低于周圍的巖體［12-13］，當(dāng)開采的影響作用于斷層時，其活化會使得巖體內(nèi)部的形變傳遞變得異常。

4 總結(jié)

本文采用“理論分析—模型構(gòu)建—實(shí)例分析”的主體研究思路，針對傳統(tǒng)單軌InSAR只對視線向形變敏感，而對于東西向和南北向變形監(jiān)測精度不高的不足，分析了InSAR 技術(shù)在礦區(qū)監(jiān)測的局限，并結(jié)合礦區(qū)開采沉陷中巖層內(nèi)部的梯度關(guān)系構(gòu)建地表的三維解算模型。得到以下結(jié)論：

（1）通過SAR 側(cè)視成像的幾何關(guān)系，分析其在礦區(qū)地表形變監(jiān)測的局限性，導(dǎo)致InSAR無法獲得準(zhǔn)確的礦區(qū)地表三維形變場。

（2）提出一種聯(lián)合開采沉陷梯度關(guān)系的單視線向解算方法，通過開采引起的水平移動與傾斜的數(shù)學(xué)關(guān)系為反演三維形變場提供很好的約束條件，構(gòu)建地表形變的三維解算模型。

（3）基于開采沉陷知識的三維形變提取方法，完成萬年礦區(qū)的地表形變監(jiān)測，提取了沉陷盆地的三維形變場，有效地揭示了該區(qū)域真實(shí)的地表形變特征，分析了地質(zhì)弱面斷層受采動影響活化使得巖層內(nèi)部形變傳遞出現(xiàn)異常。