曾 濤 鄧云開* 胡 程 田衛(wèi)明

①(北京理工大學(xué)信息與電子學(xué)院雷達(dá)技術(shù)研究所 北京 100081)

②(北京理工大學(xué)衛(wèi)星導(dǎo)航電子信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室 北京 100081)

1 引言

受到自然因素(降雨、融雪、地震等)和人為因素(采礦、地下水枯竭、植被破壞等)的影響，地質(zhì)災(zāi)害在世界范圍內(nèi)頻繁發(fā)生，每年均會(huì)造成嚴(yán)重的經(jīng)濟(jì)損失和大量的人員傷亡。2018年，我國(guó)共發(fā)生地質(zhì)災(zāi)害2966起，造成105人死亡、7人失蹤、73人受傷，直接經(jīng)濟(jì)損失達(dá)到14.7億元(自然資源部)?；聻?zāi)害是地質(zhì)災(zāi)害中發(fā)生頻率最高和危害最大的一種，如2018年11月份在川藏交界發(fā)生的金沙江滑坡堰塞湖災(zāi)害，共造成10.2萬(wàn)人受災(zāi)，3400余間房屋倒塌，農(nóng)作物受災(zāi)面積3.5千公頃，沿江部分地區(qū)道路、橋梁、電力等基礎(chǔ)設(shè)施損失較為嚴(yán)重。

一般而言，在滑坡造成的重大事故中，露天礦邊坡滑坡、山體滑坡、尾礦壩邊坡滑坡和大壩滑坡的發(fā)生最為頻繁，造成的損失也最大?；碌恼T發(fā)因素很多，在邊坡宏觀失穩(wěn)之前，均會(huì)發(fā)生應(yīng)力的改變，其表面通常會(huì)出現(xiàn)形變[1，2]。因此，為了更深入地研究邊坡滑坡的觸發(fā)機(jī)理，并實(shí)現(xiàn)對(duì)滑坡災(zāi)害的預(yù)測(cè)預(yù)警，國(guó)內(nèi)外學(xué)者開展了大量邊坡表面形變測(cè)量方面的研究。

根據(jù)形變測(cè)量過(guò)程中，測(cè)量?jī)x器是否接觸被測(cè)目標(biāo)區(qū)域，形變測(cè)量技術(shù)可以分為兩大類：接觸式測(cè)量和非接觸式測(cè)量。接觸式測(cè)量主要包含水準(zhǔn)儀測(cè)量、全站儀測(cè)量、差分GPS測(cè)量等。該項(xiàng)技術(shù)的最大優(yōu)點(diǎn)是成本低，操作便捷，但是只能對(duì)場(chǎng)景中的一些離散點(diǎn)進(jìn)行測(cè)量，難以滿足對(duì)邊坡進(jìn)行全覆蓋形變監(jiān)測(cè)的需求，而且接觸式測(cè)量需要在目標(biāo)區(qū)域進(jìn)行測(cè)量點(diǎn)布設(shè)，在一些危險(xiǎn)目標(biāo)區(qū)域很難實(shí)施[3]。非接觸式測(cè)量主要包括激光掃描儀測(cè)量、星載SAR(Synthetic Aperture Radar，合成孔徑雷達(dá))干涉測(cè)量、地基差分干涉測(cè)量等。激光掃描儀測(cè)量可以獲取較高的形變測(cè)量精度，但光學(xué)遙感頻率很高，信號(hào)波長(zhǎng)短，對(duì)邊坡進(jìn)行長(zhǎng)時(shí)間、連續(xù)形變監(jiān)測(cè)時(shí)，易受雨、雪、霧等氣象條件的影響[4]。星載SAR干涉測(cè)量技術(shù)，可以實(shí)現(xiàn)全天時(shí)全天候的大范圍監(jiān)測(cè)，但現(xiàn)階段其觀測(cè)實(shí)時(shí)性容易受到衛(wèi)星重訪周期的限制，且觀測(cè)角度難以靈活選擇，對(duì)大梯度邊坡測(cè)量時(shí)，易受疊掩、陰影等問(wèn)題的影響[5]。地基差分干涉測(cè)量在近些年取得了迅速的發(fā)展。

地基差分干涉測(cè)量雷達(dá)，由于其工作平臺(tái)在地面上，可以在目標(biāo)區(qū)域幾十米到幾公里外進(jìn)行觀測(cè)，有利于靈活地選擇布設(shè)地點(diǎn)和觀測(cè)視角，且工作在微波波段，成像時(shí)不受天氣條件的影響，具有全天時(shí)全天候的優(yōu)點(diǎn)[6]。地基差分干涉雷達(dá)的圖像獲取速度很快，一般為幾分鐘，有利于對(duì)目標(biāo)區(qū)域進(jìn)行連續(xù)實(shí)時(shí)觀測(cè)，已經(jīng)在形變監(jiān)測(cè)領(lǐng)域得到了廣泛的應(yīng)用，如建筑物、露天礦坑、山體邊坡、水壩、冰川等的監(jiān)測(cè)[7]。

2 典型系統(tǒng)

地基差分干涉測(cè)量雷達(dá)，目前多工作在X波段或者Ku波段，系統(tǒng)組成部分主要包括收發(fā)天線、供電模塊、數(shù)據(jù)采集和存儲(chǔ)單元、數(shù)據(jù)處理模塊等。按照雷達(dá)成像模式的不同，可以分為兩種類型：地基RAR(Real Aperture Radar，實(shí)孔徑雷達(dá))和地基SAR。

圖1 典型地基RARFig.1 Typical GB-RAR systems

2.1 地基RAR

地基RAR中，代表性的系統(tǒng)有澳大利亞Ground Probe公司的SSR(Slope Stability Radar)系統(tǒng)、瑞士Gamma公司的GPRI(Gamma Portable Radar Instrument)系統(tǒng)和南非REUTECH MINING公司的MSR(Movement and Surveying Radar)系統(tǒng)，參見圖1。

