李一帆 杜柏林 付鈺涵 楊紅磊

（1.中國地質(zhì)大學(xué)（北京）土地科學(xué)技術(shù)學(xué)院；2.中國礦業(yè)大學(xué)（北京）地球科學(xué)與測繪學(xué)院）

獨特的地域和地質(zhì)環(huán)境使我國成為世界上地質(zhì)災(zāi)害最嚴重、受威脅人口最多的國家之一，其中滑坡是地質(zhì)災(zāi)害中發(fā)生頻率和危害最大的一種，特別是近年來由于人類活動及極端氣候條件等因素，滑坡災(zāi)害發(fā)生的頻率越來越高［1］。我國已將地質(zhì)災(zāi)害防治視作生態(tài)文明和美麗中國建設(shè)的重要內(nèi)容［2-3］，因此加強滑坡災(zāi)害的防災(zāi)減災(zāi)研究工作十分必要。

滑坡的發(fā)生、發(fā)展和演化過程伴隨著大量的宏觀可觀測數(shù)據(jù)，如滑坡形變、物理場變化和地下水位變化等，在眾多可觀測的數(shù)據(jù)或物理量中，形變是表現(xiàn)滑坡變化最顯著、最直接的變量，其形變趨勢又與滑坡所處階段存在良好的映射或函數(shù)關(guān)系，因此研究滑坡變形監(jiān)測技術(shù)對于滑坡預(yù)報和邊穩(wěn)定性分析具有重要意義?，F(xiàn)有的基于點位監(jiān)測技術(shù)，對于高陡順層邊坡而言，存在監(jiān)測點布設(shè)困難的問題，即使是布設(shè)了密集的監(jiān)測設(shè)備也很難獲取邊坡的空間連續(xù)和關(guān)鍵部位全過程形變信息，而這些信息對于高風(fēng)險的滑坡來講又非常關(guān)鍵，近些年發(fā)展起來的地基雷達干涉測量技術(shù)（Ground-based Interferometry Synthetic Aperture Radar，GBInSAR）［4］不受霧、雨、雪、粉塵等影響，可以實現(xiàn)全天時、全天候、高精度和高時空分辨率的觀測，是對局部區(qū)域形變非接觸監(jiān)測的一種新技術(shù)手段，獲取的形變結(jié)果更有助于滑坡災(zāi)害力學(xué)模型的理解與穩(wěn)定性評價。當(dāng)前存在各種體制的地基干涉雷達設(shè)備，其中GAMMA Remote Sensing公司的研制的GPRI-II設(shè)備與同類設(shè)備相比，能夠?qū)崿F(xiàn)長距離360°監(jiān)測，具有監(jiān)測范圍大、精度高的優(yōu)勢，同時該設(shè)備與專業(yè)數(shù)據(jù)處理軟件GAMMA相適配，功能強大，市場占有率較高［5］。

與GPRI-II 設(shè)備配套的GAMMA 軟件模塊內(nèi)的程序需在控制臺上運行，無可視化界面，雖然程序的組合十分方便，但需要用戶具有一定的編程基礎(chǔ)來實現(xiàn)批處理，多采用事后的方式對數(shù)據(jù)進行處理，這種處理方式具有延后性，不利于滑坡災(zāi)害的監(jiān)測預(yù)警。市面上具有可視化界面且能夠?qū)崿F(xiàn)實時處理的軟件多于地基雷達硬件匹配，由于數(shù)據(jù)格式不兼容的原因，無法實現(xiàn)對GPRI-II 設(shè)備數(shù)據(jù)的處理。因此迫切需要在GAMMA 軟件的基礎(chǔ)上進行二次開發(fā)，編制一套具有可視化界面的地基干涉雷達實時處理軟件。為解決上述問題，本文提出一種地基干涉雷達數(shù)據(jù)實時處理流程，在GAMMA 軟件基礎(chǔ)上，采用Python語言設(shè)計與開發(fā)地基干涉雷達實時處理軟件，對于滑坡災(zāi)害監(jiān)測具有重要的實際意義。

1 原 理

