利用雷達干涉時序分析方法研究地面沉降

吳文豪

(武漢大學衛(wèi)星定位系統(tǒng)工程技術研究中心，湖北 武漢 430079)

利用雷達干涉時序分析方法研究地面沉降

吳文豪

(武漢大學衛(wèi)星定位系統(tǒng)工程技術研究中心，湖北 武漢 430079)

天津地區(qū)因過度開采地下水導致地面沉降現(xiàn)象嚴重，對城市建設和發(fā)展產(chǎn)生了巨大的負面影響，因此科學分析其地面沉降機理及相關因素非常必要。本文以天津南郊為試驗區(qū)，基于短基線集干涉技術實施地面沉降監(jiān)測，并根據(jù)誤差源的空間相關性采用干涉圖的半變異值檢測誤差是否有效剔除，實現(xiàn)了短基線集處理結(jié)果的內(nèi)部驗證，減少了對水準結(jié)果檢校的依賴。結(jié)果表明，相對于比較穩(wěn)定的天津城區(qū)，天津南郊出現(xiàn)了大面積的沉降漏斗。最后結(jié)合水文地質(zhì)數(shù)據(jù)分析了地面沉降區(qū)域與地下水系關聯(lián)。

地面沉降；短基線；半變異函數(shù)

一、引 言

地面沉降是城市的主要地質(zhì)災害之一。隨著我國城市化進程的加快，地面沉降規(guī)模擴大，危害加劇。為有效控制地面沉降的速度，減小其危害性，沉降監(jiān)測是極為重要的一個環(huán)節(jié)，也是研究地面沉降機理與發(fā)展趨勢必不可少的基礎性工作。傳統(tǒng)的地面沉降測量方法包括水準測量、基巖標測量和分層標測量。隨著現(xiàn)代測繪技術的發(fā)展，布設GPS地面沉降觀測站也成為沉降監(jiān)測的數(shù)據(jù)來源，但這些技術都是獲取地面沉降的點信息。合成孔徑雷達干涉技術是近年來發(fā)展的新技術，它能夠在廣域范圍內(nèi)監(jiān)測地面沉降。

差分雷達干涉技術對微小的地形變化敏感度極高，可以監(jiān)測雷達視線方向厘米級的微小位移變化。然而進行長時間的地表微小形變監(jiān)測中，時間和幾何去相關，以及大氣擾動能因素限制了該技術的應用。為實現(xiàn)長時間和高精度監(jiān)測地表形變，許多學者提出了相應的方法來解決去相干問題。1999年Ferretti提出了一種“永久散射體”的時間序列處理方法[1]，隨后Berardino提出了短基線集方法，研究低分辨率、大尺度上的形變[2]。與永久散射體方法相比，短基線集方法降低了幾何和時間去相關效應。2003年Mora結(jié)合PS和短基線集方法的優(yōu)點，提出了基于高相干點的形變分析方法[3]，其優(yōu)點是采用少量的SAR圖像也可以得到毫米級的結(jié)果。

與傳統(tǒng)的沉降監(jiān)測相比，差分干涉雷達廣域監(jiān)測地表變化有著不可比擬的優(yōu)勢。伴隨著雷達干涉技術的發(fā)展，相關團隊利用該技術進行廣域范圍內(nèi)的城市沉降監(jiān)測研究。國外陸續(xù)開展了采用雷達干涉技術對熱點城市進行地面沉降的研究，這些城市包括洛杉磯、舊金山、拉斯維加斯、墨西哥的墨西哥城、挪威首都奧斯陸、意大利Naples和北部Bologna地區(qū)等。國內(nèi)諸多學者也利用該技術開展了地面沉降監(jiān)測實驗，主要集中在蘇州、上海、天津地區(qū)。如王超利用歐空局ERS1/2獲取蘇州市地面沉降測量值[4]；廖明生等探討了相干點目標分析技術原理和方法，進行了上海城區(qū)的地表沉降監(jiān)測實驗[5]；張詩玉、范景輝等分別采用不同的雷達干涉時序分析方法獲取天津地面沉降速率[6-7]；葛大慶以河北廊坊城區(qū)為實驗區(qū)，運用了短基線集方法獲取了形變時間序列和形變速率[8]；伍吉倉利用27景ENVISAT ASAR數(shù)據(jù)和18景ALOS PALSAR數(shù)據(jù)進行短基線集時間序列分析方法處理后得到上海磁懸浮列車線路周邊地區(qū)的地面沉降的年沉降速率都超過了20 mm/a[9]。無論在國內(nèi)還是在國外，雷達干涉時序分析方法已經(jīng)成為研究城市沉降的一個技術熱點[10]。

