基于光學衛星影像的同震形變提取技術及形變梯度的地學意義

黃 宇 劉小利,2 莫昕宇 鄧德貝爾 阮巧喆 賈治革,2

1 中國地震局地震研究所,武漢市洪山側路40號,430071

2 湖北省地震局,武漢市洪山側路48號,430071

同震地表形變量及其在空間上的分布特征,是理解地殼彈性應變轉換為不可逆的永久性構造變形和構造地震破裂行為的關鍵,也是構造活動區各種重要基礎設施抗震設防的重要依據[1]。目前,獲取同震形變的方法主要有GNSS、InSAR、精密水準、地震波反演、遙感及多源數據聯合反演等。受限于儀器成本和信號傳輸等因素,GNSS、測震臺站多布設于重點監測區,無法精細刻畫地震同震形變的空間分布與變化特征[2]。自1993年首個InSAR同震形變場獲取以來,InSAR技術因其低成本、面狀覆蓋等優勢迅速成為研究同震形變的最重要方法之一。但該技術最大只能監測出相鄰像元間小于半波長的形變[3],且近斷層強震動作用會導致相位失相干,造成近場形變缺失[2]。此外,中強地震多發生于構造活動強烈的造山帶,地形效應對相位解纏影響較大,有望借助深度學習等方法突破瓶頸,但需要海量訓練集的支撐。為此,近年來一些學者利用光學影像中的像素點偏移來快速獲取大范圍中強地震同震形變,以彌補GNSS和InSAR技術的不足[4-5]。本文將重點討論基于光學影像的同震形變處理技術及具可操作性的優化流程,以2021年瑪多地震為例對比兩款主流軟件的差異,基于同震形變場的一階和二階導數推導具有地學意義的形變算子,旨在為光學同震形變技術的深入應用提供參考。

1 光學同震形變處理技術

理論上,同一區域、不同時相(時間基線較短)的光學影像對中相同坐標的像素點應當完全重合。但由于大氣輻射、幾何畸變及地震引起的地表變化等因素,使得不同時相影像對存在像素偏移。因此,在剔除大氣輻射、幾何畸變等影響后,影像對產生的像素偏移可用于描述同震地表形變,利用亞像素相關性匹配技術甚至可檢測1/10亞像素級別的形變[6]。

1.1 亞像素相關性匹配技術

亞像素相關性匹配技術指利用經過輻射校正、正射校正和控制點匹配等預處理后的光學影像對進行亞像素級別的匹配,計算不同時相影像對偏差帶來的位移量[6]。相位相關法[7]能精準計算偏移量估計,利用傅里葉變換將加窗后的光學影像對包含的空間域信息變換到頻率域,獲取圖像的變換關系并計算歸一化互功率譜函數相位信息,經過傅里葉逆變換后的互功率譜相位矩陣最大值即為圖像間的相對偏移量,并采用內插或者奇異值(singular value decomposition, SVD)分解等方法進行擬合估計,從而獲得亞像素級別偏移檢測精度[8]。具體計算如下:

假設fpre(x,y)和fpost(x,y)分別為震前、震后影像的像元點,(Δx,Δy)為同震地表形變量,有:

fpost(x,y)=fpre(x-Δx,y-Δy)

(1)

式中,(x,y)為像元坐標。設Fpost(u,v)為fpost(x,y)的傅里葉變換,Fpre(u,v)為fpre(x,y)的傅里葉變換,將兩幅影像加窗后經過傅里葉變換轉化為頻率域:

Fpost(u,v)=

Fpre(u,v)*exp(-j2π(uΔx+vΔy))

(2)

式中,u、v分別為行、列的頻率變化,j為復數單位。以Fpre(u,v)為參考得到歸一化互功率譜函數H(u,v):

exp(-j2π(uΔx+vΔy))

(3)

式中,上標*為共軛。理論上,歸一化互功率譜矩陣H是一個秩為1的矩陣,因此可將歸一化互功率譜函數H(u,v)分解為:

H(u,v)=exp(-j2πuΔx)exp(-j2πuΔy)=

(4)

式中,上標G為共軛轉置。進行SVD分解,用最大奇異值及相應的奇異向量對歸一化互功率譜矩陣H進行秩1估計:

