毛嘉騏,李素敏*,崔文東,廖園歡,成 睿

(1.昆明理工大學(xué) 國土資源工程學(xué)院,云南 昆明 650093;2.甘肅交建項目管理有限公司,甘肅 蘭州 730000)

滑坡具備高突發(fā)性、高破壞性的特點,給人民群眾的生命安全以及生產(chǎn)建設(shè)活動造成了極大的危害,嚴重制約著社會經(jīng)濟的發(fā)展[1]。利用傳統(tǒng)的監(jiān)測方法很難對滑坡進行大面積的監(jiān)測預(yù)警,而且監(jiān)測成本較高,效率低,易受到外界因素的影響,難以實現(xiàn)滑坡預(yù)警[2]。合成孔徑雷達干涉測量(interferometric synthetic aperture radar,InSAR),具有全天時、全天候和覆蓋范圍廣等特點,已經(jīng)成功運用于滑坡災(zāi)害的監(jiān)測中,特別是在滑坡災(zāi)害的早期隱患識別中得到了廣泛的應(yīng)用[3]。戴可人等[4]利用時序InSAR技術(shù),對雅礱江中段地區(qū)進行高山峽谷地區(qū)的早期識別工作,成功探測出8處具有滑坡隱患區(qū)域。蔡杰華等[5]基于InSAR技術(shù),對九寨溝震后滑坡隱患區(qū)域進行早期的探測識別,共發(fā)現(xiàn)7處靠近居民區(qū)的隱患區(qū)域,潛在威脅性較大。ZHU等[6]利用InSAR技術(shù)對青海省劍溝滑坡進行探測,發(fā)現(xiàn)InSAR技術(shù)對分析滑坡前的蠕變變形具備優(yōu)勢。張毅[7]利用InSAR技術(shù)識別出白龍江流域內(nèi)的133處活動斜坡。在山區(qū)地形中,現(xiàn)有研究可以從一定程度上識別出潛在滑坡區(qū)域,卻無法準確評估潛在滑坡區(qū)域的穩(wěn)定性,同時,單一的潛在滑坡隱患識別也難以具體反映滑坡體失穩(wěn)特征。

FLAC3D軟件中的有限元強度折減法,可以通過彈塑性計算得出坡體的不穩(wěn)定面,并得到相應(yīng)的安全系數(shù),無需事先假設(shè),且考慮土體的本構(gòu)關(guān)系及巖土體本身的變形,可以對坡體的穩(wěn)定性以及在坡體發(fā)生破壞前后的位移情況進行判定[8]。如蔡捷等[9]利用FLAC3D中的強度折減法,對普陽煤礦1號排土場邊坡失穩(wěn)進行模擬分析,預(yù)測可能出現(xiàn)牽引式滑坡。劉宏等[10]針對緩傾角層狀巖質(zhì)邊坡小危巖體失穩(wěn)破壞模式進行分析,指出緩傾內(nèi)層狀巖質(zhì)邊坡分別為巖層較厚、巖層較薄、巖層傾向坡外等幾種情況時對應(yīng)的破壞模式,同時對破壞模式提出了相應(yīng)的穩(wěn)定性評價理論及方法。余忠祥等[8]運用FLAC3D對邊坡失穩(wěn)破壞進行模擬,同時采用強度折減法得出邊坡在不同條件下的相應(yīng)安全系數(shù),以及影響坡體穩(wěn)定性的主要因素為坡高和坡角。FLAC3D有限元強度折減法可以對坡體穩(wěn)定性實現(xiàn)量化,同時反映坡體失穩(wěn)前后的位移形變情況,提高時序InSAR監(jiān)測結(jié)果的可靠性。

因此,本文采用時序InSAR技術(shù)獲取天水市中部地區(qū)地表形變數(shù)據(jù),識別出潛在滑坡區(qū)域;同時結(jié)合FLAC3D軟件,運用強度折減法得出隱患坡體的安全系數(shù),分析識別潛在滑坡區(qū)域坡體的穩(wěn)定性情況,為滑坡的預(yù)測和坡體的穩(wěn)定性評價提供了有力的支持,并為當?shù)氐幕骂A(yù)警監(jiān)測提供了理論支撐。

1 研究區(qū)概況及數(shù)據(jù)處理

1.1 研究區(qū)概況

