毛亞純， 曹 旺， 趙占國， 徐茂林

(1. 東北大學(xué) 資源與土木工程學(xué)院， 遼寧 沈陽 110819； 2. 中國黃金集團， 北京 100000；3. 遼寧科技大學(xué) 土木工程學(xué)院， 遼寧 鞍山 114053)

近年來，我國各類地質(zhì)災(zāi)害事件時有發(fā)生，如何為易發(fā)生地質(zhì)災(zāi)害區(qū)域的實時預(yù)警提供準確形變監(jiān)測數(shù)據(jù)已成為一個亟待解決的問題.地基合成孔徑雷達(GB-SAR)是集成了步進連續(xù)波、合成孔徑雷達干涉測量、PS-InSAR等技術(shù)為一體的新型地面形變監(jiān)測設(shè)備.由于GB-SAR具有全天候、高精度等形變監(jiān)測的技術(shù)優(yōu)勢，已被廣泛應(yīng)用于露天礦邊坡、大壩、橋梁、公路、冰川等形變監(jiān)測領(lǐng)域[1-2].

雷達數(shù)據(jù)處理結(jié)果的準確性是露天礦邊坡滑坡危害預(yù)警的重要保障.在雷達數(shù)據(jù)處理的研究中， Iglésias等[3]針對陡峭山區(qū)的雷達數(shù)據(jù)大氣相位誤差改正提出了多元回歸模型且模型高度相關(guān)，但該模型所用氣象數(shù)據(jù)源是沿山路近地表環(huán)境獲取，不適合露天礦坑環(huán)境；曲世勃等[4]分析了GB-SAR形變監(jiān)測誤差來源并驗證了GB-SAR形變監(jiān)測數(shù)據(jù)精度；徐亞明等[5]利用永久散射體改正網(wǎng)對氣象誤差進行改正.雖然上述方法對GB-SAR形變監(jiān)測誤差都有不同程度改正效果，但均忽略了方位向的軌道誤差[6]對GB-SAR形變監(jiān)測精度的影響.本文以馬蘭莊露天鐵礦GB-SAR原始影像為數(shù)據(jù)源，基于小周期處理方法首先對雷達原始影像進行干涉處理生成干涉圖，然后對干涉圖進行小波去噪處理，并采用三重閾值算法提取了雷達影像中的PS點和高質(zhì)量PS點.其中利用PS點進行形變分析，而高質(zhì)量PS點由于其在時間序列上所表現(xiàn)出的高穩(wěn)定性特征(理論干涉相位為零)，因此可認為其干涉相位是由于大氣誤差相位和軌道誤差相位共同影響所致.本文通過對馬蘭莊GB-SAR形變監(jiān)測誤差的規(guī)律分析，提出了基于高質(zhì)量PS點的相位與像元坐標建立多元回歸模型進行相位誤差改正的方法，結(jié)果表明所建模型可準確改正 GB-SAR形變監(jiān)測誤差，為邊坡的穩(wěn)定性分析提供了準確形變數(shù)據(jù)源.

1 GB-SAR原始影像采集及處理流程

本次形變監(jiān)測使用型號為LKR-05-KU國產(chǎn)地基合成孔徑雷達，其方位分辨率在1 km處為4 m，距離分辨率為0.3 m.地基合成孔徑雷達被安置在露天采場東坡和西坡之間的穩(wěn)定區(qū)域，如圖1所示，并對采場東坡進行時長為60 h的形變監(jiān)測，數(shù)據(jù)采集周期設(shè)置為12 min，共采集了300幅原始影像.本文采用小周期處理算法，其處理流程如圖2所示.

2 GB-SAR數(shù)據(jù)處理及形變分析

2.1 GB-SAR數(shù)據(jù)干涉及濾波處理

由于本實驗GB-SAR觀測距離為250～1 200 m，且在觀測過程中GB-SAR安置位置及掃描軌道參數(shù)均未發(fā)生改變，因此用于干涉計算的兩幅原始影像數(shù)據(jù)有嚴格的配準精度，不存在空間基線解算等問題.在干涉計算過程中只需將兩幅原始影像矩陣共軛相乘，即為兩幅影像的干涉圖，設(shè)回波矩陣中任意信號為I1，對應(yīng)相鄰影像中相同位置的回波信號為I2，則干涉信號I12為

(1)

干涉相位Δφ為干涉信號I12對應(yīng)的相位角，其包含雷達視線形變相位φlos、大氣誤差相位φatm、頻率偏移誤差相位φf、軌道誤差相位φtrack及噪聲誤差相位φnoise.故干涉相位Δφ可表示為

