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

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

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

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

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

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

1.1 研究區(qū)概況

天水市地處甘肅省東南部,位于104°35′～106°44′E、34°05′～35°10′N的青藏高原和黃土高原交匯地帶,地勢(shì)西北高東南低,平均海拔高度為1 100 m,境內(nèi)有渭河及其支流流經(jīng),受地質(zhì)沉陷和紅、黃土層沉積,形成黃土丘陵地貌?！疤m州—天水”地震帶穿城而過,區(qū)域內(nèi)斷裂褶皺發(fā)育,地質(zhì)構(gòu)造較為活躍。天水降雨主要集中在夏秋兩季,且呈現(xiàn)出集中性和持續(xù)性的特點(diǎn),加之近年來工程活動(dòng)較多,導(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ū)多幅時(shí)間基線較短的SAR影像形成的干涉對(duì)。干涉對(duì)通過解纏、濾波等方式去除軌道誤差、噪音以及地形的殘余相位,保證時(shí)間基線的高相干性,然后采用奇異值分解(singular value decomposition, SVD)的方法,將多個(gè)基線集聯(lián)合求解,并對(duì)時(shí)間域和空間域的濾波進(jìn)行分析,分離出殘余相位中的大氣相位和非線性形變誤差,得到目標(biāo)區(qū)域內(nèi)覆蓋整個(gè)觀測(cè)時(shí)間的地表形變信息[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,用于計(jì)算過程中去除地形相位的影響。

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

1.3 FLAC3D有限元強(qiáng)度折減法

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

(1)

(2)

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

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

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

1.4 技術(shù)路線

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

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

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

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

2.1 總體識(shí)別結(jié)果

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

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

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

2.2 重點(diǎn)區(qū)域解譯

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

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

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

3 基于FLAC3D強(qiáng)度折減法的模擬分析

3.1 模擬初始位移分析

運(yùn)用FLAC3D對(duì)丁家門村坡體進(jìn)行計(jì)算,得到位移云圖如圖7所示。由圖7可以看出:坡體中部在水平方向上沿X軸正方向位移較大,在靠近底部處達(dá)到最大,位移值為9 mm;在豎直方向上,坡體中部沿Z軸負(fù)方向位移情況明顯,最大位移值為-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ù)分析

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

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

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

4 結(jié)果對(duì)比

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

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

5 結(jié)論

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

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

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