(5)

式中,σmax為最大奇異值,umax和vmax為對應的奇異向量。將umax和vmax相位解纏,得到相位PΔx和PΔy,擬合線性相位得到平移量:

(6)

式中,kx和ky分別為pΔx和pΔy的斜率。最終,由震前、震后光學影像經過亞像素相關性匹配后即可得到東西向(EW)和南北向(NS)的水平向形變分量。

1.2 光學同震形變提取的優化流程

通過對不同時相的影像進行正射校正、幾何校正、掩膜、地理配準等預處理,可消除幾何畸變、地形、水體等帶來的噪聲。但亞像素匹配后得到的水平向形變場還存在諸多噪聲,不利于定量化的形變分析。例如,不同時相的軌道參數變化導致軌道誤差,多光譜探測器中電耦合器件(charge coupled device,CCD)線陣列錯位排列、抖動震蕩引起波動偽影條帶誤差,地表輻射性能變化造成失相關噪聲,大氣噪聲、時間失相關(區域內隨時間發生的變化)和陰影等與場地相關的因素以及熱波動和信號量化等傳感器因素導致白噪聲和加性噪聲[9]等,但這些因素在已有研究中很少被綜合考慮[10]。為有效消除上述誤差影響,經過實驗對比,提出以下對水平向形變場的精化后處理:

1)失相關區域掩膜。對同震形變場中信噪比(signal-noise ratio,SNR)低于設定經驗閾值的像素點進行掩膜,抑制失相關噪聲點對后續誤差處理的影響。

2)條帶誤差處理。為消除CCD線陣列姿態變化引起的條帶誤差,將每列條帶區域內像元形變值疊加求平均,再由每個像元形變值減去平均值,去除條帶誤差。

3)非局部均值濾波(NL-means)。相對于傳統的中值濾波,實驗選擇非局部均值濾波(NL-means)平滑噪聲、減少突跳點。 NL-means濾波是基于圖像自相似性的鄰域濾波方法,充分利用圖像中的冗余信息,在去噪的同時最大程度地保持圖像的細節特征。在條帶誤差處理后,進行NL-means濾波,可進一步平滑噪聲、突出同震形變場細節特征。其核心是計算圖像中所有像素與當前像素間的相似性,一般會設定兩個固定大小的窗口——一個大的搜索窗口和一個小的鄰域窗口,鄰域窗口在搜索窗口中進行滑動,根據鄰域間的相似性修改權值來確定對應中心像素對當前像素的影響度。

(7)

式中,I為搜索窗口區域;權值w(a,b)為像素點a和b間的相似度,它的值由以a、b為中心的矩形鄰域M(a)、M(b)間的距離‖M(a)-M(b)‖2決定:

(8)

‖M(a)-M(b)‖2=

(9)

式中,M為鄰域窗口區域,控制高斯函數的衰減程度;c為鄰域窗口中的像元點;d為鄰域窗口尺寸;h為平滑系數,h越大高斯函數變化越平緩,去噪水平越高,但同時也會導致圖像越模糊,建議采取小平滑系數的多次濾波策略。

4)軌道誤差處理。為剔除圖像中的多項式趨勢,在非形變區域均勻選取若干個像素點,通過最小二乘平差求得4個參數和多項式軌道誤差,即將一個多項式曲面擬合到圖像上,并將其從圖像中移除[5]。

5)賦空值-中值濾波。當形變場的空值較多時,會在一定程度上影響形變表達和解讀,傳統方法中對此不加考慮。實驗通過中值濾波法將失相關掩膜空值區域賦值,只對失相關區域進行中值濾波,可以在不失真的情況下盡可能減小噪聲影響。

基于精化后處理的水平形變分量,即可進一步定量化分析,提取具有地學意義的形變算子。完整的光學同震形變提取優化流程如圖1所示。

圖1 光學同震形變提取優化流程Fig.1 Optimization process for optical coseismic deformation extraction

1.3 主流光學同震形變軟件對比分析