SSR系統(tǒng)早期的兩種型號(hào)SSR-XT/MT，均采用大孔徑的拋物面天線來(lái)發(fā)射波束寬度極窄的鉛筆狀波束，然后通過(guò)高精度的伺服系統(tǒng)的控制，實(shí)現(xiàn)方位維和俯仰維的大范圍逐點(diǎn)掃描，該系統(tǒng)可以直接將雷達(dá)圖像與3維地形相匹配。以SSR-XT為例，其工作在X波段，監(jiān)測(cè)范圍為30～3500 m，1 km處分辨率為8.7 m×8.7 m，掃描180°×60°范圍耗時(shí)26 min， 85°×20°范圍4 min， 30°×15°范圍2 min[8]。MSR系統(tǒng)同樣基于實(shí)孔徑技術(shù)，采用大孔徑的拋物面天線，通過(guò)方位維和俯仰維的大范圍掃描來(lái)獲取3維點(diǎn)云，掃描120°×45°范圍耗時(shí)少于4 min，60°×80°范圍少于3 min。MSR系統(tǒng)包括4種型號(hào)：MSR60， MSR120， MSR250和MSR400，以型號(hào)MSR400為例，其最小測(cè)量距離為30 m，最大測(cè)量距離為4 km。在1 km處，距離向、方位向和高度向的分辨率分別為0.50 m， 4.40 m和0.44 m[9]。

SSR系統(tǒng)新型的兩種型號(hào)SSR-FX/OMNI，均采用線性天線，同樣可以在方位維和俯仰維實(shí)現(xiàn)大范圍的掃描。以SSR-OMNI為例，其天線長(zhǎng)度為2.74 m，最大監(jiān)測(cè)距離可達(dá)5600 m，可以實(shí)現(xiàn)方位維360°、俯仰維60°的掃描，僅耗時(shí)2 min， 1 km處分辨率為4.30 m×0.68 m。GPRI系統(tǒng)是安裝有旋轉(zhuǎn)掃描儀的FM-CW(Frequency-Modulated Continuous-Wave)雷達(dá)干涉儀，采用2.06 m長(zhǎng)的線性波導(dǎo)天線來(lái)發(fā)射波束寬度在方位維為0.4°、高度維為35.0°的扇形波束，掃描速率為每秒0.5°到10.0°[10]。以GPRI-II系統(tǒng)為例，其通過(guò)天線在方位維的大范圍掃描來(lái)獲取2維雷達(dá)圖像，距離向分辨率約為0.75 m，方位向分辨率為6.80 m。該系統(tǒng)的工作頻率范圍為17.1～17.3 GHz，測(cè)量范圍為50 m～10 km，并配備1根發(fā)射天線和2根接收天線形成一個(gè)垂直基線為25 cm的干涉陣列，從而基于干涉測(cè)高原理來(lái)獲取觀測(cè)區(qū)域的3維地形，并可實(shí)現(xiàn)形變量的3維可視化顯示。

圖2 典型直線掃描地基SARFig.2 Typical linear-scanning GB-SAR systems

2.2 地基SAR

按照合成孔徑的實(shí)現(xiàn)方式，可以將地基SAR劃分為以下3種類型：

(1)直線掃描地基SAR。直線掃描地基SAR是通過(guò)收發(fā)天線沿著高精密滑軌的移動(dòng)來(lái)獲取方位維的大合成孔徑，從而實(shí)現(xiàn)方位維的高分辨[11]。第1款商用的直線掃描地基SAR是由意大利IDS公司和佛羅倫薩大學(xué)聯(lián)合開發(fā)的IBIS(Image By Interferometric Survey)系統(tǒng)[12]。國(guó)內(nèi)外很多研究機(jī)構(gòu)和公司，開發(fā)了多款直線掃描體制的地基SAR系統(tǒng)，代表性的有意大利IDS公司的IBIS-FM系統(tǒng)、歐盟JRC的LiSA系統(tǒng)[13]、荷蘭MetaSensing公司的FastGBSAR-S系統(tǒng)[14]、西班牙UPC大學(xué)的RiskSAR系統(tǒng)[15]、中國(guó)安全生產(chǎn)科學(xué)研究院(安科院)的S-SAR系統(tǒng)[16]、北京理工雷科電子信息技術(shù)有限公司(理工雷科)的邊坡形變監(jiān)測(cè)系統(tǒng)[17]、內(nèi)蒙古自治區(qū)方向圖公司的微變監(jiān)測(cè)雷達(dá)LSA系統(tǒng)等(參見圖2)。此外，北方工業(yè)大學(xué)、日本Tohoku大學(xué)、英國(guó)Sheffield大學(xué)等也開展了各自地基SAR系統(tǒng)的研究。各款直線掃描地基SAR的工作原理相似，但在一些工作參數(shù)上有所不同，如距離向分辨率、方位向分辨率、測(cè)量范圍、掃描時(shí)間等，如表1所示。

表1 直線掃描地基SAR參數(shù)表Tab.1 Parameters of linear-scanning GB-SAR

圖3 典型弧線掃描地基SARFig.3 Typical arc-scanning GB-SAR systems

以IBIS-FM系統(tǒng)為例，其工作在Ku波段，最遠(yuǎn)探測(cè)距離為4.5 km， 1 km處的空間分辨率為0.5 m×4.4 m，快速掃描時(shí)間約為3 min，可以獲取亞毫米量級(jí)的形變測(cè)量精度。在SAR圖像合成及處理階段，其將線性調(diào)頻連續(xù)波技術(shù)、合成孔徑雷達(dá)技術(shù)、干涉測(cè)量技術(shù)和永久散射體技術(shù)相結(jié)合，可以應(yīng)用于對(duì)礦山邊坡、水電站大壩、冰川等的實(shí)時(shí)形變監(jiān)測(cè)。

