地面沉降－回彈及地下水位波動的InSAR長時序監測——以德州市為例

葛大慶，殷躍平，王 艷，張 玲，郭小方，王 毅

(1.中國地質大學(北京)水資源與環境學院，北京 100083;2.中國國土資源航空物探遙感中心，北京 100083;3.中國地質環境監測院，北京 100081)

0 引言

長期過量開采地下水，特別是深層地下水，是目前我國一些平原和三角洲地區產生地面沉降的主要原因[1]。地下水位的持續下降改變了原有的補徑排特征，形成了區域性地下水位降落漏斗，地面沉降隨之產生并擴大。我國已形成了以華北平原、長江三角洲及汾渭斷陷盆地為主的三大地面沉降嚴重區，華北平原已出現近20個大型復合水位降落漏斗，地面沉降呈現出區域性且連片分布的特征[2－3]。

雷達干涉測量(InSAR)技術為區域性地面沉降的快速準確和連續監測提供了重要手段。以永久性散射體干涉測量(permanent scatterer InSAR，PSIn-SAR)[4]為代表的相干目標(coherent target，CT)時間序列分析技術克服了差分InSAR技術所面臨的失相干和大氣相位延遲等問題，通過對相干目標的差分干涉相位序列進行時序分析，逐個估計和分離地形相位、形變相位、大氣相位和噪聲等信息，進而解算相干目標的形變速率和累積形變量。與之類似的還有短基線集(small baseline subset，SBAS)[5]以及干涉測量點目標分析(interferometric point target analysis，IPTA)[6]等技術，這些技術已在地面沉降、滑坡、開采沉陷、地震和冰川滑移等多尺度緩慢地表形變研究中得以發展和應用[7－9]。

德州市位于山東省北部，是華北平原地面沉降和深層地下水位降落漏斗極為典型的地區之一。其地下水位降落漏斗在1978年前已經形成，是目前華北平原水位降落深度最大的漏斗，與衡水、滄州等漏斗連成一片。該區的地面沉降在1989年以前就已經出現，至2010年累積沉降量超過1 m。與地面沉降對應的是，德州市地下水位呈現出多年連續下降和年內季節性波動的特征[10]。由于地面沉降與地下水位下降相關，因此在德州市深、淺層地下水位均呈現出季節性波動的情況下，地面沉降場的時空演變特征值得研究，這包括:沉降場的空間變化特征，即沉降中心的移動及其變化幅度;沉降場的時間變化特征，即季節性波動，包括是否會出現地面回彈，回彈幅度多大，地面沉降與地下水位波動是否同步等方面。針對上述問題，本文以相干目標InSAR時序分析方法為主，利用2004年1月—2010年10月近7 a的ENVISAT衛星數據研究分析德州地區地面沉降－回彈的動態變化特征及其與地下水位波動的對應關系，以揭示該地區地面沉降季節性變化的產生原因和誘發機制。

1 相干目標InSAR時序分析方法

InSAR時序分析主要有PSInSAR和SBAS兩類方法。根據干涉像對的組合模式，前者以單一主影像構成差分干涉圖序列，后者則利用短基線原則以多個主影像構成若干個干涉紋圖子集。單一主影像的干涉紋圖序列處理主要顧及大氣相位和非線性形變的統計模型，有利于對大氣相位的準確估計，但對數據量要求過高，且受相干性影響而對大變形的監測能力不足;多個主影像則降低了對數據量的要求，將時空基線小于一定閾值的干涉像對組合，生成干涉紋圖序列，增加對同一信號的采樣密度，有利于準確求解形變速率，在數據量較小的情況下也可以實現。在選擇相干目標的方式上，前者以點目標為主，后者兼顧點目標和分布式目標[5]。本文集成了2種方法的優點，以短基線為準則構成差分干涉相位圖，利用點目標識別算法提取相干目標。以Delanay三角網連接相鄰的相干像元，用二維周期圖估計點間形變速率和高程誤差改正。以此為基礎，根據大氣相位特征，利用奇異值分解算法解算空間濾波后的殘余相位，并對解算結果進行時域濾波，以求解非線性形變，得到不同時刻的累積變形量，再利用線性回歸法求解每個相干目標的變形速率。

