高二濤， 李 豪， 雍 琦， 付波霖， 范冬林

(1.桂林理工大學(xué)測繪地理信息學(xué)院，桂林 541006； 2.廣西空間信息與測繪重點(diǎn)實(shí)驗(yàn)室，桂林 541006；3.山西省工業(yè)設(shè)備安裝集團(tuán)有限公司，太原 030000)

中國地質(zhì)環(huán)境復(fù)雜、氣候條件時(shí)空差異大，地震、滑坡、泥石流等地質(zhì)災(zāi)害類型多且分布廣，難以有效防范[1-2]。在上述自然災(zāi)害中，滑坡是一種危害程度僅次于地震，排名第二的環(huán)境地質(zhì)災(zāi)害?；聫V泛發(fā)育于河流、湖泊、海洋和溝渠的岸坡，地形高差大的山區(qū)、峽谷區(qū)，道路的斜坡地段及強(qiáng)降雨地區(qū)等，具有隱蔽性、突發(fā)性、不確定性等特點(diǎn)，且破壞力巨大[3]。蒲縣屬于晉西黃土高原地質(zhì)災(zāi)害區(qū)，地貌復(fù)雜多樣，此地區(qū)應(yīng)重點(diǎn)加強(qiáng)對(duì)居民居住區(qū)和采煤區(qū)的黃土滑坡、崩塌等地質(zhì)災(zāi)害的防范?；略谕话l(fā)前一般伴隨長期緩慢的地表變形，在此背景下，基于變形監(jiān)測為主的滑坡災(zāi)害早期識(shí)別與持續(xù)監(jiān)測可以為預(yù)報(bào)預(yù)警與防災(zāi)減災(zāi)工作提供科學(xué)依據(jù)，對(duì)于主動(dòng)防范地質(zhì)災(zāi)害意義重大[4]。

水準(zhǔn)測量、全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system，GNSS)觀測、攝影測量法和位移計(jì)等傳統(tǒng)滑坡監(jiān)測方法以其點(diǎn)觀測、高精度等特點(diǎn)，目前已被廣泛應(yīng)用于滑坡監(jiān)測中[5-7]。近年來，空間對(duì)地觀測技術(shù)迅速發(fā)展，以永久散射體雷達(dá)差分干涉測量技術(shù)(persistent scatterer interferometric synthetic aperture radar，PS-InSAR)與短基線集雷達(dá)干涉技術(shù)(small baseline subset，SBAS)為代表的時(shí)間序列合成孔徑雷達(dá)干涉(time series interferometric synthetic aperture radar，TS-InSAR)技術(shù)，可以遠(yuǎn)程監(jiān)測地球表面以獲取全球的高精度地形和地表微小形變信息，具有全天時(shí)、全天候、時(shí)空分辨率高、覆蓋范圍廣等特點(diǎn)，彌補(bǔ)了傳統(tǒng)滑坡形變監(jiān)測方法計(jì)算量大且監(jiān)測成本高的不足，在滑坡、地震、城市地面沉降等地質(zhì)災(zāi)害的普查與監(jiān)測預(yù)警中展現(xiàn)出巨大的發(fā)展?jié)摿8-12]。現(xiàn)采用PS-InSAR與SBAS技術(shù)分別處理覆蓋山西蒲縣地區(qū)24景哨兵-1A(Sentinel-1A)數(shù)據(jù)，獲取研究區(qū)域2019年年平均形變速率和累計(jì)形變位移量，通過結(jié)合已有的資料和谷歌地球(Google Earth)等數(shù)據(jù)，進(jìn)一步分析是否為滑坡隱患點(diǎn)及其可能性成因。

1 研究區(qū)概況與數(shù)據(jù)簡介

1.1 研究區(qū)概況

