牛地，吳倩，朱成林

(1.安徽省地質(zhì)礦產(chǎn)勘查局332地質(zhì)隊，安徽 黃山 245000；2.安徽省地質(zhì)實驗研究所，安徽 合肥 230041)

0 引言

地面沉降是一種主要由人類工程活動引起的地殼表面標高緩慢降低的環(huán)境地質(zhì)現(xiàn)象，其變化過程緩慢，難以被及時察覺，且不易治理[1]。該現(xiàn)象自二十世紀20年代開始在我國東部沿海出現(xiàn)，現(xiàn)已廣泛分布于100多個地級市[2]，嚴重威脅社會安全與經(jīng)濟健康發(fā)展，發(fā)展科學有效的監(jiān)測手段是預防和減輕災害損失的關(guān)鍵所在。

相較于傳統(tǒng)的地面沉降監(jiān)測方法，合成孔徑雷達差分干涉測量(differential interferometric synthetic aperture radar，D-InSAR)技術(shù)具有全天候、全天時、高精度以及快速獲取大范圍地表高程和形變信息的優(yōu)勢[3]，但時空失相干和大氣效應嚴重阻礙了其獲取真實的地面形變信息。在此基礎(chǔ)上發(fā)展的基于相干目標的時序InSAR處理技術(shù)，主要包括永久散射體(permanent scatterer-interferometric synthetic aperture radar，PS-InSAR)方法和小基線集(small baseline subsets-interferometric synthetic aperture radar，SBAS-InSAR)方法，能在很大程度上減少D-InSAR方法中失相干和大氣異常引起的相位誤差。SBAS-InSAR技術(shù)已廣泛應用于地面沉降[4-5]、滑坡[6-8]、泥石流[9]等地質(zhì)災害的監(jiān)測，其中于地面沉降方面的監(jiān)測成果尤為豐富，致災機理研究多指向地下水[10-15]和礦產(chǎn)開采[16-18]，并在地質(zhì)構(gòu)造活動[19]、工程建設(shè)[20-21]的形變監(jiān)測中取得了良好的反演結(jié)果。上述研究都表明了小基線集方法在地表形變監(jiān)測預警與防治方面有著顯著的技術(shù)的優(yōu)勢。

本文基于SBAS-InSAR技術(shù)和Sentinel-1 A升軌數(shù)據(jù)，對宿州市碭山縣2017—2021年間地表形變進行反演，得到了監(jiān)測期間該區(qū)平均沉降速率和累積地面沉降量，并對地面沉降的時空分布和變化特征進行了細致的分析。在此基礎(chǔ)上結(jié)合相關(guān)水文地質(zhì)資料對該區(qū)地面沉降的主要誘發(fā)因素進行探究，以期通過對碭山縣地面沉降的監(jiān)測，為災害預警提供決策參考依據(jù)。

1 研究區(qū)概況

碭山縣地處淮北平原，位于安徽省最北部，地勢開闊平坦，地貌特征以故黃河沖積平原為主，總面積1 193 km2，下轄13個鎮(zhèn)。縣內(nèi)交通路網(wǎng)發(fā)達，隴海鐵路和310國道貫穿全境，101省道和223省道及縣鄉(xiāng)道路構(gòu)成以縣城為中心的密集交通網(wǎng)絡(luò)。縣境內(nèi)地層屬華北地層大區(qū)晉冀魯豫地層區(qū)徐淮地層分區(qū)淮北地層小區(qū)，據(jù)古近系、新近系和第四系厚度可知，古近紀區(qū)內(nèi)南部地層相對上升，北部則相對下降接受沉積，沉積厚度達數(shù)百米；新近紀開始，全區(qū)主要表現(xiàn)為振蕩下降，接受沉積，新近紀末期地殼上升，普遍形成沉積間斷；第四紀以來以下降為主，接受沉積。碭山縣地下水類型為松散巖類孔隙水，淺層孔隙水底板埋深50 m，深層孔隙水砂層頂板埋深165～195 m，城區(qū)地下水開采以深層為主。在深層地下水大量開采的漏斗區(qū)，人工開采成為深層地下水的主要排泄形式，同時地下水由漏斗外緣向漏斗中心徑流。在碭山縣北部及良梨鎮(zhèn)東北區(qū)域，有第四系全新統(tǒng)淤泥質(zhì)黏土、含淤泥質(zhì)黏土(軟土)分布，一般埋深4 ～10 m，最大埋深18.30 m，厚度小于15.56 m，對于上部建筑構(gòu)成一定潛在威脅。