(2)弧線掃描地基SAR?；【€掃描地基SAR通過(guò)收發(fā)天線在水平面內(nèi)的圓周運(yùn)動(dòng)來(lái)進(jìn)行圓弧掃描，從而獲取大的合成孔徑。由于弧掃描地基SAR采用特殊的運(yùn)動(dòng)形式來(lái)實(shí)現(xiàn)圓弧式合成孔徑，在成像算法上，與直線掃描地基SAR有所不同[18，19]。代表性的系統(tǒng)如韓國(guó)國(guó)立江原大學(xué)開發(fā)的ArcSAR(Arc-scanning SAR)[20]、意大利IDS公司的IBIS-ArcSAR系統(tǒng)[21]、中國(guó)科學(xué)院電子學(xué)研究所的Arc FMCW-SAR系統(tǒng)[22]和內(nèi)蒙古方向圖公司的微變監(jiān)測(cè)旋轉(zhuǎn)雷達(dá)RSA系統(tǒng)等(參見圖3)。

韓國(guó)ArcSAR系統(tǒng)，其有2種成像模式：聚束模式和掃描模式。在聚束模式下，可以獲取比常規(guī)線掃描地基SAR更高的方位向分辨率，在掃描模式下，可以實(shí)現(xiàn)近360°的大范圍掃描。意大利IBISArcSAR系統(tǒng)，采用4只基于MIMO技術(shù)的天線，可以實(shí)現(xiàn)自動(dòng)地理編碼，實(shí)時(shí)獲取場(chǎng)景數(shù)字地形模型，掃描360°范圍時(shí)只需要40 s，最大測(cè)量距離可達(dá)5 km，單臺(tái)系統(tǒng)便可實(shí)現(xiàn)對(duì)大型礦坑的完全覆蓋。

(3)MIMO地基SAR。直線掃描和弧線掃描兩種工作體制，均需要伺服系統(tǒng)控制收發(fā)天線進(jìn)行特定的機(jī)械掃描，圖像獲取的速度會(huì)受限，一般為幾分鐘至十幾分鐘。為減少圖像的獲取時(shí)間，提高地基雷達(dá)在快速形變監(jiān)測(cè)領(lǐng)域的應(yīng)用，國(guó)內(nèi)外一些機(jī)構(gòu)開展了MIMO(Multiple-Input Multiple-Output，多輸入多輸出)體制地基SAR的研究。地基MIMO雷達(dá)采用多輸入多輸出技術(shù)，通過(guò)多個(gè)發(fā)射天線和接收天線的特殊排列來(lái)等效成一個(gè)大的合成孔徑。由于穩(wěn)定的正交波形設(shè)計(jì)這個(gè)難題暫未得到有效解決，地基MIMO雷達(dá)工作時(shí)，各個(gè)發(fā)射天線分時(shí)發(fā)射，各個(gè)接收天線則同時(shí)接收，一次完整的掃描時(shí)間為幾毫秒到幾秒。代表性的地基MIMO雷達(dá)，如歐盟JRC機(jī)構(gòu)的MELISSA系統(tǒng)[23]、北京理工大學(xué)(北理工)研發(fā)的MIMO-SAR系統(tǒng)[24]等(參見圖4)。

以北理工MIMO-SAR系統(tǒng)為例，其采用16個(gè)發(fā)射天線構(gòu)成兩個(gè)密集子陣列和16個(gè)接收天線構(gòu)成一個(gè)稀疏子陣列，可以等效成一個(gè)擁有256個(gè)采樣點(diǎn)的合成孔徑。該系統(tǒng)工作在Ku波段，波長(zhǎng)為λ=1.86cm ，每個(gè)密集子陣列中相鄰發(fā)射天線的距離間隔為 λ/2=0.93cm，稀疏子陣列中相鄰接收天線的距離間隔為 8 λ=7.44cm。系統(tǒng)的等效合成孔徑為1.138 m，角分辨率為0.466°，測(cè)量距離范圍為30 m～3 km。

圖4 典型MIMO地基SARFig.4 Typical MIMO GB-SAR systems

圖5 部分地基RARFig.5 Some GB-RAR systems

上述地基RAR系統(tǒng)和SAR系統(tǒng)主要應(yīng)用于對(duì)較大范圍測(cè)量(雷達(dá)測(cè)量面積在1 km2左右)的場(chǎng)景，如礦區(qū)邊坡、山體、水壩等的形變監(jiān)測(cè)。各系統(tǒng)均可以實(shí)現(xiàn)方位向和距離向的高分辨(分米或者米量級(jí))，并基于差分干涉測(cè)量技術(shù)實(shí)現(xiàn)高精度的形變測(cè)量。一些地基RAR系統(tǒng)，如荷蘭MetaSensing公司的FastGBSAR-R系統(tǒng)、意大利IDS公司的IBIS-FS系統(tǒng)等(參見圖5)，均采用將收發(fā)天線固定在三腳架上的結(jié)構(gòu)，工作時(shí)不進(jìn)行方位維和俯仰維的掃描，通過(guò)對(duì)單一目標(biāo)體進(jìn)行可達(dá)上百赫茲的高頻率觀測(cè)，可以實(shí)現(xiàn)對(duì)橋梁、高塔、大樓等的振動(dòng)測(cè)量，測(cè)振精度可達(dá)0.01 mm。一些地基SAR系統(tǒng)，如圖4所示的兩款MIMO地基SAR，由于圖像獲取速度快，耗時(shí)在毫秒量級(jí)，也可以選擇出一些像素點(diǎn)來(lái)進(jìn)行振動(dòng)分析。在可控實(shí)驗(yàn)條件下，MELISSA系統(tǒng)的圖像獲取速度優(yōu)于4 ms，基于振動(dòng)角反的測(cè)量結(jié)果表明，其測(cè)振精度優(yōu)于10 μm[25]。

3 技術(shù)現(xiàn)狀

