差分干涉技術在地觀測中的應用

1sar圖像的融合方法干預波束測量技術是近10個方面開發(fā)的一種新的空間觀測技術。它與gps、vlbi和slpm等空間技術有關。該項技術的發(fā)展和應用與合成孔徑雷達(SAR,SyntheticApertureRadar)傳感器技術的發(fā)展息息相關。合成孔徑雷達(SAR)概念的提出是相對真實孔徑雷達天線而提出的。對于真實孔徑雷達,當雷達隨載體(飛機或衛(wèi)星)飛行時,向地表發(fā)射雷達波束,然后接受地面反射信號,這樣便得到了地表雷達圖像。我們知道衛(wèi)星雷達天線越長,對地物的觀測分辨率就越高。由于受雷達天線長度的限制,真實孔徑雷達的地表分辨率往往很低,難以滿足應用要求。而合成孔徑雷達正是解決了利用有限的雷達天線長度來獲取高分辨率雷達圖像的問題。其基本原理為:衛(wèi)星在軌道飛行中,在不同位置上定時地向同一地物發(fā)射多個電磁波脈沖信號,然后接收同一地物的多個回波信號,再利用時間延遲和多普勒頻移效應在斜距向和方位向進行壓縮合成,便得到了SAR圖像。從某種意義上講,可以認為是延長了雷達天線長度,從而大大提高了分辨率(舒寧,1997)。SAR圖像具有許多光學圖像無法比擬的特點。微波遙感(或雷達遙感)主要采用主動遙感方式,即由遙感平臺發(fā)射電磁波,然后接收輻射和散射回波信號,具有全天侯、全天時作業(yè)的優(yōu)勢。雷達波長在1mm至1m之間,對地表具有一定穿透能力。除此之外,SAR數(shù)據(jù)同時記錄了地物強度(灰度)和相位兩種信息。正是SAR圖像中每個像素所包含的相位信息,為探測地形高度和地表動態(tài)變化提供了可能。早在20世紀60年代,INSAR技術就被人們注意到。1969年Rogers等人將干涉雷達技術用于探測金星和月球表面地形,1974年Graham又將干涉雷達技術用在地形制圖上,至到1986年Zebker等在飛機上安裝了相距11.1m的兩個SAR天線,獲取了10m分辨率的地形數(shù)據(jù),并對機載雷達干涉算法進行了改進,才使得數(shù)字化干涉雷達技術得以真正實現(xiàn)。目前該項技術已成為大范圍獲取數(shù)字地面高程模型(DEM)的重要技術之一。2000年2月22日,裝載于“奮進號”航天飛機的單軌雙天線雷達系統(tǒng)經(jīng)過近10d的作業(yè),成功地獲取了從北緯60°至南緯56°間80%陸地面積的高精度三維雷達數(shù)據(jù),生成了空間分辨率為30m(C波段垂直精度為10m,X波段垂直精度為6m)和90m的DEM產(chǎn)品。若采用傳統(tǒng)測繪作業(yè),其工作量將難以想象。這一成果標志著人類地形測繪與制圖技術由二維轉向三維的一次飛躍(郭華東等,2000)。目前機載干涉雷達測高精度為1～10m,星載干涉雷達測高精度為2～25m。在INSAR技術的基礎上,再進行差分算法處理,便形成了用于測量地表形變場的差分干涉測量技術(D-INSAR)。1993年在日本召開的現(xiàn)代地殼運動國際討論會(CRCM’93)上,法國學者Massonnet首次展示了利用INSAR技術獲取的1992年6月28日美國加洲Landers地震造成的近場形變場測量結果,引起了美國、加拿大、德國、日本等發(fā)達國家的高度重視。目前D-INSAR技術在地震、火山、活動構造等動態(tài)形變測量工作中得到了廣泛應用(Kurtetal.,1995;Meyeretal.,1996;Peltzeretal.,1999;Freysteinn,1999)。它與其它離散點測量技術相比,其測量結果具有連續(xù)的空間覆蓋優(yōu)勢,是一項極具發(fā)展?jié)摿Φ目臻g對地觀測新技術。2地表起伏引起的斜距向微量變化的特征目前星載INSAR處理中大多采用重軌方式,其原理是:利用衛(wèi)星在相鄰軌道上對同一地區(qū)進行兩次SAR成像為基本處理數(shù)據(jù)。單純對一幅SAR圖像來講,絕對相位信息沒有被利用,而兩幅獨立的同一場景SAR圖像之間的相對相位卻含有該場景中重要的三維信息。利用兩幅SAR圖像的相對相位進行干涉處理,獲取干涉圖像,通過相位解纏,便可從干涉圖像中的相位差獲取地形高程數(shù)據(jù)。下面以星載單天線雙成像雷達系統(tǒng)為例簡單地介紹INSAR技術原理,但其推導結果完全適用于機載雙天線雷達系統(tǒng)。圖1給出了INSAR處理中所涉及的幾何參數(shù)關系圖。其中,A1和A2為衛(wèi)星瞬間對同一地區(qū)成像的天線位置,衛(wèi)星到地物的斜距分別為r和r+δre+δrt,其中δre,δrt分別為橢球面和地形起伏引起的斜距向變化量。地物位于P處,高度為h。A1位置的雷達天線高度為H。B為兩天線的基線距,其中B‖是基線距的斜距向(視線方向)分量,B⊥是基線距垂直于斜距向分量。α為基線距與水平方向夾角。δθe,δθt分別為橢球面和地形起伏引起的視角增量。設在A1和A2兩次成像之間,P點沒有地表形變,則回波在A1,A2天線處的相位差為?=4πλ(δre+δrt)(1)?=4πλ(δre+δrt)(1)首先求出地形起伏所造成的斜距向微量變化,根據(jù)余弦定理,有(r+δre)2=r2+B2-2rBsin(θ-α)(2)(r+δre)2=r2+B2?2rBsin(θ?α)(2)考慮到B<<r,忽略二次小量(δre)2和B22rB22r,則地形起伏造成的斜距向微量變化為δrt≈-Bsin(θ-α)=-B∥(3)δrt≈?Bsin(θ?α)=?B∥(3)然后再求地形起伏造成的斜距向微量變化。取橢球面參考點P0,θ0為參考側視角,即圖像上高程值為0點的視角,則橢球面起伏造成的斜距向微量變化為δre=[Bsin(θ-α)-Bsin(θ0-α)]δre=[Bsin(θ?α)?Bsin(θ0?α)]考慮到θ0與θ之間的差異很小,對上式作數(shù)學近似得δre=δθeBcos(θ-α)=δθe?B⊥(4)δre=δθeBcos(θ?α)=δθe?B⊥(4)考慮到(3)和(4)式,則(1)式可改寫為?≈-4πλ(Bsin(θ-α)+δθeB⊥)=-4πλ(B∥+δθeB⊥)(5)?≈?4πλ(Bsin(θ?α)+δθeB⊥)=?4πλ(B∥+δθeB⊥)(5)從(5)式可以看出,干涉紋圖上的相位差包含了兩部分,一部分是橢球面起伏造成的相位差:-4πλ?δθeB⊥?4πλ?δθeB⊥;另一部分是地形起伏造成的相位差:-4πλB∥。在實際處理中,為了得到反映地形起伏的凈相位,確保計算精度,通常在相位解纏前將橢球曲面造成的相位差從干涉紋圖中除去,稱之為“去平地”(FlatEarthPhaseRemoval)方法。經(jīng)“去平地”處理后,干涉圖所表示的相位差只是地形起伏造成的。?=-4πλBsin(θ-α)=-4πλB∥(6)要獲取的P點地物高程h表示為h=Η-rcos(θ)(7)基線距B及傾角α、航高H、斜距r等參數(shù)都可從雷達系統(tǒng)參數(shù)或干涉數(shù)據(jù)中獲得。視角θ是隨地表起伏變化而變化的,取決于成像參數(shù)和圖像上每一點的地形。因此必須利用式(6)計算出每個像元的精確入射角θ,才能得到各個相元位置的高程值。從公式(6)中可知,如果從相干圖上測定了相位差?,易求出各個像元的入射角θ。但在相干圖像上,相干條紋的周期性變化代表的是相位差被2π除之后的模,取值范圍在-π到+π之間,是包纏相位(PhaseWrapped)。因此,必須利用相位解纏(PhaseUnwrapped)算法求出相位整周數(shù)(Pritt,1996)。目前常用的相位解纏算法有支切法(BranchCut)、循環(huán)圓盤掩模法(IterativeDiskMasking)和最小二乘法(LeastSquare)。其基本原理是:假設兩個相鄰像元相位之差不超過一個π。如果從圖像某一點出發(fā),兩相鄰像元相位差超過π時,則增減2π。設?i,?j分別為兩鄰近點的測定相位差,Δ?i,Δ?j分別為兩鄰近點的絕對相位差,則兩鄰近點的相位關系可表示為Δ?i=2πni+?i(8)Δ?j=2πnj+?j(9)dφ=?j-?i(10)當dφ>π,nj=ni+1;dφ<-π,nj=ni-1;否則nj=ni。但在實際SAR處理過程中,常常由于受疊置、陰影、無回波信號等因素影響,造成相干圖像上出現(xiàn)非相干熱噪點,不滿足相鄰像元相位小于π的假設。在這種情況下,進行相位解纏時,熱噪點的不連續(xù)性將產(chǎn)生誤差的傳播和累加,使相位解纏復雜化。因此,在實際處理中,通常對不連續(xù)熱噪點進行分片或掩蓋,然后再進行解纏。3地表微量變形場差分利用上述INSAR處理算法,我們可得到同一地區(qū)的兩幅干涉圖。其中一幅為事件(地震、滑坡等)前兩景SAR圖像,可通過干涉處理得到;另一幅是由事件前后兩景SAR干涉處理所得。在此基礎上,把兩幅干涉圖像進行差分干涉處理(DifferentialINSAR),便可得出該地區(qū)的斜距向微量地形變。如圖1所示,假設A1和A2′是第2次干涉處理時的衛(wèi)星成像天線,且在A2′成像前,地表P點出現(xiàn)了微量形變d,δθd為地表微量形變引起的側視角增量,δrd為地表微量形變引起的斜距向變化量。第2幅干涉紋圖的基線距為B′,基線距與水平方向的夾角為α′。假設該形變與雷達分辨單元尺寸相比很小,即認為A1,A2′的雷達信號仍是相關的。這時干涉相位差可表示為?′=4πλ(δre+δrt+δrd)=-4πλ(B∥′+δθt?B⊥′-δrd)(11)經(jīng)“去平地”處理后,干涉圖所表示的相位差包含地形起伏和地表微量形變信息,即?′=-4πλ(B∥′-δrd)(12)考慮到(6)式,不難得出??t′=B∥B∥′,其中?t′是第2次干涉中地形起伏所貢獻的相位差。由此,得出地表微量形變場差分方程?′-B∥′B∥?=4πλδrd(13)由上式求出斜距向形變量后,便可通過分解得到水平形變量和垂直形變量。在此例中,我們采用了同一地區(qū)的3景SAR圖像進行差分干涉處理,稱為三通(ThreePass)模式。第1幅干涉紋圖與第2幅干涉紋圖的主圖像均由A1位置獲取,這樣可保持r和θ不變,其余兩幅分別與A1位置獲取的SAR圖像比較即可,這樣兩幅干涉紋圖相位就可相互比較。若采用4景SAR圖像進行差分處理(FourPass),兩主圖像間需要精確配準。若只有2景SAR圖像,再加上該地區(qū)的精確DEM(TWOPass),也可進行兩通的差分干涉模式處理。4瑪尼地區(qū)三維地震資料及應用1997年11月8日,在西藏自治區(qū)那曲地區(qū)瑪尼鄉(xiāng)北約150km處發(fā)生了強烈地震。根據(jù)美國國家地震信息中心(NEIC)測定的這次地震震中位置為35.11°N,87.37°E,震級為MS7.9,震源深度35km。瑪尼地震發(fā)生在藏北高原北部全長270km的瑪爾蓋茶卡斷層附近(李建華,1998;馬瑾等,2000)。野外考察資料表明(徐錫偉等,2000),瑪尼地震的地表同震破裂帶全長約120km,以左旋剪切為主,最大左旋位移量達4.5m,其東西兩側水平位移約2～3m,最大地裂縫寬度可達300～400m,地表形變相當劇烈。根據(jù)許力生等(1999)的研究表明,瑪尼地震發(fā)震應力場的P軸(方位角29°,傾角7°)和T軸(方位角122°,傾角23°)均接近水平,斷層性質為左旋-逆走滑斷層,斷層面的傾角較陡(88°),走向為北東東-南西西(250°),破裂方式為由西向東?，斈岬貐^(qū)氣候干燥,植被稀少,地形起伏,且地表有強烈地震造成的同震位移,非常適合利用D-INSAR技術提取地表同震位移場。通過對資料實際調研,我們選購了瑪尼地區(qū)的歐洲遙感衛(wèi)星1號和2號(ERS-1/2)SAR圖像數(shù)據(jù)(參數(shù)見表1)。ERS-1雷達衛(wèi)星于1991年7月發(fā)射,其上攜帶了C波段VV極化雷達,波長為5.6cm。1995年4月,歐洲空間局(ESA)又發(fā)射了與ERS-1系統(tǒng)參數(shù)相同、軌道相同的ERS-2雷達衛(wèi)星,兩顆衛(wèi)星構成雙星追逐模式(tandemmode),對同一地區(qū)的最短覆蓋周期為1d,保證了兩景SAR圖像之間具有良好的相干性。本次實驗中選取了震前1996年5月20日和1996年5月21日兩景ERS-1/2SAR追逐模式數(shù)據(jù),用來獲取高精度地形信息。同時選取震后1998年4月21日ERS-2SAR數(shù)據(jù),用于獲取同震形變信息。圖2給出了瑪尼地震ERS-1SAR圖像(見書后彩頁)。5地震地表變形場干涉圖分析我們選定了三通模式進行差分干涉處理,干涉處理流程見圖3。(1)首先選取重疊公共區(qū)。數(shù)據(jù)選取時,都是按同一中心經(jīng)緯度選取的,但由于選取數(shù)據(jù)時截取方式不同,選取后的實際數(shù)據(jù)覆蓋范圍仍有差別,因此需要確定公共重疊區(qū),以保持每景數(shù)據(jù)處理范圍的一致性。(2)圖像配準。我們將2號圖像作為主圖像,分別與1號和3號圖像進行基于軌道參數(shù)的自動配準,配準精度達0.2像素。(3)求取相位差,生成干涉圖。把配準后的兩組圖像對進行共軛相乘,生成兩幅干涉圖。并采用“去平地”方法除去地球曲面影響。2-1干涉圖由震前-震前追逐模式SAR數(shù)據(jù)獲取(圖4,見書后彩頁),通過采用循環(huán)圓盤掩模法進行相位解纏,獲取震區(qū)高精度DEM(單新建等,2001)。2-3干涉圖由地震前后兩景SAR數(shù)據(jù)獲取,通過相位解纏,獲取包含地表形變信息的干涉圖(圖5,見書后彩頁)。(4)在2-3干涉圖中,除去2-1所獲的地形因素,得到反映同震位移地表形變場的差分干涉圖(圖6,見書后彩頁)。圖4為震前相差一天的兩景SAR的干涉結果,從圖4可以看出干涉條紋清晰,兩景圖像的相干性最好,它反映了研究區(qū)的地形信息。而震前-震后的兩景SAR圖像由于時間間隔較長,相干性較差,干涉圖出現(xiàn)熱噪點(圖5),該圖包含了地形信息、地表同震形變信息。圖6為從圖5中除去圖4的地形信息而剩余的形變信息,反映了地震造成的地表同震形變場的分布特征。從圖6可以看出,條紋沿北東向分布,F為發(fā)震斷層——瑪爾蓋茶卡斷層。圖中一個干涉條紋代表28mm的斜距向,設垂直于發(fā)震斷層F的極遠處B地表同震形變量為零,由B向B1累加形變量,即每遇一個干涉條紋增加28mm的形變量。因此,F處的形變量最大,為形變中心軸。6地表同震破裂帶的斷裂發(fā)育(1)由圖6可以看出,同震形變場沿NEE向分布,與瑪爾蓋茶卡斷層平行。通過量測,發(fā)震斷層走向為NEE-SWW(250°),這一結果與瑪尼地震的野外考察、干涉結果、震源機制解和余震分布相一致(許