傅拓 ，張書畢

(1.中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇徐州 221116;2.中國礦業(yè)大學(xué)江蘇省測繪與國土信息工程重點(diǎn)實(shí)驗(yàn)室，江蘇徐州 221116)

1 引言

InSAR技術(shù)是近年來迅速發(fā)展起來的極具應(yīng)用價值的空間對地觀測新技術(shù)，具有監(jiān)測精度高、范圍大、成本低、空間連續(xù)覆蓋等優(yōu)點(diǎn)，為高邊坡地質(zhì)災(zāi)害區(qū)監(jiān)測提供了一種新型的監(jiān)測方法。但由于地質(zhì)災(zāi)害多發(fā)生在暴雨頻發(fā)、地質(zhì)地貌復(fù)雜的區(qū)域，特殊的地理位置與氣候使得InSAR技術(shù)應(yīng)用中受大氣延遲的影響非常嚴(yán)重，大氣延遲誤差是InSAR技術(shù)的主要誤差源，嚴(yán)重時還會導(dǎo)致InSAR結(jié)果的錯誤解譯［1］，因此，必須要對InSAR觀測進(jìn)行大氣延遲改正。

2 InSAR大氣延遲改正的研究

自1994年 Massonnet等［2］在利用 InSAR研究1992年發(fā)生在美國南加州地區(qū)Laders地震時首次發(fā)現(xiàn)了InSAR干涉圖中的大氣效應(yīng)后，眾多學(xué)者對消除或減弱InSAR中的大氣延遲的方法展開了研究。近年來，InSAR大氣延遲改正的理論及方法取得了快速的發(fā)展，國內(nèi)外提出了多種改進(jìn)InSAR大氣延遲的方法。以下介紹幾種主流的方法:

2.1 利用MODIS/MERIS水汽產(chǎn)品

MODIS(Moderate Resolution Imaging Spectroradiometer)［3］的近紅外水汽產(chǎn)品比世界上最密集的GPS網(wǎng)絡(luò)——美國南加州綜合GPS網(wǎng)(SCIGN)的密度還要高10倍以上。2005年國內(nèi)許多學(xué)者提出融合MODIS和 GPS 數(shù)據(jù)的 InSAR 干涉圖改正方法［4－5］。但MODIS容易受到云的影響，且水汽產(chǎn)品存在系統(tǒng)誤差，需要進(jìn)行改正。

MERIS(Medium Resolution Imaging Spectrometer Instrument)水汽產(chǎn)品與ASAR(Advanced Synthetic Aperture Radar)影像同時獲取，且得到的PWV空間分辨率高(300 m×300 m)，精度也比MODIS 還要高［6～7］，但其缺點(diǎn)是受云的影響更大。

2.2 利用PSInSAR技術(shù)

常規(guī)DInSAR技術(shù)進(jìn)行滑坡監(jiān)測會受到相位失相關(guān)和大氣延遲影響，F(xiàn)erretti［8～9］提出了僅僅跟蹤成像區(qū)域內(nèi)雷達(dá)散射特性較為穩(wěn)定的目標(biāo)而放棄那些失相關(guān)嚴(yán)重的分辨單元的方法。這些目標(biāo)(如建筑物的墻角、屋頂和裸露的巖石等)幾乎不受失相關(guān)噪聲影響，即使在多年時間間隔的干涉對中仍然保持較高的干涉相關(guān)性，把這些穩(wěn)定的目標(biāo)稱之為永久散射體(permanentscatterers，PS)。PSInSAR 技術(shù)［10～12］能克服傳統(tǒng)InSAR技術(shù)中的去相干性和大氣效應(yīng)等困難，但利用PSInSAR技術(shù)進(jìn)行InSAR大氣改正的缺點(diǎn)在于，PS點(diǎn)通常分布于市區(qū)等人工建筑較多的地方或無植被覆蓋的山峰、山脊等有裸露巨石的地方，對于研究山體滑坡等地質(zhì)災(zāi)害受到了很大的限制，要得到準(zhǔn)確可靠的結(jié)果，對同一區(qū)域內(nèi)的SAR影像數(shù)量也有較高要求，一般需要至少30幅影像。

2.3 利用地面GPS觀測