Δφ=φlos+φatm+φf+φtrack+φnoise.

(2)

由于雷達天線工作頻率比較穩(wěn)定，且在1 km處由于頻率偏移引起的誤差量級為10-4mm，故頻率偏移誤差相位φf可以不予考慮.

噪聲誤差相位主要體現(xiàn)在干涉圖中的斑點噪聲.由小波去噪原理可知，噪聲在頻域內(nèi)多出現(xiàn)在圖像的高頻處，圖像的紋理特征多出現(xiàn)在圖像的低頻處.因此本文利用小波去噪算法[7]對干涉相位進行了三層小波分解，并基于Stein的無偏似然估計自適應(yīng)閾值算法對干涉相位進行濾波處理.通過處理前后效果對比，該去噪方法可有效剔除干涉圖中的噪聲誤差相位φnoise，且能保持干涉圖的紋理特征.因此，本研究所用干涉圖都經(jīng)過小波去噪處理.

經(jīng)過小波去噪處理后的GB-SAR形變監(jiān)測誤差主要來自于大氣誤差相位和軌道誤差相位.因此文中將對上述兩項誤差的剔除方法進行論述.

2.2 GB-SAR數(shù)據(jù)PS點選取

由于影像中的部分像元不是邊坡穩(wěn)定散射體的反射信號，故需要提取穩(wěn)定PS點來反映真實的邊坡形變量[8-9].此外，為了準確建立誤差改正模型，還需要提取少量高穩(wěn)定性的強散射點，稱之為高質(zhì)量PS點.本文利用雷達原始影像計算了平均幅度、幅度離差指數(shù)和平均相干系數(shù)三個參數(shù)，提出了基于三重閾值算法來篩選特定PS點的方法.

平均幅度M是幅度在整個時間序列的均值，其值越大，反射信號越強，其計算公式為

(3)

(4)

式中：Mt為單幅影像幅度均值；m和n分別為影像方位向和距離向尺寸；A(i,j)為邊坡反射信號；k為影像數(shù)量.

幅度離差指數(shù)D是衡量一個像元內(nèi)邊坡在時間序列內(nèi)對雷達信號反射強度的穩(wěn)定性，其值越小代表越穩(wěn)定.計算公式為

(5)

式中：σi，j為像元的幅度標準差；Mi，j為像元的幅度均值.

平均相干系數(shù)L是雷達數(shù)據(jù)中每個像元在整個時間序列內(nèi)相干系數(shù)的均值，其值越接近1，相干性越好，其計算公式為

(6)

平均幅度、幅度離差指數(shù)和平均相干系數(shù)的計算結(jié)果如圖3所示，其中平均幅度在雷達圖像上用雷達反射強度表示.通過以下步驟選取PS點與高質(zhì)量PS點.

步驟1 設(shè)置幅度均值閾值M1=1.5M，并提取像元平均幅度Mi，j>M1的像元作為初始PS點，再提取像元平均幅度Mi，j>4M1的像元作為初始高質(zhì)量PS點.

步驟2 設(shè)置提取PS點的幅度離差閾值D1=0.5，提取高質(zhì)量PS點的閾值D2=0.1，并基于初始PS點提取幅度離差Di，j

步驟3 設(shè)置相干系數(shù)低閾值L1=0.85，高閾值L2=0.98，并基于過渡PS點提取平均相干系數(shù)Li，j>L1的像元作為最終PS點，基于過渡高質(zhì)量PS點提取平均相干系數(shù)Li，j>L2的像元作為最終高質(zhì)量PS點.

2.3 大氣誤差相位和軌道誤差相位改正

由于微波中的 Ku波段在傳播過程中對空氣介質(zhì)的氣象數(shù)據(jù)比較敏感，且空氣介質(zhì)在兩次采集的時段內(nèi)會發(fā)生一定變化，故干涉相位經(jīng)過濾波后，還包含一定的大氣延遲誤差相位[10-14].由于雷達每次的運行狀態(tài)和軌道會存在一定的偏差，故干涉相位中還包含軌道誤差相位.其中軌道誤差長期存在，大氣延遲誤差只在大氣環(huán)境變化劇烈的情況下影響較大，上述兩種誤差是干涉圖誤差的主要來源.干涉圖誤差相位分布如圖4所示.

假設(shè)微波在露天礦坑中傳播的空氣介質(zhì)均勻，設(shè)主影像的平均大氣延遲系數(shù)為α1，副影像的平均大氣延遲系數(shù)為α2，雷達天線距離PS點的觀測距離為Ld，則兩次觀測的大氣延遲誤差φatm可表示為

