基于InSAR技術(shù)的城市垃圾場存量計算及應(yīng)用

侯建國，白 駒,陶紅志，劉 江，王 強

(1.黑龍江工程學(xué)院 測繪工程學(xué)院，哈爾濱 150050;2.國家基礎(chǔ)地理信息中心，北京 100830；3.哈爾濱工大正元信息技術(shù)有限公司，哈爾濱 150016)

目前，InSAR領(lǐng)域研究應(yīng)用更多的是基于多幅雷達圖像的二次差分處理來監(jiān)測地球表面形變，以揭示許多地球物理現(xiàn)象，如地震形變、火山運動、冰川漂移、地面沉降以及山體滑坡等，這種技術(shù)被稱為雷達差分干涉測量(Differential Interferometric Systhetic Aperture Radar,D-InSAR)。常規(guī)D-InSAR由于時空失相關(guān)、大氣延遲和相位噪聲大等特點，精度一般在cm量級。時序InSAR(Interferometric Synthetic Aperture Radar)是基于長時間序列的SAR影像，利用分辨單元內(nèi)散射體的后向散射特性進行建模，具有長時間尺度、高空間密度和mm級精度,典型TS-InSAR(Time Series InSAR)技術(shù)有SBAS-InSAR(Small BaselineSubset InSAR)、PS-InSAR(Permanent Scatterer InSAR)和一些演變技術(shù)。InSAR技術(shù)的另一個應(yīng)用領(lǐng)域是重建數(shù)字高程模型(DEM)，DEM作為地形模擬的重要成果，已經(jīng)成為國家空間數(shù)據(jù)基礎(chǔ)設(shè)施的基本內(nèi)容之一，它可以通過全站儀、GNSS-RTK等技術(shù)野外直接測量，或者間接從航空影像、遙感影像、既有地形圖上得到，通過InSAR技術(shù)獲取DEM，具有大范圍、快速、高效等特點，如該技術(shù)能夠獲取覆蓋全球陸地的80%中等分辨率SRTM-DEM(Shuttle Radar Topography Mission Digital Elevation Model)數(shù)據(jù)、德國的TanDEM-X衛(wèi)星計劃能夠獲取全球高分辨率DEM數(shù)據(jù)。利用InSAR技術(shù)生成DEM，并對區(qū)域的兩個DEM作差進行比較，可以發(fā)現(xiàn)地形變化量，進而計算工程上需要填埋體積大小?；诖吮疚臑橛行гu估城市垃圾場垃圾存量，通過調(diào)查雷達衛(wèi)星存檔數(shù)據(jù)，利用InSAR技術(shù)得到不同時期的DEM數(shù)據(jù)，進而計算得到垃圾場歷史垃圾填埋量。

1 InSAR高程測量精度及生成DEM基本流程

InSAR技術(shù)是利用雷達回波信號的相位值來提取地面目標(biāo)高程的。其基本原理是利用單軌雙天線同時接收雷達回波或者重軌單天線間隔一定時間兩次接收雷達回波，獲取統(tǒng)一區(qū)域的兩個雷達回波信號，然后根據(jù)這兩個雷達回波信號之間產(chǎn)生的相位差和一些軌道數(shù)據(jù)來計算地面目標(biāo)的高程信息。SAR圖像的像元分辨率通常在m級以上，而SAR系統(tǒng)采用的電磁波波長一般在1 m以內(nèi)，多數(shù)在cm級，因此,采用干涉相位測量(斜)距離差精度較高，利用這些高精度的干涉相位信息可以獲取高精度的高程和DEM信息。

1.1 高程測量精度分析

星載SAR數(shù)據(jù)一般是以重復(fù)軌道模式獲取的，雷達影像復(fù)數(shù)據(jù)記錄的相位信息可以換算成目標(biāo)點到傳感器的距離，即以一定分辨率測量方位向和距離向目標(biāo)的距離。單次飛行獲得單個SAR影像簡化成垂直于方位向的一個剖面，如圖1(a)所示，S為傳感器位置，但是僅僅知道距離并不能確定目標(biāo)的位置和相對于某水準(zhǔn)面的高程，圖中凡是在波速范圍內(nèi)且位于同一弧線上的目標(biāo)所測得距離都是相等的。如果考慮從位于不同位置的傳感器獲得同一場景的另一幅影像，就可以解決單幅影像不確定性問題。圖1(b)中，利用附加的信息就可以確定出某一目標(biāo)的位置和高程z。將兩幅影像分別稱為主影像和從影像，如圖1(b)所示。

圖1 單幅SAR和兩幅SAR成像幾何關(guān)系

傳感器S1和S2之間距離就是基線B。主影像上的一點到目標(biāo)的距離是r，從影像上可以找到對應(yīng)同名點，它到目標(biāo)的距離是(r+dr)。H是航高；z是地物目標(biāo)的高程?？梢园词?1)計算。

(r+dr)2=r2+B2-2rBsin (θ-α),

(1)

z=H-rcosθ.

(2)

