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

（1.貴州省第一測繪院；2.貴州省北斗導航位置服務中心；3.贛南科技學院）

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

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

1 InSAR視線向模糊問題分析

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

假設地面目標視線向形變遠離雷達為正，靠近雷達為負時，則dLOS與東西dE、南北dN和垂直dU向形變量之間的關系可以有：

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

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

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

2.1 開采引起的水平移動與傾斜的關系

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

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

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

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

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

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

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

2.2 SAR三維形變場解算模型

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

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

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

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

聯立式（3）、式（4），可以進一步推導EW、SN向的水平位移為

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

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

3 應用實例

3.1 研究區(qū)概況

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

3.2 試驗處理

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

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

3.3 三維變形反演

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

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

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

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

4 總結

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

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

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

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