蒲縣(36°11′32″N～36°38′13″N，110°51′09″E～111°23′36″E)位于山西省西南部，臨汾盆地中心，呂梁山西南側(cè)，東部屬于呂梁山和太行山斷裂帶，西部屬于鄂爾多斯斷裂帶。蒲縣東接洪洞，西靠大寧，南依吉縣，北與隰縣、汾西毗鄰。境內(nèi)由山地、丘陵、河谷等地形組成，地勢剖面為東高西低，北、東、南三面環(huán)山，植被茂盛。西部與中部為典型的黃土高原地貌，第四紀(jì)黃土發(fā)育良好，地表大多被不同厚度的黃土覆蓋，地形破碎、垂直節(jié)理發(fā)育、溝谷深切超過200 m，易引發(fā)崩塌、滑坡、泥石流等多種次生地質(zhì)災(zāi)害。蒲縣是山西煤炭開采的重要基地之一，含煤面積占全縣總面積的九成以上[13]。蒲縣位于北半球中緯度大陸內(nèi)部，是典型的溫帶大陸性氣候。春季多風(fēng)沙，降雨量少；夏季多暴雨，雨熱同期；秋季短暫、溫和涼爽；冬季漫長、寒冷干燥。蒲縣地理位置如圖1所示。

圖1 蒲縣地理位置

1.2 實(shí)驗(yàn)數(shù)據(jù)

哨兵1號(hào)(Sentinel-1)衛(wèi)星是歐洲航天局(ESA)發(fā)射的載有C波段合成孔徑雷達(dá)的地球監(jiān)測衛(wèi)星[14]，可全天時(shí)、全天候提供連續(xù)圖像。哨兵1號(hào)是哥白尼計(jì)劃的一個(gè)重要組成部分，具有1A和1B雙系統(tǒng)。2014年4月發(fā)射的Sentinel-1A衛(wèi)星是太陽同步極地軌道衛(wèi)星，衛(wèi)星的重訪周期為12 d，軌道高度是693 km。它的主要用途是監(jiān)測地球環(huán)境變化，存在四種成像模式、兩種極化方式，具有生產(chǎn)能力強(qiáng)、重訪周期短的優(yōu)勢?，F(xiàn)采用Sentinel-1A衛(wèi)星獲取的24景干涉寬幅模式(interferometric wide swath，IW)的單視復(fù)數(shù)圖像數(shù)據(jù)(single look complex，SLC)。獲取影像的時(shí)間范圍是2019-01—2019-12，數(shù)據(jù)的詳細(xì)參數(shù)如表1所示。

表1 Sentinel-1A TOPS SAR影像數(shù)據(jù)參數(shù)

本文研究采用日本宇宙航空研究開發(fā)機(jī)構(gòu)獲取的30 m分辨率的AW3D30(ALOS World 3D-30 m)作為外部DEM(digital elevation model)數(shù)據(jù)。DEM的范圍在35°N～38°N、108°E～112°E，運(yùn)用其生成地質(zhì)背景構(gòu)造圖和去地形相位處理，以消除地形起伏的影響[15]。

2 時(shí)序InSAR數(shù)據(jù)處理

2.1 基于PS-InSAR技術(shù)的蒲縣滑坡形變監(jiān)測

首先，利用PS-InSAR技術(shù)監(jiān)測蒲縣地區(qū)滑坡形變過程，主要處理步驟如圖2所示。實(shí)驗(yàn)程序自動(dòng)選取2019-06-30的影像數(shù)據(jù)為公共主影像，其余23景從影像逐個(gè)與之配準(zhǔn)并建立連接關(guān)系，共得到23個(gè)干涉像對(duì)。SAR數(shù)據(jù)時(shí)空基線分布如圖3所示，具體時(shí)空基線如表2所示。可知，時(shí)間基線范圍為-180～168 d，空間基線范圍為-135～110 m。

表2 Sentinel-1A影像數(shù)據(jù)獲取時(shí)間及基線參數(shù)信息

圖2 PS-InSAR處理流程

圖3 PS-InSAR時(shí)空基線分布