φatm=α2·Ld-α1·Ld.

(7)

由于雷達距離向分辨率為0.3 m，故每個像元在距離向?qū)?yīng)的長度為0.3 m.設(shè)雷達起始采樣距離為Ls，故Ld可以用距離向的像元坐標x表示為

Ld=0.3x+Ls.

(8)

由雷達影像的方位向分布特征可知，軌道誤差與方位向像元坐標y呈較好的線性分布，因此軌道誤差φtrack可表示為

φtrack=β·y+ζ.

(9)

故濾波后的GB-SAR形變監(jiān)測誤差相位φGB-SAR與像元坐標x，y可建模為

(10)

式中：α為距離向改正系數(shù)；β為方位向改正系數(shù)；ζ為回歸常數(shù).

利用上述模型對干涉圖誤差相位進行改正，得到真實的形變相位，如圖5所示.

為了驗證模型的準確性，本文針對具有較大誤差的6幅干涉圖進行了分析，分別計算了影像改正前后的平均相位誤差，以及基于多元回歸模型的擬合優(yōu)度，如表1所示.

表1 影像誤差及擬合優(yōu)度統(tǒng)計表

多元回歸模型中自變量和因變量的擬合優(yōu)度越高，則誤差修正精度越高.由表1可以看出，6幅影像中像元坐標和誤差相位的擬合優(yōu)度均大于0.9，說明兩者高度相關(guān).經(jīng)過改正后平均相位誤差明顯變小，平均相位誤差降低了83.3%.

2.4 邊坡形變計算及分析

由于實驗數(shù)據(jù)采集周期為12 min，因此可認為一個小周期內(nèi)的邊坡形變量小于λ/4,即4.5 mm，故可認為處理后的干涉相位即為真實形變相位.根據(jù)雷達參數(shù)可以推導(dǎo)出相位差和形變量之間的關(guān)系[15]為

(11)

式中：ΔL為雷達視線方向的形變量；φlos為雷達視線上的形變相位；c為光速；f為雷達信號發(fā)射頻率.

本文基于Matlab編寫了GB-SAR原始數(shù)據(jù)處理程序，實現(xiàn)了GB-SAR原始數(shù)據(jù)的自動處理.計算出了每個小周期內(nèi)邊坡在雷達視線方向上的形變量，并對形變圖中空缺的像元(距離PS點小于6 m)進行插值，最后通過在時間序列上疊加計算分析了10，20，30，40，50，60 h共6個階段的形變結(jié)果，如圖6所示.

由圖6可以看出，在60 h內(nèi)邊坡逐漸出現(xiàn)了一個明顯的形變區(qū)域.通過對圖中形變區(qū)域?qū)?yīng)的形變數(shù)據(jù)進行統(tǒng)計分析，該區(qū)域內(nèi)最大形變量為-53.243 mm，形變區(qū)域大于10 mm的點共3 099個像元，形變面積約1 860 m2，平均形變量為-21.806 mm，平均形變速率為-0.36 mm/h.

通過對GB-SAR形變監(jiān)測數(shù)據(jù)與馬蘭莊露天礦DEM數(shù)據(jù)進行三維配準，確定了形變區(qū)域位置，經(jīng)過現(xiàn)場勘查發(fā)現(xiàn)在圖6f藍色區(qū)域?qū)?yīng)的實際位置出現(xiàn)了一個寬度為7 cm左右的裂縫，如圖7所示.由此可見利用本文提出的方法處理露天礦GB-SAR監(jiān)測數(shù)據(jù)可以為邊坡的穩(wěn)定性分析提供準確形變數(shù)據(jù)源.

3 結(jié) 論

1) 基于三重閾值算法可以準確提取用于形變分析和誤差改正的PS點和高質(zhì)量PS點.

2) 露天礦GB-SAR形變監(jiān)測數(shù)據(jù)的誤差主要來源于大氣延遲誤差和軌道誤差.基于多元回歸模型的誤差改正方法可以準確改正露天礦GB-SAR形變監(jiān)測的誤差相位，依據(jù)本文實驗數(shù)據(jù)可使GB-SAR形變監(jiān)測誤差降低83.3%.

3) 基于本文對露天礦GB-SAR數(shù)據(jù)的小周期處理方法可以準確提取形變區(qū)域及形變量，為邊坡的穩(wěn)定性分析提供準確形變數(shù)據(jù)源.