一種InSAR大氣相位建模與估計(jì)方法

占文俊， 李志偉*， 韋建超， 朱建軍， 汪長(zhǎng)城

1 中南大學(xué) 地球科學(xué)與信息物理學(xué)院， 長(zhǎng)沙 410083 2 中南大學(xué)有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)沙 410083

?

一種InSAR大氣相位建模與估計(jì)方法

占文俊1,2， 李志偉1,2*， 韋建超1,2， 朱建軍1,2， 汪長(zhǎng)城1,2

1 中南大學(xué) 地球科學(xué)與信息物理學(xué)院， 長(zhǎng)沙 410083 2 中南大學(xué)有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)沙 410083

為了削弱大氣延遲對(duì)干涉結(jié)果的影響以提高InSAR的測(cè)量能力，本文在InSAR大氣相位特征分析的基礎(chǔ)上，研究了一種新的InSAR大氣相位建模與估計(jì)方法.首先采用穩(wěn)健估計(jì)確定大氣垂直分層部分的模型參數(shù)，然后利用基于Matern模型的Kriging插值估計(jì)大氣紊流部分，最后應(yīng)用估計(jì)的大氣垂直分層和紊流資料改正InSAR測(cè)量結(jié)果.利用覆蓋河南義馬地區(qū)的ASAR數(shù)據(jù)對(duì)本文提出的方法進(jìn)行了驗(yàn)證，結(jié)果表明去除大氣影響后，InSAR重建的DEM與參考DEM的高程差異的均方誤差由19.5 m降至5.3 m，精度提高了約72%.同時(shí)，改正后的干涉圖更合理地揭示了義馬礦區(qū)的沉降漏斗情況，進(jìn)一步驗(yàn)證了本文方法的有效性.

合成孔徑雷達(dá)干涉；大氣相位；建模與估計(jì)；穩(wěn)健估計(jì)；Matern模型

1 引言

雷達(dá)對(duì)地觀測(cè)時(shí)，信號(hào)會(huì)受到大氣折射的影響而發(fā)生傳播延遲，已成為制約InSAR精度的關(guān)鍵(Hanssen，2001；Li et al.，2007).這種影響主要來(lái)自對(duì)流層中水汽含量時(shí)空分布的不均勻性，目前很難用一種確定性的方法來(lái)模擬并去除.1994年，Massonnet等首次發(fā)現(xiàn)了InSAR中的大氣影響，至今已有不少學(xué)者致力于InSAR大氣改正的研究.這些研究大致可以分為兩類：一類是采用外部數(shù)據(jù)直接進(jìn)行改正，另一類是利用InSAR數(shù)據(jù)本身進(jìn)行大氣相位估計(jì)及改正.然而，受站點(diǎn)距離、空間分辨率以及觀測(cè)條件(如云層)的限制，采用GPS、MERIS、MODIS、FY-1C等外部數(shù)據(jù)進(jìn)行大氣改正時(shí)需要合理的插值策略(Li et al.，2004；Li et al.，2006a；Li et al.，2006b).其中，以O(shè)nn+Kriging模型為代表的顧及了大氣延遲隨地形變化的插值算法效果最佳(Xu et al.，2011；Li et al.，2012).同時(shí)，stacking、大氣相位跟高程的相關(guān)分析、PS-InSAR和SBAS-InSAR等利用InSAR數(shù)據(jù)本身進(jìn)行大氣相位改正的方法，均須以大氣的時(shí)空分布特性為前提假設(shè)(Zebker et al.，1997；Ferretti et al.，2001；Beaudueel et al.，2000).

基于上述現(xiàn)狀，通過(guò)研究InSAR中大氣相位特征以獲取相關(guān)的先驗(yàn)信息，然后開(kāi)展大氣改正，有助于削弱InSAR中大氣延遲影響.本文在InSAR大氣相位特征分析的基礎(chǔ)上，研究了一種新的InSAR大氣相位建模與估計(jì)方法，包括采用穩(wěn)健估計(jì)確定大氣垂直分層部分的模型參數(shù)和利用基于Matern模型的Kriging插值估計(jì)大氣紊流部分，并通過(guò)真實(shí)數(shù)據(jù)實(shí)驗(yàn)對(duì)新方法進(jìn)行了驗(yàn)證.

2 InSAR大氣相位特征分析