1.1 DInSAR技術(shù)

DInSAR 技術(shù)即差分干涉測量技術(shù)，與星載SAR不同，地基SAR觀測方程如式（1）所示［5］：

式中，d 為干涉對時間間隔內(nèi)的形變量；λ 為雷達波長；δφatmo是大氣效應(yīng)引起的相位延遲誤差；δφnoise是噪聲引起的相位誤差；k為整周數(shù)；ΔφPP'為干涉相位，對應(yīng)的值為相位主值，介于( -π,π ]，需要進行解纏才能獲取到真實的形變相位。

地基干涉雷達大多采用Ku波段，為短波范疇，易受到大氣延遲誤差的影響，研究表明溫度在20℃時，距離雷達1km 處，l%的相對濕度的變化可導(dǎo)致2mm的測量誤差。在本文中基于大氣同質(zhì)（大氣擾動引起的相位差與雷達和目標(biāo)距離存在線性關(guān)系）的前提，構(gòu)建二階多項式擬合大氣擾動特征［6］，如式（2）所示。

式中，r為目標(biāo)點到雷達中心的距離；a、b、c為該方程的系數(shù)；φatm為該點解纏后的干涉相位。

利用穩(wěn)定點擾動值和該點與雷達距離求解系數(shù)，最后根據(jù)求解的大氣延遲相位模擬整個監(jiān)測場景的大氣延遲相位。采用DInSAR 技術(shù)獲取監(jiān)測目標(biāo)形變的過程，包括影像配準(zhǔn)、干涉圖像生成、相位降噪濾波、相位解纏、去除大氣延遲、形變計算、地理編碼7 個步驟［7］。由于估算大氣延遲相位時采用大氣參數(shù)同質(zhì)的假設(shè)較為理想，因此不能完整地消除影響，獲得形變監(jiān)測精度有限；對于采樣間隔長的干涉，又易受到失相干因素的影響。因此，DInSAR技術(shù)不用于變形監(jiān)測；但是其構(gòu)建的函數(shù)模型可作為時序InSAR技術(shù)的基礎(chǔ)模型。

1.2 時序InSAR技術(shù)

為解決大氣延遲與時空失相干問題，在DInSAR技術(shù)基礎(chǔ)上發(fā)展出時序InSAR技術(shù)，采用多時相SAR影像，依據(jù)每個像元在時間維的相位值質(zhì)量的穩(wěn)定性，識別出高質(zhì)量的點，稱之為PS 點，對離散的PS 點構(gòu)網(wǎng)，建立相鄰點監(jiān)測相位增量，根據(jù)相鄰點間的形變模式，建立合適的形變模型，進而獲取形變速率和累計形變量。相應(yīng)的技術(shù)有PSInSAR、IPTA 和SBAS等［8-9］，其中SBAS 技術(shù)以高時空相干性的原則，多參考影像的原則確定干涉對組合。對于地基干涉雷達，由于不存在空間基線，因此在干涉對組合時只考慮時間基線，短時間基線一方面可以減弱大氣環(huán)境變化的影響，另一方面可以保持干涉對的相干性。大氣延遲誤差按照同質(zhì)大氣和非同質(zhì)大氣分兩步進行剔除，第一步根據(jù)式（2）分離出同質(zhì)大氣的延遲相位；第二步根據(jù)初始獲取的形變分布，確定出穩(wěn)定的點，根據(jù)穩(wěn)定點的相位，采用空間插值的方法擬合出形變的延遲相位；最終把兩部分大氣延遲相位累計相加，即為總體大氣延遲相位。大氣延遲相位剔除后，根據(jù)干涉對組合方式確定的觀測方程系數(shù)舉證，采用SVD分解方法計算出累計形變量。

1.3 形變預(yù)警原理