對(duì)于不同工作體制的地基差分干涉雷達(dá)，其處理技術(shù)上會(huì)有一定的差別。本節(jié)以應(yīng)用于對(duì)較大范圍場(chǎng)景進(jìn)行形變監(jiān)測(cè)的地基SAR系統(tǒng)為例，介紹地基差分干涉雷達(dá)的技術(shù)現(xiàn)狀。雖然3種類型的地基SAR：直線掃描、弧線掃描和MIMO，在SAR技術(shù)，即高分辨成像處理上有著較大的差別，但差分干涉處理流程相似[26，27]。

在地基SAR差分干涉處理技術(shù)上，國(guó)內(nèi)外很多學(xué)者提出了不同的形變處理算法，雖然處理流程上會(huì)有一些差別，但主要的處理技術(shù)相同，包括差分干涉、PS點(diǎn)選擇、相位解纏、大氣相位補(bǔ)償、形變量解算和地理編碼等5個(gè)步驟[28]。如果有N幅地基SAR圖像，可以先進(jìn)行PS點(diǎn)選擇，然后將第1幅作為主圖像，其他N-1幅作為輔圖像，經(jīng)過(guò)差分干涉處理，獲取N-1幅差分干涉圖?；谶x擇出的PS點(diǎn)，對(duì)這N-1幅差分干涉圖進(jìn)行相位解纏和大氣相位補(bǔ)償處理，即可以實(shí)現(xiàn)形變量的解算，處理流程如圖6所示。如果地基SAR應(yīng)用于實(shí)時(shí)形變處理，考慮到地基SAR的圖像獲取速度較快，一般是每隔幾分鐘即可獲取1幅SAR圖像，1天時(shí)間即可以獲取上百幅SAR圖像，為保證形變測(cè)量的實(shí)時(shí)性，需要采用不同的差分干涉處理流程[29]。

圖6 地基SAR差分干涉處理流程Fig.6 Differential interferometric scheme of GB-SAR

3.1 差分干涉

地基SAR可以實(shí)現(xiàn)2維高分辨成像，在距離維上，通常采用調(diào)頻連續(xù)波技術(shù)或者步進(jìn)頻技術(shù)，在方位維上，則采用合成孔徑技術(shù)。圖像中每一個(gè)像素點(diǎn)均是復(fù)數(shù)，其幅度通常用來(lái)解譯成像場(chǎng)景及研究散射特性，相位則可以用來(lái)獲取目標(biāo)區(qū)域的高程或者形變信息。地基SAR進(jìn)行形變測(cè)量時(shí)，雷達(dá)位置固定不動(dòng)，不同圖像之間的空間基線為零，對(duì)兩幅圖像進(jìn)行對(duì)應(yīng)像素的復(fù)共軛相乘，即可以實(shí)現(xiàn)差分干涉處理。圖7所示為地基SAR差分干涉測(cè)量原理示意圖。理想情況下，差分干涉相位 Δ φ與視線方向的形變量ΔR線性相關(guān)，關(guān)系可以表示為

其中，λ表示信號(hào)波長(zhǎng)。

圖7 差分干涉原理示意圖Fig.7 Schematic diagram of differential interferometry

實(shí)際中，受到各種誤差源的影響，差分干涉相位Δ φ中還包含其他分量，可以建模為

其中，φdefo為形變相位分量；φatm為兩幅圖像獲取期間由氣象條件改變所導(dǎo)致的大氣相位分量；φgeom為重軌誤差所導(dǎo)致的幾何相位分量，一般可以忽略不計(jì)； φnoi為像素點(diǎn)散射特性改變及系統(tǒng)熱噪聲等帶來(lái)的誤差相位分量，在差分干涉處理后，可以經(jīng)過(guò)干涉相位濾波處理來(lái)進(jìn)行濾除。由于相位周期性的影響，差分干涉相位 Δ φ是纏繞的，處在區(qū)間 [ -,)內(nèi)，k表示相位模糊度，且是一個(gè)整數(shù)。

3.2 PS點(diǎn)選擇

在利用像素點(diǎn)的相位信息進(jìn)行形變測(cè)量時(shí)，差分干涉相位的質(zhì)量直接影響到形變測(cè)量的精度。但受到大氣擾動(dòng)和熱噪聲等非理想因素的影響，對(duì)低相位質(zhì)量的像素點(diǎn)，一方面難以進(jìn)行正確的相位解纏，另一方面對(duì)其進(jìn)行形變分析時(shí)會(huì)出現(xiàn)較大的測(cè)量誤差。因此，地基SAR差分干涉處理時(shí)，通常需要選擇出一些高質(zhì)量的像素點(diǎn)，即為PS(Permanent Scatterer，永久散射體)，來(lái)進(jìn)行形變分析。

在地基SAR領(lǐng)域，廣泛采用幅度離差法來(lái)進(jìn)行PS點(diǎn)的選擇。該方法利用對(duì)一個(gè)像素點(diǎn)的幅度穩(wěn)定性的估計(jì)來(lái)代替對(duì)其相位穩(wěn)定性的估計(jì)，一般而言，至少需要20幅SAR圖像[30]。一個(gè)像素點(diǎn)的幅度離差值(ADI， Amplitude Dispersion Index)DA定義為

其中，σA和mA分別表示該像素點(diǎn)的時(shí)序幅值序列的標(biāo)準(zhǔn)差和均值。對(duì)DA設(shè)置一定的閾值DT，選擇準(zhǔn)則為DA≤DT，即可實(shí)現(xiàn)PS點(diǎn)的選擇。DT的取值范圍一般為0.10～0.25。圖8所示分別為一植被覆蓋山體邊坡的照片和基于30幅圖像計(jì)算出的ADI圖，可以看出，無(wú)植被覆蓋區(qū)域的ADI一般在0.2以下，有植被覆蓋區(qū)域的ADI則在0.5以上，說(shuō)明了PS點(diǎn)一般處在裸露的巖石區(qū)域，植被區(qū)域像素點(diǎn)的幅度穩(wěn)定性很低[31]。

