江培華,胡寶俊,唐 偉
(1. 華北科技學(xué)院建筑工程學(xué)院,北京 東燕郊 065201;2. 中國礦業(yè)大學(xué)(北京)地球科學(xué)與測繪工程學(xué)院,北京 100083)
地面沉降是一種在自然或人為因素下導(dǎo)致地下松散巖層固結(jié)壓縮,并引起地面一定范圍內(nèi)標(biāo)高降低的地質(zhì)現(xiàn)象[1,2]。對于這種前期難發(fā)現(xiàn)、影響范圍廣、治理難度大的地質(zhì)災(zāi)害已成為影響我國經(jīng)濟社會可持續(xù)發(fā)展的重要因素。近年來,由于地下水過度開采而帶來的地面沉降問題越來越引起人們的關(guān)注,華北平原更是此類地質(zhì)災(zāi)害的受災(zāi)區(qū),地表非均勻的沉降已經(jīng)對人民的生產(chǎn)生活帶來了嚴(yán)重的經(jīng)濟損失,特別是對城鎮(zhèn)線狀工程的破壞。為了避免此類地質(zhì)災(zāi)害繼續(xù)造成威脅,必須對地表沉降情況進行有效監(jiān)測和防治。
目前,地面沉降的監(jiān)測方法主要有分層標(biāo)監(jiān)測、高精度水準(zhǔn)測量等。這些傳統(tǒng)的監(jiān)測手段不僅耗資巨大,而且效率不高,傳統(tǒng)監(jiān)測方法中水準(zhǔn)測量復(fù)測時間跨度過大,而分層標(biāo)和GPS高程測量空間分辨率低,均不利于局部地區(qū)形變場時空演化特征分析與研究。合成孔徑雷達差分干涉測量(Differential Interferometric Synthetic Aperture Radar,D-InSAR)技術(shù)是近二十年發(fā)展起來的一種新型空間對地觀測技術(shù),可以利用重復(fù)觀測的合成孔徑雷達(Synthetic Aperture Radar,SAR)影像相位信息獲取大范圍、高精度、面狀分布的地面沉降形變信息,且監(jiān)測精度理論上可以達到mm級。同時,依據(jù)雷達衛(wèi)星的過境頻率,實現(xiàn)地面沉降的動態(tài)跟蹤監(jiān)測,對地質(zhì)災(zāi)害的預(yù)防做出了突出的貢獻。
然而,傳統(tǒng)的D-InSAR技術(shù)由于其易受大氣效應(yīng)的影響,以及精度低的局限性,為此,基于傳統(tǒng)D-InSAR技術(shù)的局限性研究分析,提出了永久散射體合成孔徑雷達(Permanent Scatterer InSAR,PS-InSAR)[3]和小基線集合成孔徑雷達干涉測量技術(shù)(Small Baseline Subset,SBAS-InSAR)[4],這兩種時序InSAR技術(shù)通過提取散射信息穩(wěn)定的點反演形變信息,克服時空失相關(guān),并通過時空濾波方法改正大氣延遲誤差問題。
廊坊地區(qū)地下水的過度開采和城鎮(zhèn)化建設(shè)發(fā)展迅速,蔣喆等利用2015—2017年Sentinel-1A數(shù)據(jù)基于SBAS-InSAR技術(shù)對廊坊市進行了監(jiān)測分析[5],得出廊坊地區(qū)的廣陽區(qū)南尖塔鎮(zhèn)存在明顯的地面沉降,最大沉降速率達60.4mm/a,累計最大沉降量達107.6mm。周呂等[6]獲取2007—2010年覆蓋廊坊地區(qū)18景Envisat ASAR影像數(shù)據(jù),基于可同時獲取點目標(biāo)和分布式目標(biāo)沉降信息的多時相InSAR(Multiple Temporal InSAR,MTInSAR)技術(shù),并結(jié)合前人研究結(jié)果、地下水開采以及工業(yè)發(fā)展?fàn)顩r等研究分析了廊坊市的地面沉降時空特征。李海君等結(jié)合振幅離差指數(shù)與相干系數(shù)識別,采用2007—2010年覆蓋廊坊北部地區(qū)的33景Envisat ASAR重軌影像數(shù)據(jù)使用PS-InSAR技術(shù),選取PS點,借助GIS平臺,對研究區(qū)進行沉降監(jiān)測、沉降差異性特征與成因進行分析[7]。周旭等利用SABS算法發(fā)展而來的InSAR TS+AEM(InSAR time series with atmospheric model)時序分析模型,分析了24幅Envisat ASAR影像數(shù)據(jù),研究表明利用時序InSAR技術(shù)進行城市地表沉降監(jiān)測具有良好的精度和穩(wěn)定性[8]。上述研究均只獲得雷達視線方向(line-of-sight,LOS)的形變信息,并假設(shè)形變只發(fā)生在垂直向,將LOS形變簡單投影為垂直方向,存在一定的誤差。
廊坊市受城市工業(yè)、生活開采深層地下水的影響,形成了以城區(qū)為中心的地下水位漏斗,漏斗斷面呈上寬下窄的“V”字形[9]。近年來,廊坊市陸續(xù)出臺了地下水限采、禁采等管理辦法,并關(guān)停了城區(qū)自備井,深層地下水位開采得到了有效控制,水位下降幅度減緩。
本文基于Sentinel-1衛(wèi)星2017年1月至2020年12月升、降兩軌數(shù)據(jù),利用SBAS-InSAR方法獲取了2個軌道的LOS形變,并將兩者融合獲取了垂直向和水平向形變。在此基礎(chǔ)上,分析了廊坊市地面沉降的時空演變格局。
廊坊市位于河北省中部偏東,由于1974年行政區(qū)劃變動,北京和和天津地域相交,使廊坊形成北部三縣(市)與境域主體不相連接的版圖,其北起燕山南麓丘陵地區(qū),和首都北京為鄰,南接壤滄州,可抵黑龍港流域。西連保定市的涿州、雄縣、高碑店,東交天津市武清區(qū)、寶坻區(qū)、薊州區(qū),地處京津冀城市群核心地帶、環(huán)渤海灣腹地[10];廊坊市轄安次區(qū)、廣陽區(qū)兩個區(qū)、三河市、霸州市兩個縣級市以及香河縣、永清縣、固安縣、文安縣、大城縣和大廠回族自治縣六個縣,總面積約6429km2,截至2021年末,常住人口為553.82萬人。
廊坊市受地質(zhì)構(gòu)造的影響,大部處于凹陷地區(qū),隨著地殼下沉,地面逐漸被第四紀(jì)沉積物填平,致使新生界地層沉降厚度較大,全市地貌比較平緩單調(diào),以平原為主,一般高程在2.5~30m之間,平均海拔13m左右。由于洪積、沖積作用和河流多次決口改道淤積,沉積物交錯分布,加上風(fēng)力及人為活動的影響,境內(nèi)地貌差異性較大,緩崗、洼地、沙丘、小型沖積堆等遍布,全市地貌呈現(xiàn)大平小不平狀態(tài)。
研究區(qū)屬于華北平原中北部,廊坊市南部整個區(qū)域,介于東經(jīng)116.06°~116.96°,北緯38.17°~39.65°之間,范圍如圖1中黑色粗線包含的區(qū)域。
圖中黑色粗線為廊坊市行政區(qū)劃。方框為Sentinel-1衛(wèi)星的覆蓋范圍(紅色方框為升軌,黃色方框為降軌)。
Sentinel-1衛(wèi)星是歐洲航天局哥白尼計劃(Global Monitoring for Environment and Security,GMES)中的地球觀測衛(wèi)星,由兩顆載有C波段合成孔徑雷達的衛(wèi)星組成。本研究收集了2017年至2020年的升軌(P142F121、P142F126)和降軌(P47F461)兩個軌道數(shù)據(jù),其中升軌數(shù)據(jù)為2017年5月至2020年12月共102幅影像,覆蓋范圍如圖1中紅色框所示。降軌數(shù)據(jù)為2017年1月至2020年12月共97幅影像,覆蓋范圍如圖1中黃色框所示。
本文在選取時空基線時,時間基線閾值選擇50天,空間基線閾值選擇100m,對于時空基線不相連的情況,特別是在降軌影像間,因此手動增加相關(guān)干涉對(圖2(b)中紅色線),干涉對組合時空基線如圖2所示。最終,升軌數(shù)據(jù)共生成290個干涉對,降軌數(shù)據(jù)共產(chǎn)生248個干涉對。SBAS-InSAR算法重在應(yīng)用短基線的差分干涉圖集,集合內(nèi)的影像對基線距小,而集合間的影像對基線距大,集合構(gòu)成相位回歸分析的序列,用均值相干系數(shù)作為相干目標(biāo)識別的指標(biāo),再根據(jù)最小范數(shù)準(zhǔn)則,利用奇異值分解算法逐個分離相位成分,實現(xiàn)每個相干目標(biāo)形變序列的獲取,減少空間和時間去相關(guān)影響,保證精度與可靠性[1]。