圖1 研究區(qū)范圍示意圖(底圖源自Google Earth)Fig.1 Schematic diagram of the study area(Base map from Google Earth)

2 研究方法與數(shù)據(jù)處理

2.1 短基線集(SBAS-InSAR)方法

SBAS-InSAR方法通過限定合適的時間、空間基線閾值，將已有的SAR影像數(shù)據(jù)集劃分為若干個小集合，使得各個子集內(nèi)時空基線較小，從而克服時空失相干現(xiàn)象[22]，并能充分利用SAR影像數(shù)據(jù)，在降低解算復雜度的同時保持一定的觀測精度。

假設(shè)在t0…ti…tN時刻獲取了N+1幅SLC SAR影像，在時空基線閾值限定條件下將SAR影像分為L組，之后將影像之間相互配對組成短基線子集，每個短基線子集中的若干影像相互配對后差分干涉處理，L組影像共生成M幅差分干涉圖，規(guī)定N為奇數(shù)，則M可以表示為

(1)

式中：N為SAR影像幅數(shù)，幅；M為差分干涉圖幅數(shù)，幅。

以t0為初始時刻，某一像元(x,r)在ti(i=1,…,N)相對于t0時刻在視線向(Line Of Sight,LOS)的累積形變量d(ti,x,r)(i=1,…,N)為待求量。假設(shè)tA和tB時刻(tA早于tB)的2幅SLC SAR影像配對生成第k(k=1,…,M)幅差分干涉圖，所得相位為δφ(tk,x,r) (k=1,…,M)。

忽略大氣效應和高程殘余等噪聲引起的相位誤差，則(x,r)處的相位為

(2)

式中：δφ(tA,x,r)為tA時刻的影像相位；δφ(tB,x,r)為tB時刻影像相位；d(tA,x,r)為tA時刻視線向形變量，mm；d(tB,x,r)為tB時刻視線向形變量，mm。

這個含有M個等式的N元方程組可以用矩陣表示為

δφ=A[d]。

(3)

(4)

但是，通常情況下影像數(shù)據(jù)并非僅有一個小基線子集(L>1)，此時的ATA就成為了一個奇異矩陣,A的秩為(N-L+1)，方程組有無限多解，基于奇異值分解(singular value decomposition,SVD)方法即可求出方程組在最小范數(shù)意義上的最小二乘解。同時，估計并去除殘余地形相位和大氣延遲相位，即可得到地表形變的時間序列。

2.2 數(shù)據(jù)簡介

Sentinel-1 A是歐空局為“哥白尼計劃”(前身是GMES計劃)發(fā)射的首顆衛(wèi)星，由2顆 C 波段SAR構(gòu)成。本次研究采用2017年5月20日至2021年12月1日間的Sentinel-1 A IW(interferometric wide swath)成像模式下的VV極化升軌數(shù)據(jù)(Ascending Path 142)，共計47景，相鄰兩幅影像之間最短時間間隔均為36 d。此外，POD精密定軌星歷數(shù)據(jù)(POD Precise Orbit Ephemerides)和SRTM1 DEM數(shù)字高程模型被用作計算軌道誤差與模擬地形相位，并從干涉相位中去除。影像參數(shù)信息如表1所示。

表1 Sentinel-1 A雷達影像數(shù)據(jù)參數(shù)Tab.1 Sentinenl-1 A data parameters

2.3 SBAS-InSAR處理流程