由于地基SAR獲取圖像時(shí)的空間基線為0，時(shí)間基線為幾分鐘到十幾分鐘，一般情況下，基于幅度離差法選擇出的PS點(diǎn)，其密度和質(zhì)量均可以滿足高精度形變測(cè)量的需求。為進(jìn)一步提高PS點(diǎn)的質(zhì)量，可以在幅度離差法的基礎(chǔ)上，基于信雜比、時(shí)序相關(guān)系數(shù)等進(jìn)行二次選擇[32]。在地基SAR應(yīng)用于對(duì)緩慢形變區(qū)域進(jìn)行非連續(xù)監(jiān)測(cè)(時(shí)間基線達(dá)幾十天)，或?qū)τ兄脖桓采w區(qū)域進(jìn)行常規(guī)連續(xù)觀測(cè)時(shí)，如果采用幅度離差法來(lái)選擇PS點(diǎn)，其密度或質(zhì)量可能無(wú)法滿足形變監(jiān)測(cè)需求，可以采用在星載SAR領(lǐng)域常用的其他方法，如StaMPS方法、SqueeSAR方法和PSP方法等[33，34]。

3.3 相位解纏

由于干涉相位圖是纏繞的，為實(shí)現(xiàn)正確的形變測(cè)量，需要對(duì)其進(jìn)行相位解纏，即從式(2)中解出相位模糊度k。常規(guī)干涉SAR的相位解纏是在相鄰像素點(diǎn)之間進(jìn)行的，包括路徑跟蹤法和最小范數(shù)法等，但PS點(diǎn)是離散、非均勻地分布在地基SAR差分干涉圖上的，常規(guī)解纏方法不再適用。為實(shí)現(xiàn)對(duì)PS點(diǎn)的相位解纏，常用非均勻網(wǎng)格下的最小費(fèi)用流算法或者最小二乘方法等[35]。以最小費(fèi)用流算法為例，首先根據(jù)PS點(diǎn)的分布，采用Delaunay三角網(wǎng)生成算法構(gòu)建不規(guī)則三角網(wǎng)，然后計(jì)算每一個(gè)三角形的殘差值，采用最小費(fèi)用流算法計(jì)算相位模糊度，最后在三角網(wǎng)中對(duì)纏繞相位積分，實(shí)現(xiàn)相位解纏。

為提高相位解纏的準(zhǔn)確度，可以在干涉圖中選擇出一些高質(zhì)量參考點(diǎn)或者人工布設(shè)若干個(gè)控制點(diǎn)來(lái)輔助解纏，也可以對(duì)最小費(fèi)用流算法做出適當(dāng)?shù)母倪M(jìn)[36]。Noferini等學(xué)者[37]在采用地基SAR系統(tǒng)進(jìn)行長(zhǎng)時(shí)間基線的非連續(xù)形變監(jiān)測(cè)時(shí)，針對(duì)最小費(fèi)用流解纏算法會(huì)出現(xiàn)較大誤差的問(wèn)題，提出了考慮PS點(diǎn)形變速率的改進(jìn)方案，極大地提高了相位解纏準(zhǔn)確度，對(duì)比結(jié)果如圖9所示。

在采用地基SAR進(jìn)行實(shí)時(shí)形變監(jiān)測(cè)時(shí)，為保證相位解纏的準(zhǔn)確度，一般還會(huì)考慮時(shí)間維的1維解纏，這樣相位解纏轉(zhuǎn)化為兼顧方位維、距離維和時(shí)間維的3維解纏問(wèn)題[38]。在1維時(shí)間維上，可以采用卡爾曼濾波或者歐拉方法來(lái)進(jìn)行解纏。

圖8 典型PS點(diǎn)選擇結(jié)果Fig.8 Typical PS selection results

3.4 大氣相位補(bǔ)償

大氣會(huì)影響到電磁波的傳輸速率及路徑，因此不同時(shí)刻大氣條件(溫度、濕度、大氣壓)的改變會(huì)造成不同的傳輸延遲，從而干涉相位中存在大氣相位分量。大氣相位可以建模為

其中， ΔN表示折射率的變化，其隨時(shí)間t和空間r變化，L表示信號(hào)的傳輸路徑。

在進(jìn)行大氣相位補(bǔ)償時(shí)，一般假設(shè)大氣在空間上是均勻的、在時(shí)間上是隨機(jī)的，則大氣相位φatm呈現(xiàn)出隨斜距而線性變化的趨勢(shì)。對(duì)于Ku波段的地基SAR，在20°C， 1 km的距離上，僅1%的濕度變化就可以帶來(lái)約2 mm的測(cè)量誤差。最基本的補(bǔ)償方法包括氣象數(shù)據(jù)法、控制點(diǎn)校正法和基于PS技術(shù)的參數(shù)模型法[39]。方法1通過(guò)在目標(biāo)區(qū)域布設(shè)氣象站，獲取氣象參數(shù)，包括溫度、濕度和大氣壓，之后根據(jù)大氣折射經(jīng)驗(yàn)?zāi)Ｐ蛯?duì)大氣相位進(jìn)行定量的分析。方法2則是在目標(biāo)區(qū)域內(nèi)人工布設(shè)或者選擇若干個(gè)強(qiáng)散射體目標(biāo)，然后對(duì)這些控制點(diǎn)的干涉相位進(jìn)行分析，通過(guò)插值來(lái)消除其他像素點(diǎn)的大氣相位。在缺少氣象參數(shù)或者外部控制點(diǎn)時(shí)，可以采用方法3，建立合理的大氣相位模型，然后基于PS點(diǎn)的解纏相位迭代估計(jì)模型參數(shù)，進(jìn)一步實(shí)現(xiàn)大氣相位的補(bǔ)償[40]。圖10所示為Huang等學(xué)者[41]采用方法2，對(duì)一干涉相位圖進(jìn)行大氣相位補(bǔ)償前后的結(jié)果，補(bǔ)償后圖像中大部分像素點(diǎn)的相位在0 rad左右，有效地減少了大氣相位對(duì)形變測(cè)量的影響。