圖2 升降軌數(shù)據(jù)干涉對組合時空基線
SBAS方法為了在提高點目標(biāo)密度的同時保持干涉圖中的相干性,采用多主影像方式,通過組合時間和空間基線閾值內(nèi)的干涉對,提取在一定時間內(nèi)保持相干的分布式散射體目標(biāo)。假設(shè)獲取同一地區(qū)、時間分別為t0,t1,t2,…,tn的N+1幅SAR影像,這些影像可以生成的干涉圖數(shù)量為M,且M滿足:
(1)
其中第j幅干涉圖為ta和tb兩個時刻的SAR影像組合而成,第j幅干涉圖在某一地面目標(biāo)點處的干涉相位δφab表示為:

(2)
式中,λ為波長;da和db為ta和tb兩個時刻相對起始時刻SAR衛(wèi)星視線向LOS的累計形變量;φtopo為地形殘差相位;φatm為大氣延遲相位;φnoise為噪聲相位。在估計并移除上述三個相位分量后,上式簡化為:
(3)
為了獲取時序形變信息,上式中的相位變化與時間成線性關(guān)系:
(4)
相應(yīng)的,第j幅干涉圖的相位值可以表示為:
(5)
進一步,上式用矩陣形式表示為:
Bv=δφ
(6)
式中,B矩陣大小為M×N,由于SBAS為多主影像策略,矩陣B容易存在秩虧,利用奇異值分解方法可以求出最小范數(shù)意義上的最小二乘解[4]。
對升降軌數(shù)據(jù)進行融合分解可以獲取地表水平東西向和垂直向形變分布信息[11]。如圖3所示,展示了SAR側(cè)視成像幾何及三維坐標(biāo)下的角度參數(shù)。

