李宗仁, 沙永蓮, 辛榮芳*, 張興, 隋嘉, 孫皓
(1.青海省地質調查院, 西寧 810012; 2.青海省遙感大數據工程技術研究中心, 西寧 810012; 3.青藏高原北部地質過程與礦產資源重點實驗室, 西寧 810012; 4.青海省地質環境監測總站, 西寧 810001)
青海省地質災害主要發育在山間谷地,具有點多面廣的特點[1],面對嚴峻復雜的防災形勢,省內先后開展了系統全面的地質災害詳查和多輪次拉網式隱患排查,每年汛期不斷巡查,并在此基礎上建立了較為完善的群測群防體系。但近年來發生的很多地質災害事件中有部分并不在已知隱患范圍內,對于隱蔽性強、變形緩慢的隱患,傳統的人工排查手段已顯無力。當前防范地質災害的核心需求是要搞清楚“隱患點在哪里”“什么時候發生”,為此,如何突破傳統人工調查的局限,快速、有效地識別地質災害隱患是青海省防災工作急需解決的問題[2-3]。
目前,青海省地災隱患主要靠人工地面摸排進行識別,這僅能從點上發現明顯的災害征兆,如地裂縫、鼓包、下錯以及建筑物變形等跡象,面對大范圍的整體位移且無明顯變形跡象時難以察覺。近年來雖然開展了重點區域的光學遙感調查,實現了地質災害的影像識別與變化監測,但也只能發現已發生或正在發生的災害體形態變化,對于緩慢變形的隱患識別能力仍有不足。合成孔徑雷達干涉測量(synthetic aperture radar interferometry,InSAR)技術具有全天候、全天時、覆蓋范圍廣、空間分辨率高、非接觸、綜合成本低等優點,適宜于開展大范圍地質災害普查與長期持續觀測,特別是其大范圍連續跟蹤微小形變的特性,對正在變形區具有獨特的識別能力[4],但是,青海省地質災害監管部門還未自行開展過地質災害InSAR早期識別業務,導致防災工作比較被動。現通過長時序SAR數據對同仁市進行持續形變觀測,配合孕災背景資料及高分辨率光學影像,識別圈定滑坡災害隱患,并對隱患點形變特征進行時序分析,判斷變形階段,評估其危險性與風險,旨在為青海省地質災害防治工作提供技術支撐。
同仁市屬青海省黃南藏族自治州管轄,位于青海省東南部、黃南州東北部,市域東西寬75 km,南北長85 km,總面積3 275 km2,地理坐標為東經101°38′~102°27′,北緯35°01′~35°47′。研究區地處青藏高原與黃土高原相接的過渡地帶,境內山巒起伏,河谷相間,地勢南高北低,高差懸殊,隆務河縱貫全境南北,形成了東-西部山區和中部河谷的兩峽-谷地地形特征(圖1)。該區屬典型的高原溫帶半干旱氣候區,具有寒冷干燥,日照充足,降水集中等特點,年內降水分配極不均勻,多集中在5—9月,且降水量隨地勢增高垂直分帶性明顯。區內已造就的現代地形其發展演變過程是由地球內、外力共同作用的結果——地形破碎、坡陡溝深、溝壑縱橫、地表裸露,地質環境較為脆弱,引發的環境地質問題主要有滑坡、崩塌、泥石流等地質現象以及修建道路引起的邊坡失穩等[5]。

圖1 研究區地理位置及地貌特征圖Fig.1 Map of geographical location and geomorphic characteristics of the study area
地質災害隱患早期識別是地質災害防治工作的重要組成部分,相較于傳統地面調查與光學遙感解譯,時序InSAR技術作為一種新型的對地觀測技術,不僅可以識別連續、緩慢的地表形變信息,同時還可以提升災害隱患的識別效率[6]。研究區由于其特殊的地理環境,地質災害發育具有數量多、類型全、規模大的特點,現針對區內的滑坡及潛在滑坡,采用基于Sentinel-1A升降軌數據的差分干涉測量短基線集時序分析技術(small baselines subset InSAR, SBAS-InSAR)進行地質災害隱患的早期識別,技術流程如圖2所示。
(1)數據選擇。以Sentinel-1A衛星SAR數據為主要數據源,選用2019年1月—2020年10月覆蓋研究區的52景升軌影像和45景降軌影像進行區域形變信息提取(表1、圖3)。除SAR數據外,還采用了ALOS衛星30 m分辨率的DEM數據用于消除干涉過程中的地形相位,降低地形誤差,同時還采用了精密定軌星歷數據(POD)作為軌道數據對軌道信息進行修正,消除因軌道誤差產生的系統性誤差;采用了2020年5—9月覆蓋研究區的14景空間分辨率為1 m的國產高分二號(GF-2)光學影像對隱患點進行篩選與確認(圖3)。高分二號衛星是中國自主研制的首顆空間分辨率優于1 m的民用光學遙感衛星,具有分辨率高、定位精度準、姿態機動快等特點,其亞米級空間分辨率能夠有效識別地質災害形態要素。