GPS具有高精度、全天候、可連續(xù)估計(jì)大氣延遲的優(yōu)點(diǎn)，GPS氣象學(xué)的發(fā)展為精確反演對流層延遲提供了保障，亦為估計(jì)與改正InSAR大氣效應(yīng)提供了一種新的方法，因而利用GPS進(jìn)行InSAR大氣延遲改正的方法得到了越來越多的應(yīng)用。但SAR影像與GPS測站的空間分辨率差異很大，很大程度上限制了大氣延遲分布圖的估計(jì)精度，為了對SAR影像進(jìn)行逐像素的大氣層延遲改正，必須對GPS獲取的大氣延遲進(jìn)行加密插值，Janssen等［13］對反距離加權(quán)平均(IDW)、克里格(Kriging)、三次樣條3種不同的插值算法用InSAR大氣改正的效果進(jìn)行了比較，結(jié)果表明反距離加權(quán)平均插值和克里格插值法比三次樣條插值算法效果更佳;李志偉［14］提出了一種新的大氣校正的方法，通過計(jì)算出平均大氣延遲，然后利用GPS探測到的大氣延遲對平均延遲進(jìn)行修正，彌補(bǔ)了SAR干涉圖水平方向上各向異性的問題。但GPS改正InSAR大氣延遲仍然存在限制，當(dāng)GPS測站附近的氣象觀測缺失時，大氣延遲估計(jì)精度會降低，GPS采集數(shù)據(jù)連續(xù)運(yùn)行站點(diǎn)之間的間隔一般為幾十千米不等，測站點(diǎn)過于分散會降低觀測精度，增加觀測站的密度會受到地理環(huán)境和運(yùn)作成本等因素的限制，即使設(shè)立臨時性、低成本的觀測點(diǎn)，效果仍然非常有限。針對以上幾種方法存在的問題，本文提出了一種新的改正InSAR大氣延遲的方法。

3 GPS和NCEP FNL數(shù)據(jù)改正InSAR大氣延遲

3.1 NCEP全球分析最終產(chǎn)品FNL(Final)

NCEP/NCAR再分析資料［15］是1991年開始實(shí)施的美國國家氣象中心(NMC)氣候資料同化系統(tǒng)(CDAS)的一個重要部分，其包含了豐富的氣象觀測資料，可用來精確估計(jì)氣象參數(shù)并提高大氣延遲分布圖的精度。該系統(tǒng)由NASA資助，以NCAR以前為航空與交通部門開發(fā)的用于美國大陸的產(chǎn)品為基礎(chǔ)，將衛(wèi)星資料、模式結(jié)果與先進(jìn)的人工智能技術(shù)結(jié)合起來，快速確定和預(yù)報(bào)風(fēng)暴與潛在擾動氣流。

NCEP FNL資料［16］比再分析資料具有更高的時間、空間分辨率，其同化了全球搜集到的幾乎所有觀測信息，是全球業(yè)務(wù)模型生成的最后產(chǎn)品，可以認(rèn)為是NCEP所有產(chǎn)品中用于模式存檔的最佳選擇。其數(shù)據(jù)水平格距為1°，水平方向共360×181個格點(diǎn)，垂直方向分為地表層和從1000～10hPa共26個等壓層，有海平面氣壓、表面溫度、海冰等9個單層變量和位勢高度、溫度、東西風(fēng)分量、南北風(fēng)分量、相對濕度共5個全層次變量。從1999年7月30日至今的NCEP FNL產(chǎn)品可以免費(fèi)獲取，其空間分辨率為1°×1°，時間分辨率為 6h(http://dss.ucar.edu/datasets/ds083.2/)。NCEP FNL數(shù)據(jù)集以GRIB格式存儲，同樣可用GrADS軟件進(jìn)行讀取。圖1為逐6h再分析資料集的分布圖。

圖1 逐6h再分析資料集的分布圖

3.2 基于“凝固流”假設(shè)的InSAR大氣改正方法

“凝固流”假設(shè)由泰勒于1938年提出［17］，該假設(shè)認(rèn)為在平均風(fēng)速的驅(qū)使下，空氣流在平移過程中處于“凝固”狀態(tài)，即在某一固定點(diǎn)上觀測到的時間上的延遲波動是由某一折射率場在平均速度為V的流的驅(qū)使下，經(jīng)過該點(diǎn)時產(chǎn)生的空間波動引起的，“凝固流”假設(shè)的主要思想是空氣流在平均風(fēng)速的驅(qū)使下進(jìn)行平流輸送。

“凝固流”假設(shè)運(yùn)用到空氣中的三維濕折射率Nwet(x，y，z，t)的數(shù)學(xué)描述可以表示為:(x，y，z)表示空間的三維坐標(biāo)，t表示時間，u，v分別表示水平方向的緯向風(fēng)和經(jīng)向風(fēng)。

“凝固流”假設(shè)研域平均風(fēng)速U→—空間上為常數(shù)，但隨著時間變化，則指定空間位置(x，y，z)上的濕折射率可以認(rèn)為是在之前時刻t，由另一空間位置(x－ut，yvt，z)上的濕折射率通過平均風(fēng)速U→—=(u，v)的平流傳輸?shù)玫降摹?/p>