圖3 SAR側(cè)視成像幾何及角度參數(shù)
InSAR干涉圖中某個像素的視線向LOS的形變實際上是該像素點三維位移分量在視線向的投影之和,如下:
rLOS=uxcosφsinθ-uysinφsinθ-uzcosθ
(7)
式中,γLOS為視線向形變量;θ為雷達信號入射角,即視線向與垂直方向的夾角;φ為衛(wèi)星航向角。地面的位移矢量u=[ux,uy,uz],ux為沿東西方向的形變分量,uy為沿南北方向的形變分量,uz為沿垂直方向的形變分量,可將視線向的沉降量r,由式(7)分解得
(8)

忽略南北方向分量,即令uy為零,且考慮升降軌兩個視線向形變則列出以下兩個方程,如下:
(9)
將式(9)列為矩陣形式:
r=Pu′
(10)

u′=[PTP]-1PTr
(11)
得到ux與uz,即東西方向與垂直方向的形變。
經(jīng)過SBAS解算后得到了研究區(qū)升、降兩軌數(shù)據(jù)的LOS地表形變速率(見圖4)。其中,形變速率的正值和負值分別表示朝向衛(wèi)星方向移動和遠離衛(wèi)星方向移動。參考點位東經(jīng)116.48°、北緯39.32°,遠離地面沉降中心(見圖4中黑色倒三角)??梢钥闯錾壓徒弟壍玫降某两悼臻g分布基本一致,證明了InSAR方法提取形變速率的可行性。從圖4中可以看出,研究區(qū)存在兩個明顯的沉降區(qū)域,如圖4中橢圓實線包含區(qū)域所示,分別位于天津武清區(qū)-廊坊霸州市(圖4右側(cè)橢圓)和廊坊固安縣-保定雄縣(圖4左側(cè)橢圓),此外研究區(qū)中存在部分小局部區(qū)域的沉降,分別位于文安縣中部和廣陽區(qū)廊坊市中心。天津武清區(qū)-廊坊霸州市的沉降影響范圍最大,沉降分布差異最明顯,LOS向平均形變速率為-76.08~-1.1mm/a(表1)。廊坊固安縣-保定雄縣的沉降影響范圍較小,LOS向平均形變速率范圍在-59.6~2.6mm/a(表2)。相比于大范圍沉降中心,研究區(qū)小局部區(qū)域沉降速率普遍較小,其中文安縣內(nèi)最大形變速率為-30mm/a。