在雷達(dá)的觀測(cè)范圍過(guò)大，或者存在較大的高程差異時(shí)，大氣在空間上不再是均勻變化時(shí)，大氣相位隨斜距而線性變化的模型會(huì)存在較大的誤差。Iglesias等人[42]采用地基SAR系統(tǒng)對(duì)一高山進(jìn)行了長(zhǎng)期觀測(cè)，雷達(dá)觀測(cè)區(qū)域的高程差達(dá)400 m，并沿山體道路布設(shè)了多處氣象站，分析得出氣象條件會(huì)隨著斜距和高程而發(fā)生變化。文章提出了兼顧斜距和高程的多參數(shù)模型，如式(5)所示，β1和β2為待估計(jì)參數(shù)，hn和rn分別表示第n個(gè)PS點(diǎn)的高程和斜距，然后建立線性方程組，迭代估計(jì)出β1和 β2。

圖9 相位解纏結(jié)果Fig.9 Phase unwrapping results

圖10 干涉相位圖Fig.10 Phase interferogram

3.5 形變量解算和地理編碼

經(jīng)過(guò)上述處理，可以獲取到形變相位 φdefo，基于式(1)即可以實(shí)現(xiàn)形變量的解算。值得注意的是，雷達(dá)測(cè)量的是1維視線方向的形變量，即為目標(biāo)區(qū)域的真實(shí)形變量在雷達(dá)視線方向的投影分量，需要將形變量在3維地形上進(jìn)行準(zhǔn)確的定位及顯示，即地理編碼問(wèn)題。對(duì)于地基RAR系統(tǒng)，除基于差分干涉進(jìn)行形變測(cè)量外，一般也同時(shí)具備3維地形測(cè)量能力，可以很方便地進(jìn)行地理編碼。對(duì)于地基SAR系統(tǒng)，成像時(shí)是將3維地形在雷達(dá)的2維成像平面內(nèi)進(jìn)行投影，其地理編碼問(wèn)題，可以視為一個(gè)由雷達(dá)2維成像坐標(biāo)系向空間3維直角坐標(biāo)系轉(zhuǎn)換的問(wèn)題。考慮到地基SAR系統(tǒng)一般不具備3維地形測(cè)量能力，為實(shí)現(xiàn)地理編碼，可以使用激光掃描儀來(lái)輔助測(cè)量[43]。經(jīng)過(guò)地理編碼，有利于確定發(fā)生形變的區(qū)域，從而進(jìn)一步開展形變量和形變速率分析，實(shí)現(xiàn)滑坡災(zāi)害的預(yù)測(cè)預(yù)警等[44]。

4 應(yīng)用實(shí)例

地基SAR已經(jīng)在形變監(jiān)測(cè)領(lǐng)域得到了廣泛的應(yīng)用，國(guó)內(nèi)外學(xué)者已經(jīng)發(fā)表了很多篇文章來(lái)闡述地基SAR的應(yīng)用案例。本文以北京理工大學(xué)雷達(dá)技術(shù)研究所開展的3次監(jiān)測(cè)實(shí)驗(yàn)，即馬蘭莊露天開采邊坡監(jiān)測(cè)、貴州納雍滑坡后續(xù)監(jiān)測(cè)和橋梁振動(dòng)測(cè)量實(shí)驗(yàn)為例，展示地基差分干涉雷達(dá)的典型應(yīng)用。

圖11 實(shí)驗(yàn)信息Fig.11 Experimental information

表2 MIMO雷達(dá)參數(shù)表Tab.2 Parameters of the MIMO radar system

4.1 露天開采邊坡監(jiān)測(cè)

馬蘭莊露天開采邊坡(E118°36′， N40°06′)位于河北省遷安市馬蘭莊鎮(zhèn)，該露天礦坑呈橢球形，整個(gè)礦坑長(zhǎng)半軸約1100 m，短半軸約900 m。礦坑邊坡為典型巖質(zhì)邊坡，無(wú)植被覆蓋，最大開采深度大于200 m，邊坡傾角為38°～47°[45]。圖11(a)所示為場(chǎng)景照片，黃色矩形代表雷達(dá)的布放位置，雷達(dá)成像時(shí)的角度范圍設(shè)置為60°，紅色橢圓A和B所示為形變區(qū)域。圖11(b)所示為MIMO雷達(dá)現(xiàn)場(chǎng)布放照片，系統(tǒng)布放在一活動(dòng)板房?jī)?nèi)，避免雨、雪等天氣的影響。

為評(píng)估該邊坡的穩(wěn)定性，采用北理工研發(fā)的MIMO體制地基SAR，對(duì)該邊坡進(jìn)行形變監(jiān)測(cè)，系統(tǒng)參數(shù)如表2所示。監(jiān)測(cè)時(shí)間段為2017年7月19日14時(shí)至2017年7月29日9時(shí)，累計(jì)獲取了2010幅MIMO雷達(dá)圖像。在形變監(jiān)測(cè)階段，每一幅圖像的數(shù)據(jù)獲取時(shí)間加形變處理時(shí)間約為7 min。

圖12(a)所示為極坐標(biāo)系下的成像結(jié)果，邊坡區(qū)域內(nèi)像素點(diǎn)的幅度值分布在-30～0 dB范圍內(nèi)[46]。圖12(b)所示為基于圖像1和圖像2獲取的差分干涉相位圖，邊坡區(qū)域內(nèi)像素點(diǎn)的干涉相位在0 rad左右，邊坡區(qū)域外像素點(diǎn)的干涉相位則隨機(jī)變化[47]。

圖13(a)所示為基于這2010幅MIMO雷達(dá)圖像獲取的累積形變量結(jié)果，在這11天的監(jiān)測(cè)周期內(nèi)，有兩部分區(qū)域A和B呈現(xiàn)出明顯的負(fù)形變量，負(fù)號(hào)代表向著靠近雷達(dá)的方向形變。圖13(b)所示為將形變信息反投到3維立體圖上的結(jié)果，結(jié)合礦坑照片，可以確定形變區(qū)域處于圖11(a)中的區(qū)域A和區(qū)域B。

