基于InSAR技術(shù)的流域庫岸地質(zhì)災(zāi)害監(jiān)測研究

王國強，陳仁琛

(1.水利部水利水電規(guī)劃設(shè)計總院，北京 100120；2.國網(wǎng)物資有限公司，北京 100120)

隨著社會經(jīng)濟的發(fā)展，水利工程發(fā)揮著越來越大的作用，然而，在實現(xiàn)水利工程對社會創(chuàng)造經(jīng)濟效益的同時，也必須注意對水庫庫岸邊坡造成的影響，避免發(fā)生滑坡、泥石流等地質(zhì)災(zāi)害。這些地質(zhì)災(zāi)害的發(fā)生不僅會對水利工程造成破壞，還會嚴(yán)重威脅周邊地區(qū)人民、移民安置點居民的生命財產(chǎn)安全，造成巨大的社會經(jīng)濟損失[1,2]。作為對這些地質(zhì)災(zāi)害進行監(jiān)測分析最直接的手段，水庫庫岸邊坡高精度形變監(jiān)測具有十分重要的意義[3-6]。

目前已有多種方法被聯(lián)合用于滑坡體形變監(jiān)測，如測地測量，基于全球定位系統(tǒng)的監(jiān)測網(wǎng)絡(luò)，以及遙感圖像的解譯。盡管這些方法可以獲得高精度的形變監(jiān)測信息，但由于人力和儀器成本的限制，它們不容易獲得高密度的測量點。合成孔徑雷達(Interferometric Synthetic Aperture Radar, InSAR)技術(shù)的出現(xiàn)為庫岸邊坡形變監(jiān)測提供了契機，InSAR技術(shù)可以大規(guī)模地獲得地表形變量，精度可達到厘米級，已在地震、火山和采礦變形的監(jiān)測等多個行業(yè)中得到了廣泛的應(yīng)用[8,9]，但在水利水電行業(yè)還沒有得到很好的應(yīng)用。

本文選取金沙江流域某大型水電站所處的某段流域作為研究區(qū)，首先收集7景ALOS PALSAR2數(shù)據(jù)，基于小基線集技術(shù)(Small Baseline Subset, SBAS)，對整個研究區(qū)的邊坡進行形變區(qū)域普查；然后對普查獲取的形變區(qū)安裝GNSS設(shè)備，進行接觸式、不間斷形變信息獲取；最后基于InSAR技術(shù)獲取的形變信息與GNSS設(shè)備獲取的形變信息進行對比分析，以證明基于InSAR技術(shù)獲取的形變信息的準(zhǔn)確性，同時構(gòu)建“星-地”一體化監(jiān)測體系，通過全面的信息采集和處理，對水庫沿岸的穩(wěn)定性進行監(jiān)測，分析庫岸邊坡的地質(zhì)災(zāi)害分布點與各災(zāi)害點區(qū)域的形變特性及趨勢，為水利樞紐工程的安全運營、航道通航及移民生命財產(chǎn)安全防范提供了全面準(zhǔn)確的預(yù)警信息。

1 實驗區(qū)的基本概況

本文試驗區(qū)選在西南某大型水利樞紐工程，針對該水庫庫岸邊坡的地質(zhì)災(zāi)害問題，進行了地質(zhì)災(zāi)害危險性評估，結(jié)合所選數(shù)據(jù)的覆蓋范圍，整個評估區(qū)的流域長度約為75 km。該流域?qū)儆诮鹕辰饔颍庸瘸省癡”形谷，海拔在2 200～3 500 m之間，垂直高程落差達到1 300 m，屬于典型的構(gòu)造侵蝕高山強烈寒凍風(fēng)化地貌和大陸性季風(fēng)高原型氣候區(qū)，極易產(chǎn)生地質(zhì)災(zāi)害。流域沿岸分布有大量人類聚居區(qū)域，移民安置點及道路、橋梁等基礎(chǔ)建筑設(shè)施，實驗區(qū)范圍及流域位置如圖1所示。從圖1中可以看出，實驗區(qū)范圍遠大于河流范圍。這是由于在基于InSAR技術(shù)的數(shù)據(jù)處理過程中，需要在監(jiān)測區(qū)內(nèi)選擇穩(wěn)定點作為參考點，擴大測區(qū)范圍能夠在確立參考點時有更多選擇。

鑒于該區(qū)域地質(zhì)災(zāi)害風(fēng)險點諸多，涉及面廣，且水庫蓄水處于較高水位運行，若庫岸滑坡體整體塌滑入庫區(qū)或出現(xiàn)崩塌及巖溶塌陷等地質(zhì)災(zāi)害，這都將會對大壩及電站運行帶來較大的安全隱患，同時對航道及庫岸兩邊的城鎮(zhèn)造成巨大的損害。因此，需對該評估區(qū)內(nèi)的地質(zhì)災(zāi)害點進行形變監(jiān)測。

2 理論及方法

2.1 數(shù)據(jù)源簡介