天水市地處甘肅省東南部,位于104°35′～106°44′E、34°05′～35°10′N的青藏高原和黃土高原交匯地帶,地勢西北高東南低,平均海拔高度為1 100 m,境內(nèi)有渭河及其支流流經(jīng),受地質(zhì)沉陷和紅、黃土層沉積,形成黃土丘陵地貌。“蘭州—天水”地震帶穿城而過,區(qū)域內(nèi)斷裂褶皺發(fā)育,地質(zhì)構(gòu)造較為活躍。天水降雨主要集中在夏秋兩季,且呈現(xiàn)出集中性和持續(xù)性的特點,加之近年來工程活動較多,導(dǎo)致滑坡、泥石流等自然災(zāi)害頻發(fā)[11]。本文主要研究區(qū)為天水市中部地區(qū),研究區(qū)位置如圖1所示。

圖1 研究區(qū)域Fig.1 Research area

1.2 SBAS-InSAR數(shù)據(jù)處理方法

SBAS-InSAR是利用同一地區(qū)多幅時間基線較短的SAR影像形成的干涉對。干涉對通過解纏、濾波等方式去除軌道誤差、噪音以及地形的殘余相位,保證時間基線的高相干性,然后采用奇異值分解(singular value decomposition, SVD)的方法,將多個基線集聯(lián)合求解,并對時間域和空間域的濾波進行分析,分離出殘余相位中的大氣相位和非線性形變誤差,得到目標區(qū)域內(nèi)覆蓋整個觀測時間的地表形變信息[12]。本文選取IW模式下的2018年10月—2021年10月的37景Sentinel-1A升軌數(shù)據(jù),數(shù)據(jù)的相關(guān)參數(shù)如表1所示。為提高SAR衛(wèi)星影像的軌道精度,引入由歐空局提供的精密軌道數(shù)據(jù),使用30 m分辨率的DEM,用于計算過程中去除地形相位的影響。

表1 Sentinel-1A數(shù)據(jù)參數(shù)Tab.1 Sentinel-1A data parameters

1.3 FLAC3D有限元強度折減法

本文采用有限元強度折減法,運用FLAC3D對坡體進行強度折減計算。其原理是將坡體的巖土體的黏聚力C和內(nèi)摩擦角φ進行調(diào)整,得到一組新的黏聚力和內(nèi)摩擦角,然后利用折減后的參數(shù)進行計算,直到坡體達到極限狀態(tài),此時的折減系數(shù)F即為坡體穩(wěn)定時的安全系數(shù)[13]。

(1)

(2)

式中:C′為折減后的黏聚力;φ′為折減后的摩擦角;F為折減系數(shù)。

通過現(xiàn)場調(diào)查,計算采用摩爾庫倫本構(gòu),巖土體抗剪強度、彈性模量、泊松比等物理力學(xué)參數(shù)均通過室內(nèi)實驗獲取,如表2所示。

表2 巖土體物理力學(xué)參數(shù)Tab.2 Physical mechanics parameters of rock and soil mass

1.4 技術(shù)路線

利用時序InSAR技術(shù)進行滑坡災(zāi)害隱患點的大規(guī)模普查,查找隱患點位置,然后對普查結(jié)果用FLAC3D有限元軟件生成網(wǎng)格單元,并建立坡體模型;對坡體賦予相應(yīng)的條件和參數(shù),獲得達到初始平衡狀態(tài)的初始模型;對模型施加一定的外力,破壞其現(xiàn)有狀態(tài),繼續(xù)進行計算,使模型處于平衡或者破壞狀態(tài);最后結(jié)合坡體模型得出安全系數(shù),對坡體安全性進行綜合判斷[14]。整體技術(shù)流程如圖2所示。

圖2 滑坡災(zāi)害隱患點識別技術(shù)流程圖Fig.2 Flow chart of landslide disaster hidden danger point identification technology

黑色圓圈為沉降速率較大區(qū)域;紅色方框為城區(qū)。圖3 2018—2021年均形變速率圖Fig.3 The average variable rate map of 2018-2021

2 InSAR滑坡隱患早期識別結(jié)果

2.1 總體識別結(jié)果

