公用主影像干涉圖加權疊加方法及其在地面沉降監測中的應用

龍四春，張詩玉，馮 濤，李 黎

1.湖南科技大學 煤炭資源清潔利用與礦山環境保護湖南省重點實驗室，湖南 湘潭 411201；2.中南大學資源與安全工程學院，湖南 長沙 410083；3.商丘師范學院 環境與規劃學院，河南 商丘 476100

1 引 言

時空去相干和大氣折射延遲是限制DInSAR地面沉降監測精度的主要因素，為了削弱大氣擾動對差分干涉相位質量的影響，文獻［1］提出永久散射體雷達干涉測量技術（permanent scatterers InSAR，PS－InSAR），較好地解決了大氣擾動問題，隨后，該技術得到迅速發展，但由于要求數據量較多（通常20幅以上），且數據處理范圍通常要求在5km2以內［2－10］，因此，該技術目前很難被廣泛應用。文獻［11—19］采用基于相干點目標的DInSAR技術進行削弱大氣擾動和時空去相干的研究與試驗分析，得到的試驗區平均沉降圖能有效揭示地面形變的空間展布，得到較可靠的監測結果。可見，采用相干點目標分析方法能降低大氣延遲影響、提高信噪比，但這些研究沒有對相干點目標影像對的質量進行合理評價，如果其中部分影像質量太差，相干性太低，必將會影響到最后的沉降監測結果；文獻［20］利用6幅ENVISAT ASAR影像，應用干涉圖疊加法對珠江三角洲的地面沉降進行研究，發現因過去20年的城市化發展導致廣州、佛山、東莞等地的地面沉降現象非常明顯，但其研究是在非單一公共主影像的基礎上進行干涉圖疊加，會減弱干涉圖形變的線性關系，增大干涉圖疊加方法（Stacking）的假定條件要求（其一，單副雷達影像中大氣延遲相位為隨機分布，各差分干涉圖中大氣延遲相位影響也為隨機分布；其二，研究區的地表形變為線性特征）［11，17］，同時，沒有考慮干涉圖質量權重問題，影響了監測結果的精度。基于此，本文提出一種基于單一公用主影像且考慮干涉圖質量權重的干涉圖疊加方法，解決了以上技術存在問題，提高了監測精度，獲得較好的地面沉降監測結果。

2 公用主影像干涉圖加權疊加方法

2.1 基本原理及流程

公用主影像干涉圖加權疊加方法是選取最優公用主影像進行配準，將配準到同一主影像幾何空間的差分干涉解纏相位圖進行加權疊加，利用大氣延遲相位的隨機分布特性和地表形變信號近似線性表現的特點，提高疊加結果的信噪比。該方法假設在干涉圖中，大氣擾動的誤差相位是隨機的，而形變為線性形變［13，18，20］，并根據相干性來定權。基于這種假設，將多幅干涉圖對應的解纏相位疊加起來，所得的形變相位信息對應著累加時間基線內加權的變形量；疊加后的大氣延遲相位，卻不是單幅干涉圖中大氣相位誤差隨干涉圖數量倍數增長的結果，而只是干涉圖數量的平方根倍增長的結果。這樣，疊加相位圖中形變信息和大氣延遲項之間的信噪比就能得到提高。

公用主影像干涉圖疊加方法數據處理的基本流程包括：主影像的優化選取［21－23］、輻射定標等數據預處理；影像配準與重采樣到主影像空間；高相干區域的選取；干涉圖生成；去平地、地形生成差分干涉圖；相位解纏；根據相干圖質量定權，干涉圖解纏相位疊加；地理編碼、形變圖生成等。具體流程見圖1。

從圖1可以看出，公用主影像干涉圖加權疊加方法基本包含了常規DInSAR的所有關鍵步驟。但要進行有效的干涉圖加權疊加，首先，基于同一公用主影像的干涉圖疊加方法，應進行公用主影像的優化選取，選擇大氣擾動等誤差影響最小的影像作為公用主影像，以便最大限度減少大氣延遲和其他噪聲的影響。對于公用主影像的優化選取，既要使時間基線最佳（盡量位于時間跨度中間），又要顧及干涉對的有效空間基線和Doppler質心頻率、大氣延遲影響因子以及季節因素等對干涉圖相關性的顯著影響，尋找綜合影響最小的影像為公用主影像［21－23］。其次，要根據干涉圖相干性質量（本文根據各影像相干點的數量成正比來定權，數量最多的相干圖疊加權重qji設定為1）進行定權，如果存在有與其他干涉圖相差明顯較大的影像（即相干圖相干點的數量≤1／2最高相干圖相干點數時），將其權重設為0，不讓它參加干涉圖疊加計算與后續的數據處理。