本文收集了2018年5月至2019年5月期間的7景ALOS PALSAR2雷達數(shù)據(jù)(見表1)，其工作模式為FBS(fine beam single polarization)，L波段(波長23.6 cm)，產(chǎn)品級別為Level1.1(單視復(fù)數(shù)影像數(shù)據(jù))，HH極化，視角為32.8°，方位向像元大小為2.129 888 m，距離向像元大小為1.430 422 m。

表1 實驗區(qū)ALOS PALSAR2數(shù)據(jù)

為了去除地形相位，在進行二軌差分法差分干涉測量時需要利用高精度DEM，本文采用美國地質(zhì)調(diào)查局發(fā)布的30 m SRTM DEM數(shù)據(jù)。

2.2 SBAS-InSAR技術(shù)

SBAS-InSAR技術(shù)的基本思想是采用多干涉集合組網(wǎng)的方式，通過設(shè)置時間和空間基線閾值，組合出短時間和空間基線的干涉對。這樣可以盡量保持干涉圖的相干性，在提取一定時間基線范圍和空間基線范圍內(nèi)具有穩(wěn)定相干性的分布式點目標(biāo)時，能夠提高相干點密度，適用于自然地表的形變監(jiān)測。

假設(shè)獲取了監(jiān)測區(qū)共有N+1幅雷達影像數(shù)據(jù)，首先選取其中一幅影像作為主影像，將其他N幅影像配準(zhǔn)到主影像的影像坐標(biāo)下。然后設(shè)置合適的時間、空間基線閾值，將符合閾值條件的干涉對進行差分干涉，生成M幅差分干涉圖，M滿足以下條件：

(1)

Berardino等人在2002年提出的SBAS-InSAR算法，對差分干涉圖首先進行了解纏，也有很多后續(xù)發(fā)展的小基線算法利用的是未解纏的干涉圖，在選完點目標(biāo)之后在點目標(biāo)上進行解纏。這里假設(shè)輸入的干涉圖為解纏后的干涉圖[10]。在tA和tB時刻獲取的干涉對生成的第j幅干涉圖中，任意一個像元的干涉相位可以表示為：

δφj=φtB-φtA≈(φdtB-φdtA)+

(2)

(3)

結(jié)合差分干涉原理，上式可寫成：

(4)

式中：Δφnonlinear為非線性形變引起的相位差；Δφatm為大氣延遲引起的相位差；Δφnoi為噪聲引起的相位差。

當(dāng)公式(4)中Δφnonlinear+Δφatm+Δφnoi的絕對值小于π時進行解纏。解纏后可通過以下公式對時間序列上每個解纏結(jié)果的進度進行評價。

(5)

式中：Δωi為殘差相位。當(dāng)評價因子γ小于預(yù)先設(shè)置的閾值時，將該條邊的解纏結(jié)果舍棄到三角網(wǎng)的合邊。當(dāng)γ取最大值時，高程誤差和形變速率的估計值為最優(yōu)估計值，將這兩個估計值從公式(4)中移除，剩下的殘余相位中包括非線性形變相位、大氣延遲相位和相干噪聲相位。根據(jù)非線性形變和大氣延遲相位之間不同的特點，可以將二者逐一分離。大氣延遲相位在空間維度中表現(xiàn)為相關(guān)的低頻信息，而在時間維度上表現(xiàn)為隨機的高頻信息；而非線性形變相位在空間和時間維度中都表現(xiàn)為相關(guān)的低頻信息。因此，通過時間域高通濾波和空間域低通濾波的方法可以將大氣延遲相位分離出來[11]。

當(dāng)去掉差分相位中的DEM誤差和大氣延遲相位后，干涉相位可以簡化為：

(6)

設(shè)相鄰時間段內(nèi)的形變是線性的，即整個時間段內(nèi)的形變是分段線性的，則第j景干涉圖監(jiān)測的形變相位值可表示為：

(7)

然后將所有的解纏后的差分干涉相位用矩陣形式表示：

Bv=δφ

(8)

式中：矩陣B為一個M×N的矩陣。SBAS-InSAR方法為了防止長時間、空間基線造成的失相干現(xiàn)象，M幅差分干涉圖通常存在多個子集，當(dāng)存在不連續(xù)的子集時，矩陣B的秩小于N，即B出現(xiàn)秩虧的現(xiàn)象。為此使用奇異值分解(SVD)方法求出矩陣B的廣義逆矩陣，并求得最終的形變速率[12]。

具體的處理流程如圖2。

圖2 SBAS-InSAR數(shù)據(jù)處理流程圖

3 實驗過程及結(jié)果驗證

3.1 實驗過程

因?qū)嶒瀰^(qū)地處我國西南部地區(qū)，山高林密，植被覆蓋量常年較高，對短波信號的影響較大。L波段對植被冠層的穿透性較高，保證了ALOS-PALSAR 2在該區(qū)域的相干性[13-15]。同時，ALOS-PALSAR 2的衛(wèi)星軌道離差較小，根據(jù)實際應(yīng)用經(jīng)驗來看，有效基線一般保持在300 m以內(nèi)。但由于本次實驗收集的SAR影響數(shù)量較少，因此將時間基線閾值設(shè)置為365 d，即所有影象兩兩組合，共計構(gòu)成21組干涉對，干涉對連接圖如圖3所示。