首先對Sentinel-1 A數(shù)據(jù)進行裁剪，后對覆蓋研究區(qū)范圍的圖像進行SBAS處理(圖2)，主要包括：選取2017年11月16日的Sentinel-1A影像作為超級主影像，其余所有從影像與主影像進行配準；設(shè)定空間基線閾值130 m(2%)，時間基線閾值200 d，數(shù)據(jù)多視比4∶1(距離向∶方位向)，共生成了220幅干涉圖(圖3);基于SRTM1 DEM計算模擬地形相位對干涉圖進行去平，消除地形相位得到差分干涉圖;采用最小費用流(minimum cost flow,MCF)方法實現(xiàn)相位解纏，去除相干性低和解纏結(jié)果不理想的相對；創(chuàng)建地面控制點對數(shù)據(jù)對進行軌道精煉和重去平，估算和去除殘余的恒定相位和解纏后殘存的相位斜坡；基于SVD估算平均位移速率和殘余地形相;分別設(shè)定大氣高通(365 d)和大氣低通(1 200 m)進行時間域和空間域的濾波，估算和去除大氣延遲相位，得到研究區(qū)在監(jiān)測時間段內(nèi)的地表時序形變；對SBAS反演結(jié)果進行地理編碼，得到研究區(qū)形變監(jiān)測結(jié)果。

圖2 SBAS-InSAR技術(shù)流程Fig.2 Technical process of SBAS-InSAR

圖3 時空基線連接圖Fig.3 Space-time Baseline Connection Diagram

3 結(jié)果分析

3.1 年均形變速率

通過對碭山縣近5 a來的47期Sentinel-1 A IW SLC影像做SBAS反演處理，獲取了研究區(qū)2017—2021年視線向的年平均地表形變速率。忽略水平位移的微弱影響，將視線向監(jiān)測結(jié)果除以cosβ(β為入射角)得到垂直方向形變速率(圖4)，其中負值代表地面向遠離衛(wèi)星傳感器方向運動(沉降)，而正值則代表地面向靠近衛(wèi)星傳感器方向移動(抬升)。

圖4 碭山縣2017—2021年平均垂直方向地表形變速率Fig.4 Average vertical surface deformation rate in Dangshan County from 2017 to 2021

碭山縣地面在觀測期間整體呈現(xiàn)出“北部輕微沉降，南部穩(wěn)定，碭城鎮(zhèn)局部沉降嚴重”的特點。縣城碭城鎮(zhèn)南部以東城社區(qū)為中心，半徑5 km范圍內(nèi)均屬嚴重沉降區(qū)域，東城社區(qū)附近最大沉降速率達50 mm/a，碭城火車站附近沉降速率約為30 mm/a；高鐵新區(qū)東部與經(jīng)濟開發(fā)區(qū)西部，沉降速率約為40 mm/a，這些區(qū)域的工程建設(shè)等人類活動比較活躍。官莊壩鎮(zhèn)劉堤頭村—玄廟鎮(zhèn)孫黑樓村—周寨鎮(zhèn)大孟莊村一線沉降速率也較高，為25～30 mm/a，可能與這一帶存在第四系全新統(tǒng)淤泥質(zhì)黏土、含淤泥質(zhì)軟土分布有關(guān)。綜合來看，以碭山縣境內(nèi)隴海鐵路為界，北部大片區(qū)域沉降速率基本為6～20 mm/a，約占總形變面積的78.20%；南部區(qū)域地表穩(wěn)定，零散分布微弱沉降。

3.2 地面沉降場時空分布特征