圖1 公用主影像干涉圖加權疊加方法處理流程Fig.1 Flow of weighted stacking based on common master image

2.2 干涉圖加權疊加方法的函數模型

式中，N為干涉圖的數目；與分別為第i個像元第j幅差分干涉圖的解纏形變相位與其疊加權重；ΔTj為第j幅差分干涉圖主輔影像之間的時間間隔（時間基線）；λ為雷達波的波長。

由于單幅雷達影像中的大氣延遲相位分布是隨機的，各差分干涉圖中的大氣延遲相位影響也隨機分布，根據誤差傳播定律，則所有干涉圖疊加后，高相干點受到的大氣延遲影響可表示為

式中，為第j幅差分干涉像對對應的大氣相位延遲；ΔΦi為i像元在N個干涉圖疊加后的大氣延遲影響總和。

則干涉圖疊加后，高相干點大氣延遲對線性形變速率的影響可表示為

從式（4）可以看出，疊加后平均形變速率Vi中，形變速率的信噪比得到了提高。因此，干涉圖疊加像元中平均形變速率的標準偏差為

干涉圖加權疊加最有效的方法是對相干性較高的點相位進行疊加，在選擇高相干像元時，首先計算各雷達影像對應的相干系數，然后可用函數模型（6）［4，13］進行高相干目標點的提取

式中，m、n表示選擇的移動窗口方位向、距離向的像元數為某像元i在移動窗口中相干系數的均值，是通過包含該像元的移動窗口的所有像素復數信息來進行估算的，當高于給定閾值γc的點選定為干涉圖相位解纏的高相干目標點，將所有滿足該模型的點組合成高質量相干點集。

3 公用主影像加權疊加方法在沉降監測中的應用

3.1 研究區域與數據來源

自20世紀初以來，天津市區及近郊一直存在不同程度的沉降，某些地段沉降量已經超過3m，但1986年開始對過度開采地下水采取控沉措施后，市區地段沉降速度明顯減小，但市區以外的近郊及縣區，某些地段沉降不但沒有減弱，相反沉降量增大，沉降速率有進一步加大的趨勢［24］。

考慮到公用主影像干涉圖加權疊加方法的優越性，選擇天津市主城區及近郊作為試驗區，其中心經緯度為（39.112 88°N，117.069 63°E），面積約24km×28km，利用2003—2006年間獲取的11景ENVISAT ASAR數據（Track2447，Frame22817）進行試驗。

首先從這11幅SAR影像中擇優選擇1幅SAR影像作為公用主影像，可利用綜合函數模型［3－4］式（7）來進行確定

在雷達差分干涉測量中，通常利用美國航天飛機測圖任務SRTM 3″分辨率（90m）DEM來去除地形相位。天津地處華北平原上，地勢比較平坦，地形起伏在10m以內，具體見圖2，圖中白色方框為所選雷達影像覆蓋范圍。該試驗區地形的相位貢獻不大［11］。

表1 ASAR影像數據相關參數表Tab.1 Parameters list of ASAR image data

3.2 差分干涉測量與加權疊架

根據DInSAR的基本原理，以14529作為干涉公用主影像，將其他10幅輔影像與14529影像進行粗配準和精配準，并重采樣到主影像14529幾何空間。根據式（6）高相干點的選取模型，選定0.2作為相干系數閾值，則可得各相干圖見圖3。從圖3 中可以看出，14529－10521、14529－12024、14529－13026、 14529－17034、 14529－18537 和14529－21042相干性較好，相干點多且相干系數高，相干性基本一致，它們對應的有效空間基線較短，幾何去相干較小；而14529－16032相干性最差，它對應的有效空間基線最大，為－678.37m，是其他有效空間基線的近2倍以上，造成的幾何去相干太嚴重，高相干點不到其他相干圖的一半。如果將其直接與其他相干性好的相干圖進行疊加解算，必將會影響最終的監測精度。因此，合理設置相位疊加時各影像的權重尤為重要，權重的設置，可根據相干圖3中各影像相干點的數量成正比來定權，數量最多的相干圖疊加權重設定為1，當相干點的數量明顯低于其他相干圖時（經試驗證明，相干圖相干點的數量≤1／2最高相干圖相干點數，即≤1／2時），需將此相干圖疊加權重定權為零，可見14529－16032的疊加權重應設為0，否則會降低最終沉降監測結果。

