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

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

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

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

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

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

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

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

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

(1)

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

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

(2)

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

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

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

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

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

平均幅度M是幅度在整個(gè)時(shí)間序列的均值，其值越大，反射信號(hào)越強(qiáng)，其計(jì)算公式為

(3)

(4)

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

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

(5)

式中：σi，j為像元的幅度標(biāo)準(zhǔn)差；Mi，j為像元的幅度均值.

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

(6)

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

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

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

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

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

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

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

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

(7)

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

Ld=0.3x+Ls.

(8)

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

φtrack=β·y+ζ.

(9)

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

(10)

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

利用上述模型對(duì)干涉圖誤差相位進(jìn)行改正，得到真實(shí)的形變相位，如圖5所示.

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

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

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

2.4 邊坡形變計(jì)算及分析

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

(11)

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

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

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

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

3 結(jié) 論

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

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

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