傅拓 ,張書畢
(1.中國礦業大學環境與測繪學院,江蘇徐州 221116;2.中國礦業大學江蘇省測繪與國土信息工程重點實驗室,江蘇徐州 221116)
InSAR技術是近年來迅速發展起來的極具應用價值的空間對地觀測新技術,具有監測精度高、范圍大、成本低、空間連續覆蓋等優點,為高邊坡地質災害區監測提供了一種新型的監測方法。但由于地質災害多發生在暴雨頻發、地質地貌復雜的區域,特殊的地理位置與氣候使得InSAR技術應用中受大氣延遲的影響非常嚴重,大氣延遲誤差是InSAR技術的主要誤差源,嚴重時還會導致InSAR結果的錯誤解譯[1],因此,必須要對InSAR觀測進行大氣延遲改正。
自1994年 Massonnet等[2]在利用 InSAR研究1992年發生在美國南加州地區Laders地震時首次發現了InSAR干涉圖中的大氣效應后,眾多學者對消除或減弱InSAR中的大氣延遲的方法展開了研究。近年來,InSAR大氣延遲改正的理論及方法取得了快速的發展,國內外提出了多種改進InSAR大氣延遲的方法。以下介紹幾種主流的方法:
MODIS(Moderate Resolution Imaging Spectroradiometer)[3]的近紅外水汽產品比世界上最密集的GPS網絡——美國南加州綜合GPS網(SCIGN)的密度還要高10倍以上。2005年國內許多學者提出融合MODIS和 GPS 數據的 InSAR 干涉圖改正方法[4-5]。但MODIS容易受到云的影響,且水汽產品存在系統誤差,需要進行改正。
MERIS(Medium Resolution Imaging Spectrometer Instrument)水汽產品與ASAR(Advanced Synthetic Aperture Radar)影像同時獲取,且得到的PWV空間分辨率高(300 m×300 m),精度也比MODIS 還要高[6~7],但其缺點是受云的影響更大。
常規DInSAR技術進行滑坡監測會受到相位失相關和大氣延遲影響,Ferretti[8~9]提出了僅僅跟蹤成像區域內雷達散射特性較為穩定的目標而放棄那些失相關嚴重的分辨單元的方法。這些目標(如建筑物的墻角、屋頂和裸露的巖石等)幾乎不受失相關噪聲影響,即使在多年時間間隔的干涉對中仍然保持較高的干涉相關性,把這些穩定的目標稱之為永久散射體(permanentscatterers,PS)。PSInSAR 技術[10~12]能克服傳統InSAR技術中的去相干性和大氣效應等困難,但利用PSInSAR技術進行InSAR大氣改正的缺點在于,PS點通常分布于市區等人工建筑較多的地方或無植被覆蓋的山峰、山脊等有裸露巨石的地方,對于研究山體滑坡等地質災害受到了很大的限制,要得到準確可靠的結果,對同一區域內的SAR影像數量也有較高要求,一般需要至少30幅影像。
GPS具有高精度、全天候、可連續估計大氣延遲的優點,GPS氣象學的發展為精確反演對流層延遲提供了保障,亦為估計與改正InSAR大氣效應提供了一種新的方法,因而利用GPS進行InSAR大氣延遲改正的方法得到了越來越多的應用。但SAR影像與GPS測站的空間分辨率差異很大,很大程度上限制了大氣延遲分布圖的估計精度,為了對SAR影像進行逐像素的大氣層延遲改正,必須對GPS獲取的大氣延遲進行加密插值,Janssen等[13]對反距離加權平均(IDW)、克里格(Kriging)、三次樣條3種不同的插值算法用InSAR大氣改正的效果進行了比較,結果表明反距離加權平均插值和克里格插值法比三次樣條插值算法效果更佳;李志偉[14]提出了一種新的大氣校正的方法,通過計算出平均大氣延遲,然后利用GPS探測到的大氣延遲對平均延遲進行修正,彌補了SAR干涉圖水平方向上各向異性的問題。但GPS改正InSAR大氣延遲仍然存在限制,當GPS測站附近的氣象觀測缺失時,大氣延遲估計精度會降低,GPS采集數據連續運行站點之間的間隔一般為幾十千米不等,測站點過于分散會降低觀測精度,增加觀測站的密度會受到地理環境和運作成本等因素的限制,即使設立臨時性、低成本的觀測點,效果仍然非常有限。針對以上幾種方法存在的問題,本文提出了一種新的改正InSAR大氣延遲的方法。
NCEP/NCAR再分析資料[15]是1991年開始實施的美國國家氣象中心(NMC)氣候資料同化系統(CDAS)的一個重要部分,其包含了豐富的氣象觀測資料,可用來精確估計氣象參數并提高大氣延遲分布圖的精度。該系統由NASA資助,以NCAR以前為航空與交通部門開發的用于美國大陸的產品為基礎,將衛星資料、模式結果與先進的人工智能技術結合起來,快速確定和預報風暴與潛在擾動氣流。
NCEP FNL資料[16]比再分析資料具有更高的時間、空間分辨率,其同化了全球搜集到的幾乎所有觀測信息,是全球業務模型生成的最后產品,可以認為是NCEP所有產品中用于模式存檔的最佳選擇。其數據水平格距為1°,水平方向共360×181個格點,垂直方向分為地表層和從1000~10hPa共26個等壓層,有海平面氣壓、表面溫度、海冰等9個單層變量和位勢高度、溫度、東西風分量、南北風分量、相對濕度共5個全層次變量。從1999年7月30日至今的NCEP FNL產品可以免費獲取,其空間分辨率為1°×1°,時間分辨率為 6h(http://dss.ucar.edu/datasets/ds083.2/)。NCEP FNL數據集以GRIB格式存儲,同樣可用GrADS軟件進行讀取。圖1為逐6h再分析資料集的分布圖。

圖1 逐6h再分析資料集的分布圖
“凝固流”假設由泰勒于1938年提出[17],該假設認為在平均風速的驅使下,空氣流在平移過程中處于“凝固”狀態,即在某一固定點上觀測到的時間上的延遲波動是由某一折射率場在平均速度為V的流的驅使下,經過該點時產生的空間波動引起的,“凝固流”假設的主要思想是空氣流在平均風速的驅使下進行平流輸送。
“凝固流”假設運用到空氣中的三維濕折射率Nwet(x,y,z,t)的數學描述可以表示為:(x,y,z)表示空間的三維坐標,t表示時間,u,v分別表示水平方向的緯向風和經向風。

“凝固流”假設研域平均風速U→—空間上為常數,但隨著時間變化,則指定空間位置(x,y,z)上的濕折射率可以認為是在之前時刻t,由另一空間位置(x-ut,yvt,z)上的濕折射率通過平均風速U→—=(u,v)的平流傳輸得到的。
為了利用GPS天頂濕延遲(ZWD)時間序列估計平均風速,假設研究區域內的GPS天頂濕延遲(ZWD)在空間上的統計特性為各向同性,且在時間上服從廣義穩態隨機過程[18],則兩個GPS測站上時間連續的ZWD間的互相關函數可表示為:

R=[RxRy]T表示 GPS 測站間的空間距離,ρ、ρ+R分別表示GPS測站的空間位置,△t表示ZWD獲取的時間間隔,t、t+△t分別表示ZWD的獲取時刻,V為平均風速。ZWD間的自相關函數可表示為:

連續時間上的二維空間ZWD場的互相關函數可以通過平均風速與其自相關函數進行轉換,將式(2)線性化有:

?表示Hamilton算子。
利用最小二乘原理,建立GPS ZWD序列的自相關函數模型與互相關之間的關系后,即根據式(4)解算出平均風速V。ZWD在連續觀測時刻的互相關函數可表示為:

式中,S(Rk)={(ρa,ρb):ρa-ρb≈Rk};k 表示 GPS測站個數,Rk表示 ZWD 間的空間距離;ρa、ρb分別表示時刻的 ZWD 位置。
利用式(4),平均風速可以用自相關函數進行估計:

式中,tm、ts分別表示主副SAR影像的獲取時刻。
由于地面GPS測站的分辨率低,利用其插值得到的濕延遲分布不能真實反映濕延遲空間分布狀況。基于泰勒提出的“凝固流”假設,Onn[19]將這種思想用于InSAR大氣改正,并指出這種方法比僅用GPS觀測時更為有效。其中,準確地估計研究區域的平均風速是至關重要的。平均風速的大小及方向決定了“擴展”的ZWD控制點的位置,決定了基于這些控制點插值得到的濕延遲分布圖的精度,從而決定了該方法對In-SAR大氣延遲效應的改正效果。
但是,將研究區域內ZWD的變化全部歸因于平均風速的變化是不嚴密的,GPS ZWD的變化是多方面因素的綜合結果,除了與風速有關外,還與研究區域的溫度等其他氣象因素有關。
針對上述問題,本文提出了融合GPS和NCEP全球分析最終產品FNL(Final)改正InSAR大氣延遲的方法。本方法的思想就是:“以時間換空間”。
利用NCEP FNL對研究區域的平均風速矢量進行估計時,首先獲取UTC時間00h,06h,12h和18h上的平均風速的u分量和v分量序列,然后將u分量和v分量序列插值到SAR影像過境時刻。在利用NCEP FNL資料獲取了研究區域內的平均風速后,即可利用“凝固流”假設,通過式(1)來獲取空間分辨率更高的濕延遲控制點。假定GPS天線的位置為(xg,yg),則在時刻tni該GPS測站獲取的ZWD可以平流輸送至SAR影像獲取時刻ti的新位置(x,y)。式(1)可表示為:

針對GPS改正InSAR大氣延遲時,地面GPS測站的低空間分辨率影響了其生成的大氣濕延遲分布的精度,從而影響了GPS對InSAR大氣的改正效果這一問題,本文利用NCEP全球操作分析最終產品FNL,通過時間和空間插值,獲取研究區域平均風速大小和方向的方法,提高濕延遲控制點的空間分辨率和差分濕延遲分布圖的精度。改進方法與在理論上證明了減少InSAR干涉圖中的殘余大氣相位的可行性,這為以后InSAR在地質災害區域進行山體滑坡等的監測起到了很重要的推進作用。
[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,Frascati,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.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.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].氣象與環境科學,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.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.Modeling water vapor using GPS with application to mitigating InSAR atmospheric distortions[D].Ph.D dissertation,Stanford University,2006;P.176