基于時序InSAR獲取天水市中部地區(qū)2018年10月—2021年10月的形變數(shù)據(jù),其年均形變速率圖如圖3所示。從圖3可以看出:整個研究區(qū)形變情況分布不均勻,位于中部的主城區(qū)附近較為穩(wěn)定,形變速率為10 mm/a,未見有較大的形變情況發(fā)生;形變嚴重的區(qū)域主要集中在研究區(qū)北部,區(qū)域內(nèi)地形起伏較大,年均形變速率變化范圍為-92～68 mm/a;研究區(qū)南部地形也有較大起伏,形變速率在-40～35 mm/a之間。結(jié)合研究區(qū)光學(xué)遙感影像及實地調(diào)查,發(fā)現(xiàn)研究區(qū)南部植被覆蓋率較高,會導(dǎo)致嚴重的失相干情況;研究區(qū)北部植被覆蓋率較低,形變監(jiān)測結(jié)果良好,故在研究區(qū)北部區(qū)域篩選滑坡隱患區(qū)域。

結(jié)合坡度(圖4)、高程(圖1)、形變速率等分析,位于研究區(qū)北部的秦州區(qū)太京鎮(zhèn)丁家門村和麥積區(qū)能干村的坡體存在隱患,其最大形變速率為92 mm/a,坡度均大于35°,且所處位置海拔較高。其中,位于丁家門村坡體靠近斷裂帶,能干村坡體靠近渭河、隴海鐵路,存在造成次生災(zāi)害的風(fēng)險。

圖4 研究區(qū)坡度Fig.4 Slope of research area

2.2 重點區(qū)域解譯

第一處滑坡隱患較大的區(qū)域位于秦州區(qū)太京鎮(zhèn)丁家門村附近坡體。該處坡體最大沉降值為-120 mm,從上至下選取4個點位進行時序分析,如圖5所示。從圖5可以看出:坡體整體沉降速率較大,最大沉降速率為44 mm/a;2018年10月—2020年4月,坡體變化趨勢較為穩(wěn)定,之后以較大的速率持續(xù)沉降。巖土體裸露情況較為嚴重,植被種類單一且稀少,周邊有部分耕地,在強降雨及持續(xù)性集中降雨條件下,坡體易受侵蝕,存在較高滑坡成災(zāi)風(fēng)險。

第二處隱患區(qū)位于麥積區(qū)能干村,坡體周邊植被較少,計算結(jié)果整體相干性較好。該坡體最大沉降值為-130 mm,坡體的中部和下部都保持著持續(xù)且較大沉降速率,取該坡體4個點位進行時序監(jiān)測,如圖6所示。從圖6可以看出:自監(jiān)測之日起,該坡體除頂部點位沉降速率較為平緩?fù)?其余監(jiān)測點位持續(xù)保持一定速率沉降,最大沉降速率為54 mm/a。坡體頂部附近為農(nóng)田,周邊溝壑落差較大,坡體側(cè)臂較陡,毗鄰渭河、隴海鐵路,如若發(fā)生滑坡,易造成次生災(zāi)害的發(fā)生,影響渭河流域及隴海鐵路的安全性。

紅色框線內(nèi)為農(nóng)田;黑色線內(nèi)為坡體邊界;白色框內(nèi)1～4為選取的時序點位。圖5 丁家門村2018—2021年均沉降速率圖Fig.5 Average settlement rate map of Dingjiamen Village from 2018 to 2021

紅色框線內(nèi)為農(nóng)田;黑色線內(nèi)為坡體邊界;白色框內(nèi)1～4為選取的時序點位。圖6 能干村2018—2021年均沉降速率圖Fig.6 Average settlement rate map of Nenggan Village from 2018 to 2021

3 基于FLAC3D強度折減法的模擬分析

3.1 模擬初始位移分析

運用FLAC3D對丁家門村坡體進行計算,得到位移云圖如圖7所示。由圖7可以看出:坡體中部在水平方向上沿X軸正方向位移較大,在靠近底部處達到最大,位移值為9 mm;在豎直方向上,坡體中部沿Z軸負方向位移情況明顯,最大位移值為-7 mm。

能干村坡體的整體位移情況較為明顯,位移云圖如圖8所示。由圖8可以看出:在水平方向上,坡體中下部沿X軸正方向位移情況突出,最大位移值為5 mm;在豎直方向上,整體位移較大,在坡體中上部尤為顯著,最大位移值為-5 mm。

(a)水平方向 (b)豎直方向圖7 丁家門村坡體位移云圖Fig.7 Cloud map of Dingjiamen Village slope displacement