以監(jiān)測起始時間2017年11月16日的影像為參考基準計算其余時相的累計形變，繪制得到碭山縣2017—2021年地表累計時序形變圖(圖5)。從圖中可以看出，監(jiān)測期間內(nèi)碭山縣的沉降區(qū)域不斷擴大，累計沉降量逐年增大，沉降區(qū)域主要集中在以碭城鎮(zhèn)為中心的縣城周圍。碭城鎮(zhèn)沉降中心從2018年初顯形態(tài)，此時沉降量為55 mm，2019年之后逐漸明顯，至2019年下半年，沉降范圍基本確定，即以縣城為中心，向外輻射半徑5 km左右，沉降面積約78.50 km2。隨著時間的推移，該區(qū)域沉降量日益增加，到2021年12月形變中心沉降量達186 mm。此外還在縣城與李莊鎮(zhèn)之間監(jiān)測到一條較明顯的條帶狀沉降區(qū)，其沉降進程與縣城沉降中心基本一致，隴海鐵路、310國道正位于這條沉降帶上。縣城以北的官莊壩鎮(zhèn)、玄廟鎮(zhèn)，在2017—2021年期間，累計沉降超過120 mm，沉降速率約為30 mm/a。碭山縣東北部分的葛集鎮(zhèn)、唐寨鎮(zhèn)沉降也較為明顯，葛集鎮(zhèn)南部與唐寨鎮(zhèn)北部相鄰，此處沉降量最大超過100 mm，沉降速率約為20 mm/a。碭山縣西南地區(qū)，趙屯鎮(zhèn)至關(guān)帝廟鎮(zhèn)一帶，地表有輕微抬升，推測可能與鄭徐客運專線工程建設(shè)有關(guān)。

圖5 碭山縣2017—2021年地表累計時序形變量Fig.5 The cumulative displacement of Dangshan County from 2017 to 2021

3.3 時間序列分析

為更直觀地表達地表位移隨時間的變化，在嚴重沉降區(qū)、一般沉降區(qū)和穩(wěn)定區(qū)分別選定Ⅰ、Ⅱ、Ⅲ三點進行時序分析(三點位置見圖4)。圖6為Ⅰ、Ⅱ、Ⅲ三點的時序累計形變曲線，其中Ⅰ點位于地面沉降漏斗中心，沉降速率高達55.66 mm/a；Ⅱ點靠近沉降中心，沉降速率約為33.15 mm/a；Ⅲ點處于整體穩(wěn)定區(qū)。監(jiān)測期間各點的沉降量與時間基本呈線性關(guān)系。

圖6 時間序列形變量Fig.6 Time series deformation

值得注意的是，Ⅰ、Ⅱ兩點累計形變在2018年7—8月出現(xiàn)小規(guī)?；貜?，可能與該時段的異常強降雨有關(guān)(圖7)。受“溫比亞”臺風的影響，宿州市 8月16日開始降小到中雨，17日普降大到暴雨，18日普降大暴雨，部分地區(qū)特大暴雨，具有暴雨強度大、范圍廣等特點。18日，暴雨中心區(qū)域位于碭山縣全境，碭山縣1 d、2 d、3 d最大降雨量分別為216.80 mm、255 mm、264.40 mm，均居有資料記錄以來第一位，重現(xiàn)期分別為 60 a、超百年、90 a。2018年碭山縣地下水資源量比2017 年增長12.40%，比2019年增長24%，比多年平均值增長28.10%。

嚴重沉降區(qū)域分布在人口密集的城區(qū)，且未采取相應防治措施，沉降量級和范圍仍在快速增大，應盡快開展防治預警工作，合理限制地下水開采，開展第四系結(jié)構(gòu)分布詳細調(diào)查，條件允許的情況下采取人工回灌措施。

圖7 2018年每月降水量統(tǒng)計(蕭縣氣象站)Fig.7 Statistics of Monthly Precipitation in 2018(Meteorological Station in Xiao County)

3.4 沉降機理分析

前人大量研究證實，過量開采地下水是導致地面沉降的外因，而中等、高壓縮性黏土層和承壓含水層的壓實則是地面沉降的內(nèi)因[23]。

碭山縣城區(qū)有開采井20余眼，井深多為150～300 m，主要開采中、下更新統(tǒng)及新近系含水層。2016—2020年，碭山縣深層地下水年平均開采量約2 268萬 m3，其中城區(qū)約占50%。據(jù)碭山縣2008年水文調(diào)查資料，縣城區(qū)漏斗中心水位埋深約為44 m，最大動水位埋深達61.20 m。與已發(fā)生地面沉降且地質(zhì)背景相似的阜陽市類比[24]，在降落漏斗區(qū)，當區(qū)域地下水位埋深下降至20 m時，即可引發(fā)地面沉降。