圖2 本次研究技術流程圖Fig.2 Technical flow chart of the study

表1 研究所用Sentinel-1A數據參數表Table 1 Main parameters of Sentinel-1A data

圖3 本次研究所用衛星數據覆蓋分布圖Fig.3 Satellite data coverage map used in this study
(2)SAR數據處理。對覆蓋研究區的總計97景升降軌數據進行SBAS-InSAR技術處理,包括連接圖生成、差分干涉、相位解纏、軌道精煉與重去平、SBAS反演、地理編碼等一系列工作流[7-8]。首先進行干涉像對的配對,對升、降軌數據分別設置最大時間基線為60 d和48 d,通過參數設置及經驗分析,剔除質量差的干涉對,最終得到有效干涉對183組和105組,影像配對情況良好;其次進行濾波與干涉處理,對升降軌數據均設置多視比為10∶2,提高數據信噪比,降低方位向與距離向分辨率之間的差異。同時采用高斯濾波(Goldstein)方法進行濾波處理,提高干涉條紋的清晰度和減少時空基線引起的失相干,從而提升差分干涉的穩定與可靠性;利用狄洛尼最小費用流(Delaunay minimum cost flow, DMCF)方法進行相位解纏處理[9-10],實現大范圍的全局最優解纏結果,處理過程中對相干性低或存在明顯大氣相位的像對進行了移除,升降軌數據分別移除干涉對18組和35組;通過選取地面控制點(GCP)實現干涉結果的軌道精煉與重去平,利用相位解纏圖和相干性圖輔助完成GCP點的選取,選擇的控制點均勻分布在高相干區域;SBAS形變反演分兩次進行,第一次反演是在忽略大氣效應的基礎上進行地表形變速率和高程改正數估算,同時通過二次解纏進行干涉圖優化。第二次反演是使用定制濾波去除低頻形變信息、高程誤差以及大氣誤差,得到最終的時間序列位移結果;最后將SBAS反演結果轉化成統一的地理坐標,實現SAR數據從斜距標坐標到地理坐標系的轉換。
(3)滑坡隱患識別與特征解析。在SBAS-InSAR反演的地表形變信息基礎上,以“青海省同仁縣1∶5萬地質災害詳細調查(2018)”成果資料中地質災害分布情況為參考,根據GF-2影像中判識的地形地貌、植被覆蓋等特征,圈定具有孕災條件且有威脅對象斜坡處的形變集中區,根據光學影像特征確定滑坡隱患的邊界范圍,經與收集的“青海省地質災害隱患數據庫(2020年)”套合比對,完成“是否新增隱患”屬性的判識。并選取典型滑坡隱患從地表形變的時空分布、光學影像特征、實地變形特征等方面進行特征分析。
基于SBAS-InSAR技術獲得了同仁市在2019年1月—2020年10月的區域地表形變信息。從升軌數據反演結果中可以看出,監測周期內大部分區域處于穩定狀態,其中遠離衛星視線方向的最大累積形變量為130 mm,靠近衛星視線方向的最大累積形變量為56 mm,平均累積形變量為5 mm,強形變區域多分布于隆務河西岸的陡峭斜坡上[圖4(a)]。從中識別出8個形變特征點A01~A08,經時序特征分析,以2019年1月9日為基準,各特征點均表現為隨著時間的變化其累積形變量不斷增加。從量級上看,監測時段內A01與A02的累積形變量最大,達105 mm。根據研究區降雨特征,在2020年汛期,各特征點形變量呈明顯加速趨勢(圖5)。
基于降軌數據的反演結果顯示,研究區在監測周期內遠離衛星視線方向的最大累積形變量為95 mm,靠近衛星視線方向的最大累積形變量為73 mm,平均累積形變量為1 mm[圖4(b)],共識別出11個形變特征點,全部特征點均呈持續變形趨勢,不同特征點的累積形變量級與形變特征存在差異,其中,D08的累積形變量最大,達85 mm;D06、D07、D10、D11的累積形變量較小,且變化相對較緩。受降雨影響,雨季時(曲線圖中橙色部分)各特征點均出現形變加速或波動現象,此外,D10和D11在2020年3月出現波動是因冰雪消融產生的影響(圖6)。