(a)水平方向 (b)豎直方向圖8 能干村坡體位移云圖Fig.8 Cloud map of Nenggan Village slope displacement

3.2 模擬安全系數(shù)分析

運用強度折減法對邊坡進行折減計算,對坡體在不同折減系數(shù)下的最大豎直方向位移變化情況進行分析,從而確定坡體的最終安全系數(shù)。位于秦州區(qū)丁家門村的坡體在不同折減系數(shù)下位移變化情況如圖9(a)所示。由圖9(a)可知:當折減系數(shù)為0.900～1.000時,坡體最大位移變化平穩(wěn);當折減系數(shù)大于1.003時,最大位移出現(xiàn)陡增,坡體失穩(wěn)。因此,坡體安全系數(shù)為1.003,小于坡體失穩(wěn)時安全系數(shù)需大于1.1[14]的要求,且坡體最大位移值為231 mm,坡體形變量較大,發(fā)生位移突變,坡體處于不穩(wěn)定狀態(tài),存在發(fā)生滑坡的風(fēng)險。

圖9 不同折減系數(shù)下的位移變化Fig.9 Reduction coefficient diagram of slope numerical simulation

位于麥積區(qū)能干村的坡體在不同折減系數(shù)下位移變化情況如圖9(b)所示。由圖9(b)可知:坡體以0.80為起始折減值,折減系數(shù)為0.800～1.040時,最大位移處于平穩(wěn)上升狀態(tài),之后隨著折減系數(shù)不斷增大,位移值攀升。因此,坡體安全系數(shù)為1.040,小于1.1,最大位移值為135 mm,之后坡體失穩(wěn),坡體屬于不穩(wěn)定坡體,存在發(fā)生滑坡的風(fēng)險。

4 結(jié)果對比

將SBAS-InSAR計算的沉降結(jié)果與數(shù)值模擬的位移結(jié)果進行對比,如圖10所示。圖10左圖顯示,丁家門村坡體沉降值為-120～-90 mm,能干村坡體沉降值為-130～-40 mm;圖10右顯示,丁家門村坡體最大合位移值為-105 mm,能干村坡體最大合位移值為-137 mm。由圖10可知:通過時序InSAR技術(shù)測量得到的沉降量較大區(qū)域與數(shù)值模擬得出的坡體合位移較大區(qū)域基本一致,坡體存在較大形變情況。丁家門村坡體安全系數(shù)為1.003,小于1.1,屬于潛在滑坡易發(fā)區(qū);坡體中部形變速率較大,與數(shù)值模擬得出的初始位移情況區(qū)域一致。能干村坡體最終安全系數(shù)為1.040,屬于潛在滑坡易發(fā)區(qū);坡體中部及下部形變速率較大,與數(shù)值模擬得出初始模型存在位移隱患區(qū)域基本一致,且由于該處坡體靠近渭河及隴海線,存在安全隱患。

圖10 SBAS-InSAR計算沉降結(jié)果與數(shù)值模擬位移結(jié)果對比Fig.10 SBAS-INSAR calculation settlement results and numerical simulation displacement results comparison

5 結(jié)論

本文使用Sentinel-1A數(shù)據(jù),基于時序InSAR技術(shù)對天水市中心的山區(qū)地帶滑坡隱患開展早期的探測識別工作,共識別出2處隱患較大的典型區(qū)域作為重點進行探究。得出如下結(jié)論:

1)建立數(shù)值模擬的計算模型,對模型材料的抗剪強度進行折減,得出坡體水平和豎直方向位移變化和安全系數(shù),其中,丁家門村與能干村坡體安全系數(shù)分別為1.003和1.040。綜合研判確定丁家門村坡體與能干村坡體均屬于“潛在滑坡易發(fā)區(qū)”。

2)丁家門村與能干村坡體形變區(qū)域的沉降值較大區(qū)域與數(shù)值模擬得出的相對位移較大區(qū)域基本一致。因此,利用時序InSAR技術(shù)結(jié)合FLAC3D有限元軟件,不僅可以兼顧快速普查滑坡隱患區(qū)域的需求,還能實現(xiàn)對于重點區(qū)域的精準查詢,從而對隱患區(qū)域的安全性進行判斷,為滑坡隱患識別、滑坡預(yù)警等提供數(shù)據(jù)與技術(shù)支持。