表1 天津武清區(qū)-廊坊霸州市特征點LOS向速率

表2 廊坊固安縣-保定雄縣特征點LOS向速率

圖4 LOS向平均形變速率
InSAR獲得的形變是地表三維形變的LOS向的一維投影,僅以LOS向作為觀測結(jié)果難以反映真實的地面沉降漏斗。多角度InSAR觀測數(shù)據(jù)融合通過組合來自兩個或者多個獨立InSAR觀測幾何的LOS測量值,可以分離垂直和水平分量。由于升、降軌數(shù)據(jù)提取的像素點不一致,分離前需對升降兩軌數(shù)據(jù)進行重采樣,以使二者具有相同的規(guī)則網(wǎng)格。LOS 測量值隨入射角變化很大。詳見2.2節(jié)中所述,當(dāng)舍棄式 (8) 中的uy,即忽略LOS觀測中的南北分量時,估計的東西方向和垂直方向分量會更加準(zhǔn)確。這僅適用于南北分量與東西 分量具有相同數(shù)量級的地表變形現(xiàn)象。如果南北分量顯著大于東西分量時(如地震和滑坡變形的情況下),則不應(yīng)忽略南北分量。本文研究區(qū)屬于典型地下水開采引起的形變,以垂直沉降為主,水平變形一般具有相似量級,且朝著沉降中心移動。因此,在融合升、降軌數(shù)據(jù)時,本文忽略南北方向形變,僅分析垂直向和東西向形變。圖5展示了研究區(qū)垂直方向和東西方向的形變速率,垂直向沉降速率最大達到-100.2mm/a(圖5(a)),東西向形變速率最大為-29.6mm/a(圖5(b)),所以區(qū)域內(nèi)的形變原因主要以垂直方向為主。廊坊固安縣-保定雄縣沉降區(qū)(圖5左側(cè)橢圓)沒有明顯的水平形變,而天津武清區(qū)-廊坊霸州市沉降區(qū)(圖5右側(cè)橢圓)水平變形較為顯著,整體上向沉降中心移動。在圖5中,對于東西向水平形變,負值表示向西移動,正值表示向東移動。