3.3 平均風(fēng)速的估計(jì)

為了利用GPS天頂濕延遲(ZWD)時間序列估計(jì)平均風(fēng)速，假設(shè)研究區(qū)域內(nèi)的GPS天頂濕延遲(ZWD)在空間上的統(tǒng)計(jì)特性為各向同性，且在時間上服從廣義穩(wěn)態(tài)隨機(jī)過程［18］，則兩個GPS測站上時間連續(xù)的ZWD間的互相關(guān)函數(shù)可表示為:

R=［RxRy］T表示 GPS 測站間的空間距離，ρ、ρ+R分別表示GPS測站的空間位置，△t表示ZWD獲取的時間間隔，t、t+△t分別表示ZWD的獲取時刻，V為平均風(fēng)速。ZWD間的自相關(guān)函數(shù)可表示為:

連續(xù)時間上的二維空間ZWD場的互相關(guān)函數(shù)可以通過平均風(fēng)速與其自相關(guān)函數(shù)進(jìn)行轉(zhuǎn)換，將式(2)線性化有:

?表示Hamilton算子。

利用最小二乘原理，建立GPS ZWD序列的自相關(guān)函數(shù)模型與互相關(guān)之間的關(guān)系后，即根據(jù)式(4)解算出平均風(fēng)速V。ZWD在連續(xù)觀測時刻的互相關(guān)函數(shù)可表示為:

式中，S(Rk)={(ρa(bǔ)，ρb):ρa(bǔ)－ρb≈Rk};k 表示 GPS測站個數(shù)，Rk表示 ZWD 間的空間距離;ρa(bǔ)、ρb分別表示時刻的 ZWD 位置。

利用式(4)，平均風(fēng)速可以用自相關(guān)函數(shù)進(jìn)行估計(jì):

式中，tm、ts分別表示主副SAR影像的獲取時刻。

由于地面GPS測站的分辨率低，利用其插值得到的濕延遲分布不能真實(shí)反映濕延遲空間分布狀況。基于泰勒提出的“凝固流”假設(shè)，Onn［19］將這種思想用于InSAR大氣改正，并指出這種方法比僅用GPS觀測時更為有效。其中，準(zhǔn)確地估計(jì)研究區(qū)域的平均風(fēng)速是至關(guān)重要的。平均風(fēng)速的大小及方向決定了“擴(kuò)展”的ZWD控制點(diǎn)的位置，決定了基于這些控制點(diǎn)插值得到的濕延遲分布圖的精度，從而決定了該方法對In-SAR大氣延遲效應(yīng)的改正效果。

但是，將研究區(qū)域內(nèi)ZWD的變化全部歸因于平均風(fēng)速的變化是不嚴(yán)密的，GPS ZWD的變化是多方面因素的綜合結(jié)果，除了與風(fēng)速有關(guān)外，還與研究區(qū)域的溫度等其他氣象因素有關(guān)。

3.4 改進(jìn)的InSAR大氣改正方法

針對上述問題，本文提出了融合GPS和NCEP全球分析最終產(chǎn)品FNL(Final)改正InSAR大氣延遲的方法。本方法的思想就是:“以時間換空間”。

利用NCEP FNL對研究區(qū)域的平均風(fēng)速矢量進(jìn)行估計(jì)時，首先獲取UTC時間00h，06h，12h和18h上的平均風(fēng)速的u分量和v分量序列，然后將u分量和v分量序列插值到SAR影像過境時刻。在利用NCEP FNL資料獲取了研究區(qū)域內(nèi)的平均風(fēng)速后，即可利用“凝固流”假設(shè)，通過式(1)來獲取空間分辨率更高的濕延遲控制點(diǎn)。假定GPS天線的位置為(xg，yg)，則在時刻tni該GPS測站獲取的ZWD可以平流輸送至SAR影像獲取時刻ti的新位置(x，y)。式(1)可表示為:

4 結(jié)語