1.1 干涉相位時序分析的基本模型

對于給定的M景雷達數據，在短基線條件下生成N個差分干涉圖(一般情況下，N＞M)，對于其中的任一差分干涉圖k(k=1，2，…，N)，其相位模型為

顯然，待求解量d與InSAR參數不直接相關，但受其他分量的影響。對于大氣分量，其本質上是2次大氣延遲影響的疊加，在空間上具有一定的相關性[4]。大氣波動在時間上的隨機特性決定了2次互差也為隨機量。需要重點區分的是大氣相位與軌道誤差，雖然二者都具有空間低頻的特性，但各自的影響范圍有所不同?；€誤差多是全局性的，而大氣影響半徑有限，對于大范圍處理，二者應區別對待。高程誤差與垂直基線相關，垂直基線是變量。變形量與干涉圖的時間基線相關，時間間隔是變量，相對變形的速率與時間相關。

1.2 地表形變參數的提取算法

1.2.1 線性分量的估計

考慮到大氣的空間相關性，如果對相鄰2點(小于大氣相關距離)求差，則可削弱大氣影響。相鄰相干目標i和j差分干涉相位的互差△為

由于相鄰目標間的相對變形速率與時間間隔相關，干涉圖的時間基線為自變量，因而可將形變量d分解為線性形變速率和非線性形變量2部分。如此，時序分析過程則轉為參數估計，即

式中:CB為與垂直基線相關的系數;T為時間基線;△ε為相對高程誤差;△υ為相對形變速率;μNL為非線性形變量;α為大氣相位;n為噪聲。

對于纏繞相位，其周期數由相對高程誤差和線性形變速率共同決定，采用二維參數估計方法可實現對相位周期的估算。由此，構建目標函數

并將其從相位互差中減去，得到殘余相位，即

顯然，在準確估計目標函數的參數△ε和△υ時，殘余相位將最小化。

由于相鄰2點存在多個干涉圖，即存在N維B和T，因而，對△ε和△v的估計使模型相關系數最大化，即

在纏繞相位條件下，△ε和△v為周期函數的二維頻率，可利用二維周期圖估計使模型相關系數最大化;在解纏相位條件下，則轉換為二維線性函數，利用參數估計方法可實現對△ε和△v的求解。

在相鄰2點解算時需對所有相干目標建立連接關系，整體求解形變場的速率和高程誤差改正量。這一目標的實現可通過建立Delanay三角網的方法完成，也可以利用冗余網構建更為復雜的連接關系，強化對待解算方程組的約束。利用鄰近法則將所有距離滿足大氣相關距離的相干目標連接起來，在求解完成相鄰點間的互差后，通過最小二乘或加權平均的方法求解每個目標相對于參考點的變形速率，即總體形變速率場。

1.2.2 非線性分量的估計

對于具有顯著非線性形變過程，仍需對殘余相位進行更為復雜的處理，以提取非線性形變量。從差分相位中去除高程誤差估計值和線性速率后，殘余相位為

由于已解算出整周相位，因而殘余相位為相位主值，其大小

這里，首先要實現對單個差分干涉圖中殘余相位的濾波處理。由于大氣表現出空間低頻特性，而非線性變形的空間變化范圍較小，相對大氣而言表現為高通特性。因而，對殘余相位圖進行空間低通濾波可以進一步弱化大氣影響。在短基線集條件下求解的非線性形變量，不同于PSInSAR經典處理策略，不直接對大氣進行估計，而是通過濾波減弱大氣的影響。殘余相位經空間高通濾波后，大氣分量已經減弱。此時，進一步的處理是求解與雷達數據對應的不同時刻的非線性變形量，包含噪聲等影響。根據主輔影像的關系，可將解纏后的殘余相位分解為