鑒于短基線集技術具有克服時間和空間失相關的優(yōu)勢，本文采用短基線方法獲取天津南部郊區(qū)沉降結(jié)果，并根據(jù)差分干涉影像的空間相關性采用半變異函數(shù)檢驗結(jié)果的正確性，最后根據(jù)天津水文地質(zhì)資料分析發(fā)現(xiàn)沉降漏斗空間分布與地下水系的分布存在相似性。

二、短基線集數(shù)據(jù)處理

1.短基線集理論

短基線集技術采用基線較短的干涉對組合進行干涉測量，不但繼承了常規(guī)差分干涉的優(yōu)點，并且避免了長空間基線的幾何去相干問題和長時間基線的時間去相干問題，能夠獲取長時間緩慢地表形變的演變規(guī)律，提高了形變監(jiān)測的時間分辨率。在短基線集中，高相干點像元依然是通過它們的相位屬性來提取的[11]。差分相位干涉圖中的兩個高相干點(p，q)的觀測模型可表示為

式中，ψx，i是第i副干涉圖中像元x的纏繞相位；N為干涉圖的數(shù)量；x，i是空間相關部分的估計值，如地形沉降信息，大氣擾動，軌道誤差；是空間非相關部分視角誤差，如地形誤差；γ即為像元在時間維上的相干系數(shù)，它的x值越大，越有可能成為高相干點。在數(shù)據(jù)處理過程中，大氣擾動、軌道誤差及地表沉降都被認為在小區(qū)域內(nèi)是強相關的。

2.數(shù)據(jù)處理

本文選取天津南郊作為試驗區(qū)(圖1中灰色框圖區(qū)域)，采用該區(qū)域降軌ASAR影像22景，時間跨度為2007年5月至2009年6月。

剪切后的數(shù)據(jù)基于短基線集理論進行處理，具體處理流程如下：

(1)差分干涉圖生成

1)配準：選取2008年3月14日的ASAR數(shù)據(jù)作為主圖像，其他影像與主影像進行配準，對于多普勒中心頻率差異超過100 Hz的影像對，配準前進行多普勒中心頻率濾波。配準使主影像和輔影像具有統(tǒng)一的幾何投影，是地面沉降干涉處理的關鍵環(huán)節(jié)，配準精度影響地面沉降的精度。

2)干涉對組合：一般認為時間間隔為35 d的干涉對沉降量很小，干涉圖中只含有噪聲，因此選取時間基線大于35 d且小于600 d，軌道基線小于450 m的影像對進行組合。干涉對組合如圖2所示。

圖1 研究區(qū)域與SAR影像覆蓋范圍

圖2 短基線干涉對組合

3)差分干涉處理：首先將干涉對影像共軛相乘獲得干涉相位，干涉相位由地形相位和形變相位組成；然后利用DEM數(shù)據(jù)去除地形相位信息，剩余的形變相位即為干涉條紋圖。差分干涉圖如圖3所示。

(2)高相干點及其時間序列的獲取

高相干點就是多景差分干涉圖中保持穩(wěn)定相干特征的目標點。通過監(jiān)測相干點的相位變化來分析地表沉降，因此相干點提取是差分干涉技術的前提條件。獲取高相干點的主要任務是完成部分誤差的分離，評估干涉圖相位在時間維的穩(wěn)定性。具體步驟如下：

1)相干系數(shù)計算：理論上高相干點是基于相位信息定義的。為了減少數(shù)據(jù)運算量，首先根據(jù)影像強度采用振幅離差法獲取高相干候選點。以某一個候選點為中心，指定半徑圓形區(qū)域內(nèi)使用自適應濾波器消除大氣擾動和軌道誤差，以及沉降信息，即x，i。地形誤差產(chǎn)生的相位與基線呈線性關系，可采用最小二乘法進行估計，據(jù)此可以將不同的相位成分進行分離，然后根據(jù)式(2)計算γx值?；谝堰x擇的高相干點和改正后的地形誤差，重復進行式(2)的計算，直至計算γx值的均方根變化很小。

圖3 短基線集差分干涉圖

