朱嘉鑫，李永建，肖 剛，余東洋，付振華，張志家，朱晨浩，戴穎超

（1.北京理工大學(xué) 信息與電子學(xué)院，北京 100081；2.四川省國土空間生態(tài)修復(fù)與地質(zhì)災(zāi)害防治研究院，四川 成都 610000；3.重慶市地質(zhì)礦產(chǎn)勘查開發(fā)局107 地質(zhì)隊，重慶 401120；4.重慶市一零七市政建設(shè)工程有限公司，重慶 401120；5.蘇州理工雷科傳感技術(shù)有限公司，江蘇 蘇州 215000）

我國處于歐亞大陸板塊和太平洋板塊的交接地帶，各種地質(zhì)災(zāi)害發(fā)生頻繁，尤其在廣大西南地區(qū)，每年因地質(zhì)災(zāi)害造成的死亡和失蹤人數(shù)約占自然災(zāi)害的1/3[1]。據(jù)自然資源部《全國地質(zhì)災(zāi)害通報》，2019年全國共發(fā)生地質(zhì)災(zāi)害6 181 起，其中滑坡4 220起、崩塌1 238 起，災(zāi)害共造成211 人死亡，直接經(jīng)濟損失達27.7 億元，因此加強重大滑坡災(zāi)害的防災(zāi)減災(zāi)工作已經(jīng)成為促進經(jīng)濟建設(shè)和構(gòu)建和諧社會的國家需求[2-4]。

地基干涉雷達是一種新型的形變測量工具，具有采樣周期短、空間分辨率高、精度高、全天時全天候作業(yè)的優(yōu)勢[5]，可無接觸的獲取被測區(qū)域的表面形變信息，已經(jīng)在形變監(jiān)測領(lǐng)域得到了廣泛的應(yīng)用，常應(yīng)用于監(jiān)測山體邊坡、露天礦坑、水壩、冰川、建筑物等[6-7]。根據(jù)雷達數(shù)據(jù)獲取方式的不同，地基干涉雷達可以分為連續(xù)測量模式和間斷測量模式。在連續(xù)測量模式下，雷達圖像的獲取周期在分鐘量級，可對目標(biāo)區(qū)域?qū)崿F(xiàn)近實時的形變測量，該模式常應(yīng)用于高滑坡風(fēng)險的礦區(qū)邊坡監(jiān)測、山體二次滑坡監(jiān)測等[8-10]。但實際中存在著大量的滑坡隱患點，其形變速率很慢，每年僅發(fā)生幾厘米的形變。如果采用連續(xù)測量模式，將帶來很高的觀測成本和很大的數(shù)據(jù)冗余。在間斷測量模式下，可以根據(jù)形變區(qū)域的運動學(xué)規(guī)律設(shè)定采集周期，每隔數(shù)周或者數(shù)月，重新安裝設(shè)備并進行監(jiān)測，然后基于間斷測量技術(shù)，獲取目標(biāo)區(qū)域的表面形變信息[11]。

1 地基干涉雷達基本原理和流程

1.1 干涉測量原理

地基干涉雷達把同一目標(biāo)區(qū)域在不同時間獲取的SAR 復(fù)圖像結(jié)合起來，比較不同時刻目標(biāo)的相位差，獲得監(jiān)測目標(biāo)的毫米級精度位移信息[12]。地基干涉雷達監(jiān)測周期較短，相鄰2 次測量時重軌誤差可以忽略不計，適用于重復(fù)軌道干涉測量模型。地基干涉雷達差分干涉處理流程是通過對SAR 圖像進行高相干點選擇、干涉相位圖生成、相位解纏和誤差相位補償后，即可以實現(xiàn)形變分析[13-14]。

1.2 間斷測量的技術(shù)難點

間斷測量與連續(xù)測量最大的不同在于，需要定期的將設(shè)備拆下重裝，在實測環(huán)境下，無法保證重復(fù)安裝時，設(shè)備位于完全相同的位置，因此雷達中心與合成孔徑方向會發(fā)生偏移。由于地基干涉雷達的方位角分辨率在若干mrad（毫弧度）量級，軌道的輕微偏移就會導(dǎo)致圖像之間發(fā)生嚴(yán)重的失相干。在進行差分干涉測量前，需要保證2 幅圖像中相同位置的像素點，完全對應(yīng)著地面的同一個目標(biāo)點，因此對間斷測量圖像處理前，必須進行圖像配準(zhǔn)處理。

1.3 圖像合成

在間斷測量中，利用地基雷達對同一個場景進行多次間斷性觀測，每次獲取到數(shù)十幅雷達圖像，利用圖像的相干合成，可以提高信噪比，降低干涉相位誤差。對單次測量的數(shù)十幅雷達圖像進行幅度合成和相位合成，最終獲取1 幅合成圖。通過圖像合成，具有強反射特性的巖體、弱反射特性的植被和無反射區(qū)域之間的區(qū)分度會十分明顯。相位合成后，干涉條紋的變化更加均勻，噪聲點也得到了抑制。圖像合成流程圖如圖1。