針對GPS改正InSAR大氣延遲時，地面GPS測站的低空間分辨率影響了其生成的大氣濕延遲分布的精度，從而影響了GPS對InSAR大氣的改正效果這一問題，本文利用NCEP全球操作分析最終產(chǎn)品FNL，通過時間和空間插值，獲取研究區(qū)域平均風(fēng)速大小和方向的方法，提高濕延遲控制點(diǎn)的空間分辨率和差分濕延遲分布圖的精度。改進(jìn)方法與在理論上證明了減少InSAR干涉圖中的殘余大氣相位的可行性，這為以后InSAR在地質(zhì)災(zāi)害區(qū)域進(jìn)行山體滑坡等的監(jiān)測起到了很重要的推進(jìn)作用。

［1］Hanssen，R.F.，Radar interferometry:data interpretation and error analysis［M］.xviii，308pp.，Kluwer Academic，Dordrecht;Boston，2001

［2］Massonnet，D.，K.Feigl，M.Rossi，and F.Adragna，Radar interferometric mapping of deformation in the year after the Landers earthquake［J］.Nature，369，227 ～230，1994

［3］Gao，B.C.;Kaufman，Y.J.Water vapor retrievals using Mod-erate Resolution Imaging Spectro－radiometer(MODIS)near－infrared channels［J］.J.Geophys.Res.2003，108，4389～4398

［4］Li，Z.H.;Muller，J.－P.;Cross，P.;Fielding，E.J.Interferometric synthetic aperture radar(InSAR)atmosphericcorrection:GPS，ModerateResolutionImagingSpectroradiometer(MODIS)，and InSAR integration［J］.J.Geophys.Res.2005，110，B03410

［5］Li，Z.W.，2005.Modeling atmospheric effects on repeatpass InSAR measurements［D］.PhD Dissertation.The Hong Kong Polytechnic University，Hong Kong

［6］Zhenhong LI.Correction of Atmospheric Water Vapour Effects on Repeat－Pass SAR Interferometry Using GPS，MODIS and MERIS Data［D］.University College London，2005.5

［7］Li，Z.H.;Muller，J.P.;Cross，P.Tropospheric correction techniques in repeat－pass SAR interferometry［C］.Proceedings of the FRINGE 2003 workshop，ESA ESRIN，F(xiàn)rascati，Italy，1 ～5，December 2003

［8］FerrettiA，PratiC and RoccaF.Perm anentscatterersin SAR interferometry［J］.IEEE Transactions on Geoscience and Rem ote Sensing，2000，38(5):2202 ～ 2212

［9］FerrettiA，PratiC and Rocca F.Perm anentscatterers inSAR interferom etry［J］.IEEE Transactions on Geoscience and Rem ote Sensing，2001，39(1):8 ～20

［10］Ferretti，A.;Prati，C.;Rocca，F(xiàn).Nonlinear Subsidence Rate Estimation Using Permanent Scatters in Differential SAR Interferometry［J］.IEEE T.Geosci.Remote Sens.2000，38，2202 ～2212

［11］Ferretti，A.;Prati，C.;Rocca，F(xiàn).Permanent Scatters in SAR interferometry［J］.IEEE T.Geosci.Remote Sens.2001，39，8 ～20

［12］Hooper，A.;Zebker，H.;Segall，P.;Kampes，B.A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers［J］.Geophys.Res.Lett.2004，31，L23611

［13］Janssen，V.;Ge，L.L.;Rizos，C.Tropospheric correction to SAR interferometry from GPS observations［J］.GPS Solut.2004，8，140 ～ 151.China，November.1997［C］.［s.l.］:［s.n.］，1997

［14］Li，Z.W.;Ding，X.L，Huang，C.;Wadge，G.;Zheng，D.W.Modeling of atmospheric effects on InSAR measurements by incorporating terrain elevation information［J］.J.Atmos.Terr.Phys.2006，66，1189 ～1194

［15］Kalnay E，M Kanamit su，R Kistler，et al.The NECP/NCAR 40 year reanalysis project［J］.Bull.Amer.Meteor.Soc.，1996，77(3)，437 ～470

［16］鄧偉，陳海波，馬振升等.NCEP FNL全球分析資料的解碼及其圖形顯示［J］.氣象與環(huán)境科學(xué)，2009，32(3):79～82

［17］Taylor，G.I.，1938.The Spectrum of Turbulence［J］.P.Roy.Soc.Lond.A.Mat.，164，476 ～490

［18］Emardson，T.R.andWebb，F(xiàn).H.Estimating the motion of atmospheric water vapor using the global positioning system.GPS Solutions［J］.2002，6(1－2):58 ～64

［19］Onn，F(xiàn).Modeling water vapor using GPS with application to mitigating InSAR atmospheric distortions［D］.Ph.D dissertation，Stanford University，2006;P.176