基于雷達遙感衛星數據的礦區形變三維測量研究

馬 飛，張垚杰

（長治學院 計算機系，山西 長治 046011）

目前礦區沉降監測的常規方法有水準測量、全站儀測量、GPS 測量及三維激光掃描技術等測量方法。傳統的方法是在采空區上方設置監測點，對監測點進行連續不斷的觀測，得到測量值，經過平差測量后分析其形變趨勢。在我國黃土高原區的煤礦儲藏區地質結構復雜，地表塬梁縱橫，地表沉陷形態多樣化，若長期采用常規監測方法，觀測站的離散沉降值無法反映出礦區的真實沉降全貌。礦區的常規沉降監測方法僅限于單點的沉降監測。隨著技術的發展，沉降監測要求不僅要進行單點觀測，還提出了快速獲取地面大面積的沉降圖的要求。隨著雷達衛星技術的發展，InSAR技術被逐漸應用到礦區地表沉降監測的領域[1]。

InSAR 技術應用于沉降監測源于20 世紀90年代，國內外許多學者應用InSAR 技術并取得很多成果。該技術與常規測量技術相比有很大的優勢，其應用也愈來愈廣泛，實用價值愈來愈高。

劉國祥教授在香港赤鱲角國際機場進行干涉測量技術實驗，發現SAR 影像可以成功解譯到連續一年的沉降量，且與傳統監測方法得到的數據一致，精度優于1 cm[2]。

秦姍蘭等人為分析甘肅平涼華亭煤礦沉陷區的特點，將2007 至2010 年連續4 年的SAR 數據進行干涉組合，通過干涉測量技術獲取了差分干涉圖，對干涉圖中的數據進行分析，反演獲取了礦區地表下沉速率，為礦區地質災害評估和安全生產提供了數據依據[3]。

張景發教授等人利用多源SAR 數據，對主輔影像進行優化組合，通過干涉測量獲取了河北武安礦區的差分干涉圖，定量地分析了該礦不同時期的幅度信息和地面沉降分布[4]。

張學冬等人將連續兩年多的ENVISAT 影像進行二軌差分處理，發現在2004 年11 月至2006年3 月間荊各莊煤礦出現嚴重的下沉，且下沉速率逐年遞增[5]。

鄧喀中、祝傳廣等學者研究了城市地下水開采導致的地表下沉，以天津市塘沽地區為例，處理了連續時序的19 景ERS-1/2 影像，反演獲取了該地區1992～1998 年間城市地表沉降歷史[6]。

獨知行等人通過對InSAR 技術和GPS 基本概念，模型算法研究，發現這兩項技術的融合研究應用優勢很大，GPS 可以在精度評定，質量控制等方面對InSAR 定量分析作補充，提出一套GPS 和InSAR 數據融合的研究框架[7]。

阮崇籍等人深入分析了InSAR 高程測量模型，對基線估算、測高模型、仿真處理等模型進行分析和實驗，得到多種SAR 衛星平臺的基線測高模型和誤差模型[8]。

近幾年，廖明生、李志偉、吳立新、汪云甲等許多學者也進行大量的D-InSAR 研究與應用工作，根據我國煤礦區地質結構特性提出了不同的InSAR 技術在礦區開采沉陷監測中的應用思路。然而，專門針對我國黃土高原地質特性的InSAR技術應用研究成果還較少，尤其是西北地區黃土高原地下采空區沉陷監測的研究工作還處于探索階段，許多技術問題有待進一步研究。

1 合成孔徑雷達干涉測量基本原理

1.1 InSAR 技術的一般原理

InSAR 測量是傳統的微波遙感技術與光波干涉技術結合的產物。微波遙感技術分為機載或星載雷達兩種。該技術利用平臺搭載衛星的側視成像技術，和常規監測技術相比，具有空間分辨率高，全天候、全天時大面積地獲取地面精確三維信息的優勢。該技術通過向地球表面主動發射電磁波輻射，傳感器捕獲其后向散射強度信息和相位信息，再對覆蓋同一地區的不同時域的兩幅圖像配準，干涉處理后提取相位差圖，即相干圖。假設當地在這時間段內沒有發生形變，利用該技術可建立數字高程模型。正如我們所知，為了得到影像的目標點的高程，需通過采取立體攝影測量的基本方式，在星載SAR 衛星上只有一個傳感器。為了對同一觀測點進行兩次觀測，星載SAR采用重復觀測模式，雷達干涉測量的一般原理如圖1 所示[9]。

圖1 雷達干涉測量成像一般原理圖

上圖中S1 和S2 代表衛星兩次觀測位置，H是S1 到地面的高度，B 是量測觀測點間的距離（基線），為S1 的入射角，為S1 到目標點P 的距離，為S2 到目標點P 的距離，h 為我們需要計算的測量值，即目標點到地面的高程。

當S1、S2、P 之間的幾何關系穩定或滿足一定精度可以用來解算，則目標點的高程可按下式進行計算：

1.2 二軌法D-InSAR 技術原理