圖1 圖像合成流程圖

1.4 圖像配準(zhǔn)

1.4.1 基本流程

圖像配準(zhǔn)是指將不同時間、不同視角及不同成像條件下獲取的雷達圖像進行匹配和疊加的過程，是干涉圖產(chǎn)生的基礎(chǔ)，對形變測量結(jié)果精準(zhǔn)與否起著重要作用。通過配準(zhǔn)處理，每個場景目標(biāo)點在主輔圖像中對應(yīng)著相同的像素點，減小了重軌誤差對干涉相位的影響。由于間斷測量的時間間隔較大，為減小時間去相干對圖像配準(zhǔn)的影響，經(jīng)過組內(nèi)圖像合成后，采用對相鄰的2 幅圖像進行配準(zhǔn)。配準(zhǔn)流程分為粗配準(zhǔn)和精配準(zhǔn)，圖像配準(zhǔn)方法的流程圖如圖2。

圖2 圖像配準(zhǔn)方法流程圖

分析間斷測量時雷達采集的2 組數(shù)據(jù)，對這2組數(shù)據(jù)進行圖像合成。配準(zhǔn)前以每一個像素點為中心構(gòu)建l1×l2的矩形框計算相干系數(shù)。第k 個像素點的相關(guān)系數(shù)γk計算公式為：

式中：l1×l2為矩形窗的大?。粸閳D像1、圖像2 中以像素點k 為中心的矩形窗。

1.4.2 粗配準(zhǔn)

圖像粗配準(zhǔn)包括參考點選擇、粗偏移量計算、圖像裁切，以雷達圖像1 為基準(zhǔn)，選擇信噪比較高的像素點作為參考點。分別以每一個參考點為中心點，在雷達圖像1 和雷達圖像2 中分別選擇不同大小的匹配窗和搜索窗，匹配窗、搜索窗示意圖如圖3。

圖3 匹配窗、搜索窗示意圖

圖3 中小矩形代表雷達圖像1 中的匹配窗，大矩形代表雷達圖像2 中的搜索窗。沿著行向和列向，在搜索窗中順序移動匹配窗，每個位置計算1 次相關(guān)系數(shù)，然后根據(jù)相干系數(shù)最大值所出現(xiàn)的位置，確定匹配窗的粗偏移量。

對于每個參考點，進行相同處理。在實現(xiàn)所有參考點的粗偏移量的計算后，以出現(xiàn)概率最高的行、列偏移量為基準(zhǔn)，對圖像2 進行裁切。經(jīng)過粗配準(zhǔn)后，提取相干系數(shù)較大的像素點作為參考點，以每個參考點為中心點，在2 個圖像中設(shè)定矩形窗，確定相干系數(shù)最大時的行和列偏移量。

1.4.3 精配準(zhǔn)

圖像中所有像素點的偏移量為：

式中：△xk為第k 個高相干點在行上的配準(zhǔn)偏移量；△yk為第k 個高相干點在列上的配準(zhǔn)偏移量；xk為行坐標(biāo)；yk為列坐標(biāo)；a1、a2、a3、b1、b2、b3為待估計的配準(zhǔn)系數(shù)，可采用最小二乘法估算獲得。

根據(jù)圖像中所有像素點的行坐標(biāo)和列坐標(biāo)，能夠得到實現(xiàn)偏移量的估計值。基于高相干點的偏移量對圖像2 進行插值，獲取精配準(zhǔn)后的圖像。

1.5 差分干涉處理

將配準(zhǔn)后的雷達圖像進行差分干涉處理，即可以實現(xiàn)形變分析。2020 年3 月至2020 年11 月，利用地基雷達在監(jiān)測點進行了10 次間斷測量并對雷達圖像進行圖像合成和配準(zhǔn)，對10 幅配準(zhǔn)后的圖像進行差分干涉處理，平均相干系數(shù)圖如圖4。

圖4 平均相干系數(shù)圖

2 監(jiān)測案例

重慶市彭水縣聯(lián)合鄉(xiāng)馬巖為1 處危巖體區(qū)域。該危巖成帶狀，總長約3.9 km。陡崖帶上發(fā)育有規(guī)模不等的25 個危巖單體，危巖體總體積68.77 萬m3。根據(jù)陡崖帶延伸方向、危巖體發(fā)育、崩落情況、保護對象等將馬巖危巖區(qū)域分為：A 區(qū)、B 區(qū)、C 區(qū)。使用1 臺雷達對3 個區(qū)域進行循環(huán)間斷測量。

通過對監(jiān)測區(qū)域進行實地勘察，確定A、B、C 3個監(jiān)測區(qū)域雷達的部署位置，依次命名：INSAR01、INSAR02、INSAR03。INSAR01坐標(biāo)為北緯29°35′40.46"，東經(jīng)108°25′52.29"；INSAR02 坐標(biāo)為北緯29°35′54.18"，東經(jīng)108°25′45.77"；INSAR03 坐標(biāo)為北緯29°36′27.47"，東經(jīng)108°25′43.54"。