干涉圖中的大氣延遲是兩次SAR成像時(shí)刻大氣折射引起的路徑延遲之差.由于每次雷達(dá)成像時(shí)大氣分布狀況都不相同，因此，大氣延遲影響不可避免.一般地，地表任一點(diǎn)s(x,y)在雷達(dá)成像時(shí)刻t的大氣延遲相位τ(s)可以表示為折射率N(s,z)在湍流有效高度H內(nèi)沿雷達(dá)視線方向(LOS)積分的結(jié)果(Hanssen，2001)(如圖1)：

圖1 雷達(dá)成像時(shí)大氣延遲示意圖Fig.1 Sketch map of atmospheric delay at radar acquisition time

(1)

式中，λ為雷達(dá)波長(zhǎng)，θ為入射角.圖1中虛線方格代表湍流有效高度H內(nèi)折射率N(s,z)的三維分布，其中z為湍流高度值.

(2)式中，hp、hr和sp、sr分別為點(diǎn)r、p對(duì)應(yīng)的高程和平面位置.第一項(xiàng)為垂直分層部分Δτtopography，第二項(xiàng)為水平紊流部分Δτturbulence.后續(xù)表達(dá)中，假設(shè)τ為干涉圖中差分大氣延遲，即τ=τ(tmaster)-τ(tslave).

為了分析不同地形條件下干涉圖中大氣相位特性，以上海、南加州和Etna火山三個(gè)典型區(qū)域的ERSTandem干涉對(duì)為例進(jìn)行分析，SAR數(shù)據(jù)基本參數(shù)見(jiàn)表1.其中，上海和Etna火山代表典型的沖積平原和高山地形，且臨近海域，大氣影響復(fù)雜；而南加州地區(qū)地勢(shì)略有起伏、構(gòu)造和人類活動(dòng)復(fù)雜，是InSAR研究的典型區(qū)域.由于各干涉對(duì)時(shí)間間隔均僅一天，可以忽略形變信號(hào)對(duì)干涉結(jié)果的影響，則差分干涉結(jié)果主要反映差分大氣相位的變化，如圖2.其中地形相位采用90m分辨率的SRTMDEM進(jìn)行模擬并去除.圖2c中可以看出，Etna火山地區(qū)的大氣延遲幅度最大可以達(dá)到約16rad，換算成LOS向的形變值約為7.1cm，這可能會(huì)完全掩蓋火山運(yùn)動(dòng)引起的地表形變.

表1 SAR影像基本參數(shù)Table 1 Basic parameters of SAR images

2.1 垂直分層效應(yīng)

垂直分層效應(yīng)是由于雷達(dá)成像時(shí)刻沿垂直方向上大氣折射率變化引起.由于低層大氣特別是水汽密度隨海拔高度增加而遞減，因此大氣延遲中垂直分層效應(yīng)并不是隨機(jī)的，可以認(rèn)為是高程的函數(shù)，這一點(diǎn)也可由公式(2)推出.Taylor和Peltzer等研究指出，干涉SAR中差分大氣延遲相位與高程之間呈線性或指數(shù)遞減關(guān)系(Taylor and Peltzer，2006；Peltzer et al.，2006).

圖3利用差分干涉相位與對(duì)應(yīng)的SRTM DEM分別進(jìn)行了線性和指數(shù)模型回歸：

(3)

(4)

其中，τ(h)為高程h處回歸的大氣延遲，a0、a1、b0、b1和b2為待確定的模型參數(shù).從圖3中可以看出，地形起伏越大，垂直分層現(xiàn)象越明顯，相位與高程的擬合度越高，且指數(shù)模型相對(duì)于傳統(tǒng)的線性模型能以更高的相關(guān)度擬合差分干涉相位.比如，上海地區(qū)地形平坦，水汽分布呈現(xiàn)隨機(jī)性，大氣紊流效應(yīng)表現(xiàn)明顯(如圖3a)，所以不論是線性模型還是指數(shù)模型，其與高程的擬合系數(shù)均很低；而Etna火山地區(qū)地形起伏較大，水汽分布隨地形顯著變化，大氣垂直分層效應(yīng)表現(xiàn)明顯(如圖3c)，所以相位與高程的擬合程度最好，擬合系數(shù)最高；而南加州地區(qū)地形介于二者之間，既有較明顯的大氣紊流效應(yīng)，也有顯著的垂直分層效應(yīng)，擬合系數(shù)也處于二者之間(圖3b).

2.2 紊流效應(yīng)

紊流效應(yīng)是由大氣中的湍流過(guò)程引起大氣折射分布時(shí)空變化引起的.Tatarski研究表明，假設(shè)大氣紊流各向同性，大氣折射率N的空間分布服從Kolmogorov冪次紊流定律，其結(jié)構(gòu)函數(shù)DN(ρ)為(Tatarski，1961)：

(5)