圖3 雷達干涉對連接圖

為了避免低相干點帶來的誤差，在差分干涉相位圖中對高相干點目標(biāo)進行選取，將本次相干性閾值設(shè)置為0.3。對干涉對進行查分干涉處理、濾波和基于最小費用流的相位解纏，獲取到21幅相位解纏圖，如圖4所示。

圖4 相位解纏圖集

在SBAS-InSAR時序監(jiān)測過程中，對殘余相位進行時間域上的高通濾波，計算結(jié)果中包含有大氣相位以及相干噪聲。在此基礎(chǔ)上，聯(lián)合所有相干圖像元在同一差分干涉圖中進行空間域上的低通濾波，以濾除相干噪聲，剩余相位即大氣相位。經(jīng)過不斷地迭代反演，對地形誤差相位、基線誤差相位以及大氣延遲相位進行改正去除，獲取了最終時序形變相位結(jié)果，將結(jié)果進行相位到形變量的轉(zhuǎn)換，并以第一景影像獲取時間為起始形變時間，獲取實驗區(qū)金沙江流域5 km緩沖區(qū)內(nèi)的地表累積形變量圖如圖5所示。根據(jù)實驗結(jié)果，通過目視解譯，共發(fā)現(xiàn)3處形變坡體。

3.2 精度驗證

選擇圖5中3號形變區(qū)域作為驗證坡體，圖6展示的為3號形變區(qū)域多期形變監(jiān)測結(jié)果。由監(jiān)測結(jié)果顯示，3號邊坡的主要形變發(fā)生在該坡體的東北坡，而西北坡在監(jiān)測結(jié)果上顯示為疊掩區(qū)域，這是由于SAR影像采用側(cè)視成像的方式，傳感器視角接近坡度角是，坡頂與坡底的信號將集中在一起，形成疊掩，而無法獲取有效數(shù)據(jù)。因此，在坡體東北面上共計安裝4組GNSS設(shè)備，其中1組作為控制基站(圖7中基站A)，3組作為測量基站(圖7中基站B～D)?；静贾眉捌麦w情況如圖7所示。

圖5 實驗區(qū)金沙江流域5 km緩沖區(qū)地表累積形變量

由于D基站在運行過程中出現(xiàn)故障，我們采集了B和C兩個基站在2018年5月28日-2019年5月27日之間豎直方向的累計形變信息，監(jiān)測頻率為每5 d一次。同時，由于InSAR監(jiān)測獲取的形變?yōu)橐暰€向(Line of Sight, LOS)的形變，因此，我們將InSAR監(jiān)測獲取的累計形變量投影到豎直方向，與GNSS監(jiān)測獲取的結(jié)果進行對比，對比曲線如圖8(a)和圖8(b)所示。

圖8 兩種監(jiān)測結(jié)果對比統(tǒng)計圖

根據(jù)兩個基站的監(jiān)測結(jié)果可以看出兩種監(jiān)測方式獲得的監(jiān)測結(jié)果的變化趨勢幾乎相同，但GNSS監(jiān)測獲取的形變量要高于InSAR監(jiān)測獲取的形變量。這是由于GNSS獲取的是該點在豎直方向上實際的形變量，而InSAR僅能獲取該點雷達LOS方向上形變的豎直分量。通過兩種監(jiān)測結(jié)果的對比分析，證明了基于L波段雷達影像的星載InSAR能夠用于流域邊坡的形變普查，同時，結(jié)合GNSS監(jiān)測手段，可以提高邊坡監(jiān)測結(jié)果的時間分辨率，提高監(jiān)測精度。

4 結(jié) 語

周期性的對水利工程庫岸邊坡進行形變監(jiān)測，積極防治水庫沿岸地質(zhì)災(zāi)害的發(fā)生，這對保障周邊居民生命財產(chǎn)安全和確保水利樞紐正常運行具有重要意義。本文利用星載InSAR技術(shù)對實驗區(qū)進行邊坡形變普查，再利用GNSS對普查發(fā)現(xiàn)的高危坡體進行重點觀測，同時將GNSS的監(jiān)測數(shù)據(jù)與InSAR監(jiān)測結(jié)果進行了驗證，提高了監(jiān)測結(jié)果的時間分辨率，構(gòu)建了“星-地”一體化的庫岸邊坡監(jiān)測體系。本文提出的監(jiān)測方案可以解決傳統(tǒng)的庫岸邊坡監(jiān)測手段難以獲取大范圍高精度形變監(jiān)測信息的局限性，保障了庫區(qū)移民安置居民和水庫的運行安全，同時可以大幅降低成本，具有較強的推廣應(yīng)用價值，可以向路橋等行業(yè)進行推廣應(yīng)用。

□