為了更好地說(shuō)明區(qū)域A和區(qū)域B的形變情況，結(jié)合形變量曲線和形變速率曲線來(lái)進(jìn)行進(jìn)一步的分析。形變量測(cè)量曲線如圖14(a)所示，區(qū)域A和區(qū)域B的最大形變量分別達(dá)到了-15.45 mm和-30.13 mm，且均隨時(shí)間連續(xù)變化。對(duì)形變量測(cè)量曲線進(jìn)行時(shí)域差分及濾波處理，即可以得到形變速率測(cè)量曲線，如圖14(b)所示。區(qū)域A和區(qū)域B的最大形變速率分別達(dá)到了-0.519 mm/h和-2.232 mm/h，且均出現(xiàn)在2017年7月21日上午6時(shí)左右。原因可能是在該時(shí)間點(diǎn)，施工人員對(duì)該礦坑開展了爆破工作，導(dǎo)致這兩個(gè)區(qū)域出現(xiàn)了明顯的形變加速。

圖12 MIMO雷達(dá)圖像與干涉相位圖Fig.12 MIMO radar image and phase interferogram

圖13 形變測(cè)量結(jié)果Fig.13 Deformation measurement results

圖14 形變分析結(jié)果Fig.14 Deformation analysis results

4.2 山體滑坡后續(xù)監(jiān)測(cè)

2017年8月28日上午10點(diǎn)40分，貴州省畢節(jié)市納雍縣張家灣鎮(zhèn)發(fā)生一起較大規(guī)模的山體垮塌，且后續(xù)發(fā)生了若干次小規(guī)?；隆榱藢?duì)滑坡后的山體邊坡開展形變監(jiān)測(cè)，采用理工雷科公司開發(fā)的邊坡形變監(jiān)測(cè)雷達(dá)(直線掃描地基SAR，見圖2(e))，于2017年8月31日21時(shí)至2017年9月8日16時(shí)，對(duì)垮塌殘余體進(jìn)行了不間斷實(shí)時(shí)監(jiān)測(cè)?，F(xiàn)場(chǎng)監(jiān)測(cè)照片如圖15所示[48]。

圖16所示分別為從2017年8月31日21時(shí)開始形變監(jiān)測(cè)，到2017年9月2日16時(shí)、2017年9月4日16時(shí)、2017年9月6日16時(shí)和2017年9月8日16時(shí)的積累形變測(cè)量結(jié)果?？梢钥闯?，隨著監(jiān)測(cè)時(shí)間的增加，部分區(qū)域的形變量明顯增大。

將2017年9月8日累積形變測(cè)量結(jié)果投影到3維地形圖上，如圖17(a)所示，可以看出該山體邊坡上有兩塊區(qū)域(區(qū)域A和區(qū)域B)的形變量較大，結(jié)合形變量曲線和形變速率曲線來(lái)分析區(qū)域A和區(qū)域B的形變情況。可以看出，區(qū)域A和區(qū)域B的形變變化趨勢(shì)很接近，均在2017年9月2日5時(shí)左右和2017年9月5日9時(shí)左右出現(xiàn)明顯的變形加速。結(jié)合現(xiàn)場(chǎng)的天氣條件來(lái)看，很可能是由于降雨導(dǎo)致的變形加速。2017年9月5日上午，區(qū)域B的形變速率達(dá)到了最大值，約-8.2 mm/h，隨后區(qū)域B發(fā)生了崩塌。這也說(shuō)明了基于地基差分干涉雷達(dá)來(lái)進(jìn)行滑坡災(zāi)害預(yù)報(bào)預(yù)警的可行性。

圖15 納雍滑坡災(zāi)后現(xiàn)場(chǎng)監(jiān)測(cè)照片F(xiàn)ig.15 On-site monitoring photo after the Nayong landslide

圖16 累積形變測(cè)量結(jié)果Fig.16 Cumulative deformation measurements

4.3 橋梁振動(dòng)測(cè)量實(shí)驗(yàn)

基于微波干涉原理，地基雷達(dá)可以對(duì)橋梁、高樓、高塔等結(jié)構(gòu)體的振動(dòng)進(jìn)行測(cè)量，其具有快速、高精度、非接觸等優(yōu)點(diǎn)[49]。但振動(dòng)測(cè)量技術(shù)，與較大范圍場(chǎng)景的形變測(cè)量技術(shù)相比，有著較大差別。一般來(lái)說(shuō)，振動(dòng)測(cè)量是對(duì)單一像素點(diǎn)的時(shí)序變化序列進(jìn)行分析，形變測(cè)量則側(cè)重于對(duì)整幅雷達(dá)圖像進(jìn)行分析。振動(dòng)測(cè)量的基本處理流程，主要包括回波信號(hào)處理、強(qiáng)散射點(diǎn)提取、雜波抑制、振動(dòng)點(diǎn)檢測(cè)、振動(dòng)參數(shù)估計(jì)等步驟[50]。

圖17 形變分析結(jié)果Fig.17 Deformation analysis results

現(xiàn)階段國(guó)內(nèi)外應(yīng)用于振動(dòng)測(cè)量領(lǐng)域的雷達(dá)系統(tǒng)主要是采用RAR體制，地基RAR系統(tǒng)僅具有1維距離向分辨率，對(duì)干擾信號(hào)的抑制能力弱，且不利于對(duì)振動(dòng)區(qū)域進(jìn)行準(zhǔn)確定位。地基MIMO雷達(dá)的圖像獲取速度較快，具有2維高分辨能力，有利于對(duì)振動(dòng)區(qū)域進(jìn)行準(zhǔn)確地識(shí)別，且從雷達(dá)體制上提高了干擾信號(hào)的抑制能力，展現(xiàn)了其在振動(dòng)測(cè)量領(lǐng)域進(jìn)行應(yīng)用的潛力。本文采用北理工自主研發(fā)的MIMO雷達(dá)，開展了目標(biāo)對(duì)象分別為角反和橋梁的兩次測(cè)振實(shí)驗(yàn)，初步探討了其應(yīng)用于振動(dòng)測(cè)量的可行性。