2)高相干點提取：由于試驗數(shù)據(jù)量有限，時間維相干系數(shù)γx的直方圖分布與大樣本容量的直方圖相比會出現(xiàn)偏差。為了消除這種偏差，文中選用模擬隨機噪聲的方法修正試驗數(shù)據(jù)中的相干系數(shù)分布規(guī)律。然后設定γx大于0.7的像元為高相干點。

3)相位解纏：差分圖中的相位是纏繞的，為了獲取真實的相位信息，需要將纏繞的相位信息還原。相位解纏的基本假設是臨近相位差分絕對值小于π，但是非空間相關相位噪聲很可能導致相位差分出現(xiàn)絕對值大于π的情況，因此解纏前應盡量消除非空間相關誤差。高相干點在空間維是稀疏點，而開源解纏軟件snaphu只能對規(guī)則格網(wǎng)進行解纏。為解決這一問題，首先針對高相干點在空間維建立狄洛尼三角網(wǎng)，在時間維建立方格網(wǎng)；然后采用最鄰近插值法將三角網(wǎng)重采樣成為規(guī)則格網(wǎng)，這樣就將三角網(wǎng)解纏轉(zhuǎn)換為規(guī)則格網(wǎng)解纏；最后采用snaphu軟件完成短基線集的三維解纏。

4)大氣效應、軌道誤差的消除?？臻g非相關噪聲在沉降模型中可以直接視為隨機噪聲，但大氣效應和軌道誤差因其空間相關性會嚴重影響地面沉降計算結(jié)果。因此本文通過在時間維上的高通濾波和空間維上的低通濾波消除大氣效應和軌道誤差。與其他雷達干涉時間序列分析方法相比，該方法不需要提前預知地面沉降的時間序列模型，直接利用信號的空間相關性來估計像元相位穩(wěn)定性[12]，如圖4所示。

圖4 短基線集方法獲得的研究區(qū)域平均沉降速率圖

3.結(jié)果驗證

如式(1)所述，大氣效應和軌道誤差是雷達干涉的主要誤差源，具有較強的空間相關性。因此分析差分干涉圖的空間相關性可以判斷大氣擾動和軌道誤差是否被有效剔除。半變異函數(shù)作為表示區(qū)域化變量空間相關性的重要函數(shù)，可以定量化描述差分干涉圖相位的空間相關性[13－14]。因此本文采用半變異函數(shù)分析差分干涉圖在空間相關性誤差消除前后的差異。半變異函數(shù)假設某種屬性的的空間變異既不是完全隨機的，也不是完全確定的。一般認為空間變異可能包含3種影響因素：表征區(qū)域變量變異的空間相關因素；表征趨勢的“漂移”或結(jié)構，以及隨機誤差[15]。半變異的數(shù)學表達式為

式中，γ(h)為已知點xi和xj的半變異；h為這兩個點的空間距離；z為空間區(qū)域變量，在差分干涉圖中z為相位值。如果差分干涉相位不存在空間相關，那么近距離范圍內(nèi)點之間的半變異很??；相反，較遠距離的已知點的半變異變大?？紤]到大氣擾動的各向異性，相同的距離值可能出現(xiàn)多個半變異。因此在文中采用平均半變異和平均距離來繪制半變異圖。為了簡化計算半變異函數(shù)的數(shù)據(jù)運算量，又進一步對高相干點在空間維上作了降采樣處理[16]。

考慮到地面沉降有著很強的空間相關性，本文選取時間間隔較近的4幅干涉圖進行半變異分析，半變異圖如圖5所示。當差分干涉圖因軌道和大氣擾動的存在，解纏后相位值存在斜坡效應時，半變異表達式可變化為

式中，{E[z(xi)－z(xj)]}2稱為偏項，如果相位值不是平穩(wěn)隨機分布，區(qū)域變量存在趨勢項，則偏項值不為0，半變異值γ(h)呈現(xiàn)開口朝上的拋物線分布。因此圖中深色曲線則預示了差分干涉圖存在軌道誤差和大氣擾動，淺色曲線則表示偏項很小或者為0，大氣擾動和軌道誤差消除效果較好。所以利用半變異值可以對短基線干涉處理結(jié)果進行內(nèi)部檢驗，驗證誤差項是否有效剔除，彌補因缺少GPS數(shù)據(jù)或者水準數(shù)據(jù)無法進行結(jié)果驗證的不足。

圖5 距離分組后的半變異圖

三、結(jié)果分析