使用AW3D30 DEM對(duì)干涉相位進(jìn)行去平處理，實(shí)驗(yàn)采用的DEM精度和分辨率越高，去除地形相位的效果越理想。打開Goldstein濾波器，將距離向多視數(shù)設(shè)置為5，方位向多視數(shù)設(shè)置為1，輸出差分干涉圖。選擇永久散射體(PS點(diǎn))時(shí)，應(yīng)選取有穩(wěn)定散射特性的點(diǎn)，如橋墩、大壩、房頂?shù)热嗽斓匚锘蛘呗銕r等自然地物。通過密集分布的PS點(diǎn)消除由于信號(hào)在大氣中傳播延遲造成的波動(dòng)，從而有效改善實(shí)驗(yàn)精度。為提高穩(wěn)定散射體點(diǎn)的準(zhǔn)確性，采用振幅離差指數(shù)法確定PS點(diǎn)，將相干系數(shù)閾值設(shè)置為0.75，對(duì)研究區(qū)域內(nèi)的穩(wěn)定散射目標(biāo)進(jìn)行有效識(shí)別[16]，如圖4所示，在圖4中用綠點(diǎn)標(biāo)注，PS點(diǎn)目標(biāo)主要分布在裸露的巖石、公路上，而植被覆蓋率高的地區(qū)幾乎沒有PS點(diǎn)。

圖4 PS-InSAR技術(shù)永久散射體(PS點(diǎn))分布

利用重去平的干涉圖建立線性模型，估算變形率、殘余高度等。大氣效應(yīng)的估算和清除通過兩種濾波程序進(jìn)行，分別是大氣高通濾波和大氣低通濾波。對(duì)大氣校正后的結(jié)果進(jìn)行地理編碼，即可得到蒲縣地區(qū)2019-01—2019-12年平均沉降速率(圖5)。對(duì)PS-InSAR技術(shù)所得結(jié)果圖進(jìn)行觀察分析可知，研究區(qū)域內(nèi)形變速率范圍為-53～48 mm/a。研究區(qū)域絕大部分地區(qū)處于穩(wěn)定狀態(tài)，形變速率范圍為-15～10 mm/a。

沿雷達(dá)視線方向LOS (line of sight)地表抬升的地區(qū)用紅色表示，發(fā)生地面沉降的區(qū)域顯示為藍(lán)紫色

2.2 基于SBAS技術(shù)的蒲縣滑坡形變監(jiān)測

為驗(yàn)證PS-InSAR方法監(jiān)測的準(zhǔn)確度與可靠性，采用SBAS方法在同一數(shù)據(jù)源的基礎(chǔ)上對(duì)山西蒲縣地區(qū)進(jìn)行滑坡監(jiān)測，SBAS處理流程如圖6所示。

圖6 SBAS處理流程

采用蒲縣地區(qū)2019-01—2019-12的24景Sentinel-1A數(shù)據(jù)，設(shè)置臨界基線最小、最大百分比分別為0和2，時(shí)間基線閾值設(shè)定為60 d，系統(tǒng)自主選取2019-06-18的影像為超級(jí)主影像。根據(jù)設(shè)定得到閾值，共生成170個(gè)干涉像對(duì)，圖7為SBAS技術(shù)時(shí)空基線分布圖。

圖7 SBAS時(shí)空基線分布

為使實(shí)驗(yàn)結(jié)果在距離向和方位向的分辨率盡可能保持相同的效果，實(shí)驗(yàn)將方位向和距離向多視數(shù)分別設(shè)置為1和4，與超級(jí)主影像保持一致。選用常被應(yīng)用于植被密集或環(huán)境潮濕地區(qū)的最小費(fèi)用流法(minimum cost flow， MCF)進(jìn)行相位解纏，將解纏相關(guān)系數(shù)閾值設(shè)置為0.23，利用Goldstein濾波方法以提高研究區(qū)中低相干位置的濾波強(qiáng)度。

