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

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

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

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

1 光學(xué)同震形變處理技術(shù)

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

1.1 亞像素相關(guān)性匹配技術(shù)

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

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

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

(1)

式中,(x,y)為像元坐標(biāo)。設(shè)Fpost(u,v)為fpost(x,y)的傅里葉變換,Fpre(u,v)為fpre(x,y)的傅里葉變換,將兩幅影像加窗后經(jīng)過傅里葉變換轉(zhuǎn)化為頻率域:

Fpost(u,v)=

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

(2)

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

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

(3)

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

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

(4)

式中,上標(biāo)G為共軛轉(zhuǎn)置。進(jìn)行SVD分解,用最大奇異值及相應(yīng)的奇異向量對歸一化互功率譜矩陣H進(jìn)行秩1估計:

(5)

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

(6)

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

1.2 光學(xué)同震形變提取的優(yōu)化流程

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

1)失相關(guān)區(qū)域掩膜。對同震形變場中信噪比(signal-noise ratio,SNR)低于設(shè)定經(jīng)驗(yàn)閾值的像素點(diǎn)進(jìn)行掩膜,抑制失相關(guān)噪聲點(diǎn)對后續(xù)誤差處理的影響。

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

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

(7)

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

(8)

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

(9)

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

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

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

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

圖1 光學(xué)同震形變提取優(yōu)化流程Fig.1 Optimization process for optical coseismic deformation extraction

1.3 主流光學(xué)同震形變軟件對比分析

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

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

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

Cost=C1*(1-C3)

(10)

C1=1-Cmin

(11)

(12)

(13)

式中,Cor為歸一化相關(guān)性系數(shù),Cmin為相關(guān)性最小閾值,γ為相關(guān)影響程度。

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

表1 COSI-Corr與MicMac結(jié)果差異

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

在數(shù)學(xué)上,梯度作為矢量表示為函數(shù)在某點(diǎn)沿法向的變化。將光學(xué)同震形變場看作二維離散函數(shù),同理可進(jìn)行梯度計算?；谙袼仄频奶荻扔嬎憧傻玫轿灰圃诳臻g某一位置沿某一方向的變化量,即位移矢量場。梯度變化最大的方向和變化量指示同震破裂在地表的投影。

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

?u(x,y)/?x=

(14)

?u(x,y)/?y=

(15)

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

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

(16)

(17)

(18)

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

綜上,位移梯度的旋度、膨脹量、剪切應(yīng)力算子可以直觀地刻畫同震地表破裂形跡、指示發(fā)震斷層的運(yùn)動學(xué)性質(zhì),還可以基于其數(shù)值開展統(tǒng)計分析。

3 2021年瑪多MW7.4地震實(shí)例分析

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

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

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

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

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

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

4 結(jié) 語

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