圖5 研究區(qū)垂直方向和東西水平方向平均形變速率
武清區(qū)和霸州市由于開采大量深層地下水用于工業(yè)和農(nóng)業(yè)發(fā)展,早在1992—2000年就出現(xiàn)了沉降現(xiàn)象,沉降速率為30~60mm/a[12]。之后在2003—2004年兩個區(qū)域的沉降速率進一步增大,其中武清區(qū)王慶坨鎮(zhèn)和霸州市勝芳鎮(zhèn)為兩個區(qū)域的沉降中心,霸州市勝芳鎮(zhèn)的沉降速率一度高達192mm/a。2003—2010年和2012—2014年武清區(qū)和霸州市的沉降影響范圍分別以王慶坨鎮(zhèn)和勝芳鎮(zhèn)為中心進一步擴大,其中武清區(qū)王慶坨鎮(zhèn)沉降速率從100mm/a增大為153mm/a[12]。如圖6所示,2017—2020年,武清區(qū)和霸州市沉降區(qū)不僅連接在一起,沉降范圍還進一步擴大至安次區(qū)南部、文安縣東部、大成縣北部、北辰區(qū)西部和靜海區(qū)西北部,沉降影響面積達到約1470km2。天津武清區(qū)-廊坊霸州市主要包含兩個顯著的沉降中心,分別位于王慶坨鎮(zhèn)在內(nèi)的整個武清區(qū)南部和北辰區(qū)西部地區(qū)以及勝芳鎮(zhèn)在內(nèi)的霸州市中部地區(qū)。
為量化展示霸州市勝芳鎮(zhèn)范圍內(nèi)的累積沉降量,選取橫跨沉降中心的一條剖面線A-B(圖6(a)中黑色直線),繪制了剖面線上地面沉降隨時間的變化累積(圖7)。A-B剖面上最大累積沉降量為187.2mm。
將2017—2020年以年為單位分別做累積沉降分析,如圖8所示,發(fā)現(xiàn)年度累積沉降量在逐年減小,2017年累積最大沉降量89mm,2018年累積最大沉降量67mm,2019年累積最大沉降量64mm,2020年累積最大沉降量65mm,2018—2020年穩(wěn)定控制在65mm/a。其中由于2020年受新冠疫情和降水偏豐影響,地下水供水量相比2019年減少0.9×108m3,降幅超23%,為2016年以來地下水供水量最大降幅[13],因此2020年累積沉降量略有所減緩。從側(cè)面反映出2014年以來政府陸續(xù)出臺的各項地下水限采、禁采政策以及南水北調(diào)工程的建成通水使研究區(qū)深層地下水明顯回升,有效緩解了地面沉降。蔣喆等人通過對廊坊市沉降時空分布特征進行分析和提取監(jiān)測,得到的平均沉降形變速率區(qū)間是-17.9~8.5mm/a,在沉降空間分布和幅度上與本文的結(jié)果基本一致[5]。

圖8 天津武清區(qū)-廊坊霸州市沉降區(qū)2017-2020年逐年累積沉降量
此外,提取和了沉降區(qū)內(nèi)6個點(圖6中P1~P6點)的沉降時間序列,繪制于圖9,可以看出這6個點均呈線性沉降趨勢,最大沉降速率位于P6點,速率達到68.5mm/a。P1、P2點位于沉降區(qū)邊緣,其累積沉降值最小,分別為195mm和200mm。P3、P4、P5、P6點累積沉降量分別為210mm、215mm、249mm、274mm。其中P2點在勝芳汽車站附近,P4點在勝芳火車站附近,P5點靠近勝芳鎮(zhèn)政府,三者均在一定程度上受該沉降漏斗的影響。

固安縣是華北地?zé)豳Y源開發(fā)利用條件較好的地區(qū)。地?zé)豳Y源在上世紀(jì)七十年代石油地質(zhì)勘查中被發(fā)現(xiàn),有利的開采地段主要分布在牛駝鎮(zhèn)地區(qū),面積約120平方公里內(nèi)。雄縣地?zé)豳Y源豐富,境內(nèi)60%以上的區(qū)域儲藏著優(yōu)質(zhì)溫泉,基巖熱儲面積320平方千米,占牛駝鎮(zhèn)地?zé)崽锟偯娣e的50%。地?zé)豳Y源的開采成為兩地交界處資源開采發(fā)展而引起的地面沉降的首要因素。此外,西侯留工業(yè)區(qū)、趙崗工業(yè)區(qū)、仇小王工業(yè)區(qū)、陳臺工業(yè)區(qū)、米東工業(yè)區(qū)在內(nèi)的17個工業(yè)園區(qū)在距離沉降區(qū)20km范圍內(nèi),工業(yè)區(qū)的大量抽取水資源成為引起地面沉降的另一大因素。
如圖10所示,該沉降區(qū)域東起固安縣羊駝鎮(zhèn),貫穿固安縣馬莊鎮(zhèn)直達保定雄縣金三角經(jīng)濟開發(fā)區(qū)內(nèi)的板西村小陽村,西至大營鎮(zhèn)全長29.6km。水平方向變形不顯著。

圖10 廊坊固安-保定雄縣沉降區(qū)垂直沉降和水平變形
本文獲取了廊坊范圍內(nèi)的三個特征沉降點Q1、Q2、Q3(點的位置如圖10(a)所示),繪制了它們的沉降時間序列,如圖11所示。其中Q1和Q3在沉降中心邊緣,累積沉降值Q1點最小(175.4mm),Q2點最大 (206.8mm)。沉降范圍內(nèi)有眾多交通要道,榮烏高速、大廣高速、霸州立交橋、京環(huán)線等交通線路都穿過沉降區(qū)。