圖4 基于升降軌數據的累積形變量反演結果Fig.4 Inversion results of cumulative deformation based on ascending and descending data

圖6 降軌數據中形變特征點時序曲線圖Fig.6 Time series curve of deformation characteristic points in descending data
通過時序SBAS-InSAR技術反演的區域地表形變結果是進行滑坡隱患早期識別的重要依據[11],在形變速率、累計形變量分析基礎上,以形變特征點的形變數量級、坡體是否具備孕災條件、形變速率是否突變、是否具有威脅對象等綜合因素為參考,結合GF-2光學影像進行篩查,剔除明顯受地形或大氣影響、不滿足滑坡發育條件、長時間處于冰雪覆蓋、農田耕作區等形變特征區域(點),在區內共識別出19處滑坡隱患,集中分布于隆務河兩側的斜坡上,其中升軌數據中8處,降軌數據中11處,主要涉及黃乃亥鄉、年都乎鄉、蘭采鄉、隆務鎮、加吾鄉、曲庫乎鄉、扎毛鄉等7個鄉鎮,配合GF-2光學影像明確了各滑坡隱患坡體的位置及其邊界范圍,滑坡隱患早期識別詳情見表2,各隱患分布情況如圖7所示。
基于升軌數據識別出的8處隱患中A-02、A-03、A-04、A-05和A-08與歷史地質災害面重合,A-01、A-06和A-07為新識別隱患,降軌數據中D-02、D-05、D-09與已知災害點重合,各隱患點形變特征詳如圖8所示。
3.3.1 形變速率精度分析
經形變速率平均精度計算公式統計,升軌數據中共獲取有效觀測值16 346 546個,其中有效形變點的年均形變速率集中于-10~10 mm/a,形變速率精度范圍為0~8 mm,平均形變速率精度2 mm,標準差為4 mm;降軌數據中共獲取有效觀測值15 748 139個,其中有效形變點的年均形變速率集中于-20~20 mm/a,形變速率精度范圍為0~11 mm,平均形變速率精度為0,標準差為5 mm(圖9)。總體而言,研究區內相干性較高,基于升、降數據反演的形變速率精度可達毫米級,反演結果精度優于8~11 mm,升軌形變速率精度整體優于降軌,監測結果有效可靠。

表2 基于升降軌數據的滑坡隱患早期識別詳情表Table 2 Details of Early Identification of Landslide Hazards Based on ascending and descending data

圖8 基于升降軌數據的形變速率及隱患點形變特征圖Fig.8 Deformation feature map of deformation rate and hidden danger points based on ascending and descending data
3.3.2 結合幾何畸變的識別結果差異性分析
星載SAR數據的幾何畸變與雷達成像側視角、地形坡度坡向有著直接的聯系,受幾何畸變影響而導致的相位丟失或重疊是雷達衛星側視成像的先天不足。從所用SAR數據幾何畸變分析結果可知,畸變類型主要分為主動疊掩、被動疊掩,主動陰影、被動陰影以及疊掩陰影,主、被動根據觀測地物本身產生畸變或因周圍地物影響產生畸變而劃分。升降軌數據中幾何畸變主要發生在隆務河東西兩岸,區域內受幾何畸變影響較降軌數據更為明顯。經統計,升降軌數據主要受到疊掩畸變的影響,且升軌受影響范圍比降軌更廣(圖10)。

圖9 形變速率精度反演結果與統計圖Fig.9 Inversion results and statistical chart of deformation rate accuracy