斜坡的穩(wěn)定性與斜坡當(dāng)前變形階段聯(lián)系緊密。松散土質(zhì)斜坡和具有時效變形特點的巖質(zhì)斜坡，變形往往呈現(xiàn)出蠕變的特點［10］，其變形過程一般可分為初始變形、等速變形和加速變形3 個階段［11］，在加速變形階段又包含初加速、中加速和臨滑3 個亞階段［12］。斜坡一旦進入加速階段，遲早會發(fā)生滑坡，當(dāng)處于臨滑階段時，斜坡整體失穩(wěn)即將滑坡［13］。閾值預(yù)警法是目前國際上最常用的方法［14］，即通過統(tǒng)計分析等給定預(yù)警閾值，當(dāng)監(jiān)測數(shù)據(jù)達到閾值時發(fā)出警示信息。由于不同物質(zhì)組成、不同規(guī)模、不同成因類型的滑坡變形閾值差異非常大，不存在統(tǒng)一的閾值［14］。本文使用速率曲線和加速度曲線的雙重閾值法，達到閾值進行警報。

1.4 實時處理流程

本文基于DInSAR 技術(shù)與時序InSAR 技術(shù)提出一種實時處理的流程，如圖1所示。基于地基SAR零空間基線和短時間基線的特點［15］，假設(shè)短時間內(nèi)監(jiān)測目標(biāo)形變呈現(xiàn)線性趨勢來對小基線集方法進行簡化，每次選取連續(xù)的一定數(shù)量（20 景左右）的監(jiān)測數(shù)據(jù)組合干涉對，根據(jù)時序干涉相干性的平均值和標(biāo)注差雙閾值確定高相干點。對每個干涉對進行大氣延遲相位分離，剩下的相位作為觀測值，采用最小二乘線性擬合［16］，即可得最終的形變量。

式中，A為時間間隔系數(shù)矩陣；b為形變值矩陣。

以一個像元點( i,j )為例，設(shè)有n副干涉影像參與解算，實際計算形式如式（4）所示，其中vˉ(i,j)是擬合出的該點的平均形變速率（一階項），c(i,j)是擬合出的該點的常數(shù)項，d(i,j)n表示第n景干涉圖像( i,j )點的形變值。

2 地基干涉雷達軟件系統(tǒng)開發(fā)

地基干涉雷達實時處理軟件是一個在GAMMA軟件的基礎(chǔ)上進行二次開發(fā)的軟件。該軟件使用并整合GAMMA 命令集，以InSAR 數(shù)據(jù)處理算法為核心，以合理的體系架構(gòu)模塊設(shè)計為基礎(chǔ)，來為用戶提供友好的界面和良好的使用體驗。

2.1 開發(fā)環(huán)境與配置

軟件開發(fā)采用Python 語言，使用Pyside2 庫（Qt for Python）并結(jié)合如Numpy 等Python 眾多第三方庫在GAMMA軟件的基礎(chǔ)上進行二次開發(fā)。

Qt為Python編程環(huán)境提供支持，開發(fā)可在Python上使用的庫PySide2，該庫包含絕大多數(shù)Qt 的類與函數(shù)，用法與Qt 內(nèi)的類與函數(shù)十分接近。安裝PySide2及其它第三方庫僅需通過python的包安裝組件pip下載即可。

GAMMA 腳 本 文 件 使 用 如sh、csh、tcsh、bash、perl、python等多種腳本語言完成編寫，要想成功運行GAMMA軟件，除要含有GAMMA軟件本體外，還應(yīng)當(dāng)含有對這些腳本語言的環(huán)境支持。由于本文軟件開發(fā)使用Python語言，在文件打包時便已經(jīng)包含Python語言環(huán)境支持，而且sh、csh、tcsh、bash 等都是Unix 或Linux 系統(tǒng)的命令交互語言，因此，在win 系統(tǒng)上需安裝對Unix/Linux 命令及Perl 語言的支持，在Unix 和Linux僅需安裝Perl語言支持即可。

同時，GAMMA 公司已經(jīng)推出Python 模塊文件py_gamma.py，在完成環(huán)境變量配置并正確導(dǎo)入該文件后，可實現(xiàn)在Python 中運行GAMMA 命令。原始的py_gamma.py文件需要存放在GAMMA 軟件根目錄下才能使用，這將在軟件打包后產(chǎn)生找不到文件的錯誤。查看相關(guān)源代碼后發(fā)現(xiàn)，該文件在原理上是通過添加環(huán)境變量的方式來查找同目錄下GAMMA 軟件的可執(zhí)行文件與腳本文件，故可對該部分代碼稍作修改，使得在打包后的軟件根目錄內(nèi)也能向環(huán)境變量添加正確的文件路徑。

