侯建國,白 駒,陶紅志,劉 江,王 強
(1.黑龍江工程學院 測繪工程學院,哈爾濱 150050;2.國家基礎地理信息中心,北京 100830;3.哈爾濱工大正元信息技術有限公司,哈爾濱 150016)
目前,InSAR領域研究應用更多的是基于多幅雷達圖像的二次差分處理來監測地球表面形變,以揭示許多地球物理現象,如地震形變、火山運動、冰川漂移、地面沉降以及山體滑坡等,這種技術被稱為雷達差分干涉測量(Differential Interferometric Systhetic Aperture Radar,D-InSAR)。常規D-InSAR由于時空失相關、大氣延遲和相位噪聲大等特點,精度一般在cm量級。時序InSAR(Interferometric Synthetic Aperture Radar)是基于長時間序列的SAR影像,利用分辨單元內散射體的后向散射特性進行建模,具有長時間尺度、高空間密度和mm級精度,典型TS-InSAR(Time Series InSAR)技術有SBAS-InSAR(Small BaselineSubset InSAR)、PS-InSAR(Permanent Scatterer InSAR)和一些演變技術。InSAR技術的另一個應用領域是重建數字高程模型(DEM),DEM作為地形模擬的重要成果,已經成為國家空間數據基礎設施的基本內容之一,它可以通過全站儀、GNSS-RTK等技術野外直接測量,或者間接從航空影像、遙感影像、既有地形圖上得到,通過InSAR技術獲取DEM,具有大范圍、快速、高效等特點,如該技術能夠獲取覆蓋全球陸地的80%中等分辨率SRTM-DEM(Shuttle Radar Topography Mission Digital Elevation Model)數據、德國的TanDEM-X衛星計劃能夠獲取全球高分辨率DEM數據。利用InSAR技術生成DEM,并對區域的兩個DEM作差進行比較,可以發現地形變化量,進而計算工程上需要填埋體積大小。基于此本文為有效評估城市垃圾場垃圾存量,通過調查雷達衛星存檔數據,利用InSAR技術得到不同時期的DEM數據,進而計算得到垃圾場歷史垃圾填埋量。
InSAR技術是利用雷達回波信號的相位值來提取地面目標高程的。其基本原理是利用單軌雙天線同時接收雷達回波或者重軌單天線間隔一定時間兩次接收雷達回波,獲取統一區域的兩個雷達回波信號,然后根據這兩個雷達回波信號之間產生的相位差和一些軌道數據來計算地面目標的高程信息。SAR圖像的像元分辨率通常在m級以上,而SAR系統采用的電磁波波長一般在1 m以內,多數在cm級,因此,采用干涉相位測量(斜)距離差精度較高,利用這些高精度的干涉相位信息可以獲取高精度的高程和DEM信息。
星載SAR數據一般是以重復軌道模式獲取的,雷達影像復數據記錄的相位信息可以換算成目標點到傳感器的距離,即以一定分辨率測量方位向和距離向目標的距離。單次飛行獲得單個SAR影像簡化成垂直于方位向的一個剖面,如圖1(a)所示,S為傳感器位置,但是僅僅知道距離并不能確定目標的位置和相對于某水準面的高程,圖中凡是在波速范圍內且位于同一弧線上的目標所測得距離都是相等的。如果考慮從位于不同位置的傳感器獲得同一場景的另一幅影像,就可以解決單幅影像不確定性問題。圖1(b)中,利用附加的信息就可以確定出某一目標的位置和高程z。將兩幅影像分別稱為主影像和從影像,如圖1(b)所示。

圖1 單幅SAR和兩幅SAR成像幾何關系
傳感器S1和S2之間距離就是基線B。主影像上的一點到目標的距離是r,從影像上可以找到對應同名點,它到目標的距離是(r+dr)。H是航高;z是地物目標的高程。可以按式(1)計算。
(r+dr)2=r2+B2-2rBsin (θ-α),
(1)
z=H-rcosθ.
(2)
假設式中r,dr,B,α已知,α是基線與水平方向夾角,θ是視線與垂直方向夾角,由式(1)求出θ,由式(2)求出高程z。考慮r?B,r?dr,可以認為S1和S2的視線是平行的,從式(1)可得:
dr=-Bsin(θ-α).
(3)
根據成像幾何關系和提取高程精度之間聯系,假設準確知道傳感器參數和軌道參數B,H,θ,α等,r即可確定,dr精度主要取決于斜距方向的分辨率,斜距方向的分辨率越高,dr測量精度越高。從式(2)可知,最終高程測量精度與dr誤差之間的關系:
(4)
寫成中誤差形式為:
(5)
式中:mz為高程測量誤差,mr為斜距測量誤差。從式(5)可以看出,高程測量精度取決于斜距和基線比r/B。在SAR成像系統中,這個比值很大,斜距r的微小誤差會經過放大后傳遞給高程。想要改善高程精度就需要基線拉長,但是如果基線拉長后,兩次成像視角相差太大,以至后向散射波發生很大變化,加上斑點噪聲影響,主從影像之間的相似性蕩然無存,也就無法找到同名點。因此,增加基線不能解決根本問題。雷達干涉測量中斜距之差dr是由影像之間的相位差得到。
(6)
相位測量誤差mφ算式為:
(7)
設C波段波長 5.6 cm,θ=23°,相位φ的噪聲標準差為30°,則距離測量誤差為2.3 mm。若基線110 m,其高程精度5 m。
InSAR通過重復軌道模式,獲取地面同一區域的復數影像對,由于目標與天線位置的幾何關系,在復圖像上產生相位差,形成干涉條紋圖。干涉條紋圖中包含了斜距向上的點與兩天線位置之差的精確信息。因此,利用傳感器高度、雷達波長、波束視向及天線基線距之間的幾何關系,可以精確測量出圖像中每一點的三維位置,生成該地區的DEM。InSAR影像提取DEM的算法主要包括:復圖象的配準、干涉圖的生成、相位解纏、數字高程模型重建。圖2是雷達信號到DEM產品的處理流程。

