基于SAR 圖像配準的地基干涉雷達間斷測量技術(shù)

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

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

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

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

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

1.1 干涉測量原理

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

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

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

1.3 圖像合成

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

圖1 圖像合成流程圖

1.4 圖像配準

1.4.1 基本流程

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

圖2 圖像配準方法流程圖

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

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

1.4.2 粗配準

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

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

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

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

1.4.3 精配準

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

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

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

1.5 差分干涉處理

將配準后的雷達圖像進行差分干涉處理，即可以實現(xiàn)形變分析。2020 年3 月至2020 年11 月，利用地基雷達在監(jiān)測點進行了10 次間斷測量并對雷達圖像進行圖像合成和配準，對10 幅配準后的圖像進行差分干涉處理，平均相干系數(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坐標為北緯29°35′40.46"，東經(jīng)108°25′52.29"；INSAR02 坐標為北緯29°35′54.18"，東經(jīng)108°25′45.77"；INSAR03 坐標為北緯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)過圖像配準CBO1 角反射器的形變曲線如圖5，A 區(qū)經(jīng)圖像配準CBO1 角反射器的形變曲線如圖6，A 區(qū)全站儀采集的CBO1 位移數(shù)據(jù)曲線如圖7。

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

圖6 A 區(qū)經(jīng)過圖像配準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)圖像配準CBO7 角反射器的形變曲線如圖8，B 區(qū)經(jīng)過圖像配準CBO7 角反射器的形變曲線如圖9，B 區(qū)全站儀采集的CBO7 位移數(shù)據(jù)曲線如圖10。

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

圖9 B 區(qū)經(jīng)過圖像配準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)圖像配準CBO12 角反射器的形變曲線如圖11，C 區(qū)經(jīng)過圖像配準CBO12 角反射器的形變曲線如圖12，C 區(qū)全站儀采集的CBO12 位移數(shù)據(jù)曲線如圖13。

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

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

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

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

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

3 結(jié)語

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