2.2 軟件組成

本文軟件模塊組成與關(guān)系如圖2所示。

（1）格式轉(zhuǎn)換模塊。本軟件不僅支持GPRI-II 設(shè)備，還可以將其他品牌的地基干涉雷達數(shù)據(jù)轉(zhuǎn)換為GAMMA格式，如Fast GBSAR設(shè)備等。

（2）差分干涉處理模塊。該模塊對不同時間獲取的兩景SAR 影像進行處理，包括配準(zhǔn)、干涉、濾波和相位解纏等功能，最終獲得兩景影像之間的形變圖像。

（3）時序處理模塊。該模塊每次對連續(xù)的19 個干涉對間形變圖進行處理。該模塊流程含有生成累積形變圖和形變速率2 個步驟。模塊一次運行結(jié)束后獲得19 個累積形變圖和1 個形變速率圖。其中19個累計形變圖分別記錄對應(yīng)影像與第一景影像之間的累積形變，與之前數(shù)據(jù)有關(guān)聯(lián)；那一個形變速率圖記錄該20景影像內(nèi)的平均速率，與之前數(shù)據(jù)無關(guān)聯(lián)。

（4）地理編碼和成果輸出模塊。該模塊將雷達坐標(biāo)系下的成果數(shù)據(jù)轉(zhuǎn)為地理坐標(biāo)系下的成果，并轉(zhuǎn)換輸出可被如ArcGIS、QGIS 和Google Earth 等軟件使用的文件。

（5）預(yù)警處理模塊。該模塊對所有累積形變圖和形變速率圖上的某一特定區(qū)域（人為設(shè)置的預(yù)警區(qū)域）進行處理。獲得統(tǒng)計曲線并對標(biāo)記形變速率和加速度高于預(yù)警閾值的時間段。

（6）項目文件與日志模塊。該模塊與軟件設(shè)置和軟件運行信息有關(guān)，可分為兩部分：項目文件以及日志。項目文件儲存軟件相關(guān)的所有設(shè)置，能夠保證軟件在打開時恢復(fù)之前運行進度。日志部分將軟件運行時的報錯信息與運行信息輸出到文件中，能夠更好地發(fā)現(xiàn)錯誤的原因。

（7）顯示模塊。該模塊將數(shù)據(jù)及相關(guān)信息顯示在軟件上，主要為影像的顯示和統(tǒng)計曲線的顯示兩部分。影像的顯示包含影像數(shù)據(jù)轉(zhuǎn)為彩色圖像、彩色圖像的移動與縮放、彩色圖顏色欄的顯示等。統(tǒng)計曲線的顯示包含曲線的展繪、曲線上某數(shù)據(jù)點信息的顯示、曲線自適應(yīng)標(biāo)簽等。此外，還有數(shù)據(jù)信息的顯示如影像某一點的數(shù)據(jù)信息顯示等。

（8） 自動監(jiān)測與警報模塊。在結(jié)果達到閾值時向外發(fā)送聲音、短信、郵件等警報信息。

3 案例分析

3.1 數(shù)據(jù)來源

本文使用河北馬蘭莊礦區(qū)邊坡監(jiān)測數(shù)據(jù)進行分析，數(shù)據(jù)獲取所使用的設(shè)備為中國地質(zhì)大學(xué)（北京）與北京理工大學(xué)共同研制的地基MIMO 成像雷達系統(tǒng)，共計382景影像，采樣間隔3 min。

3.2 結(jié)果分析

3.2.1 軟件運行時間