D-InSAR 技術是從InSAR 技術發展而來的，它考慮到了InSAR 技術在測量地表高程時目標區域產生形變信息時，通過已知的高程信息去除高程相位，解譯提取地面的微小形變相位信息，該技術理論觀測精度可達亞厘米級。D-InSAR 具有全天時，不受天氣的影響，覆蓋范圍廣，分辨率高及地表垂直靈敏度高等優點。目前常用的D-InSAR 技術有二軌法測量技術、三軌法測量技術及四軌法測量技術，這三種方法的幾何原理基本相同。隨著D-InSAR 技術在礦區沉降監測的應用日益增加，也出現了一些相應的D-InSAR改進技術[10]。

二軌法差分干涉測量技術是利用已知的外部DEM 消除雷達干涉圖中相位的地形相位。該方法先對兩景SAR 數據進行影像聚焦處理，配準后進行相位干涉處理，利用外部DEM 仿真地面高程干涉相位，在獲取的SAR 干涉圖中去除仿真的高程相位，最終獲取到只包含有地表信息的形變相位。在二軌法的處理影像的過程中，一般先去除地形相位，再去除平地相位。理論上相位圖中只包含形變相位，但由于噪聲相位的誤差和大氣相位的干擾，需進一步進行濾波等操作。具體的流程圖如圖2 所示。其優勢是數據處理流程中跳過了InSAR 相位解纏步驟，工作量少，操作步驟簡單，運行效率高，但可能存在將外部DEM 的誤差帶入干涉相位差中的問題[11]。

圖2 二軌法數據處理流程

數據配準過程可分為去噪、控制點自動搜索、重采樣與配準質量評價等。數據配準的算法較多，包括領域分析、基于幅度信息的配準、SIFT 算子配準等等。配準精度要求達到子像素級，每個像素點的配準誤差與均值差不得超過三倍標準差時。經過該處理，影像噪聲則可被抑制或消除[12-13]。

濾波是指對SAR 圖像或者干涉圖中的噪音信息濾除，一般會在干涉測量之前對原始SAR 影像進行濾波，去除雜質微波信號；干涉處理之后對干涉圖中的誤差相位進行濾波，去除干涉測量系統誤差。濾波處理算法包括頻率域濾波和空間域濾波，可以實現幅度信息和相位信息的濾波。

干涉圖生成需進行相位干涉，即對兩景影像的匹配點的數值做共軛相乘，得到地面任意點的干涉后影像圖的數值。通常計算得到的相位差由干涉圖的灰度表示，并以條紋形狀表現。

在平地效應過程中，地面某兩點高程不同，其相位差可能有兩部分造成的，其一為平地兩點的相位差；其二為兩點間存在的高差而造成相位差。通過去平地效應，使得到干涉圖的相位差只由地形引起，從而得到更為準確的地形數據。

針對相位解纏，由于干涉圖中的相位值取值范圍為（-π,π]，即從干涉圖得到的相位缺少了相位整數主值，為了得到真實相位差需在這個干涉圖相位值的基礎上加上2π 的整數倍。

由干涉圖得到高程值后，只是得到各點在SAR 坐標系下的高程值。為了得到對應坐標系下的數字高程模型，需將SAR 坐標系下的每一點高程值進行坐標轉換，得到最終的地面形變值。上述過程為地理編碼。

2 實驗研究與結果分析

2.1 數據來源及研究區概況

SAR 影像數據為日本ALOS 衛星的兩景PALSAR 影像數據，其衛星搭載L 波為工作波，波長為23.5cm。該影像覆蓋我國陜西省咸陽市彬縣亭南礦區。該兩景數據獲取時間分別為2009 年8 月24 日和2009 年10 月09 日。根據該兩景影像數據，選取2009 年8 月24 日獲取的影像為主影像；選取2009 年10 月09 日獲取的影像為輔影像，具體參數見表1.

表1 SAR 影像數據參數

利用二軌法D-InSAR 技術需要研究區域的外部DEM，其精度對干涉影像差分的精度有較大的影響。文章選取的DEM 為SRTM 獲取的研究區地面雷達數據。SRTM（全稱為Shuttle Radar Topography Misson）由美國航空航天局、地理空間情報局等航天機構組成，該機構于2002 年開始，獲取了全球絕大部分的地面雷達數據。其數據的平面參考系為WGS84 坐標系，其垂直參考系為WGS84 EGM96 大地水準面[14]。

當利用經緯度對SRTM 數據劃分時，其精度可分為90 m 和30 m 兩類。根據已有條件，選用精度為90 m 的DEM 數據。該數據每隔90 m 取一個高程點，該高程點的高程由9 個精度為30 m的數值平均計算得到，滿足精度要求[15]。

本次研究區域的緯度范圍為34 40 至35 50 ，經度范圍為107 30 至108 40 ，需先對DEM 進行拼接，形成新的覆蓋研究區域的DEM。本次研究沉降監測的礦區為亭南井田。表2 為亭南井田的煤炭開發數據表。