(6)

圖2 地理編碼后的差分干涉圖(疊加地形等高線) (a) 上海; (b) 南加州; (c) Etna火山.Fig.2 Differential interferogram (overlaying topographic contours) (a) Shanghai; (b) South California; (c) Etna volcano.

圖3 回歸結(jié)果圖(黑線：線性模型，紅線：指數(shù)模型) (a) 上海；(b) 南加州；(c) Etna火山.Fig.3 The regression results (black line: linear model; red line: exponential model) (a) Shanghai; (b) South California; (c) Etna volcano.

(7)采用公式(7)和大氣紊流樣本可計(jì)算實(shí)驗(yàn)變差函數(shù)樣本值，同時(shí)根據(jù)大氣紊流部分結(jié)構(gòu)函數(shù)服從Kolmogorov冪次紊流定律，可以對(duì)實(shí)驗(yàn)變差函數(shù)樣本值進(jìn)行模型擬合.由Kriging插值原理可知，實(shí)驗(yàn)變差函數(shù)模型的最佳擬合是Kriging插值的關(guān)鍵.常用的理論變差函數(shù)模型如球狀模型、指數(shù)模型等已被廣泛使用，然而考慮到Matern模型可以很好地描述Kolmogorov紊流定律，且可以實(shí)現(xiàn)自適應(yīng)平滑，因此，本文利用Matern模型擬合實(shí)驗(yàn)變差函數(shù)樣本(KnospeandJonsson，2009)：

γ(ρ)=

(8)

其中，K是第二類修正Bessel函數(shù)，Γ是Gamma函數(shù)，d是平滑因子，a是變程，C是拱高，N是塊金值.由于干涉圖中大氣信號(hào)呈空間自相關(guān)，連續(xù)性較好，因此塊金值N通常趨于0.

在2.1節(jié)中，將回歸的參數(shù)值b0、b1和b2回代到指數(shù)模型中，根據(jù)影像范圍內(nèi)各像素點(diǎn)的高程值可計(jì)算垂直分層部分，并從差分干涉相位中減去，得到殘余的相位可認(rèn)為主要包含大氣紊流部分.對(duì)紊流部分根據(jù)公式(7)計(jì)算相應(yīng)的實(shí)驗(yàn)變差函數(shù)樣本值，并用球狀模型和Matern模型進(jìn)行擬合，結(jié)果如圖4.從圖中可以看出，Matern模型較傳統(tǒng)的球狀模型能更好的貼合樣本點(diǎn)，擬合相關(guān)度也更高，因此更好地描述大氣紊流狀態(tài).其中，Etna火山地區(qū)的變程相對(duì)較小(如圖4c)，這是因?yàn)镋tna火山地區(qū)大氣垂直分層效應(yīng)明顯，大氣紊流效應(yīng)較弱，大氣紊流部分的空間自相關(guān)性較差.

圖4 實(shí)驗(yàn)變差函數(shù)擬合(紅線：Matern模型，黑線：球狀模型) (a)上海；(b)南加州；(c)Etna火山地區(qū).Fig.4 Variogram fitting (red line: Matern model; black line: spherical model) (a) Shanghai; (b) South California; (c): Etna volcano.

3 大氣相位建模與估計(jì)

本文擬對(duì)SAR干涉圖中不含形變信號(hào)的區(qū)域進(jìn)行InSAR大氣相位的建模與估計(jì).采用“二軌法”對(duì)干涉對(duì)進(jìn)行差分干涉處理，則干涉圖中不含形變區(qū)域的差分干涉相位由大氣相位、地形殘差相位和失相干噪聲組成.其中，大氣相位與其他相位分量的空間分布特性不同，因此可以采取一定方法進(jìn)行估計(jì).

前文分析指出，對(duì)InSAR大氣相位垂直分層部分的建模，指數(shù)模型要優(yōu)于線性模型；同時(shí)，在描述InSAR大氣相位紊流部分的結(jié)構(gòu)函數(shù)時(shí)，Matern模型較傳統(tǒng)的球狀模型更為理想.本文將基于這些分析結(jié)果建立更高精度的重軌InSAR大氣相位模型.為了便于程序?qū)崿F(xiàn)，且二次多項(xiàng)式可以很好地逼近指數(shù)函數(shù)曲線，對(duì)公式(4)在h=0處進(jìn)行Taylor級(jí)數(shù)展開(kāi)并取前三項(xiàng)，建立大氣相位數(shù)學(xué)模型：

(9)

式中，a0、a1、a2為大氣垂直分層部分待確定的參數(shù)；τturbulence為大氣紊流部分，利用Matern模型描述其結(jié)構(gòu)函數(shù).

