基于SAR 圖像配準(zhǔn)的地基干涉雷達(dá)間斷測(cè)量技術(shù)

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

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

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

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

1 地基干涉雷達(dá)基本原理和流程

1.1 干涉測(cè)量原理

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

1.2 間斷測(cè)量的技術(shù)難點(diǎn)

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

1.3 圖像合成

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

圖1 圖像合成流程圖

1.4 圖像配準(zhǔn)

1.4.1 基本流程

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

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

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

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

1.4.2 粗配準(zhǔn)

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

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

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

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

1.4.3 精配準(zhǔn)

圖像中所有像素點(diǎn)的偏移量為：

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

根據(jù)圖像中所有像素點(diǎn)的行坐標(biāo)和列坐標(biāo)，能夠得到實(shí)現(xiàn)偏移量的估計(jì)值?；诟呦喔牲c(diǎn)的偏移量對(duì)圖像2 進(jìn)行插值，獲取精配準(zhǔn)后的圖像。

1.5 差分干涉處理

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

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

2 監(jiān)測(cè)案例

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

通過(guò)對(duì)監(jiān)測(cè)區(qū)域進(jìn)行實(shí)地勘察，確定A、B、C 3個(gè)監(jiān)測(cè)區(qū)域雷達(dá)的部署位置，依次命名：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"。

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)語(yǔ)

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