表2 亭南井田煤炭開采數據表

2.2 SAR 數據干涉測量處理

文章采用重復軌道干涉測量軟件包（repeat orbit interferometry package, ROI_PAC）完成兩景影像的處理工作。ROI_PAC 軟件基于C 語言開發，主要運行在Unix 操作系統和Linux 操作系統，是較為成熟的SAR 影像處理軟件。在初步處理中，利用該軟件的處理器功能完成方位向和距離向的壓縮、對焦等步驟。具體的處理流程為如圖3。

至于幾種證明方法，課前都有準備，但不都是自己想出來的，很多是以前的學生想出來的．課前設想只有思路3和思路4必講，尤其是思路4，用向量的方法解決問題是一種意識，教材中雖有所涉及，但不夠系統，需教師自己去總結．至于其它方法如果學生不提出來，可能略講或不講，一切取決于學生的需求．解題方法的呈現要有適當的時機，要有充分的理由，也就是說要講理．

圖3 SAR 初步處理流程圖

在初步處理中，由衛星的數據信號得到具有復數結構的影像圖，即得到兩景數據的單視復影像和參數文件，為后面的差分干涉提供基本的影像數據。圖3 為兩景影像初步處理后的單視復影像。

圖3 兩景影像的單視復影像圖

主輔影像的差分干涉測量處理是基于軟件干涉處理模塊，最終得到的結果為經坐標轉換和地理編碼的差分干涉圖。具體的處理流程如圖4 所示。

圖4 兩景影像的差分干涉處理流程圖

圖5 地理編碼后的差分干涉相位圖

得到經坐標轉換和地理編碼的差分干涉圖，該干涉圖的坐標系統為WGS84 坐標系，而礦區已得水準數據的坐標系統為1954 北京坐標系統。因此需要先利用Sufer軟件得到其柵格矢量文件，其儲存格式為點號、經度、緯度、形變量，需對SAR 影像處理結果進行坐標轉換處理，具體的流程如圖6 所示。

圖6 基于MATLAB 與Sufer 軟件的后期處理流程圖

根據前面介紹的影像裁剪方法，按經緯度對原影像圖進行裁剪，得到亭南井田沉降圖。為了更形象地表現亭南礦區的沉降效果，利用Sufer8.0 軟件得到該礦區將WGS84 坐標轉換北京54 坐標的點文件。該過程的實現依賴于坐標系統轉換。本次坐標轉換的方法是七參數轉換，即利用研究區域的已知北京54 坐標和WGS84 坐標兩套坐標系統的4 個控制點，先得到各自的空間坐標系的坐標，再用對應點的空間坐標系下的坐標得到7 個轉換參數，最后利用該轉換參數，將亭南礦區的WGS84 坐標轉換到北京54 坐標系下，得到54 坐標點文件，實現了坐標轉換的目的。其中7 參數指坐標轉換模型中，采用三個3 個平移參數，3 個旋轉參數及1 個尺度縮放因子的布爾莎模型[16]。坐標轉換的流程如圖7 所示。

圖7 七參數坐標轉換的流程圖

因亭南礦區范圍內的沉降量的變化范圍較大，裁剪得到礦區的沉降圖，如圖8 所示。

圖8 裁剪后的亭南礦區沉降圖

利用得到的北京54 坐標點文件模擬亭南礦區的沉降區域圖，如圖9 所示。

圖9 亭南礦區沉降區三維視圖

上圖中，由矩形經緯度區域轉換成北京54坐標系呈梯形區域，沉降區域的形狀幾乎為漏斗，越是靠近沉降中心的點，其沉降量越大，代表其對應的地表被破壞越嚴重。由沉降圖可得到該礦區在兩景影像獲取時間段內地表的沉降量。再結合該時間段內的地下采煤，工作面的推進情況，可作為預測該礦區未來某時間段內由地下開采引起的地表沉降災害的依據。

3 結束語

（1）基于雷達遙感數據，采用二軌法干涉測量數據處理可以快速研判煤礦區地表沉降區域位置，實現遙感定量獲取地表沉降數據。

（2）利用二軌法干涉測量地表形變數據進行三維可視化分析，可以繪制形變區等值線圖，得到其在研究時間段內的呈現漏斗狀三維圖，其最大沉降量為4.89 cm。

（3）在亭南礦區沉降監測的后續研究中，可以從以下幾個方面進行改進：獲取現今服役的新型高分辨率衛星的SAR 影像，獲取的充足的影像數據，獲得的成果不僅有沉降圖，還可獲取沉降速度等因素用于影像處理的評價因子；可根據近幾年礦區的沉降范圍及沉降量，礦區的開采量，事先布設角反射（CR）點，利用CR 觀測模型來觀測礦區地面的較大的沉降；在獲取影像的時間段內，測定研究區域一定數量點的GPS 數據，用其來矯正SAR 影像的相對大氣延遲誤差，采用常規技術與D-InSAR 融合技術進行影像處理。