在每個監(jiān)測點布設(shè)一定數(shù)量的角反射器。每次間斷測量后，利用全站儀獲取角反射器的位移數(shù)據(jù)。最后將全站儀獲取的位移數(shù)據(jù)與雷達的監(jiān)測數(shù)據(jù)進行對比，判斷SAR 間斷測量方法的有效性。

1）A 區(qū)監(jiān)測數(shù)據(jù)。2020 年3 月24 日至2020 年9 月23 日，雷達非連續(xù)獲取4 組累計180 幅雷達圖像。A 區(qū)的圖像信息和全站儀采集的位移數(shù)據(jù)見表1。A 區(qū)未經(jīng)過圖像配準(zhǔn)CBO1 角反射器的形變曲線如圖5，A 區(qū)經(jīng)圖像配準(zhǔn)CBO1 角反射器的形變曲線如圖6，A 區(qū)全站儀采集的CBO1 位移數(shù)據(jù)曲線如圖7。

圖5 A 區(qū)未經(jīng)圖像配準(zhǔn)CBO1 角反射器的形變曲線

圖6 A 區(qū)經(jīng)過圖像配準(zhǔn)CBO1 角反射器的形變曲線

圖7 A 區(qū)全站儀采集的CBO1 位移數(shù)據(jù)曲線

表1 A 區(qū)的圖像信息和全站儀采集的位移數(shù)據(jù)

2）B 區(qū)監(jiān)測數(shù)據(jù)。2020 年3 月18 日至2020 年9 月24 日，雷達非連續(xù)獲取了6 組累計238 幅雷達圖像。B 區(qū)的圖像信息和全站儀采集的位移數(shù)據(jù)見表2。B 區(qū)未經(jīng)圖像配準(zhǔn)CBO7 角反射器的形變曲線如圖8，B 區(qū)經(jīng)過圖像配準(zhǔn)CBO7 角反射器的形變曲線如圖9，B 區(qū)全站儀采集的CBO7 位移數(shù)據(jù)曲線如圖10。

圖8 B 區(qū)未經(jīng)圖像配準(zhǔn)CBO7 角反射器的形變曲線

圖9 B 區(qū)經(jīng)過圖像配準(zhǔn)CBO7 角反射器的形變曲線

圖10 B 區(qū)全站儀采集的CBO7 位移數(shù)據(jù)曲線

表2 B 區(qū)的圖像信息和全站儀采集的位移數(shù)據(jù)

3）C 區(qū)監(jiān)測數(shù)據(jù)。從2020 年3 月18 日至2020年9 月24 日，雷達非連續(xù)獲取了6 組、累計238 幅雷達圖像。C 區(qū)的圖像信息和全站儀采集的位移數(shù)據(jù)見表3。C 區(qū)未經(jīng)圖像配準(zhǔn)CBO12 角反射器的形變曲線如圖11，C 區(qū)經(jīng)過圖像配準(zhǔn)CBO12 角反射器的形變曲線如圖12，C 區(qū)全站儀采集的CBO12 位移數(shù)據(jù)曲線如圖13。

表3 C 區(qū)的圖像信息和全站儀采集的位移數(shù)據(jù)

圖11 C 區(qū)未經(jīng)圖像配準(zhǔn)CBO12 角反射器的形變曲線

圖12 C 區(qū)經(jīng)過圖像配準(zhǔn)CBO12 角反射器的形變曲線

圖13 C 區(qū)全站儀采集的CBO12 位移數(shù)據(jù)曲線

綜合A 區(qū)、B 區(qū)和C 區(qū)的形變圖和形變曲線可知，未經(jīng)圖像配準(zhǔn)的形變區(qū)域沒有延續(xù)性，經(jīng)過圖像配準(zhǔn)后的形變區(qū)域具有延續(xù)性；從曲線圖像對比后可以看出，受到間斷測量的影響，未經(jīng)圖像配準(zhǔn)時形變結(jié)果會有大的突變，經(jīng)過圖像配準(zhǔn)后的形變曲線有著較好的連續(xù)性。對比雷達和全站儀的監(jiān)測結(jié)果發(fā)現(xiàn)，經(jīng)過圖像配準(zhǔn)后的形變數(shù)據(jù)結(jié)果和趨勢與全站儀的結(jié)果更為接近。

3 結(jié)語

基于地基差分干涉測量原理，提出了一種適用于間斷測量模式的形變處理方法。利用地基干涉雷達對一處危巖體區(qū)域進行了長時間觀測，對不同時段獲取的雷達圖像進行圖像合成與圖像配準(zhǔn)，并對配準(zhǔn)后的圖像做差分干涉處理。實測數(shù)據(jù)結(jié)果表明，此方法能夠有效地解決雷達圖像之間的失相干問題，處理后的形變區(qū)域具有較好的延續(xù)性。通過對比全站儀與地基干涉雷達的測量結(jié)果，經(jīng)過圖像配準(zhǔn)后的形變趨勢與全站儀的結(jié)果更為接近。