針對(duì)這兩部分大氣相位的分布特性，分別采取相應(yīng)的方法進(jìn)行估計(jì).詳細(xì)流程如圖5，通過(guò)對(duì)比干涉圖中不含形變區(qū)域大氣改正前后InSAR重建的DEM的精度來(lái)間接評(píng)定大氣相位估計(jì)結(jié)果的精度和可靠性.其中，針對(duì)InSAR重建DEM過(guò)程中相位解纏難的問(wèn)題，采用“高程補(bǔ)償”方法降低干涉條紋率，保證相鄰像素的干涉相位之差在(-π,π)之間(吳宏安等，2009).

3.1 估計(jì)垂直分層部分

采用穩(wěn)健估計(jì)中的選權(quán)迭代法，實(shí)現(xiàn)參數(shù)ak(k=0,1,2)的穩(wěn)健估計(jì)，以求得大氣相位垂直分層部分.其中心思想是采用加權(quán)最小二乘法(WLSE)估計(jì)模型參數(shù)，根據(jù)殘差的大小反復(fù)迭代選權(quán)，實(shí)現(xiàn)參數(shù)估計(jì)的穩(wěn)健性(周江文，1989).設(shè)定優(yōu)化的目標(biāo)函數(shù)如下：

min =Q(a0,a1,a2)

(10)

(11)

(12)

圖5 大氣相位估計(jì)與改正流程Fig.5 The flow chart of atmospheric phase estimation and correction

以認(rèn)為迭代收斂.

由此確定模型參數(shù)ak(k=0,1,2)后，將研究區(qū)域內(nèi)各像素點(diǎn)的高程值代入公式(9)，估計(jì)干涉圖中大氣相位垂直分層部分.

3.2 估計(jì)紊流部分

前已述及，在描述InSAR大氣相位紊流部分的結(jié)構(gòu)函數(shù)時(shí)，Matern模型較傳統(tǒng)的球狀模型更為理想.由于在估計(jì)InSAR大氣相位模型時(shí)，一般只選取一些有代表性的高相干點(diǎn)，這樣還需把他們插值到整個(gè)干涉圖.當(dāng)樣本數(shù)n足夠大時(shí)，采用基于Matern模型的Kriging插值算法將扣除垂直分層部分后的殘余相位插值到整個(gè)SAR干涉圖空間.由于地形殘差和失相關(guān)噪聲在空間上都表現(xiàn)為高頻信號(hào)，而Kriging插值等效加權(quán)平均過(guò)程，因此可以認(rèn)為插值結(jié)果為大氣相位紊流部分(萬(wàn)青等，2012).

4 實(shí)例研究與結(jié)果分析

4.1 研究區(qū)域

以覆蓋河南部分地區(qū)的ASAR影像為干涉數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，影像參數(shù)如表2，覆蓋范圍如圖6中虛線方框所示，圖中地形由90 m分辨率的SRTM DEM提供.該干涉對(duì)垂直基線為160 m，對(duì)應(yīng)模糊高約58 m，半個(gè)周期的大氣傳播延遲將產(chǎn)生約29 m的DEM偏差.而且，由于研究區(qū)域高程差最大達(dá)到1 km以上，地形起伏較大造成水汽分布不均，大氣延遲對(duì)該區(qū)域開(kāi)展干涉測(cè)量的影響較為嚴(yán)重，因此，有必要對(duì)干涉結(jié)果進(jìn)行大氣改正.

表2 ASAR影像基本參數(shù)Table 2 Basic parameters of the ASAR images

4.2 數(shù)據(jù)處理

對(duì)該干涉對(duì)進(jìn)行“二軌法”差分干涉處理，為了抑制干涉圖中失相干噪聲的影響，方位向(azimuth)和距離向(range)分別做10∶2的多視平均，處理后的地面分辨率約為40 m×40 m，并采改進(jìn)的Goldstein濾波算法進(jìn)一步降低干涉圖中的噪聲(Li et al.，2008).干涉圖中的地形相位采用90m分辨率的SRTM DEM進(jìn)行模擬并去除，圖7a為相位解纏后的差分干涉圖.其中，為了確保干涉圖與地形的對(duì)應(yīng)關(guān)系，本文利用模擬的SAR幅度圖與真實(shí)的SAR幅度圖之間的配準(zhǔn)多項(xiàng)式對(duì)初始地理編碼表進(jìn)行了精化.

