馮曉松 楊成生 趙超英 牛玉芬 高 涵
1 長安大學(xué)地質(zhì)工程與測繪學(xué)院,西安市雁塔路126號,710054 2 河北工程大學(xué)礦業(yè)與測繪工程學(xué)院,河北省邯鄲市太極路19號,056038 3 邯鄲市自然資源空間信息重點(diǎn)實(shí)驗(yàn)室,河北省邯鄲市太極路19號,056038 4 云南省地震局信息中心,昆明市茨壩藍(lán)桉路2號,650201
地裂縫是一種由內(nèi)部和外部地質(zhì)動力過程所引起的淺表破裂[1],通常表現(xiàn)為土壤地表層在各種人為因素(采礦、澆灌、抽水等)[2-4]及其他自然因素(構(gòu)造斷裂[4]、地震[5])影響下,土層的連續(xù)性被力的相互作用打斷,并以一種重大地裂縫和地表裂縫破壞現(xiàn)象呈現(xiàn)出來[6]。地裂縫直接或間接破壞了人類和野生動物的生存環(huán)境[7],并且地裂縫周邊地區(qū)的地表往往會出現(xiàn)沉降現(xiàn)象[8]。
隆堯地裂縫位于華北平原隆堯縣內(nèi),形成于1966年邢臺地震[9],震后裂縫進(jìn)入非活躍期,2003年裂縫再次暴露[10],2006年以來裂縫規(guī)模及活躍度迅速增加[11]?,F(xiàn)今,隆堯地裂縫是華北平原規(guī)模最大、災(zāi)害程度最為嚴(yán)重的構(gòu)造地裂縫[12],分布在近EW向的20多個村莊里[13]。地表裂縫最大寬度為20 cm,可見深度達(dá)數(shù)m[14]。Yang等[15]研究了2007~2011年間隆堯地裂縫的形變特征,并進(jìn)行建模處理,但隨著《全國地面沉降防治規(guī)劃(2011~2020年)》、《京津冀平原地面沉降綜合防治總體規(guī)劃》等文件的發(fā)布及地下水控采方案和南水北調(diào)工程的實(shí)施,隆堯地裂縫的特征可能發(fā)生了改變,但目前缺乏對其現(xiàn)階段活動性的研究。
有學(xué)者使用全球?qū)Ш叫l(wèi)星系統(tǒng)GNSS、全站儀等手段對地裂縫進(jìn)行了大量的形變監(jiān)測研究,但傳統(tǒng)測量技術(shù)成本昂貴且數(shù)據(jù)獲取不連續(xù),嚴(yán)重阻礙了對地裂縫的進(jìn)一步研究[16]。干涉合成孔徑雷達(dá)InSAR技術(shù)是一項(xiàng)基于地球物理學(xué)、信號學(xué)、計算機(jī)等學(xué)科發(fā)展而來的新型技術(shù),被廣泛應(yīng)用于火山監(jiān)測、地震監(jiān)測、滑坡監(jiān)測、地裂縫監(jiān)測等地質(zhì)災(zāi)害監(jiān)測領(lǐng)域。本文采用時序InSAR技術(shù)探測隆堯地裂縫及其周邊地區(qū)的形變,并進(jìn)一步分析現(xiàn)階段隆堯地裂縫的活動特征。
中國大陸位于印度、亞歐和太平洋等板塊大量匯聚形成的沖擊區(qū)和擠壓區(qū),華北平原長期遭受拉張應(yīng)力[17],因此1 000多條裂縫遍布華北平原[18]。隆堯地裂縫位于華北平原的邢臺市,處于華北地塊的核心位置,擁有3個地質(zhì)塊:臨清地塹、邢-衡隆起、寧晉-衡水?dāng)鄬?,形成了一個“三叉戟”的形狀。該區(qū)域分布有NE向和EW向斷層,其中主要裂隙是太行山山前斷裂和巨鹿EW向斷裂。
隆堯地裂縫是華北平原最長、最活躍的地裂縫,位于隆堯縣城與魏家鎮(zhèn)之間(圖1 紅色實(shí)線),全長約36 km,地表裂縫最大寬度為20 cm,可見深度達(dá)數(shù)m[19],沿途的道路、河流、農(nóng)田、房屋等都被其切斷,造成了巨大的經(jīng)濟(jì)損失。在隆堯地裂縫周邊分布有太行山山前斷裂、臨城斷裂、巨鹿西部斷裂、巨鹿東部斷裂、雞澤斷裂和明華斷裂,這些斷裂與隆堯地裂縫相互作用,使得隆堯地裂縫周邊地區(qū)的地表沉降尤為復(fù)雜。由于研究區(qū)域較為廣闊,本文采用InSAR技術(shù)對河北隆堯地裂縫及其周邊的地表形變進(jìn)行監(jiān)測。
圖1 隆堯地裂縫及其周圍地區(qū)地形Fig.1 Topographic of Longyao ground fissure and its surrounding area
哨兵一號Sentinel-1采用的是C波段,其中心頻率為5.405 GHz,C波段在X波段和S波段之間,兼具兩者的特性。本文采用31景VV極化的Sentinel-1數(shù)據(jù),所選時段為2019-01-03~2020-12-23,如表1所示。數(shù)據(jù)覆蓋范圍如圖1所示,圖中藍(lán)色實(shí)線代表裁剪后的雷達(dá)影像范圍,即本文的研究區(qū)域。本文使用的外部DEM數(shù)據(jù)是由美國地質(zhì)勘探局USGS提供的SRTM DEM,分辨率為30 m。
由于研究區(qū)域植被覆蓋度較高,失相干現(xiàn)象較為嚴(yán)重,為了通過不同時間基線干涉圖提供的多余觀測來降低失相干的影響,并去除干涉圖中的各種誤差(如大氣誤差、DEM殘差等),本文首先采用SBAS-InSAR技術(shù)獲得隆堯地裂縫及周邊地區(qū)的時間形變序列。SBAS方法是在利用小基線干涉圖集進(jìn)行形變測量的前提下,合成符合條件的小基線干涉圖對,然后基于形變速率的最小范數(shù)準(zhǔn)則,通過SVD方法獲得系數(shù)矩陣的偽逆,最終獲得相干目標(biāo)的形變速率及時間序列[20]。計算公式如下:
(1)
表1 實(shí)驗(yàn)數(shù)據(jù)
為獲取研究區(qū)域含有較小誤差的形變速率,使用Stacking InSAR技術(shù)對解纏后的差分干涉圖進(jìn)行處理。Stacking InSAR技術(shù)假設(shè)大氣在時間上的擾動是隨機(jī)的,且地表形變?yōu)榫€性形變,因此可以將所有符合運(yùn)算條件、經(jīng)過解纏的干涉圖進(jìn)行線性疊加,從而減少信號穿過對流層、電離層所產(chǎn)生的大氣誤差對監(jiān)測結(jié)果的影響,提高監(jiān)測精度[21]。公式如下:
(2)
式中,Vi為第i幅干涉圖的相位變化速率,φi為第i幅干涉圖的解纏相位,ΔTi為第i幅干涉圖的成像時間間隔。
數(shù)據(jù)處理流程如下:1)首先對原始數(shù)據(jù)進(jìn)行預(yù)處理,轉(zhuǎn)化為SLC格式,然后進(jìn)行基線估計,空間基線設(shè)置為200 m,時間基線設(shè)置為100 d,時空基線如圖2所示;2)將得到的干涉對進(jìn)行共軛相乘,得到差分干涉圖,并對干涉圖進(jìn)行濾波以提高影像的信噪比;3)通過最小費(fèi)用流法對濾波后的干涉圖進(jìn)行相位解纏,得到地表像元的真實(shí)相位;4)分別通過SBAS-InSAR技術(shù)和Stacking InSAR技術(shù)對研究區(qū)域的數(shù)據(jù)進(jìn)行時序處理。經(jīng)過上述數(shù)據(jù)處理得到研究區(qū)域的形變時間序列和形變速率,其中形變速率如圖3所示,2019-01-03~2020-12-23的累積形變量如圖4所示。
圖2 時空基線Fig.2 Spatial-temporal baseline
圖3 2019~2020年隆堯地裂縫周邊地區(qū)形變速率Fig.3 Deformation rate of Longyao ground fissure from 2019 to 2020
圖4 2019~2020年隆堯地裂縫周邊地區(qū)累積沉降量Fig.4 Cumulative subsidence around Longyao groundfissure from 2019 to 2020
3.2.1 隆堯地裂縫
由圖3和4可見,隆堯地裂縫的南北兩側(cè)有一定的形變速率差異,為了更加直觀形象地觀察隆堯地裂縫南北兩側(cè)的速率差,對隆堯地裂縫南北兩側(cè)進(jìn)行剖線分析。為與Yang等[15]2007~2011年的ALOS PALSAR結(jié)果進(jìn)行比較,本文選擇相近位置的剖線(圖3黑色虛線EF)。隆堯地區(qū)形變速率剖線所得結(jié)果如圖5所示,圖中淺灰色實(shí)線為隆堯地裂縫,實(shí)線左側(cè)為地裂縫北側(cè)、右側(cè)為地裂縫南側(cè)。由圖可見,相對于隆堯地裂縫南側(cè),隆堯地裂縫北側(cè)呈抬升趨勢,二者的形變速率差達(dá)4 cm/a。
圖5 沿剖線EF隆堯地裂縫南北兩側(cè)形變速率差異Fig.5 Deformation rate of both sides of Longyao ground fissure along profile EF
將本文采用Sentinel-1衛(wèi)星影像獲取的2019~2020年隆堯地裂縫南北側(cè)形變速率剖面,與采用日本地球觀測衛(wèi)星合成孔徑雷達(dá)ALOS PALSAR-1獲取的2007~2011年形變速率剖面[14]進(jìn)行對比分析發(fā)現(xiàn),前者的隆堯地裂縫南側(cè)沉降速率比后者的相對沉降速率約低1 cm/a。搜集資料后發(fā)現(xiàn),國務(wù)院在2012年發(fā)布了《全國地面沉降防治規(guī)劃(2011~2020年)》,隨后自然資源部又組織編制了《京津冀平原地面沉降綜合防治總體規(guī)劃》。隨著地下水控采措施和南水北調(diào)工程的實(shí)施,京津冀地區(qū)的沉降狀況在2016年以后得到了大幅度的緩解,沉降范圍和沉降中心形變量級的降低可能影響了隆堯地裂縫的活動性。此外,由剖線圖可以明顯看出,隆堯斷層兩側(cè)形變特征存在明顯差異,這一現(xiàn)象與Yang等[15]的研究中L波段獲取的形變特征類似。
3.2.2 地裂縫周邊地區(qū)
由圖3和4可見,除巨鹿縣縣城所在位置有輕微的抬升外,境內(nèi)其他地區(qū)均處于沉降狀態(tài)。其中,固城店鎮(zhèn)、魏家鎮(zhèn)、官莊鎮(zhèn)是明顯的沉降中心,3個區(qū)域年沉降速率均達(dá)6~8 cm/a。對固城店鎮(zhèn)、魏家莊鎮(zhèn)、官莊鎮(zhèn)的形變時間序列進(jìn)行統(tǒng)計分析(圖6)可知,固城店鎮(zhèn)、魏家莊鎮(zhèn)及官莊鎮(zhèn)處于周期沉降狀態(tài),即當(dāng)年2~9月處于快速沉降狀態(tài),當(dāng)年9月至次年2月又存在一定程度的抬升,但其總體形變趨勢仍是沉降。
圖6 固城店鎮(zhèn)、魏家莊鎮(zhèn)及官莊鎮(zhèn)單點(diǎn)形變時序Fig.6 Time series subsidence of Guchengdian, Weijiazhuang and Guanzhuang
調(diào)查發(fā)現(xiàn),固城店鎮(zhèn)、魏家莊鎮(zhèn)、官莊鎮(zhèn)都種植了大面積的農(nóng)作物,包括冬小麥、玉米、棉花、蔬菜等,每年需要抽取大量的地下水進(jìn)行灌溉。由于河北獨(dú)特的氣候條件,當(dāng)?shù)氐闹饕r(nóng)作物(小麥、玉米)為一年兩收模式,冬小麥通常在當(dāng)年10月播種并進(jìn)行初次灌溉,次年開春時進(jìn)行第2次灌溉。由此推測,當(dāng)?shù)貙Φ叵滤拇罅砍槿〖爸苓厰鄬拥挠绊懺斐闪?~10月的沉降狀態(tài),之后隨著灌溉抽水用量的急劇減少及冬季降雪和周圍地下水回流,區(qū)域表現(xiàn)為緩慢的地表抬升[22]。此外,隆堯地裂縫東南端有一個較大的形變區(qū)域(圖3和4中區(qū)域D),查閱資料發(fā)現(xiàn),當(dāng)?shù)赜?003-11成立了隆堯工業(yè)園,主要生產(chǎn)方便面和飲料,其工業(yè)過程消耗大量的地下水,可能是造成當(dāng)?shù)氐孛娉两档闹饕颉?/p>
研究表明,彈性均勻空間中的蠕變位錯[6,23-24]可用來解釋活動斷層的大地觀測。一個活動斷層的地表變形可以被建模為反正切,為斷層滑動速率和鎖定深度提供了估計方法。垂直走滑斷層運(yùn)動的一個簡單模型是均質(zhì)彈性半空間的,可以描述為:
(3)
式中,V(x)為沿斷層垂直剖面估計的各點(diǎn)平行斷層速度,x為離斷層的距離,S為斷層滑動速率,D為位錯深度,α為靜態(tài)偏移量。
為更好地了解隆堯地裂縫的地面變形,使用上述模型對InSAR形變剖面進(jìn)行建模。通過計算得知,地裂縫深度為5 m,幾乎破裂至地表,滑移速率為27 mm/a,相比于2007~2011年的48 mm/a,滑移速率減緩約21 mm/a,說明現(xiàn)階段隆堯地裂縫雖仍處于活躍狀態(tài),但隨著對京津冀地區(qū)地下用水的宏觀調(diào)控,地裂縫活躍度降低。
本文采用SBAS-InSAR技術(shù),基于31景Sentinel-1衛(wèi)星影像對邢臺市隆堯地裂縫及其周邊區(qū)域進(jìn)行時間序列分析和物理模型構(gòu)建。結(jié)果表明,隨著對京津冀地區(qū)地下水使用的宏觀調(diào)控,隆堯地裂縫及其周邊地區(qū)形變趨勢有所減緩,活躍度降低。具體結(jié)論如下:
1)現(xiàn)階段隆堯地裂縫兩側(cè)相對形變有所降低,相較于2007~2011年降低約1 cm/a;
2)隆堯地裂縫周邊地區(qū)如固城店鎮(zhèn)、魏家莊鎮(zhèn)、官莊鎮(zhèn)是明顯的沉降中心,3個區(qū)域年沉降速率均達(dá)6~8 cm/a,且呈現(xiàn)周期性的形變特征;
3)物理模型構(gòu)建結(jié)果表明,隆堯地裂縫幾乎破裂至地表,滑移速率約為27 mm/a,相較于2007~2011年減緩約21 mm/a。