打開存放差分干涉成果的文件夾，仔細(xì)查看每一個(gè)像對(duì)的相干性圖片以及解纏后的圖片，用連接圖編輯工具剔除相干性低和含有明顯大氣相位的像對(duì)。軌道精煉與重去平是為了估算差分后依舊殘存的恒定相位和解纏后還留有的相位坡道，利用GCP點(diǎn)對(duì)所有的數(shù)據(jù)對(duì)執(zhí)行重去平以消除地形相位。在實(shí)驗(yàn)過程中選取了遠(yuǎn)離形變區(qū)域、孤立相位、殘余地形條紋等位置，位于相對(duì)穩(wěn)定地點(diǎn)的70余個(gè)GCP點(diǎn)。在第一次反演中，對(duì)相干點(diǎn)建模并構(gòu)建方程組矩陣，由于線性模型的穩(wěn)定性最佳，因此，實(shí)驗(yàn)選用線性模型求解所有影像對(duì)的變形率和新的DEM等。基于第一次反演所得形變速率，第二次反演是通過設(shè)置大氣高通和大氣低通來估計(jì)和清除實(shí)驗(yàn)中的大氣影響，獲得更為準(zhǔn)確純凈的時(shí)間序列中的位移結(jié)果并進(jìn)行地理編碼，將其轉(zhuǎn)換到地理坐標(biāo)系，獲得形變信息和年平均位移速率，如圖8所示。研究區(qū)域形變速率范圍為-70～52 mm/a，有4處地表沉降較為明顯的疑似滑坡隱患地點(diǎn)，大部分地表形變?cè)?10～10 mm/a，趨于穩(wěn)定狀態(tài)。通過SBAS方法采集的蒲縣地區(qū)2019年清晰明了的時(shí)序地表變化過程如圖9所示。

紅色表示雷達(dá)視線方向存在抬升，藍(lán)色表示下沉

圖9 基于SBAS技術(shù)的蒲縣地區(qū)2019年地表累計(jì)形變

3 時(shí)序形變結(jié)果分析

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

基于PS-InSAR技術(shù)和SBAS方法，利用哨兵升軌SAR數(shù)據(jù)分別獲取蒲縣地區(qū)年平均形變速率，通過目視解譯選取出兩種技術(shù)所得結(jié)果中共同存在的4處較為明顯和典型的形變位置進(jìn)行對(duì)比分析，所選位置分布如圖10所示，包括太林鄉(xiāng)半溝村羅克線沿線、水泉窊015鄉(xiāng)道兩側(cè)、郭家山周邊、井子洼村趙克路附近，具體地理位置如表3所示。圖11為4處形變速率放大圖，從放大圖中可以清楚地看到4個(gè)疑似隱患點(diǎn)的沉降信息。除圖11所示4處較明顯的疑似滑坡隱患點(diǎn)，研究區(qū)域仍存在部分升降不一的形變點(diǎn)，為使結(jié)果更加精準(zhǔn)可靠，只對(duì)上述4處共同存在的疑似滑坡點(diǎn)做進(jìn)一步對(duì)比分析。

圖11 4處疑似滑坡點(diǎn)SBAS技術(shù)形變速率放大圖

表3 4處疑似滑坡點(diǎn)詳細(xì)地理位置

圖10 4處疑似滑坡點(diǎn)位置

為了更加精準(zhǔn)地分析PS-InSAR技術(shù)與SBAS技術(shù)對(duì)于同一數(shù)據(jù)源的結(jié)果精度對(duì)比情況，基于4處疑似滑坡點(diǎn)累計(jì)形變量的詳細(xì)信息，計(jì)算每一期數(shù)據(jù)的差值，并用兩種數(shù)據(jù)所得結(jié)果的平均值代替真值計(jì)算其均方根誤差(RMSE)，RMSE代表每個(gè)數(shù)據(jù)與真值間的偏離程度，均方根誤差越大，表示兩種技術(shù)獲取的結(jié)果差異越大，可靠性越低；反之，則表示測量的精度越高，進(jìn)一步說明時(shí)序InSAR技術(shù)識(shí)別并監(jiān)測滑坡的有效性，如表4所示。

由表4可得，一號(hào)疑似滑坡點(diǎn)PS-InSAR技術(shù)與SBAS技術(shù)所得最終累計(jì)沉降量分別為-51.23 mm 和-52.32 mm，差值為1.09 mm，RMSE范圍在0～2.25 mm，RMSE均值為0.75 mm；二號(hào)疑似滑坡點(diǎn)PS-InSAR技術(shù)與SBAS技術(shù)所得最終累計(jì)沉降量分別為-48.03 mm和-49.51 mm，差值為1.48 mm，RMSE范圍在0～2.06 mm，RMSE均值為0.77 mm；三號(hào)疑似滑坡點(diǎn)PS-InSAR技術(shù)與SBAS技術(shù)所得最終累計(jì)沉降量分別為-40.87 mm和-42.42 mm，差值為1.55 mm，RMSE范圍在0～3.25 mm，RMSE均值為1.22 mm；四號(hào)疑似滑坡點(diǎn)PS-InSAR技術(shù)與SBAS技術(shù)所得最終累計(jì)沉降量分別為-46.91 mm和-39.80 mm，差值為7.11 mm，RMSE范圍在0～3.55 mm，RMSE均值為1.05 mm。由上述數(shù)據(jù)可得PS-InSAR技術(shù)與SBAS技術(shù)監(jiān)測結(jié)果的均方根誤差在0～4 mm，均值在0.5～1.5 mm，充分說明了PS-InSAR技術(shù)與SBAS技術(shù)的一致性較高。