再將配準的干涉圖與SRTM3DEM進行差分干涉，得到對應的差分干涉圖，見圖4。

從圖4中可以明顯地看出，14529－16032的差分干涉圖具有明顯的誤差。主要是由于其干涉圖相干點分布稀少，解纏相位不連續，很難形成正確的干涉影像。

高相干目標在空間分布上是離散的、間斷的，采用Delaunay法則能將所有離散相干點用若干個沒有重疊的三角形連接起來［20］，保證每一個點至少是一個三角形的頂點，假定每一個三角形的重心是一個相位解纏的節點，沿該三角形3條邊做3個頂點的相位閉合線積分，如果任意兩頂點間相位差的絕對值不大于π，則其積分值為0，若積分值為－2π或＋2π，則表示該節點有一個負留數或正留數，計算每一個三角形對應的留數，為了防止解纏誤差擴散，將正負留數用弧段成對地連接起來，在解纏時，避免積分路徑穿過這些連有正負留數的弧段［18］，再利用最小費用流法［6］來進行相位解纏，可得各解纏相位，見圖5。根據前面定權的理論，計算各影像的權重，其中14529－16032干涉對解纏相位的權重為0。通過以上10幅相位解纏圖疊加解算后，得到的年平均沉降速率見圖6，從圖6中可以看出，天津市中心城區的沉降趨勢基本一致，大部分沉降速率在2cm／a左右，但北辰區和西南角存在明顯的沉降中心，其最大沉降速率達8cm／a，可見天津郊區沉降速率要比市中心城區大，其主要原因是2003—2006年間天津工業向郊區發展和轉移，導致郊區地下水的抽取嚴重，地面下沉加速，與官方監測［24］（水準測量結果）一致。

根據式（5），對解算的平均形變速率計算標準偏差，可得如下形變速率標準偏差，見圖7。由圖7可知，市區線性形變標準偏差較小且比較均勻，最大不超過0.5cm／a，表明該地區在2003—2006年內地表沉降速率的線性特征明顯，監測結果可靠。

3.3 與高精度水準數據的對比

天津沉降辦公室提供了47個水準監測點2003年、2004年和2005年數據［24］，對獲得的這些離散監測水準數據，求取2003—2006年度內平均沉降速率，基于 Matlab軟件，采用Kriging方法進行內插計算，其內插結果見圖8。

從圖8中，可以看出北辰區的水準監測內插結果與雷達干涉圖加權疊加監測結果一致，其最大沉降速率都為8cm／a，沉降漏斗明顯。表明公用主影像干涉圖加權疊加方法與高精度水準監測精度相當。

3.4 與DInSAR及等權Stacking方法的對比

根據DInSAR的基本原理與流程，對以上11景SAR影像數據進行優化組合試驗，得到該地區最優監測結果見圖9。與公用主影像加權疊加方法圖6相比，DInSAR處理結果明顯存在噪聲的影響，且與水準監測數據最大沉降速率8cm／a相比，存在2cm／a的差異。同時，進行了隨機主影像等權Stacking方法，得到了最后形變監測結果如圖10所示，與本文公用主影像干涉圖加權疊加方法監測結果圖6相比，其連續性較差且偏差稍大，且其沉降漏斗北辰區宜興埠的最大形變監測值為7cm／a，與水準監測結果和新加權Stacking方法的監測結果有1cm／a的差異。基于以上試驗分析，初步推論公用主影像加權疊加方法具有更優的地表形變監測能力。

圖2 天津地表高程變化圖Fig.2 Tendency chart of topographic elevation in Tianjin area

圖3 相干圖集Fig.3 Coherent images

圖4 差分干涉圖集Fig.4 Differential interferograms

圖5 解纏相位圖集Fig.5 Unwrapping phase diagrams

圖6 平均形變速率Fig.6 Mean rate of deformation

圖7 形變速率標準差Fig.7 Standard deviation of deformation rate