式中:p和q分別表示生成第k(k=1，…，N)景差分干涉圖的主輔影像的獲取時間。由于采用短基線原則，在求解每個時刻對應的殘余相位時出現秩虧方程組，解決這一問題的辦法就是奇異值分解法(singular value decomposition，SVD)[3]。在此基礎上，對M景雷達數據的殘余相位進行時域低通濾波處理，以提取最終的非線性形變量，可表示為

在求解出非線性形變量后，每個相干目標對應的相變序列為

式中vest為線性形變速率。

1.2.3 平均形變速率的重新估計

對于顯著的非線性形變過程，利用前述二維線性模型估計出來的線性速率往往偏小。在恢復整個形變序列的基礎上，可根據對應的每個目標的形變序列，利用最小二乘法重新估計形變速率，即

式中d和T分別為相對起始時刻的形變量和時間間隔。

2 雷達數據及InSAR時序分析處理

研究中共獲取2004年1月—2010年10月近7 a的ENVISAT衛星ASAR數據。該數據覆蓋了德州地區主要沉降區，總數據量為46景(Track－447，Frame－2853)，其基線分布如圖1所示。

圖1 差分干涉圖序列的短基線關系Fig.1 Small baseline of the differential interferograms series

需要說明的是，利用幅度離散指數[1]、相干系數[5]以及子視相關[6]均可以很好地識別出相干目標候選點。對于候選點，只有在滿足模型相關函數最大化約束條件下才能稱為時序分析中的相干目標。實際上，一些目標后向散射特性穩定，但非線性變化強烈，不滿足該條件。因而，可通過降低相關系數的閾值，將其納入到形變解算網絡中，以增加監測點的密度。

時空基線閾值的設置取決于SAR系統參數。空間基線主要考慮雷達系統臨界基線，而時間基線則應根據變形幅度大小，短時間基線相當于降低了形變場的梯度，便于相位解纏和變形參數估計。因而，本研究中設定的短基線條件為:時間基線＜450 d，空間基線＜300 m。

圍繞德州地區地面沉降和地下水位漏斗的基本分布特征，分別針對區域性大范圍沉降和重點沉降中心進行了2類處理:利用2007—2010年4 a間的數據提取了德州地區100 km×100 km平均地面沉降速率;分析全區地面沉降和地下水位變化關系，并確定時序分析過程中穩定的地面參考點。在此基礎上，對德州城區內20 km×20 km范圍的主要沉降中心的7 a數據進行了時序分析處理，提取每個相干目標的沉降率和累積沉降量序列，用以詳細分析多年和年度內沉降中心的時空變化特征。

3 區域性地面沉降與地下水位變化

3.1 區域地面沉降分布

圖2示出德州地區約50 km×60 km范圍內2007年初—2010年底4 a間平均地面沉降速率。

圖2 德州地區2007—2010年CTInSAR監測平均地面沉降率Fig.2 Average subsidence velocity of Dezhou area from 2007 to 2010 derived by CTInSAR data

從區域分布上來看，德州市區沉降較為顯著，沉降區向西北方向延伸，與河北省景縣和故城相連，構成連片沉降區。區內景縣沉降向西北方向擴展，是目前衡水—德州沉降區的重要組成部分，沉降速率均達到40～50 mm/a。德州沉降區向東南方向，分布著3個明顯的沉降中心，分別為陵縣、臨邑縣和平原縣，其中以陵縣最為突出，中心位于縣城西北菜園村，最大沉降速率達75 mm/a。臨邑縣沉降中心位于臨盤采油廠，中心處沉降速率達50～60 mm/a。

3.2 與深層地下水位漏斗的比較分析