本文收集了碭山縣12眼監(jiān)測井2017—2022年的年平均水位變化情況，12眼監(jiān)測井可按成井深度分為50 m、150 m和350 m(兩眼200 m左右)3個層級，結(jié)尾編號分別為A、B和C，各監(jiān)測井的歷年水位變化情況如圖8所示。DS01A、DS01B、DS01C這3眼井位于碭城鎮(zhèn)李洼村，緊鄰300 mm厚度軟土層，2017—2022年間，3眼監(jiān)測井地下水位不斷下降，DS01A地下水位從42.06 m下降至45.96 m，下降了3.9 m；DS01B地下水位從40.11 m下降至43.6 m，下降了3.49 m；DS01C成井深度350.3 m，6 a內(nèi)地下水位下降了16.28 m，周邊村莊沉降速率為15～16 mm/a。210A、210B、210C位于玄廟鎮(zhèn)花園村花園小學附近，210A地下水位基本穩(wěn)定，210B和210C水位略有下降，下降值3～4 m，周邊村莊基本穩(wěn)定，未見明顯沉降。610A位于良梨鎮(zhèn)派出所，水位下降了3.04 m；610B、610C位于良梨鎮(zhèn)禮河集污水處理廠，610B年均水位下降1.93 m，610C缺失2017年監(jiān)測數(shù)據(jù)，后4 a年均水位抬升0.55 m，周邊地區(qū)基本穩(wěn)定。614A、614B、614C位于關(guān)帝廟鎮(zhèn)黃樓農(nóng)場，614A監(jiān)測井缺失2018年和2019年數(shù)據(jù)，其余4 a水位穩(wěn)定；614B水位有波動，波動范圍在2 m以內(nèi)；614C缺失2017、2018年監(jiān)測數(shù)據(jù)，后4 a年均水位下降了5.56 m。

(a)DS01A 成井深度50.2m (b)210A 成井深度46.19m (c)614A 成井深度48.56m

碭山縣城區(qū)開采井過于集中，而地表巖性為第四系全新統(tǒng)粉土，地下松散沉積物(N+Q)厚度近500 m，極易形成局部降落漏斗。碭山縣的軟土區(qū)主要分布在縣城北部地區(qū)，占全縣總面積的40.71%，地層巖性主要為第四系全新統(tǒng)淤泥質(zhì)黏性土?？h域北部2條軟土埋深等值線顯示軟土埋深為10 m，2條等值線之間正是官莊壩、玄廟鎮(zhèn)和周寨鎮(zhèn)(圖9)，與黃河故道以北玄廟鎮(zhèn)EW向沉降區(qū)空間位置基本吻合?？梢酝茢啵叵滤_采和軟土變形是造成碭山縣地面沉降的主要因素。

圖9 碭山縣軟土埋深等值線Fig.9 Soft soil depth isoline in Dang Shan County

4 結(jié)論

本文基于SBAS-InSAR技術(shù)和Sentinel-1 A數(shù)據(jù)反演了碭山縣2017—2021年地表形變，探究了沉降的分布規(guī)律和演化樣式。進一步結(jié)合水文地質(zhì)與巖土資料，得出以下結(jié)論：

(1)監(jiān)測時間段內(nèi)碭山縣地區(qū)地面沉降顯著，中心沉降速率高達50 mm/a。

(2)地面沉降中心主要位于人口密集的碭城鎮(zhèn)，沉降漏斗已形成規(guī)模。

(3)碭山縣北部地區(qū)的沉降與淤泥質(zhì)黏土(軟土)廣泛分布有關(guān)，軟土上部荷載超過其承載力時，易發(fā)生差異性沉降。

(4)地下水開采過程中含水層與隔水層壓實效應是導致沉降發(fā)生的直接原因。

為避免或減少碭山縣地面沉降可能引發(fā)的災害，應對重點沉降的城區(qū)進行實時地面監(jiān)測與防治，并在規(guī)劃地下水開采與城區(qū)建設(shè)選址中充分考慮地面沉降的變化趨勢和潛在危險，緩解工程性地面沉降。

致謝：南京大學常豐年博士為本研究與論文寫作提供了非常有益的幫助，匿名審稿人的建設(shè)性意見極大地改善了本文的質(zhì)量，在此深表感謝！