圖10 升降軌數據幾何畸變分布及統計圖Fig.10 Geometric distortion distribution and statistical chart of ascending and descending data
為直觀反映升降軌幾何畸變與坡度坡向的關系,以黃乃亥鄉為例,升軌數據中疊掩區域主要分布于隆務河東西兩側東南朝向的斜坡上,陰影區集中分布于坡度較大的區域;降軌數據中疊掩區域主要分布在靠近隆務河西側的東北朝向斜坡上,陰影區則分布在背向雷達信號的陡峭斜坡上。此外,疊掩區主要為坡向復雜區域,陰影區則為坡度較大區域,難以獲取有效的監測信號,幾何畸變的類型及其分布在升降軌數據中表現差異較大(圖11)。
為進一步說明幾何畸變對升降軌識別結果的影響,以D-04號隱患為例,該坡體在升軌數據中受到疊掩影響,形變速率反演結果中存在明顯的信號丟失,形變速率密度及精度均較低;而在降軌數據中雖受到輕微陰影影響,但得到了有效且完整的形變速率結果,反映出了明顯的形變特征(圖12),這表明升降軌數據在形變解算時具有互補性,聯合不同入射角的SAR影像獲取不同成像幾何下目標地物的形變特征,在一定程度上彌補了單一成像方式帶來的監測盲區,識別結果更加全面,從而提高了滑坡隱患早期識別結果的準確性和可靠性。
在滑坡監測中,地表形變的時間和空間變化是評價滑坡階段的重要指標,通過SBAS-InSAR時序方法獲取到滑坡隱患不同時期的形變量級與速率參數,通過形變特征分析為滑坡監測與預警提供可借鑒的參考依據。現選取西山滑坡群、羊智滑坡群、尕沙日滑坡作為典型滑坡隱患進行特征解析。

圖11 黃乃亥鄉升降軌數據幾何畸變分布與坡度坡向關系圖Fig.11 Mapping of geometric distortion distribution and slope aspect of ascending and descending data in Huangnaihai Township

圖12 D-04號隱患升降軌數據幾何畸變與形變速率關系圖Fig.12 Relationship between geometric distortion and deformation rate of ascending and descending data at D-04 hidden danger
西山滑坡群位于同仁市隆務鎮夏瓊南路西側斜坡地帶,二郎寺西側,隆務河西岸,中心經緯度坐標為102°00′E,35°31′N,由Ⅰ~Ⅻ號12個滑坡體組成,為典型的大型土質滑坡群。該滑坡群發育于臨夏組地層上部、殘坡積層及隆務河四級階地上,以西山南段古滑坡為主,嚴重威脅坡下居民點及學校等市政設施,隱患風險等級較高[12]。
由覆蓋西山滑坡群的斜坡向形變速率分布圖可以看出,整個坡面中存在2處強形變區域,呈環形、漏斗狀形態分布于Ⅱ號和Ⅴ號滑坡體,表現出持續變形趨勢。在兩處強形變分布區沿斜坡方向(西北向東南)繪制了形變速率剖面,其中,Ⅱ號滑坡中部的形變速率最大,達-13 m/a;Ⅴ號滑坡的最大形變速率為-12 mm/a,位于坡體后緣,該滑坡后緣和中部的形變速率明顯大于前緣,穩定性較差(圖13)。經實地調查發現,Ⅴ號滑坡體后緣有3條雁形排列的拉張裂縫,中部發育6、7條裂縫,呈橫向斷續或連續發育,局部有下錯現象(圖14),這與InSAR形變監測結果相一致。
監測周期內西山滑坡群最大形變量為24 mm,最大形變中心位于Ⅴ號滑坡。在Ⅱ號和Ⅴ號滑坡分別從上部、中部、下部選取特征點進行時序形變分析,結果顯示,在2019年6月23日之前形變呈線性勻速下降趨勢,至2019年12月5日,整個坡體形變出現線性波動,表現出多次極快的小幅度下降和抬升,這與夏季降雨集中導致降雨入滲裂縫,表層土壤擾動有關,坡體穩定性變差。而在2020年1月—2020年5月18日,形變恢復平穩,呈緩慢下降,隨著6月雨季的到來,形變明顯加快,并在2020年9月之后進入新一輪的波動變化(圖15)。