由于兩景影像僅間隔70天，差分干涉結(jié)果可不計(jì)地表形變的影響.然而，義馬市境內(nèi)礦區(qū)較多，有代表性的如圖7a中虛線方框標(biāo)記的義馬礦區(qū)帶，造成短時(shí)間內(nèi)局部地表的較大沉降.此外，盡管研究區(qū)域植被茂盛，但由于影像獲取時(shí)間在冬季，干涉圖仍保持非常好的相干性(圖7b).

為了盡可能消除形變信號(hào)以及相位信息不可靠性的低相干點(diǎn)的影響，文中掩膜掉義馬礦區(qū)帶后，依據(jù)相干圖，定義9×9的窗口，選取搜索窗口中相干值最大的像素點(diǎn)作為大氣相位估計(jì)的樣本點(diǎn)，共選取的樣本數(shù)為4.9284×104.樣本點(diǎn)上的差分干涉相位由大氣相位、地形殘差相位和失相干噪聲組成.建立如公式(10)的優(yōu)化目標(biāo)函數(shù)，利用穩(wěn)健估計(jì)的方法得到模型參數(shù)a0=9.4123，a1=-0.1074×10-3，a2=2.8705×10-6.然后，根據(jù)干涉圖中各像素的高程值可以估計(jì)大氣相位垂直分層部分(如圖8a).

圖6 研究區(qū)域的地形(虛線矩形框?yàn)橛跋窀采w范圍)Fig.6 Topography of study area (The dashed rectangle indicates the coverage of the ASAR image)

對(duì)殘余相位進(jìn)行地學(xué)統(tǒng)計(jì)分析，估計(jì)大氣相位紊流部分.首先，對(duì)殘余相位計(jì)算各樣本點(diǎn)對(duì)(ui,uj)的實(shí)驗(yàn)變差函數(shù)值γ*(ρij)(其中，ρij=range(ui-uj)),并利用Matern模型擬合實(shí)驗(yàn)變差函數(shù)值以求得模型參數(shù)N、C、d和a(如圖8b).然后，將參數(shù)回代到Matern模型中估計(jì)出任一點(diǎn)對(duì)(ui,uj)的變差函數(shù)值γ(ρij).最后，采用Kriging插值算法將上述選取的高相干點(diǎn)上的殘余相位重采樣到SAR格網(wǎng)空間，插值結(jié)果可以認(rèn)為是大氣相位紊流部分(如圖8c).Kriging插值算法公式如下：

圖7 差分干涉處理結(jié)果(雷達(dá)坐標(biāo)系下) (a) 差分干涉圖；(b) 相干圖.Fig.7 Results of Differential interferometry (Radar-geocoded). (a) Unwrapped differential interferogram; (b) Coherence map.

圖8 大氣相位估計(jì) (a)垂直分層部分；(b)Matern模型擬合實(shí)驗(yàn)變差函數(shù)；(c)紊流部分.Fig.8 Estimation of atmospheric phase (a) The stratified atmospheric component; (b) Variogram fitting with Matern model; (c) The turbulent atmospheric component.

(13)

(14)

其中，μ為拉格朗日常數(shù)，γ(ρij)為Matern模型估計(jì)的變異函數(shù)值.

至此，已經(jīng)估計(jì)出了研究區(qū)域各像素點(diǎn)的大氣相位垂直分層部分和紊流部分.將兩者疊加在一起，即為估計(jì)出的大氣相位，如圖9.

4.3 結(jié)果與分析

圖9 估計(jì)出的大氣相位圖Fig.9 Map of estimated atmospheric phase

(15)

式中估計(jì)的高程殘差hres實(shí)際上反映了InSARDEM與SRTM高程值之間的差值，因此，將高程殘差hres補(bǔ)償?shù)絊RTMDEM中，即為大氣改正后InSAR重建的DEM.圖10a為地理編碼到WGS-84坐標(biāo)系的InSARDEM，圖9c為圖10a中虛線方框所示的局部放大圖.圖10b為大氣改正前InSAR重建的DEM的局部放大圖，由包含大氣影響的差分干涉相位經(jīng)相位-高程轉(zhuǎn)換和地理編碼得到.

圖10 大氣改正前后InSAR重建的DEM比較(單位：m) (a)大氣改正后的InSAR DEM；(b)大氣改正前的InSAR DEM局部放大圖；(c)大氣改正后的InSAR DEM局部放大圖； (d)Aster GDEM局部放大圖.圖a中的虛線即為圖b,c和d所示DEM的范圍，紅線為截取的高程剖面.Fig.10 Comparisons between the InSAR re-constructed DEM before and after atmospheric correction (unit:m) (a) InSAR DEM after atmospheric correction; (b) Enlarged map of the InSAR DEM before atmospheric correction; (c) Enlarged map of the InSAR DEM after atmospheric correction; (d) Aster GDEM. The dashed line in Fig.10a is the spatial extent of the DEM shown in Fig.10b, c and d, and the red solid line is the location where elevation profile extracted.