表4 4處疑似滑坡點(diǎn)累計(jì)形變量的詳細(xì)信息

為了使兩種時(shí)序InSAR技術(shù)所得監(jiān)測結(jié)果對(duì)比的更為直觀簡明，基于上述4處疑似滑坡點(diǎn)的累計(jì)形變量數(shù)據(jù)，對(duì)選取出的4處點(diǎn)繪制累計(jì)形變折線圖進(jìn)行對(duì)比分析。圖12為識(shí)別出的4處疑似滑坡點(diǎn)2019-01-01—2019-12-15的LOS向累積形變位移量對(duì)比。總體看來，兩種技術(shù)所得結(jié)果的沉降趨勢較為一致，說明2019年蒲縣地區(qū)部分地區(qū)發(fā)生地表沉降，有疑似滑坡的可能，更進(jìn)一步說明了時(shí)序InSAR識(shí)別并監(jiān)測滑坡的可行性與準(zhǔn)確性。1號(hào)疑似滑坡點(diǎn)在1—5月地表相對(duì)穩(wěn)定，6月出現(xiàn)小幅沉降至20 mm，7月地面出現(xiàn)略微抬升，8—10月緩慢沉降，并在11、12月急劇下沉至-50 mm左右。2號(hào)疑似滑坡點(diǎn)1—6月緩慢沉降至-21 mm左右，并在7、8月出現(xiàn)小幅抬升，隨后逐漸沉降至-49 mm左右。3號(hào)疑似滑坡點(diǎn)沉降成線性變化，4號(hào)疑似滑坡點(diǎn)1—4月地表處于穩(wěn)定狀態(tài)，隨后呈線性下沉模式。

圖12 4處疑似滑坡點(diǎn)累計(jì)形變位移量PS-InSAR與SBAS監(jiān)測結(jié)果對(duì)比

3.2 可能性誘因分析

蒲縣是典型的溫帶大陸性氣候，雨熱同期。由上述分析可知，疑似滑坡點(diǎn)7、8月地面出現(xiàn)略微抬升，可能是由于雨季降雨量驟增的緣故，地表水滲入地下，導(dǎo)致地下含水層增厚。11、12月地表下沉明顯，可能是由于冬季降水量過少，當(dāng)?shù)鼐用翊罅渴褂玫叵滤木壒蕦?dǎo)致地下水位下降。3號(hào)疑似滑坡隱患點(diǎn)位于郭家山周邊，郭家山地區(qū)存在采煤區(qū)，該點(diǎn)的沉降疑似受到煤炭開采等人類活動(dòng)的影響。

4 結(jié)論

利用時(shí)序InSAR技術(shù)探測研究蒲縣地區(qū)的滑坡形變，獲取研究區(qū)域內(nèi)2019年年平均形變速率和累計(jì)形變位移量，確定太林鄉(xiāng)半溝村羅克線沿線、水泉窊015鄉(xiāng)道兩側(cè)、郭家山周邊、井子洼村趙克路附近等4處位置為疑似滑坡點(diǎn)，并結(jié)合煤炭開采、人類活動(dòng)、降雨量信息等，對(duì)該地區(qū)分布的疑似滑坡隱患進(jìn)行早期識(shí)別與監(jiān)測預(yù)警，為該地區(qū)滑坡誘因提供科學(xué)的依據(jù)，主動(dòng)規(guī)避自然地質(zhì)災(zāi)害風(fēng)險(xiǎn)，對(duì)未來滑坡災(zāi)害的監(jiān)測與預(yù)警提供思路與方法。