圖13 西山滑坡群形變速率分布及剖面圖Fig.13 Distribution and profile of deformation rate of Xishan landslide group

圖14 Ⅴ號滑坡中部(左)和后緣(右)拉張裂縫實地照片Fig.14 Field photos of tension cracks in the middle (left) and rear edge (right) of No.V landslide

圖15 西山滑坡群累積形變量分布及特征點時序曲線圖Fig.15 Time series curve of cumulative deformation variables and characteristic points of Xishan Landslide Group
羊智滑坡群位于黃乃亥鄉東北部,地處隆務河西岸,海拔2 275~2 864 m,平均坡度42.5°,中心地理坐標為102°12′E,35°40′N,為典型土質滑坡區,所在坡體植被較少,黃土覆蓋,南側存在歷史泥石流沖溝,坡下為谷地平原,分布有居民點和農田,嚴重威脅坡頂的羊智村,若發生滑動也會對隆務河產生影響,隱患風險等級較高[13]。
由形變速率反演結果可知,該斜坡范圍內存在4個明顯的強形變中心,自北向南編號為YZ_1#、YZ_2#、YZ_3#、YZ_4#,最大形變速率為-64 mm/yr,最大形變中心位于YZ_1#,平均形變速率為-13 mm/a。選取4條剖面線沿斜坡向下方向繪制了形變速率變化曲線,結果顯示,各滑坡強形變區均位于斜坡中部或后緣位置,其中YZ_1#滑坡形變核心位于坡頂往下450 m處,形變速率達-64 mm/a,在滑坡體前緣古滑坡面范圍內形變速率較小,變化范圍在-20~-10 mm/a;YZ_2#與YZ_1#滑坡形變速率變化相似;YZ_3#滑坡形變速率相對較小,最大形變速率為-27 mm/a,形變核心位于坡體中部;YZ_4#滑坡形變核心位于坡體后緣。其中YZ_3#靠近村莊居民點,人類活動跡象明顯(圖16)。

圖16 羊智滑坡群形變速率分布及剖面圖Fig.16 Distribution and profile of deformation rate of Yangzhi landslide group
結合GF-2光學影像圈定了該滑坡群內各滑坡的邊界范圍,明確了強形變中心對應的滑坡體,各形變中心均位于滑坡體的中上部區域,坡體陡峭,植被稀疏,黃土覆蓋,為滑坡災害的發生提供了環境條件及物質來源,滑坡群區域內無斷層發育,斷裂運動影響較小。各滑坡在GF-2影像中特征明顯,可清晰識別出YZ_1#滑坡前緣及中部有次級滑坡痕跡;YZ_2#中部發育有次級滑坡;YZ_3#滑坡后緣及中部有次級滑坡,沿斜坡方向自上而下發育有歷史泥石流沖溝;YZ_4#后緣呈典型的“圈椅狀”地貌,且發育有次級小型滑坡(圖17)。
羊智滑坡群范圍內最大累積形變量達-122 mm,位于YZ_1#滑坡。整體上累積形變量的分布特征與形變速率分布具有很好的一致性。針對4個滑坡體共選取了14個特征點繪制時間序列曲線圖,其中,YZ_1#滑坡最大累積形變量接近120 mm,在監測時段內滑坡前緣形變較小,后緣則處于快速形變階段,在2020年5月之后,因降雨增多形變量呈加速遞增趨勢,而斜坡中部處于緩慢變形階段;YZ_2#滑坡整體形變特征與YZ_1#相似,形變集中于斜坡中部偏后緣,其形變量隨著時間的變化逐漸增加;YZ_3#滑坡中部累積形變量明顯大于前緣和后緣,整體形變量隨著時間變化不斷遞增;YZ_4#滑坡后緣形變明顯,在2020年之后其形變量變化趨勢明顯加快。整體而言,羊智滑坡群受夏季降雨影響存在較大隱患(圖18)。
沙日村位于同仁盆地中的河谷階地區,處于隆務河西岸,地理坐標為102°01′33″E,35°36′6″N。村莊后側所靠山體為南北走向,海拔2381~3005 m,相對高差約為620 m,其山體中下部由水平產狀的第三系紅色砂泥巖構成,尕沙日滑坡上部覆蓋有風積黃土,縱向坡度約65°。該滑坡發生于2009年9月13日,滑體前緣在河谷階地上形成碎屑流沖擊,直接損壞公路、房屋、引水渠等,并有人員傷亡,造成直接經濟損失1 256.4萬元。
尕沙日滑坡范圍內最大形變速率達-41 mm/a,形變核心集中于坡體中部及后緣位置,監測時段內后壁處于不穩定狀態。由于受到幾何畸變的影響,后壁靠近滑坡邊界位置失相干,存在監測空值。沿斜坡向下方向(G-G′)繪制的形變速率剖面圖中顯示,坡體中后部的形變速率明顯大于斜坡前緣,即斜坡體后緣處于不穩定狀態(圖19)。
經GF-2影像解譯發現,滑坡體后壁“圈椅狀”特征明顯,且發育有多處橫縱向拉裂縫,最長裂隙約300 m,坡體后壁失穩存在滑動的可能性,這與時序InSAR的形變監測結果相吻合(圖20)。

