侯建國，白 駒,陶紅志，劉 江，王 強(qiáng)

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

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

1 InSAR高程測(cè)量精度及生成DEM基本流程

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

1.1 高程測(cè)量精度分析

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

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

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

(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測(cè)量精度越高。從式(2)可知，最終高程測(cè)量精度與dr誤差之間的關(guān)系:

(4)

寫成中誤差形式為：

(5)

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

(6)

相位測(cè)量誤差mφ算式為：

(7)

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

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

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

圖2 DEM重建流程

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

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

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

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

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

圖3 InSAR 生成的DEM

3 結(jié)果及分析

3.1 不同時(shí)間段的地形變化

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

不同時(shí)期DEM在坐標(biāo)統(tǒng)一的基礎(chǔ)上，兩個(gè)DEM作差進(jìn)行比較，圖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)象，為了對(duì)2004年1月與1997年10月期間DEM進(jìn)行比較，試驗(yàn)中對(duì)數(shù)據(jù)缺失部分的地形數(shù)據(jù)采用克里金插值方法進(jìn)行了填充。在圖4中還存在大范圍的趨勢(shì)性變化，可能是由SAR信號(hào)傳播過程中的大氣延遲誤差引起的誤差，可以忽略。

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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