本文以日本空間局的AsterGDEM作為參考值(圖10d)，對(duì)InSAR重建的DEM的精度進(jìn)行評(píng)價(jià).AsterGDEM的地面分辨率約30m，因此，需將其采樣到InSARDEM格網(wǎng)大小.從圖10(b—d)可以看出，三者的趨勢(shì)大致相同，經(jīng)大氣改正后，InSARDEM的細(xì)節(jié)信息更接近AsterGDEM.圖11a為截取的高程剖面的對(duì)比圖，大氣改正前的InSARDEM平均下降了31.6m左右；而4.2節(jié)中估計(jì)的大氣相位在虛線方框內(nèi)的均值約為2.8rad，將產(chǎn)生約為25.8m的DEM偏差.兩者的分析結(jié)果較為接近.圖11b為大氣改正前后InSAR重建的DEM與AsterGDEM高程差的統(tǒng)計(jì)直方圖擬合曲線.從圖中可以看出，經(jīng)過(guò)大氣改正后，重建的DEM與AsterGDEM的標(biāo)準(zhǔn)差由19.5m減至5.3m，精度提高了約72%.

同時(shí)，我們也對(duì)圖7a方框所示的義馬礦區(qū)帶地表沉降情況進(jìn)行了大氣改正前后的對(duì)比分析.圖12(a,b)為義馬礦區(qū)帶地表沉降經(jīng)大氣改正前后的對(duì)比圖.圖12(c,d)為截取其中兩個(gè)沉降漏斗的剖面圖(剖面位置見(jiàn)圖a)，可以看出，大氣改正前沉降漏斗的形變量大部分呈正值, 顯然不符合沉降的物理規(guī)律，大氣改正后正值基本消失，較好地揭示了該時(shí)間段礦區(qū)沉降情況.同時(shí)，大氣改正后a—b和c—d剖面的沉降量分別增加了17.1mm和21.4mm左右，而4.2節(jié)中估計(jì)的大氣相位在圖7a虛線方框內(nèi)的均值約為2.9rad，將產(chǎn)生約為25.8mm的形變誤差.以上兩者的分析結(jié)果較為接近.

5 結(jié)論

本文在InSAR大氣相位特性分析的基礎(chǔ)上，研究了一種新的InSAR大氣相位建模與估計(jì)方法.該方法首先采用穩(wěn)健估計(jì)確定大氣垂直分層部分的模型參數(shù)，然后運(yùn)用基于Matern模型的Kriging插值估計(jì)大氣紊流部分，最后利用估計(jì)的大氣垂直分層和紊流資料改正InSAR測(cè)量結(jié)果.實(shí)驗(yàn)結(jié)果顯示該方法準(zhǔn)確估計(jì)出了干涉結(jié)果中大氣相位分量，進(jìn)而大大提高了InSAR重建的DEM的精度.相比于Beauducel、Remy等學(xué)者通過(guò)相位與高程的相關(guān)分析去除大氣垂直分層部分的方法(Beaudueeletal.，2000；Remyetal.，2003)，本文方法同時(shí)去除了垂直分層部分和隨機(jī)的紊流部分，并且在紊流部分估計(jì)中采用Matern模型代替了傳統(tǒng)的變異函數(shù)模型.

圖11 (a)高程剖面對(duì)比圖；(b)大氣改正前后InSAR重建的DEM與Aster GDEM的高程差分布Fig.11 (a) Comparisons of different elevation profiles; (b) Distribution of the elevation differences between the Aster GDEM and the InSAR reconstructed DEM with/without atmospheric phase correction

圖12 (a)大氣改正前義馬礦區(qū)帶地表沉降；(b)大氣改正后義馬礦區(qū)帶地表沉降； (c)a—b剖面大氣改正前后的形變對(duì)比；(d)c—d剖面大氣改正前后的形變對(duì)比Fig.12 (a) Displacement of Yima mining area without atmospheric correction; (b) Displacement of Yima mining area with atmospheric correction; (c) Comparison of the displacements along profile a—b with/without atmospheric correction;(d) Comparison of the displacements along profile c—d with/without atmospheric correction