圖8 平均沉降速率Kriging內插結果Fig.8 Mean rate of subsidence by Kriging interpolation

圖9 DInSAR監測結果Fig.9 Monitoring results of conventional DInSAR

圖10 傳統Stacking方法監測結果Fig.10 Monitoring results of conventional Stacking

4 結 論

在現有Stacking方法的基礎上，提出一種削弱大氣延遲影響的公用主影像干涉圖加權疊加方法。該方法主要貢獻：① 考慮時間基線、空間垂直基線、多普勒質心頻率差異的綜合影響，選取最優公用主影像，消弱主影像大氣延遲影響；② 顧及相干圖質量（即高相干點數量）進行定權，解決低相干影像圖對監測結果的負面影響，通過加權疊加后，相位圖的形變信息和大氣噪聲之間的信噪比得到明顯提高。該方法要求干涉影像數不少于3幅，且研究區域的形變為緩慢形變，它能在數據量較少的情況下（與PS方法相比），得到較好的監測結果。通過對天津市區2003—2006年11景ENVISAT ASAR影像進行試驗，經比較分析，監測結果優于傳統DInSAR和Stacking方法，其精度和高精度水準測量相當。可見，公用主影像干涉圖加權疊加方法解決了PS－InSAR數據量要求較多的問題，削除了傳統Stacking方法中低相干影像對監測結果的影響，提高了信噪比。

［1］ FERRETTI A，PRATI C，ROEEA F.Permanent Scatterers InSAR Interefromety［J］.IEEE Transactions on Geoscience and Remote Sensing，2001，39（l）：8－19.

［2］ COLESANTI C，FERRETTI A.SAR Monitoring of Progressive and Seasonal Ground Deformation Using the Permanent Scatterers Technique［J］.IEEE Transactions on Geoscience and Remote Sensing，2003，41（7）：1685－1701.

［3］ KAMPES B M.Displacement Parameter Estimation Using Permanent Scatterer Interferometry ［D］.Delft：Delft University of Technology，2005.

［4］ HOOPER.Persistent Scatterer Radar Interferometry for Crustal Deformation Studies and Modeling of Volcanic Deformation［D］.Stanford：Stanford University，2006.

［5］ PERISSIN D，ROCCA F.High－Accuracy Urban DEM Using Permanent Scatterers［J］.IEEE Transactions on Geoscience and Remote Sensing，2006：44（11）：3338－3347.

［6］ KAMPES B M.Radar Interferometry Persistent Scatterer Technique［M］.Berlin：Springer，2006.

［7］ FERRETTI A，SAVIO.Submillimeter Accuracy of InSAR Time Series：Experimental Validation［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45（5）：1142－1153.

［8］ JUNG H，KIMB S，JUNGA H，et al.Satellite Observation of Coal Mining Subsidence by Persistent Scatterer Analysis［J］.Engineering Geology，2007，92：1－13.

［9］ BELL J，AMELUNG F，FERRETTI A，et al.Permanent Scatterer InSAR Reveals Seasonal and Long－term Aquifersystem Response to Ground Water Pumping and Artificial Recharge［J］.Water Resources Research，2008，44（2）：402－407.

［10］ LONG Sichun，LI Tao，LIU Jingnan.Application Study of PS－DInSAR Technique Fusing Multi－metadata in Urban Ground Deformation Survey［C］∥Proceedings of International Conference on Environmental Science and Information Application Technology.［S.l.］：ESIAT，2009：326－330.

［11］ LONG Sichun.Advanced DInSAR Techniques Based on Coherent Targets and Its Application in Ground Subsidence Monitoring［D］.Wuhan：Wuhan University，2009.（龍四春.基于相干目標的DInSAR高級技術及其在地面沉降監測中的應用［D］.武漢：武漢大學，2009.）

［12］ STROZZI T，WEGMüLLER U，WERNER C，et al.Measurement of Slow Uniform Surface Displacement with mm／year Accuracy ［C］∥ Proceedings of IGARSS.Hawaii：IEEE，2000.

［13］ WERNER C，WEGMULLER U，WIESMANN A，et al.Interferometric Point Target Analysis with JERS－1L－band SAR Data［C］∥Proceedings of IGARSS.Muri：IEEE，2003.

［14］ BERARDINO P，CASU F，FOMARO G，et al.Small Baseline DInSAR Techniques for Earth Surface Deformation Analysis［C］∥Proceedings of FRING 2003Workshop.Franscati：［s.n.］，2003.