圖11 廊坊固安-保定雄縣沉降區(qū)Q1~Q3點時序形變圖
自然或人為地下水補給和排放造成的水頭變化會引起季節(jié)性地表沉降和抬升,本文通過對研究區(qū)沉降產(chǎn)生原因分析主要以人為因素為主,忽略了自然因素。廊坊市地表水分布稀少,淺層地下水埋藏淺,主要接受大氣降水補給,其次為側(cè)向徑流補給、河渠滲流補給、地表水灌溉和井灌回歸補給[8]。以人工開采消耗為主,其次為蒸發(fā)及側(cè)向排泄。深層地下水主要為側(cè)向徑流補給和少量越流補給,主要是人工開采產(chǎn)生的消耗。如圖12所示,廊坊市降水季節(jié)分布不均,主要集中在夏季,6~8月三個月降水量一般可達全年總降水量的70%~80%。研究區(qū)地表水屬海河流域系統(tǒng)的永定河水系,主要河流有永定河、龍河、鳳河等,均為季節(jié)性河流。如圖12所示,通過對廊坊固安-保定雄縣區(qū)域特征點時序形變?nèi)ゾ€性化分析發(fā)現(xiàn)沉降存在周期性波動,周期1年。相對于降雨期延后三個月,周期性抬升在該年11月至次年1月,周期性下降在次年1月至4月,研究區(qū)內(nèi)地下水漏斗區(qū)長系列動態(tài)變化為年際間有小的起伏,但總的趨勢成波狀形下降。

圖12 2017—2020年廊坊固安-保定雄縣形變?nèi)ゾ€性化和逐月降水量
本文主要獲取2017—2020年覆蓋研究區(qū)范圍的Sentinel-1A數(shù)據(jù),即使有升降軌數(shù)據(jù)互相驗證,但缺少地面監(jiān)測數(shù)據(jù)進行對比分析。后續(xù)還可獲取GPS監(jiān)測數(shù)據(jù)進而對比分析InSAR獲取地表形變的可靠性。
由于廊坊市城區(qū)水源以深層地下水為主,本研究僅簡單分析地表水流向,缺少地下水位相關(guān)數(shù)據(jù),僅僅依賴地表水和大氣降水無法完整的解釋地表沉降產(chǎn)生原因,故將地下水概況加入沉降形成原因也必不可少。
地下水過度開采只是城市地面沉降的主要因素,該研究希望將地面荷載、城市地下建設(shè)等因素也考慮進去,以增加形變成因的分析研究。
(1) 廊坊市地面不均勻沉降明顯,存在天津武清-廊坊霸州市和廊坊固安縣—保定雄縣2 個明顯沉降區(qū)。其中,天津武清—廊坊霸州市是最嚴(yán)重的沉降區(qū),最大沉降速率達到了 100.2mm/a。此外,研究區(qū)內(nèi)散布著許多小型沉降區(qū)。
(2) 基于升降軌融合技術(shù),對LOS向平均沉降速率分別提取垂直向和東西向沉降速率,研究區(qū)內(nèi)的形變主要來源于垂直向沉降,印證了地下水過度開采,造成深層水位下降,導(dǎo)致地面沉降增加。研究區(qū)內(nèi)沉降漏斗的出現(xiàn)伴隨著工業(yè)區(qū)而出現(xiàn),沉降漏斗的空間分布與工業(yè)區(qū)和城市建設(shè)區(qū)存在較高的相關(guān)性,而沉降區(qū)周圍交通發(fā)達,榮烏高速、大廣高速、霸州立交橋以及勝芳火車站等交通相關(guān)部門應(yīng)引起重視。
(3) 通過對研究區(qū)內(nèi)兩個明顯沉降區(qū)以及年降水?dāng)?shù)據(jù)進行時序分析,廊坊固安縣—保定雄縣沉降區(qū)的沉降發(fā)育具有1年周期的季節(jié)性變動。