MIMO雷達(dá)系統(tǒng)的參數(shù)如表2所示，其發(fā)射信號(hào)為調(diào)頻連續(xù)波，則發(fā)射信號(hào)時(shí)寬等價(jià)于PRT(Pulse Repetition Time， 脈沖重復(fù)周期)。MIMO雷達(dá)工作時(shí)，16個(gè)發(fā)射天線依次分時(shí)發(fā)射，16個(gè)接收天線分為4組，每4個(gè)接收天線為1組，各組依次分時(shí)接收，則MIMO雷達(dá)的一次完整采樣周期對(duì)應(yīng)64個(gè)PRT，即一次振動(dòng)采樣周期為64個(gè)PRT。2次實(shí)驗(yàn)中，均將PRT設(shè)置為0.25 ms，則振動(dòng)采樣周期為16 ms，系統(tǒng)可測(cè)量的最大振動(dòng)頻率為31.25 Hz。

首先利用可裝置角反的振動(dòng)校準(zhǔn)儀來(lái)開展驗(yàn)證實(shí)驗(yàn)。實(shí)驗(yàn)中，調(diào)整角反的振動(dòng)方向，使其沿雷達(dá)視線方向，連續(xù)采集100個(gè)振動(dòng)采樣周期作為1組數(shù)據(jù)，則每組數(shù)據(jù)的采集時(shí)長(zhǎng)為1.6 s。通過(guò)調(diào)整校準(zhǔn)儀的振動(dòng)幅度和振動(dòng)頻率(見表3)，采集了7組實(shí)驗(yàn)數(shù)據(jù)。圖18所示分別為實(shí)驗(yàn)場(chǎng)景和校準(zhǔn)儀照片，角反的棱長(zhǎng)為8 cm，校準(zhǔn)儀與雷達(dá)之間的距離約為9 m[51]。

圖19(a)所示為MIMO雷達(dá)成像結(jié)果。在進(jìn)行振動(dòng)測(cè)量時(shí)，基于每組數(shù)據(jù)的100幅MIMO圖像，首先選擇出一些強(qiáng)散射點(diǎn)，圖中幅值最強(qiáng)點(diǎn)對(duì)應(yīng)角反，然后提取每一個(gè)強(qiáng)散射點(diǎn)的相位序列，進(jìn)行雜波抑制和振動(dòng)參數(shù)提取。圖19(b)所示基于第1組數(shù)據(jù)，獲取的部分強(qiáng)散射點(diǎn)的振動(dòng)頻率測(cè)量結(jié)果。

這7組數(shù)據(jù)的振動(dòng)測(cè)量結(jié)果如表3所示。分析結(jié)果表明，MIMO雷達(dá)的振動(dòng)頻率測(cè)量精度很高，角反的振動(dòng)頻率在10～20 Hz范圍內(nèi)時(shí)，測(cè)量誤差在5% 以內(nèi)。

之后采用MIMO雷達(dá)對(duì)一橋梁進(jìn)行了振動(dòng)測(cè)量，實(shí)驗(yàn)地點(diǎn)選擇在北理工校園西北側(cè)的蘇州橋，實(shí)驗(yàn)時(shí)間為2018年7月20日。圖20所示分別為實(shí)驗(yàn)場(chǎng)景照片和成像結(jié)果。

實(shí)驗(yàn)時(shí)連續(xù)獲取了100幅雷達(dá)圖像，并選擇橋墩進(jìn)行振動(dòng)分析。圖21(a)所示為橋墩處像素點(diǎn)的時(shí)序相位變化序列，在這100幅圖像中呈現(xiàn)出明顯的正弦變化形式。對(duì)該相位序列進(jìn)行相位濾波和傅里葉變換，振動(dòng)頻率提取結(jié)果如圖21(b)所示，圖中峰值點(diǎn)的頻率約為0.63 Hz。本次實(shí)驗(yàn)初步驗(yàn)證了MIMO雷達(dá)對(duì)橋梁振動(dòng)的測(cè)量能力，還需要結(jié)合其他設(shè)備來(lái)驗(yàn)證結(jié)果的可靠性，有待對(duì)數(shù)據(jù)進(jìn)行進(jìn)一步處理。

表3 振動(dòng)測(cè)量結(jié)果Tab.3 Vibrating measurement results

圖18 實(shí)驗(yàn)信息Fig.18 Experimental information

圖19 振動(dòng)分析結(jié)果Fig.19 Vibration analysis results

5 總結(jié)

本文綜述了地基差分干涉雷達(dá)的發(fā)展現(xiàn)狀及應(yīng)用實(shí)例。首先介紹了地基差分干涉雷達(dá)的典型系統(tǒng)，按照成像模式的不同，將其劃分為地基RAR和地基SAR兩種類型，然后按照合成孔徑的實(shí)現(xiàn)方式，將地基SAR劃分為3種類型：直線掃描、弧線掃描和MIMO。之后以地基SAR為例，詳細(xì)介紹了差分干涉形變測(cè)量技術(shù)，包括差分干涉、PS點(diǎn)選擇、相位解纏、大氣相位補(bǔ)償、形變量解算和地理編碼等5個(gè)關(guān)鍵技術(shù)。最后以3次地基差分干涉雷達(dá)監(jiān)測(cè)實(shí)驗(yàn)為例，展現(xiàn)了其在露天邊坡監(jiān)測(cè)、山體滑坡監(jiān)測(cè)和橋梁振動(dòng)測(cè)量方面的應(yīng)用。

圖20 實(shí)驗(yàn)信息Fig.20 Experimental information

圖21 實(shí)驗(yàn)分析結(jié)果Fig.21 Experimental analysis results