InSAR監測結果準確揭示了德州地區區域性沉降和單個漏斗的分布。連續多年的地下水位變化監測結果表明，德州市深層地下水位降落漏斗是目前魯北地區最大的漏斗，與周邊的小漏斗相連接，引起區域性地面沉降[8－9]。深層地下水開采(第Ⅳ，Ⅴ含水層組)是該地區地面沉降的主要原因。以臨邑為例，該沉降中心深層地下水降落漏斗的－30 m等水壓線與德城區深層漏斗相連，已成為德州漏斗的一部分;－40 m等水壓線包圍臨邑縣城和臨盤鎮，圈閉面積達192 km2。漏斗中心在臨盤采油廠一帶，水位埋深超過59.70 m，水位標高－41.82 m左右。圖2中的沉降中心準確地反映了目前該漏斗的分布位置和影響范圍。由于本區地下水水位動態屬于較穩定的連續開采消耗型，年水位變化趨勢與開采量大小密切相關。臨邑深層地下水開采量達3.56×104m3/a。其中，第Ⅳ和第Ⅴ含水層組地下水開采量約占90%，地下水位多年平均降速2.0 m/a，從而使臨邑深層地下水降落漏斗迅速擴展，成為較大的區域性降落漏斗。

4 德城區沉降中心的沉降與回彈

4.1 德城區主要沉降中心

2010年之前連續7 a的InSAR時序分析結果(圖3)表明，德城區有3個典型沉降中心，分別位于德城區陳莊(德棉一廠)、長莊(華源紡織廠)和肖何莊。其中陳莊沉降中心呈南北向分布，長約2.2 km，寬1.5 km，中心最大沉降速率達到55 mm/a;長莊沉降中心呈北東向，長為2 km，寬為1.5 km，中心最大速率也超過50 mm/a。連片分布沉降出現在德城區東部的付莊和宋官屯鎮，這一地區為近年來快速發展的新區。根據1991—2010年水準測量結果[11]，德城區沉降中心累計沉降量為－1 186.9～－636.9 mm，位于德棉一廠(即陳莊漏斗)位置，中心多年平均速率為59.35 mm/a。

圖3 德城區2004—2010年CTInSAR監測平均地面沉降率Fig.3 Average subsidence velocity of Decheng district from 2004 to 2010 derived by CTInSAR data

其中，2005—2006年沉降量為－56 mm;2006—2007年沉降量為－89.0 mm;2007—2010年總沉降量為93.2 mm，年平均沉降量為31.1 mm，與 InSAR監測反映的德城區陳莊沉降中心的變化完全一致。

4.2 地面沉降與回彈的時間效應

由2004—2010年InSAR監測結果發現:德城區多個沉降中心年內均表現出季節性沉降與回彈特征(圖4)。陳莊、長莊、肖何莊以及德棉二廠4個典型沉降中心均出現回彈，回彈最顯著的位于陳莊、長莊和德棉二廠沉降中心。其中陳莊的沉降和回彈幅度均為最大，對應的年份為2007年，2—8月間沉降幅度超過80 mm，至翌年2月回彈15～20 mm。除德棉二廠外，其他3個中心2008年沉降幅度均有所減緩，范圍縮小，2009年繼續減緩，而市區東部在8月份之后沉降加快，2010年繼續增加。

圖4 德城區典型沉降中心的季節性沉降與回彈Fig.4 Seasonal subsidence and rebound of the major center in Decheng district

圖4 所示5個年度的變化特征表明，以陳莊為代表的沉降中心具有明顯的季節性變化特征，表現為每年的2—3月開始快速沉降，至8—9月達到最大，此后逐步回彈，至翌年2—3月重新開始沉降，沉降與回彈之比約為4∶1。

5 地面沉降及回彈與地下水位變化的關系

5.1 持續沉降與深層地下水位的關系

