地基干涉雷達(dá)實(shí)時(shí)處理軟件設(shè)計(jì)與開發(fā)*

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

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

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

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

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

1 原 理

1.1 DInSAR技術(shù)

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

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

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

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

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

1.2 時(shí)序InSAR技術(shù)

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

1.3 形變預(yù)警原理

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

1.4 實(shí)時(shí)處理流程

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

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

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

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

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

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

軟件開發(fā)采用Python 語言，使用Pyside2 庫（Qt for Python）并結(jié)合如Numpy 等Python 眾多第三方庫在GAMMA軟件的基礎(chǔ)上進(jìn)行二次開發(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等多種腳本語言完成編寫，要想成功運(yùn)行GAMMA軟件，除要含有GAMMA軟件本體外，還應(yīng)當(dāng)含有對(duì)這些腳本語言的環(huán)境支持。由于本文軟件開發(fā)使用Python語言，在文件打包時(shí)便已經(jīng)包含Python語言環(huán)境支持，而且sh、csh、tcsh、bash 等都是Unix 或Linux 系統(tǒng)的命令交互語言，因此，在win 系統(tǒng)上需安裝對(duì)Unix/Linux 命令及Perl 語言的支持，在Unix 和Linux僅需安裝Perl語言支持即可。

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

2.2 軟件組成

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

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

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

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

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

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

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

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

（8） 自動(dòng)監(jiān)測(cè)與警報(bào)模塊。在結(jié)果達(dá)到閾值時(shí)向外發(fā)送聲音、短信、郵件等警報(bào)信息。

3 案例分析

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

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

3.2 結(jié)果分析

3.2.1 軟件運(yùn)行時(shí)間

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

3.2.2 軟件運(yùn)行結(jié)果

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

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

3.2.3 精度分析

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

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

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

4 結(jié) 論

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