在InSAR形變監(jiān)測(cè)中，可以利用非形變區(qū)域的干涉相位進(jìn)行大氣相位建模與估計(jì)，從而對(duì)研究區(qū)域內(nèi)的像素逐個(gè)進(jìn)行大氣改正.本文實(shí)驗(yàn)選用的干涉圖中，形變區(qū)域集中在礦區(qū)，具有一定的特殊性.然而，在一般情況下，形變區(qū)域是大范圍的，且形變和大氣信號(hào)是相互混淆的，在無(wú)先驗(yàn)知識(shí)的情況下，如何有效地選取非形變區(qū)域，這需要進(jìn)一步利用統(tǒng)計(jì)學(xué)等相關(guān)知識(shí)進(jìn)行分析.

致謝 感謝歐空局提供的ENVISATASAR(AO-4458)和ERS-Tandem(F-2979)數(shù)據(jù).

Beaudueel B, Briole P, Froger J L. 2000. Volcano-wide fringes in ERS synthetic Aperture radar interferograms of ETNA (1992-1998): Deformation or tropospheric effect?.JournalofGeophysicalResearch, 105(B7): 16391-16402.Ferretti A, Prati C, Rocca F. 2001. Permanent scatterers in SAR interferometry.IEEETransactionsonGeoscienceandRemoteSensing, 2001, 39(1): 8-20.

Hanssen R H. 2001. Radar Interferometry-Data Interpretation and Error Analysis. Netherlands: Springer Kluwer Academic Publishers.Knospe S, S Jonsson. 2009. Covariance estimation for dInSAR surface deformation measurements in the presence of anisotropic atmospheric noise.IEEETransactionsonGeoscienceandRemoteSensing, 48(4): 2057-2065.Li Z W, Ding X L, Liu G X. 2004. Modeling atmospheric effects on InSAR with meteorological and continuous GPS observations: algorithms and some test results.JournalofAtmosphericandSolar-TerrestrialPhysics, 66(11): 907-917.

Li Z W, Ding X L, Huang C, et al. 2006a. Modeling of atmospheric effects on InSAR measurements by incorporating terrain elevation information.JournalofAtmosphericandSolar-TerrestrialPhysics, 68(11): 1189-1194.Li Z H, Fielding E J, Cross P, et al. 2006b. Interferometric synthetic aperture radar atmospheric correction: GPS topography-dependent turbulence model.JournalofGeophysicalResearch, 111

(B2): B024048113-8127.

Li Z W, Ding X L, Huang C, et al. 2007. Atmospheric effects on repeat-pass InSAR measurements over Shanghai region.JournalofAtmosphericandSolar-TerrestrialPhysics, 69(12): 1344-1356.Li Z W, Ding X L, Huang C, et al. 2008. Improved filtering parameter determination for the goldstein radar iInterferogram filter.ISPRSJournalofPhotogrammetryandRemoteSensing, 63(6): 621-634.Li Z W, Xu, W B, Feng, G C, et al. 2012. Correcting atmospheric effects on InSAR with MERIS water vapour data and elevation-dependent interpolation model.GeophysicalJournalInternational, 189(2): 898-910.Peltzer G, Socquet A, Lasserre C, et al. 2006. InSAR observations of interseimsmic strain along the central Altyn Tagh fault consistent with Holocene slip-rate. ∥EOS trans., AGU, 87(52), Fall Meet. Suppl., Abstract T21E-02P.,

Remy D, Bonvalot S, Briole P, et al. 2003. Accurate measurements of tropospheric Effects in volcanic areas from SAR interferometry data: application to Sakurajima volceano (Japan).EarthandPlanetaryScienceLetters, 213(3-4): 299-310.Tatarski V I. 1961. Wave propagation in a turbulent medium. Silverman R A Trans. New York: Mc Graw-Hill. Translated from Russian,.

Taylor M, Peltzer G. 2006. Current slip rates on conjugate strike-slip faults in central Tibet using synthetic aperture radar interferometry.JournalofGeophysicalResearch, 111(B12): doi: 10.1029/2005JB004014.402-411.

Wan Q, Zhang L, Jiang H J, et al. 2012. Estimation and removal of atmospheric effects in InSAR topographic mapping.JournalofRemoteSensing(in Chinese), 16(5): 1075-1089.

Wu H A, Zhang H, Wang C, et al. 2009. A high resolution InSAR topographic height reconstruction algorithm in rugged terrain based on SRTM DEM.JournalofRemoteSensing(in Chinese), 13(1): 145-151.Xu B, Li Z W, Wang Q J, et al. 2013. A rRefined sStrategy for rRemoving cComposite eErrors of SAR iInterferogram.IEEEGeoscienceandRemoteSensingLetters, 11 (1): 143-147.