根據(jù)短基線集處理結(jié)果，天津東南郊區(qū)均出現(xiàn)明顯的沉降現(xiàn)象，其中大寺鎮(zhèn)和南河鎮(zhèn)出現(xiàn)兩個典型的沉降漏斗，并且與團泊洼水庫組成一個更大的沉降漏斗。天津地區(qū)在氣候上屬于半干旱地區(qū)，水資源一直是制約該地區(qū)經(jīng)濟發(fā)展的重要因素，過度抽取地下水已經(jīng)造成地下水位大幅度下降，出現(xiàn)大面積的地面沉降。近年來，盡管天津市區(qū)的地面沉降已經(jīng)得到了基本控制，但是隨著天津鄉(xiāng)鎮(zhèn)企業(yè)的迅速發(fā)展，郊縣和農(nóng)村地區(qū)的地下水開采量大幅度增長，天津市區(qū)的外圍郊區(qū)出現(xiàn)了新的地面沉降。根據(jù)天津水務局所提供的資料，天津地區(qū)主要由堆積平原區(qū)組成，第四系覆蓋層厚度巨大，含水層組分帶比較復雜，埋深400 m范圍自上而下可劃分為5個含水層組：從上到下依次為第Ⅰ、Ⅱ、Ⅲ、Ⅳ含水層組和咸水含水層組，受地質(zhì)條件影響，咸水含水層組一般位于第一、二層含水組中[17]。天津的主要采水層為第Ⅱ含水層組，以及第Ⅲ含水層組一部分，水位埋深在100 m左右，地表水很難進行有效補充。短基線集結(jié)果顯示即使在團泊洼水庫周邊地區(qū)同樣出現(xiàn)了大范圍的沉降，這也證明了天津地下水的抽取主要是深層地下水，僅僅依靠土壤滲透，即使地表水充足，地表補給量遠遠小于地下開采量導致地面沉降，地面同樣會出現(xiàn)沉降。

一般來說，地下水開采量、水位降深在時空分布上和地面沉降幅度上有明顯的一致性，即開采量和水位降深越大的地區(qū)和層位，地面沉降量越大，反之亦然。但是短基線集結(jié)果與水務局提供的水文資料對比表明津南區(qū)含水組漏斗與地面沉降區(qū)域并不重合。水文資料顯示第Ⅱ、Ⅲ含水組的漏斗中心在津南區(qū)的辛莊鎮(zhèn)、咸水沽鎮(zhèn)以及西青區(qū)的楊柳青和中北鎮(zhèn)，然而這些區(qū)域地面沉降并不嚴重。并且第Ⅱ、Ⅲ含水組水位埋深一致處于增加的狀態(tài)(見表1)，不存在地下水回灌的可能。我們推斷含水組漏斗和地面沉降區(qū)域不一致現(xiàn)象的出現(xiàn)與沉降區(qū)主抽水層的厚度、固結(jié)度和土體壓縮差異性有關。

天津地下水流動系統(tǒng)可分為多個地下水子系統(tǒng)[18]，圖4中黑色虛線為海河干流沖海積地下水系統(tǒng)子區(qū)與子牙河古河道帶地下水系統(tǒng)子區(qū)的分界線。本文研究的大寺鎮(zhèn)、南河鎮(zhèn)、團泊洼沉降區(qū)域以及整個東南郊區(qū)均屬于海河干流沖海積地下水系統(tǒng)子區(qū)，而左側(cè)的子牙河古河道帶地下水系統(tǒng)子區(qū)并未出現(xiàn)明顯的沉降。一般情況下，沉降區(qū)域應該是中心向外輻射分布，而該處沉降區(qū)域卻是沿著地下水子系統(tǒng)分布的。而圖中的沉降區(qū)域輪廓恰恰與地下水子系統(tǒng)的分界線是一致的，這說明地下水子系統(tǒng)之間的側(cè)向補給與各個子系統(tǒng)內(nèi)部補給還是有差異的，地面沉降區(qū)域往往與水文地質(zhì)的構造有關。

四、結(jié)束語

本文采用短基線集方法監(jiān)測天津南郊地表變化，結(jié)果表明包含大寺鎮(zhèn)、南河鎮(zhèn)及團泊洼庫區(qū)的城區(qū)東南郊出現(xiàn)范圍的沉降。為保證結(jié)果的準確性，基于差分干涉圖相位的空間相關性，對差分干涉圖進行半變異分析，驗證了大氣擾動和軌道誤差的有效剔除。但是在進行半變異參數(shù)的計算時，降采樣處理沒有考慮數(shù)據(jù)的空間相關性，在沉降區(qū)域沒有作掩膜處理，這有可能會影響到半變異函數(shù)的結(jié)果，需要作進一步的改進。最后結(jié)合水文地質(zhì)資料對津南地區(qū)的地面沉降作了分析，發(fā)現(xiàn)地面沉降區(qū)域輪廓與地下水系子系統(tǒng)分布是相似的，本文推斷這是地下水子系統(tǒng)之間的側(cè)向補給障礙造成的。這一結(jié)論對治理地面沉降，合理開采地下水有著指導意義。