目前,用于光學地表形變處理的亞像素相關性匹配軟件包括COSI-Corr、MicMac和Medicis等[10]。其中,COSI-Corr和Medicis均采用頻域相關器,但后者去噪效果較差[11],這里不作討論。COSI-Corr、MicMac的差異主要在于像素對偏移檢測相關器的影響、檢測窗口設置、噪聲抑制等方面。

COSI-Corr是一款依賴ENVI平臺的非開源免費軟件。COSI-Corr軟件采用頻域相關器完成像素偏移檢測。首先對固定尺寸的初始窗口內逐像素偏移進行粗略估計,并以此設定截止窗口進行亞像素相關性檢測。當影像噪聲較大或預期形變較大時,初始窗口尺寸設置較大為宜,而截止窗口尺寸應設計成能容納一定噪聲的最小尺寸。COSI-Corr軟件根據對數-交叉頻譜振幅大小對噪聲頻率進行屏蔽,篩選噪聲值,修改每次重新計算的頻率屏蔽次數,減少測量中的噪聲值。

MicMac是完全開源的數字攝影測量和三維重建軟件,采用正則化相關器,通過檢測窗口尺寸和移動步長來調整檢測精度和速度,從而完成亞像素相關性匹配[10];采用非線性相關性系數(式(10))在像素對偏移檢測的同時限制噪聲對形變檢測結果的影響:

Cost=C1*(1-C3)

(10)

C1=1-Cmin

(11)

(12)

(13)

式中,Cor為歸一化相關性系數,Cmin為相關性最小閾值,γ為相關影響程度。

針對相同的影像對數據,兩款軟件相關器窗口大小對噪聲的敏感度稍有差異,如表1所示。對比結果表明,MicMac始終保持較好的去噪效果。

表1 COSI-Corr與MicMac結果差異

2 基于光學同震形變場的定量化形變算子

在數學上,梯度作為矢量表示為函數在某點沿法向的變化。將光學同震形變場看作二維離散函數,同理可進行梯度計算。基于像素偏移的梯度計算可得到位移在空間某一位置沿某一方向的變化量,即位移矢量場。梯度變化最大的方向和變化量指示同震破裂在地表的投影。

圖像梯度算子有多種,與Milliner[11]采用的二階精準中心差分近似,本文選擇具有較強抗噪聲能力的Sobel算子,它通過對空間鄰域加權計算位移梯度,能充分考慮在整個空間鄰域對中心形變像素點求取梯度的影響,還可以降低場地因素和傳感器因素帶來的復合噪聲對梯度的直接影響,對形變反映更加敏感,如式(14)、(15):

?u(x,y)/?x=

(14)

?u(x,y)/?y=

(15)

式中,u為同震位移場,px為同震位移場分辨率,下文同理。

場是指物體在空間的分布情況。根據場的定義,可由位移梯度進一步計算形變場的旋度(curl,ω)、膨脹量(dilatation,δ)和剪切應變(shear-strain,γ):

(16)

(17)

(18)

式中,uEW和uNS為同震位移場的東西(EW)和南北(NS)形變分量。其中旋度(ω)是指矢量作旋轉運動的方向和強度,其正負符號指示了斷層的運動性質,旋度為正表示逆時針旋轉,為左旋走滑;旋度為負表示順時針旋轉,為右旋走滑。膨脹量(δ)是指矢量場中某點周圍的矢量向該點聚集或發散,指示斷層的擠壓或拉張運動性質,膨脹量不為0表示該點處存在垂向運動,為正表示拉張性,為負表示擠壓性。剪切應變(γ)是指南北、東西向位移分量分別在正東、正北方向的梯度之和,指示剪切應力的集中程度,剪切應變為正(夾角減小)表示左旋走滑為主,為負(夾角增大)表示右旋走滑為主。

綜上,位移梯度的旋度、膨脹量、剪切應力算子可以直觀地刻畫同震地表破裂形跡、指示發震斷層的運動學性質,還可以基于其數值開展統計分析。

3 2021年瑪多MW7.4地震實例分析