圖18 羊智滑坡群累積形變量分布及特征點時序曲線圖Fig.18 Time series curve of cumulative deformation variables and characteristic points of Yangzhi Landslide Group

圖19 尕沙日滑坡形變速率分布及剖面圖Fig.19 Distribution and profile of deformation rate of Gashari landslide

圖20 尕沙日滑坡光學影像解譯圖Fig.20 Optical image interpretation map of Gasari landslide
該滑坡在監測周期內最大累積形變量達-79 mm,整體處于蠕變狀態。沿斜坡方向選取5個特征點繪制了時序曲線,可以看出,除了靠近斜坡前緣的G4和G5兩個點從2019年7月—2020年1月存在一個緩慢抬升的趨勢,位于斜坡中部和后緣部位的G1、G2和G3的形變量隨著時間的變化而不斷增加。整體而言,斜坡中部處于不穩定的緩慢蠕動變形階段,而后緣則處于劇烈變形狀態,由于滑坡體后壁拉張裂縫的發育,在降雨的影響下,極易再次發生滑動(圖21)。

圖21 尕沙日滑坡累積形變量及時序位移曲線圖Fig.21 Cumulative deformation and time sequence displacement curve of Gasari Landslide
基于時序SBAS-InSAR技術,聯合升降軌SAR影像從不同成像幾何下反演獲取了研究區地表形變信息,結合孕災背景資料,輔以高分辨率光學影像篩選,對區內滑坡隱患進行了早期識別,得到如下結論。
(1)在2019年1月—2020年8月監測周期內,區內共識別出滑坡隱患19處,其中歷史滑坡8處,新發現隱患11處。所有隱患集中分布于隆務河兩岸和隆烏古曲河畔,涉及黃乃亥鄉、年都乎鄉、蘭采鄉、隆務鎮、加吾鄉、曲庫乎鄉、扎毛鄉等7個鄉鎮。
(2)經形變速率精度分析,本次形變監測結果精度優于8~11 mm,監測結果有效;結合區內升降軌SAR影像幾何畸變類型及分布,對識別出的19處隱患進行了相干性和幾何畸變復核檢驗,所有隱患均處在非幾何畸變區,形變監測結果可靠。同時結合幾何畸變對升降軌識別結果進行了差異性分析,19處隱患中,升軌識別8處,降軌識別11處,聯合升降軌數據從不同成像幾何下對目標地物的形變特征進行提取,在一定程度上彌補了單一成像幾何帶來的監測盲區,識別結果更加全面,提高了滑坡隱患早期識別的準確性。
(3)通過典型滑坡特征分析發現,本次研究識別出的隱患均處于坡度陡、高差大的斜坡上,沖溝發育、溯源侵蝕現象嚴重,為其滑動提供了良好的地形條件;滑坡類型均為土質滑坡,所處黃土區,坡體結構松散、空隙大,具有較強的透水性,獨特的巖土體條件極易發育形成滑坡災害;坡體植被稀疏,多分布有地表裂縫,為地表水的入滲提供了通道,致使坡體穩定性變差;從典型隱患形變的時間序列分析結果可知,各隱患坡體均呈緩慢變形趨勢,穩定性較差,在研究區5—10月,尤其是6—8月集中降雨期間,形變呈加速趨勢,降雨是該區誘發滑坡災害的主要因素。