表1 含水組漏斗中心水位埋深m

[1] FERRETTI A，PRATI C，ROCCA F.Permanent Scatterers in SAR Interferometry[J].IEEE Transactions on Geoscience and Remote Sensing，2001，39(1)：8-20.

[2] BERARDINO P，F(xiàn)ORNARO G.A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J].IEEE Transactions on Geoscience and Remote Sensing，2002，40：2375-2383.

[3] MORA O，MALLRQUI J J，BROQUETAS A.Linear and Nonlinear Terrain Deformation Maps from a Reduced Set of Interferometric SAR Images[J].IEEE Transactions on Geoscience and Remote Sensing，2003，40(10)：2243-2253.

[4] 王超，張紅，劉智，等.星載合成孔徑雷達干涉測量[M].北京：科學出版社，2002：198-205.

[5] 廖明生，盧麗君，王艷，等.基于點目標分析的InSAR技術檢測地表微小形變的研究[J].城市地質(zhì)，2006，1(2)：38-41.

[6] 張詩玉，李陶，夏耶.基于InSAR技術的城市地面沉降災害監(jiān)測研究[J].武漢大學學報：信息科學版，2008，33(8)：850-853.

[7] 范景輝，郭華東，郭小方.基于相干目標的干涉圖疊加方法監(jiān)測天津地區(qū)地面沉降[J].遙感學報，2008，12(1)：111-116.

[8] GE Daqing，WANG Yan，ZHANG Ling，et al.Mapping Urban Subsidence with TerraSAR-X Data by PSI Analysis[J].IGARSS，2010：3322-3325.

[9] WU Jicang，ZHANG Lina，LI Tao.Monitoring Ground Subsisdence in Shanghai Maglev Area Using PALSAR and ASAR Data[C]∥Proceeding of the Joint Interantional Symposium on Deformation Monitoring.Hongkong：[s. n.]，2011.

[10] 張永紅，吳宏安，孫廣通.時間序列InSAR技術中的形變模型研究[J].測繪學報，2012，41(6)：865-876.

[11] 張紅，王超，吳濤.基于相干目標的DInSAR方法研究[M].北京：科學出版社，2009：30-39.

[12] SOUSA J J，HOOPER A，HASSEN R F.Comparative Study of Two Different PS_INSAR Approaches：DePSI VS StaMPS[C]∥Proceedings of Fringe 2009 Workshop，F(xiàn)rascati，Italy：[s.n.]，2009.

[13] BONCORI J P M，MOHR J J.Statistical Description of Tropospheric Delay for InSAR：Overview and a New Model[J].IGARSS，2007：4483-4486.

[14] HASSEN R F.Radar Interferometry：Data Interpretation and Error Analysis[M].[S.l.]：Kluwer Academic Publishers，2001：154-159.

[15] 陳健飛，張筱林.地理信息系統(tǒng)導論[M].北京：科學出版社，2010：341-346.

[16] LOHMAN R B，SIMONS M.Some Thoughts on the Use of InSAR Data to Constrain Models of Surface Deformation：Noise Structure and Data Downsampling[J].Geochemistry Geophysics Geosystems，2005，6(1)：1-12.

[17] 袁宏利，董民，賈國臣，等.天津平原水利工程地質(zhì)環(huán)境概論[M].鄭州：黃河水利出版社，2008：158-161.

[18] 董克剛，王威，于強，等.土水比指標在天津沉降區(qū)地下水資源管理中的應用[J].水資源保護，2009，25(6)：51-55.

Research on Ground Subsidence by InSAR Time Series

WU Wenhao

P237

B

0494-0911(2014)11-0011-05

2013-11-18

國家自然科學基金(41274048)

吳文豪(1987—)，男，河南駐馬店人，博士生，主要從事雷達干涉測量的研究和應用。

吳文豪.利用雷達干涉時序分析方法研究地面沉降[J].測繪通報，2014(11)：11-15.

10.13474/j.cnki.11-2246.2014.0352