假設(shè)式中r,dr,B,α已知，α是基線與水平方向夾角，θ是視線與垂直方向夾角，由式(1)求出θ，由式(2)求出高程z?？紤]r?B,r?dr，可以認(rèn)為S1和S2的視線是平行的，從式(1)可得：

dr=-Bsin(θ-α).

(3)

根據(jù)成像幾何關(guān)系和提取高程精度之間聯(lián)系，假設(shè)準(zhǔn)確知道傳感器參數(shù)和軌道參數(shù)B,H,θ,α等，r即可確定，dr精度主要取決于斜距方向的分辨率，斜距方向的分辨率越高，dr測量精度越高。從式(2)可知，最終高程測量精度與dr誤差之間的關(guān)系:

(4)

寫成中誤差形式為：

(5)

式中:mz為高程測量誤差，mr為斜距測量誤差。從式(5)可以看出，高程測量精度取決于斜距和基線比r/B。在SAR成像系統(tǒng)中，這個比值很大，斜距r的微小誤差會經(jīng)過放大后傳遞給高程。想要改善高程精度就需要基線拉長，但是如果基線拉長后，兩次成像視角相差太大，以至后向散射波發(fā)生很大變化，加上斑點噪聲影響，主從影像之間的相似性蕩然無存，也就無法找到同名點。因此，增加基線不能解決根本問題。雷達干涉測量中斜距之差dr是由影像之間的相位差得到。

(6)

相位測量誤差mφ算式為：

(7)

設(shè)C波段波長 5.6 cm，θ=23°，相位φ的噪聲標(biāo)準(zhǔn)差為30°，則距離測量誤差為2.3 mm。若基線110 m，其高程精度5 m。

1.2 SAR干涉數(shù)據(jù)生成DEM處理流程

InSAR通過重復(fù)軌道模式，獲取地面同一區(qū)域的復(fù)數(shù)影像對，由于目標(biāo)與天線位置的幾何關(guān)系，在復(fù)圖像上產(chǎn)生相位差，形成干涉條紋圖。干涉條紋圖中包含了斜距向上的點與兩天線位置之差的精確信息。因此，利用傳感器高度、雷達波長、波束視向及天線基線距之間的幾何關(guān)系，可以精確測量出圖像中每一點的三維位置，生成該地區(qū)的DEM。InSAR影像提取DEM的算法主要包括：復(fù)圖象的配準(zhǔn)、干涉圖的生成、相位解纏、數(shù)字高程模型重建。圖2是雷達信號到DEM產(chǎn)品的處理流程。

圖2 DEM重建流程

2 研究區(qū)使用SAR數(shù)據(jù)

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

表1 研究區(qū)SAR數(shù)據(jù)列表

表2 雷達數(shù)據(jù)干涉處理生成的干涉對

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

圖3 InSAR 生成的DEM

3 結(jié)果及分析

3.1 不同時間段的地形變化

3.1.1 地形變化最大區(qū)域干涉處理結(jié)果

不同時期DEM在坐標(biāo)統(tǒng)一的基礎(chǔ)上，兩個DEM作差進行比較，圖4中(a)表示1997年10月與1996年1月的地形變化，圖4中(b)表示2004年1月與1997年10月的地形變化，從結(jié)果圖中可以明顯看出，在1997年10月與1996年1月期間，研究區(qū)內(nèi)有明顯地形變化區(qū)，由于2004年的干涉結(jié)果受到失相干因素的影響存在數(shù)據(jù)缺失和空值現(xiàn)象，為了對2004年1月與1997年10月期間DEM進行比較，試驗中對數(shù)據(jù)缺失部分的地形數(shù)據(jù)采用克里金插值方法進行了填充。在圖4中還存在大范圍的趨勢性變化，可能是由SAR信號傳播過程中的大氣延遲誤差引起的誤差，可以忽略。

圖4 不同時間段的地形變化檢測結(jié)果(°)

3.1.2 地形變化統(tǒng)計分析

對圖4中1996年1月與1997年10月的地形重點變化區(qū)域進行統(tǒng)計，結(jié)果如圖5所示。地形抬升變化主要分布在兩個高度段，4 m以內(nèi)和6～14 m，沿著圖4(a)中a和b提取了地形變化曲線，如圖6和圖7所示。

圖5 1997年與1996年地形變化統(tǒng)計結(jié)果

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

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

按照實際邊界范圍裁剪1996年1月的DEM與1997年10月的DEM，求差得到區(qū)域變化值展繪到光學(xué)影像上，如圖8所示，同樣裁剪1997年10月的DEM與2004年1月的DEM，變化后地形如圖9所示。根據(jù)體積變化的像元進行統(tǒng)計得到，地形變化區(qū)的面積約為18 800 m2；根據(jù)地形變化區(qū)域的面積在不同像元上高程的變化值，計算得到1996—2004年該地區(qū)地形變化的體積約為146.25×104m3。

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

圖9 1997年10月至2004年1月地形變化圖

4 結(jié) 論

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