朱嘉鑫,李永建,肖 剛,余東洋,付振華,張志家,朱晨浩,戴穎超
(1.北京理工大學 信息與電子學院,北京 100081;2.四川省國土空間生態修復與地質災害防治研究院,四川 成都 610000;3.重慶市地質礦產勘查開發局107 地質隊,重慶 401120;4.重慶市一零七市政建設工程有限公司,重慶 401120;5.蘇州理工雷科傳感技術有限公司,江蘇 蘇州 215000)
我國處于歐亞大陸板塊和太平洋板塊的交接地帶,各種地質災害發生頻繁,尤其在廣大西南地區,每年因地質災害造成的死亡和失蹤人數約占自然災害的1/3[1]。據自然資源部《全國地質災害通報》,2019年全國共發生地質災害6 181 起,其中滑坡4 220起、崩塌1 238 起,災害共造成211 人死亡,直接經濟損失達27.7 億元,因此加強重大滑坡災害的防災減災工作已經成為促進經濟建設和構建和諧社會的國家需求[2-4]。
地基干涉雷達是一種新型的形變測量工具,具有采樣周期短、空間分辨率高、精度高、全天時全天候作業的優勢[5],可無接觸的獲取被測區域的表面形變信息,已經在形變監測領域得到了廣泛的應用,常應用于監測山體邊坡、露天礦坑、水壩、冰川、建筑物等[6-7]。根據雷達數據獲取方式的不同,地基干涉雷達可以分為連續測量模式和間斷測量模式。在連續測量模式下,雷達圖像的獲取周期在分鐘量級,可對目標區域實現近實時的形變測量,該模式常應用于高滑坡風險的礦區邊坡監測、山體二次滑坡監測等[8-10]。但實際中存在著大量的滑坡隱患點,其形變速率很慢,每年僅發生幾厘米的形變。如果采用連續測量模式,將帶來很高的觀測成本和很大的數據冗余。在間斷測量模式下,可以根據形變區域的運動學規律設定采集周期,每隔數周或者數月,重新安裝設備并進行監測,然后基于間斷測量技術,獲取目標區域的表面形變信息[11]。
地基干涉雷達把同一目標區域在不同時間獲取的SAR 復圖像結合起來,比較不同時刻目標的相位差,獲得監測目標的毫米級精度位移信息[12]。地基干涉雷達監測周期較短,相鄰2 次測量時重軌誤差可以忽略不計,適用于重復軌道干涉測量模型。地基干涉雷達差分干涉處理流程是通過對SAR 圖像進行高相干點選擇、干涉相位圖生成、相位解纏和誤差相位補償后,即可以實現形變分析[13-14]。
間斷測量與連續測量最大的不同在于,需要定期的將設備拆下重裝,在實測環境下,無法保證重復安裝時,設備位于完全相同的位置,因此雷達中心與合成孔徑方向會發生偏移。由于地基干涉雷達的方位角分辨率在若干mrad(毫弧度)量級,軌道的輕微偏移就會導致圖像之間發生嚴重的失相干。在進行差分干涉測量前,需要保證2 幅圖像中相同位置的像素點,完全對應著地面的同一個目標點,因此對間斷測量圖像處理前,必須進行圖像配準處理。
在間斷測量中,利用地基雷達對同一個場景進行多次間斷性觀測,每次獲取到數十幅雷達圖像,利用圖像的相干合成,可以提高信噪比,降低干涉相位誤差。對單次測量的數十幅雷達圖像進行幅度合成和相位合成,最終獲取1 幅合成圖。通過圖像合成,具有強反射特性的巖體、弱反射特性的植被和無反射區域之間的區分度會十分明顯。相位合成后,干涉條紋的變化更加均勻,噪聲點也得到了抑制。圖像合成流程圖如圖1。

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

圖2 圖像配準方法流程圖
分析間斷測量時雷達采集的2 組數據,對這2組數據進行圖像合成。配準前以每一個像素點為中心構建l1×l2的矩形框計算相干系數。第k 個像素點的相關系數γk計算公式為:

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

圖3 匹配窗、搜索窗示意圖
圖3 中小矩形代表雷達圖像1 中的匹配窗,大矩形代表雷達圖像2 中的搜索窗。沿著行向和列向,在搜索窗中順序移動匹配窗,每個位置計算1 次相關系數,然后根據相干系數最大值所出現的位置,確定匹配窗的粗偏移量。
對于每個參考點,進行相同處理。在實現所有參考點的粗偏移量的計算后,以出現概率最高的行、列偏移量為基準,對圖像2 進行裁切。經過粗配準后,提取相干系數較大的像素點作為參考點,以每個參考點為中心點,在2 個圖像中設定矩形窗,確定相干系數最大時的行和列偏移量。
1.4.3 精配準
圖像中所有像素點的偏移量為:

式中:△xk為第k 個高相干點在行上的配準偏移量;△yk為第k 個高相干點在列上的配準偏移量;xk為行坐標;yk為列坐標;a1、a2、a3、b1、b2、b3為待估計的配準系數,可采用最小二乘法估算獲得。
根據圖像中所有像素點的行坐標和列坐標,能夠得到實現偏移量的估計值。基于高相干點的偏移量對圖像2 進行插值,獲取精配準后的圖像。
將配準后的雷達圖像進行差分干涉處理,即可以實現形變分析。2020 年3 月至2020 年11 月,利用地基雷達在監測點進行了10 次間斷測量并對雷達圖像進行圖像合成和配準,對10 幅配準后的圖像進行差分干涉處理,平均相干系數圖如圖4。

圖4 平均相干系數圖
重慶市彭水縣聯合鄉馬巖為1 處危巖體區域。該危巖成帶狀,總長約3.9 km。陡崖帶上發育有規模不等的25 個危巖單體,危巖體總體積68.77 萬m3。根據陡崖帶延伸方向、危巖體發育、崩落情況、保護對象等將馬巖危巖區域分為:A 區、B 區、C 區。使用1 臺雷達對3 個區域進行循環間斷測量。
通過對監測區域進行實地勘察,確定A、B、C 3個監測區域雷達的部署位置,依次命名:INSAR01、INSAR02、INSAR03。INSAR01坐標為北緯29°35′40.46",東經108°25′52.29";INSAR02 坐標為北緯29°35′54.18",東經108°25′45.77";INSAR03 坐標為北緯29°36′27.47",東經108°25′43.54"。
在每個監測點布設一定數量的角反射器。每次間斷測量后,利用全站儀獲取角反射器的位移數據。最后將全站儀獲取的位移數據與雷達的監測數據進行對比,判斷SAR 間斷測量方法的有效性。
1)A 區監測數據。2020 年3 月24 日至2020 年9 月23 日,雷達非連續獲取4 組累計180 幅雷達圖像。A 區的圖像信息和全站儀采集的位移數據見表1。A 區未經過圖像配準CBO1 角反射器的形變曲線如圖5,A 區經圖像配準CBO1 角反射器的形變曲線如圖6,A 區全站儀采集的CBO1 位移數據曲線如圖7。

圖5 A 區未經圖像配準CBO1 角反射器的形變曲線

圖6 A 區經過圖像配準CBO1 角反射器的形變曲線

圖7 A 區全站儀采集的CBO1 位移數據曲線

表1 A 區的圖像信息和全站儀采集的位移數據
2)B 區監測數據。2020 年3 月18 日至2020 年9 月24 日,雷達非連續獲取了6 組累計238 幅雷達圖像。B 區的圖像信息和全站儀采集的位移數據見表2。B 區未經圖像配準CBO7 角反射器的形變曲線如圖8,B 區經過圖像配準CBO7 角反射器的形變曲線如圖9,B 區全站儀采集的CBO7 位移數據曲線如圖10。

圖8 B 區未經圖像配準CBO7 角反射器的形變曲線

圖9 B 區經過圖像配準CBO7 角反射器的形變曲線

圖10 B 區全站儀采集的CBO7 位移數據曲線

表2 B 區的圖像信息和全站儀采集的位移數據
3)C 區監測數據。從2020 年3 月18 日至2020年9 月24 日,雷達非連續獲取了6 組、累計238 幅雷達圖像。C 區的圖像信息和全站儀采集的位移數據見表3。C 區未經圖像配準CBO12 角反射器的形變曲線如圖11,C 區經過圖像配準CBO12 角反射器的形變曲線如圖12,C 區全站儀采集的CBO12 位移數據曲線如圖13。

表3 C 區的圖像信息和全站儀采集的位移數據

圖11 C 區未經圖像配準CBO12 角反射器的形變曲線

圖12 C 區經過圖像配準CBO12 角反射器的形變曲線

圖13 C 區全站儀采集的CBO12 位移數據曲線
綜合A 區、B 區和C 區的形變圖和形變曲線可知,未經圖像配準的形變區域沒有延續性,經過圖像配準后的形變區域具有延續性;從曲線圖像對比后可以看出,受到間斷測量的影響,未經圖像配準時形變結果會有大的突變,經過圖像配準后的形變曲線有著較好的連續性。對比雷達和全站儀的監測結果發現,經過圖像配準后的形變數據結果和趨勢與全站儀的結果更為接近。
基于地基差分干涉測量原理,提出了一種適用于間斷測量模式的形變處理方法。利用地基干涉雷達對一處危巖體區域進行了長時間觀測,對不同時段獲取的雷達圖像進行圖像合成與圖像配準,并對配準后的圖像做差分干涉處理。實測數據結果表明,此方法能夠有效地解決雷達圖像之間的失相干問題,處理后的形變區域具有較好的延續性。通過對比全站儀與地基干涉雷達的測量結果,經過圖像配準后的形變趨勢與全站儀的結果更為接近。