［15］ BERARDINO P，CASU F.A Quantitative Analysis of the SBAS Algorithm Performance.Geoscience and Remote Sensing Symposium ［C］∥ Proceedings of IGARSS.［S.l.］：IEEE，2004.

［16］ WEGMüLLER U.GAMMA IPTA Processing Example Luxemburg，GAMMA Technical Report［R］.Luxeburg：GAMMA，2005.

［17］ CASU F，MANZO M，LANARI R.A Quantitative Assessment of the SBAS Algorithm Performance for Surface Deformation Retrieval from DInSAR Data［J］.Remote Sensing of Environment，2006，102：195－210.

［18］ GE Daqing，WANG Yan，GUO Xiaofang，et al.Surface Deformation Monitoring with Multi－Baseline D－InSAR Based on Coherent Point Target［J］.Journal of Remote Sensing，2007，11（4）：574－580.（葛大慶，王艷，郭小方，等.基于相干點目標的多基線D－InSAR技術與地表形變監測［J］.遙感學報，2007，11（4）：574－580.）

［19］ ZHANG Yonghong，ZHANG Jixian，GONG Wenyu.Monito ring Urban Subsidence Based on SAR Interferometric Point Target Analysis［J］.Acta Geodaetica et Cartographica Sinica，2009，38（6）：482－493.（張永紅，張繼賢，龔文瑜.基于SAR干涉點目標分析技術的城市地表形變監測［J］.測繪學報，2009，38（6）：482－493.）

［20］ ZHAO Qing，LIN Hui，JIANG Liming.Ground Deformation Monitoring in Pearl River Delta Region with Stacking DInSAR Technique［J］.Geoinformatics and Joint Conference on GIS and Built Enviroment，2008，7145：1－9.

［21］ LONG Sichun，LIU Jingnan，LI Tao，et al.Method for Optimum Selection of Common Master Acquisition for PS－DInSAR Fusing GPS Data［J］.Journal of Tongji University（Natural Science），2010，38（3）：453－458.（龍四春，劉經南，李陶，等.融合GPS數據的PS－DInSAR公用主影像的優化選取［J］.同濟大學學報：自然科學版，2010，38（3）：453－458.）

［22］ CHEN Qiang，DING Xiaoli，LIU Guoxiang.Optimum Selection of Common Master Acquisition for PS－DInSAR［J］.Acta Geodaetica et Cartographica Sinica，2007，36（4）：395－399.（陳強，丁曉利，劉國祥.PS－DInSAR公共主影像的優化選取［J］.測繪學報，2007，36（4）：395－399.）

［23］ ZHANG Hua，ZENG Qiming，LIU Yihua，et al.The Optimum Selection of Common Master Image for Series of Differential SAR Processing to Estimate Long and Slow Ground Deformation ［C］∥Proceedings of IGARSS.Seoul：IEEE，2005：4586－4589.

［24］ Tianjin Office of Ground Subsidence Monitoring.Annual Report on Tianjin Ground Subsidence［R］.Tianjin：Tianjin Water Resources Bureau，2006.（天津市控制地面沉降工作辦公室.天津市地面沉降年報［R］.天津：天津市水利局，2006.）

［25］ YU Xiaoxin，YANG Honglei，PENG Junhuan.A Modified Goldstein Algorithm for InSAR Interferogram Filtering［J］.Geomatics and Information Science of Wuhan University，2011，36（9）：1051－1054.（于曉歆，楊紅磊，彭軍還.一種改進的Goldstein InSAR干涉圖濾波算法［J］.武漢大學學報：信息科學版，2011，36（9）：1051－1054.）

［26］ XIA Y，KAUFMANN H，GUO X F.Differential SAR Interferometry Using Corner Reflectors［C］∥Proceedings of IGARSS.Toronto：IEEE，2002：1243－1246.

［27］ XIA Y，KAUFMANN H，GUO X F.Landslide Monitoring in the Three Georges Area Using D－InSAR and Corner Reflectors［J］.Photogrammetric Engineering and Remote Sensing，2004，7（10）：1167－1172.

［28］ XIA Y，Bam Earthquake：Surface Deformation Measurement Using Radar Interferometry［J］.Acta Seismologica Sinica，2005，18（4）：451－459.