測試軟件的電腦配置為Intel Core i7-4810MQ CPU@2.80 GHz，內(nèi)存32 GB。對382 景SLC 影像數(shù)據(jù)處理時間進行記錄并成圖，如圖3 所示?？梢钥闯?，單張數(shù)據(jù)處理時間全都在25 s 以內(nèi)，小于3min 采樣間隔，速度較快。該軟件處理效率與監(jiān)測時刻大氣環(huán)境的復(fù)雜性有關(guān)，對于多變的大氣環(huán)境需要多次迭代，耗費時間較長，反之在20s內(nèi)即可完成計算。

3.2.2 軟件運行結(jié)果

圖4 為獲取的累計形變圖，形變區(qū)域明顯，深色部分相對于周圍淺色部分有著明顯不同，淺色部分數(shù)值表示形變值大概在0 mm 附近，而深色部分對應(yīng)-30 mm 以下，能夠直觀地看出較大形變區(qū)域大概范圍和所在的位置。

圖5 位單點統(tǒng)計曲線。累積形變曲線前期變化較快，但變化率不斷減小，后期趨近于直線，比較平緩；形變速率曲線能夠直觀地看出該區(qū)域前期有滑坡現(xiàn)象但后期逐漸穩(wěn)定的態(tài)勢，速率趨勢與累積形變趨勢對應(yīng)，算法在時間線上無明顯錯誤。

3.2.3 精度分析

在分析結(jié)果精度時，由于第三方同步監(jiān)測數(shù)據(jù)的缺失，無法進行外符合精度評定。假設(shè)區(qū)域穩(wěn)定點在監(jiān)測時段內(nèi)計算得到的累積形變量理想情況下應(yīng)當(dāng)為零。因此可以使用監(jiān)測區(qū)域穩(wěn)定點的累積形變曲線是否穩(wěn)定，即殘差的標(biāo)準(zhǔn)差是否較小，來判斷處理結(jié)果精度的高低。殘差為干涉值減去模型值的剩余部分，模型值由兩部分構(gòu)成，基于距離函數(shù)生成的大氣模型與對形變進行時間擬合的形變模型。由于地基SAR 系統(tǒng)監(jiān)測精度為mm 級，可以認為實時處理中數(shù)據(jù)標(biāo)準(zhǔn)差小于1 mm的點精度尚可。

在穩(wěn)定點區(qū)域4個角落（圖4中A、B、C和D）選取5×5大小共100個點進行數(shù)據(jù)殘差的標(biāo)準(zhǔn)差計算，如圖6所示。

由圖6 可知，殘差的標(biāo)準(zhǔn)差小于1 mm 的較高精度的點所占百分比即為77%。與后處理相比，實時處理不能對整體數(shù)據(jù)在處理時進行參數(shù)的調(diào)節(jié)，較高精度穩(wěn)定點占比未能達到接近100% 的水準(zhǔn)，精度有所損失，但在該精度下能夠滿足對累積形變曲線走勢和滑坡預(yù)警的判斷。

4 結(jié) 論

針對當(dāng)前GPRI-II 設(shè)備配套的GAMMA 軟件不具備實時處理及可視化的問題，在滑坡災(zāi)害監(jiān)測預(yù)警方面應(yīng)用存在缺陷，提出一種地基干涉雷達影像數(shù)據(jù)的實時處理流程，并在GAMMA 軟件的基礎(chǔ)上，結(jié)合監(jiān)測預(yù)警的實際需要，采用Python語言進行設(shè)計和二次開發(fā)，編制了地基干涉雷達實時處理軟件，該軟件具有良好的人機交互界面，可以實現(xiàn)地基干涉雷達數(shù)據(jù)預(yù)處理到形變監(jiān)測產(chǎn)品制作全過程。以河北馬蘭莊礦區(qū)邊坡監(jiān)測數(shù)據(jù)為例進行分析，實時處理時單景數(shù)據(jù)處理時間優(yōu)于25s，小于數(shù)據(jù)采集時間；監(jiān)測結(jié)果精度優(yōu)于1mm，滿足形變監(jiān)測預(yù)警的需求。該軟件的研制能夠彌補GPRI-II 在監(jiān)測應(yīng)用中的不足，同時為拓展地基干涉雷達技術(shù)提供了技術(shù)支持。