2021-05-22青海瑪多發生MW7.4地震,是繼2008年汶川地震之后中國大陸發生的震級最高的一次地震,造成沿昆侖山口-江錯斷層158 km長的地表破裂[12],并在東、西尾端發生分叉破裂[13]。瑪多地震發生后,我們選用完全覆蓋震區的10 m分辨率的Sentinel-2影像提取同震形變場。考慮到氣候、時間、拍攝角度等影響,選擇地震前、后同一時間段(震前:2020-10-12,震后:2021-10-17)、太陽方位角相近(159°～ 163°)、含云量低的各3景數據。考慮到短波紅外波段(Band8)形變值標準差最小,選擇Band8作為輸入進行亞像素相關性匹配處理。震區水系分布較廣,為最大程度抑制水體噪聲,對其進行掩膜和空值賦值處理。其中將COSI-Corr中的初始窗口和截止窗口分別設置32×32 px和16×16 px、迭代次數為2、掩膜閾值為0.9;MicMac中的滑動窗口設置為3×3px、相關性最小閾值為0.5、相關影響程度為2、正則化系數為0.3。

基于COSI-Corr和MicMac軟件得到的亞像素相關性匹配檢測結果整體趨勢接近,一致地刻畫了發震斷層的空間位置,與基于cm級無人機航片解譯的地表破裂帶形跡(圖2(a)中黑色實線)吻合度較高。比較分析可見,MicMac軟件降噪抑制效果更佳(圖2(d))。亞像素相關性匹配后得到的同震形變場含有大量的白噪聲和加性噪聲(圖3(a))。經過條帶誤差處理后,能有效去除條帶誤差和衛星姿態角誤差帶來的影響(圖3(b))。經過非局部均值濾波后的形變場更加平滑(圖3(c)),有利于快速判定跨斷層的位錯量和斷層的精準位置。經過軌道誤差處理后整體趨勢不變,但跨斷層的位錯量有約5%的減小(圖3(d))。精化后處理(圖3(e))的結果在有效降低干擾噪聲和抑制誤差的同時細節特征得到保留(圖2(b)、2(d)),圖2(c)跨斷層形變剖線拐點指示破裂在局部出現次級分支,與野外調查結果完全吻合[14]。

圖2 EW方向同震形變場Fig.2 EW coseismic displacement field

圖3 瑪多地震精細化后處理Fig.3 Refined post-processing of the Maduo earthquake

旋度場(圖4(a))、膨脹量場(圖4(b))、剪切應變場(圖4(c))直觀刻畫了地表破裂形跡。沿斷層走向,旋度普遍大于0(圖4(d)、(e)),表示逆時針旋轉運動,指示斷層發生了左旋走滑運動;膨脹量值較小、穩定性差,局部段落出現連續正值或負值(圖4(f)、(g)),表示斷層局部存在垂向運動(隆升或沉降)。剪切應變與旋度符號一致(圖4(h)、(i)),表明地震以左旋走滑為主,局部段落兼傾滑分量,與已有結果基本一致[13]。

圖4 旋度、膨脹量、剪切應變Fig.4 Curl, dilatation, and shear strain

4 結 語

本文系統介紹光學同震形變提取的關鍵技術和優化流程,改進濾波平滑方法,提出精細化后處理能有效減小噪聲、降低誤差。分析COSI-Corr和MicMac兩款主流軟件的關鍵差異,結果顯示,COSI-Corr采用小相關窗口噪聲會比較突出且位錯振幅減少,隨著滑動窗口尺寸減小,其噪聲抑制受到弱化,而MicMac因采用正則化算法,即使采用較小的滑動窗口也不會降低噪聲抑制效果。基于同震形變場采用Sobel算子推導具有地學意義的位移梯度、旋度、膨脹量、剪切應變等算子。最后以2021年青海瑪多MW7.4地震為例,驗證了基于Sentinel-2光學影像的同震形變場提取以及精細化后處理優化流程,并推導出瑪多地震的旋度場、膨脹量場和剪切應變場,由此得出瑪多地震的發震斷層為左旋走滑,斷層兩盤存在傾滑(逆沖)運動。研究表明,對亞像素相關性匹配結果的后處理可有效消除軌道誤差、條帶誤差、場地因素和傳感器因素導致的白噪聲和加性噪聲,精化形變場,更加客觀地描述同震形變的空間分布特征。