Xu W B, Li Z W, Ding X L, et al. 2011. Interpolating atmospheric water vapor delay by incorporating terrain elevation information.JournalofGeodesy, 85 (9): 555-564.Zebker H A, Rosen P A, Hensley S. 1997. Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps.JournalofGeophysicalResearch, 102(B4): 7547-7563.Zhou J W. 1989. Classical theory of errors and robust estimation.ActaGeodaeticaetCartographicaSinica(in Chinese), 18(2): 115-120.

附中文參考文獻(xiàn)

萬(wàn)青, 張路, 蔣厚軍, 等. 2012. InSAR地形重建中大氣效應(yīng)的估計(jì)和去除. 遙感學(xué)報(bào), 16 (5): 10745-10889.

吳宏安, 張紅, 王超, 等. 2009. 基于SRTM DEM的InSAR高分辨率山區(qū)地表高程重建算法. 遙感學(xué)報(bào), 13(1): 145-151.

周江文. 1989. 經(jīng)典誤差理論與抗差估計(jì). 測(cè)繪學(xué)報(bào), 18(2): 115-120.

(本文編輯 汪海英)

A strategy for modeling and estimating atmospheric phase of SAR interferogram

ZHAN Wen-Jun1,2, LI Zhi-Wei1,2*, WEI Jian-Chao1,2, ZHU Jian-Jun1,2, WANG Chang-Cheng1,2

1SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,ChangSha410083,China2HunanKeyLaboratoryofNon-ferrousResourcesandGeologicalHazardDetection,CentralSouthUniversity,Changsha410083,China

The atmospheric effect is one of the most important limiting factors to Interferometric Synthetic Aperture Radar (InSAR) measurements. It can seriously degrade the quality of InSAR measurements or even render the technology unusable. It is of great importance to develop methods to model and mitigate this effect.A novel method for modeling and estimating the atmospheric phase in SAR interferograms is proposed. The method starts with the spatial and temporal analysis of the characteristics of the atmospheric effects in InSAR. Based on this, it constructs relevant models to describe the stratification and turbulent mixing of atmospheric effects, respectively. Then, it develops the method of robust estimation to determine the model parameters of the stratified atmospheric components, and the method of Matern variogram model-based Kriging interpolation to estimate the parameters of turbulent atmospheric components. Finally, it applies the developed method to estimate and correct the atmospheric effects in InSAR measurements.The proposed method is validated with one ASAR pair over the Yima area in Henan Province. The results show that after the atmospheric correction, the root mean square error of the differences between the InSAR-reconstructed and the reference DEM reduce from 19.5 m to 5.3 m, representing an improvement by 72%. In addition, after the correction, the sign of the line-of-sigh (LOS) range change in a mining area varies from positive to negative, indicating that subsidence rather than uplift happening in this area. The corrected interferogram much better reveals the development of the “subsidence bowl” in the mining area.This paper developed a novel method for modeling and mitigating atmospheric effects in InSAR. The method fully exploits the spatial characteristics of atmospheric effects in InSAR, and considers the stratified and the turbulent mixing atmospheric effects as well. Besides, it adopts the Matern model rather than the traditional variogram model in the turbulent mixing modeling. So it works very well. Future work will focus on validating the method in different study areas.

InSAR； Atmospheric phase； Modeling and estimation； Robust estimation； Matern model

10.6038/cjg20150710.

國(guó)家高技術(shù)研究發(fā)展計(jì)劃項(xiàng)目(2012AA121301)，國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(2012CB719903)，國(guó)家自然科學(xué)基金(41222027，41474007，41404013)，湖南省杰出青年科學(xué)基金(13JJ1006)，教育部博士點(diǎn)基金(20130162110015)資助.

占文俊，男，1987年生，碩士研究生，研究方向?yàn)镮nSAR大氣噪聲建模及大氣延遲改正.

*通訊作者 李志偉，男，1974年生，博士，博士生導(dǎo)師，測(cè)繪與遙感科學(xué)系主任.主要從事InSAR大地測(cè)量與遙感研究.E-mail:zwli@csu.edu.cn

10.6038/cjg20150710

P225

2014-02-12，2015-07-07收修定稿

占文俊，李志偉，韋建超等. 2015. 一種InSAR大氣相位建模與估計(jì)方法.地球物理學(xué)報(bào)，58(7):2320-2329,

Zhan W J, Li Z W, Wei J C, et al. 2015. A strategy for modeling and estimating atmospheric phase of SAR.ChineseJ.Geophys. (in Chinese),58(7):2320-2329,doi:10.6038/cjg20150710.