德州地區地下水系統分為4個孔隙含水巖組:第Ⅰ含水巖組潛水層，是主要開采層;第Ⅱ含水巖組屬于咸水層，基本不開采;第Ⅲ含水巖組，底板埋深450～500 m;第Ⅳ含水巖組，底板埋深800～950 m。第Ⅲ，Ⅳ含水巖組屬于承壓含水層，水頭埋深70～100 m，水質普遍較好，是城鎮居民生活和工業用水的主要開采層，德州地下水位降落漏斗屬于該層位，其中心水位(位于陳莊國棉一廠)埋深超過120 m[10]。根據德州深層地下開采引發地面沉降變化閾值的研究[11]:持續開采深層地下水是地面沉降的主要原因，地面沉降與深層水位降落呈現出顯著的線性關系，深層地下水位每下降1 m，產生的地面沉降約為17.5 mm。該經驗值與近4 a來地面測量結果一致，速率為30～40 mm/a。

圖5所示為德城區202長觀孔(位于陳莊漏斗區)1996—2010年的深層水位變化情況，表現出持續下降和年內波動的特點。德城區淺層地下水位變化曲線如圖6所示。

圖5 德城區202長觀孔深層地下水位變化曲線Fig.5 Dynamic curve of deep groundwater level of No.202 hole in Decheng district

圖6 德城區淺層地下水位變化曲線Fig.6 Dynamic curve of shallow groundwaterlevel in Decheng district

5.2 季節性沉降－回彈與地下水位波動的關系

以陳莊漏斗為例，其動態變化(圖3和圖4)與深、淺層地下水位的季節性波動密切相關。每年1—5月降水少，農業用水量大，淺層地下水位急劇下降;7—8月雨量增多、河道水量補給和開采量相對減少，地下水埋深逐漸回升，至年底達到次年最高點。深層承壓水位變化主要取決于開采量，年水位變化與開采量關系密切。每年3—7月，深層地下水連續大量開采，漏斗內壓力水頭急劇下降;8月—翌年2月由于開采量減小，加之大量引入黃河水，水位有所回升，回升幅度一般小于2.0 m。

開采量決定了地面沉降的范圍和強度。目前，在深層水位降落漏斗中心區(陳莊)，開采井密度為334眼/km2。2007年以前，由于深層地下水開采量大，漏斗內壓力水逐年下降，年下降幅度較大，一般2～4 m;2008年以后，由于采取封井措施，深層地下水開采量減少，水位呈現回升趨勢[11]。

陳莊、長莊等中心的沉降與回彈的波動與上述變化十分吻合。圖7分別給出了陳莊(國棉一廠)、長莊(華源紡織廠)、國棉二廠、肖何莊和宋官屯5個沉降中心連續7 a的沉降過程(監測日期為每年的1月30日)，均表現出季節性的非線性波動(紅色線表示初次估計沉降速率)。以陳莊為例，其年度內沉降與回彈與近年深層水位波動相對應，沉降速率在2008年后出現回落，小于2007年之前。

圖7 德城區典型沉降中心2004—2010年地面沉降及深層水位變化序列Fig.7 Subsidence history of coherent target in the typical subsidence center of Decheng district

從地下水天然補給而言，春旱夏澇的氣候特征決定了德州地區春季地下水補給量的不足[11]，而同時用水量又進一步增大，特別是陳莊、長莊和國棉二廠3個沉降中心，均為棉紡、化工等耗水企業，需要大量采水，而汛期澇季時，地表徑流補給及時，地下水開采總量減小，水位出現回升，地面沉降也隨之減緩，并在部分時段內出現回彈的現象。

6 結論

1)2004—2010年間的連續InSAR監測有效揭示了德州地區區域性地面沉降和典型漏斗季節性地面沉降－回彈的發展與變化特征。德州市及其周邊主要地下水降落漏斗與地面沉降中心相互一致，證明了深層地下水開采是本地區地面沉降的主要原因。

2)德城區陳莊、長莊等典型漏斗的季節性沉降－回彈與該地區深淺層地下水位波動相關，地面沉降與地下水位變化表現出良好的相關關系，變化時間基本同步，滯后期約1個月。地面沉降漏斗呈現出3—8月快速沉降，沉降幅度達50～80 mm，8月—翌年2月地面回彈，回彈幅度達10～20 mm的變化特點。