圖2 DEM重建流程
哈爾濱市區某垃圾場,1996年建廠,2005年關閉,為了評估現存垃圾對周圍環境的影響,計算出垃圾堆放量,通過調查雷達衛星存檔數據,根據InSAR技術進行DEM獲取的原理及數據處理流程,采用目前國際上最先進的InSAR數據處理軟件—瑞士GAMMA干涉雷達圖像處理軟件,開展了數字高程模型生成處理,得到不同時期研究區的DEM,對DEM做差求出變化值,最后計算得到垃圾歷史存量數據。本項目訂購了8景歐空局的ERS衛星數據和11景Envisat衛星數據,如表1所示。通過常規雷達數據干涉處理,共獲得11個干涉數據對,如表2所示。

表1 研究區SAR數據列表

表2 雷達數據干涉處理生成的干涉對
選擇3個干涉質量較高的干涉對,計算出1996年1月、1997年10月、2004年1月的DEM數據,用于本次地形變化檢測。1996年1月和1997年10月的DEM,由于SAR數據相隔1 d,數據保持了較好的相干性,能夠恢復出完整的DEM信息,而2004年的兩個影像相隔35 d,影像受到失相干因素的影響,導致部分地區的地面高程信息不能完全恢復,變化檢測計算時對其進行插值處理。圖3是InSAR生成的 1996年1月垃圾場位置的DEM。

圖3 InSAR 生成的DEM
3.1.1 地形變化最大區域干涉處理結果
不同時期DEM在坐標統一的基礎上,兩個DEM作差進行比較,圖4中(a)表示1997年10月與1996年1月的地形變化,圖4中(b)表示2004年1月與1997年10月的地形變化,從結果圖中可以明顯看出,在1997年10月與1996年1月期間,研究區內有明顯地形變化區,由于2004年的干涉結果受到失相干因素的影響存在數據缺失和空值現象,為了對2004年1月與1997年10月期間DEM進行比較,試驗中對數據缺失部分的地形數據采用克里金插值方法進行了填充。在圖4中還存在大范圍的趨勢性變化,可能是由SAR信號傳播過程中的大氣延遲誤差引起的誤差,可以忽略。

圖4 不同時間段的地形變化檢測結果(°)
3.1.2 地形變化統計分析
對圖4中1996年1月與1997年10月的地形重點變化區域進行統計,結果如圖5所示。地形抬升變化主要分布在兩個高度段,4 m以內和6~14 m,沿著圖4(a)中a和b提取了地形變化曲線,如圖6和圖7所示。

圖5 1997年與1996年地形變化統計結果

圖6 沿圖4a中(a)線地形變化剖線

圖7 沿圖4(a)中b線地形變化剖線
按照實際邊界范圍裁剪1996年1月的DEM與1997年10月的DEM,求差得到區域變化值展繪到光學影像上,如圖8所示,同樣裁剪1997年10月的DEM與2004年1月的DEM,變化后地形如圖9所示。根據體積變化的像元進行統計得到,地形變化區的面積約為18 800 m2;根據地形變化區域的面積在不同像元上高程的變化值,計算得到1996—2004年該地區地形變化的體積約為146.25×104m3。

圖8 1996年1月至1997年10月地形變化圖

圖9 1997年10月至2004年1月地形變化圖
根據3個不同時間段的DEM數據和垃圾場填埋垃圾范圍,計算得到該地區累計地形變化量,進而得到垃圾場建成使用至關閉期間的垃圾填埋量為146.25×104m3。垃圾場占地面積約18 800 m2,SAR處理時,影像的相干系數為0.36,根據InSAR測量高程精度分析結果,本次DEM測量精度小于5 m,計算體積誤差為9.4×104m3,存量垃圾計算誤差約為6.4%。通過對ERS和Envisat影像對的干涉處理,得到垃圾填埋區的DEM,對其進行定性與定量分析表明,該技術檢測的地形變化量,計算得到垃圾場存量數據是可靠的。本項目試驗研究中也發現幾點影響高程精度關鍵步驟,雷達像對配準的精度,會直接影響相位差的計算,最終導致高程誤差;解纏算法本身所造成的高程誤差;利用相位計算高程值的過程中,由基線模型所造成的誤差等。