志謝:感謝山東省魯北工程地質勘察院鄒祖光、陳松、楊亞賓以及河北水文地質工程地質四隊杜興明、田小偉等專家的指導和深入討論。

[1]殷躍平，張作辰，張開軍.我國地面沉降現狀及防治對策研究[J].中國地質災害與防治學報，2005，16(2):1－8.Yin Y P，Zhang Z C，Zhang K J.Land subsidence and countermeasures for its prevention in China[J].The Chinese Jounal of Geological Hazard and Control，2005，16(2):1－ 8.

[2]葛大慶，郭小方，張 玲，等.華北平原地面沉降區InSAR監測[R].北京:中國國土資源航空物探遙感中心，2010.Ge D Q，Guo X F，Zhang L，et al.Subsidence monitoring with SAR interferometry in the North China Plain[R].Beijing:China Aero Geophysical Surveying and Remote Sensing Center for Land and Resources，2010.

[3]王 艷，葛大慶，郭小方，等.長江三角洲地區地面沉降InSAR監測[R].北京:中國國土資源航空物探遙感中心，2010.Wang Y，Ge D Q，Guo X F，et al.Subsidence monitoring with SAR interferometry in the Yangtze River Delta area[R].Beijing:China Aero Geophysical Surveying and Remote Sensing Center for Land and Resources，2010.

[4]Ferretti A，Prati C，Rocca F.Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry[J].IEEE Trans Geosci Remote Sensing，2001，38(5):2202－2212.

[5]Mora O，Mallorqui J J，Broquetas A.Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images[J].IEEE Trans Geosci Remote Sensing，2003，41(10):2243－2253.

[6]Charles W，Wegmuller U，Andreas W，et al.Interferometric point target analysis with JERS－1 L－band SAR data[C]//Proceedings of the IEEE Transactfions on Geoscience and Remote Sensing.Toulouse，France，2003:4359－4361.

[7]Hoffmann J，Zebker H A，Galloway D L，et al.Seasonal subsidence and rebound in Las Vegas Valley，Nevada，observed by synethetic aperture Radar interferometry[J].Water Resources Research，2001，37(6):1551－1566.

[8]王 艷，廖明生，李德仁，等.利用長時間序列相干目標獲取地面沉降場[J].地球物理學報，2007，50(2):598－604.Wang Y，Liao M S，Li D R，et al.Subsidence velocity retrieveal from long term coherent targets in Radar interferometric stacks[J].Chinese Jounal Geophys，2007，50(2):598－604.

[9]http://www.terrafirma.eu.com/documents.htm

[10]王小剛，陳 松，王秀芹，等.德州市地面沉降成因及防治對策淺析[J].地質災害與環境保護，2006，17(3):62－66.Wang X G，Chen S，Wang X Q，et al.A Discussion on the causes of the ground subsidence and its countermeasures in Dezhou City[J].Journal of Geological Hazards and Environment Preservation，2006，17(3):62－66.

[11]楊麗芝.德州深層地下水位降落漏斗演變機制與可調控性研究[D].北京:中國地質科學院，2009.Yang L Z.Form principle and controlling－ adjusting research about deep ground water depression cone in Dezhou[D].Beijing:Chinese Academy of Geologieal Seienee，2009.

[12]王勝嶺，宋 波，王德生，等.德州市臨盤采油區地面沉降監測[J].山東國土資源，2009，25(1):25－28，32.Wang S L，Song B，Wang D S，et al.Land subsidence monitoring of Linpan oil mining area in Dezhou City[J].Land and Resources in Shandong Province，2009，25(1):25－28，32.

[13]何國榮，馮在敏，馮克印，等.水準測量網監測在德州市地面沉降中的應用[J].山東國土資源，2012，28(6):28－30.He G R，Feng Z M，Feng K Y，et al.Application of leveling surveying for subsidence monitoring in Dezhou City[J].Land and Resources in Shandong Province，2012，28(6):28－30.