基于光學(xué)衛(wèi)星影像的同震形變提取技術(shù)及形變梯度的地學(xué)意義

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

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

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

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

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

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

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

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

假設(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是一個(gè)秩為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)的奇異向量對(duì)歸一化互功率譜矩陣H進(jìn)行秩1估計(jì):

(5)

式中,σmax為最大奇異值,umax和vmax為對(duì)應(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)化流程

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

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

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

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

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

5)賦空值-中值濾波。當(dāng)形變場(chǎng)的空值較多時(shí),會(huì)在一定程度上影響形變表達(dá)和解讀,傳統(tǒng)方法中對(duì)此不加考慮。實(shí)驗(yàn)通過中值濾波法將失相關(guān)掩膜空值區(qū)域賦值,只對(duì)失相關(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é)同震形變軟件對(duì)比分析

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

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

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

Cost=C1*(1-C3)

(10)

C1=1-Cmin

(11)

(12)

(13)

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

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

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

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

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

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

?u(x,y)/?x=

(14)

?u(x,y)/?y=

(15)

式中,u為同震位移場(chǎng),px為同震位移場(chǎng)分辨率,下文同理。

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

(16)

(17)

(18)

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

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

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

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

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

圖2 EW方向同震形變場(chǎng)Fig.2 EW coseismic displacement field

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

旋度場(chǎng)(圖4(a))、膨脹量場(chǎng)(圖4(b))、剪切應(yīng)變場(chǎng)(圖4(c))直觀刻畫了地表破裂形跡。沿?cái)鄬幼呦?旋度普遍大于0(圖4(d)、(e)),表示逆時(shí)針旋轉(zhuǎn)運(yùn)動(dòng),指示斷層發(fā)生了左旋走滑運(yùn)動(dòng);膨脹量值較小、穩(wěn)定性差,局部段落出現(xiàn)連續(xù)正值或負(fù)值(圖4(f)、(g)),表示斷層局部存在垂向運(yùn)動(dòng)(隆升或沉降)。剪切應(yīng)變與旋度符號(hào)一致(圖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)窗口噪聲會(huì)比較突出且位錯(cuò)振幅減少,隨著滑動(dòng)窗口尺寸減小,其噪聲抑制受到弱化,而MicMac因采用正則化算法,即使采用較小的滑動(dòng)窗口也不會(huì)降低噪聲抑制效果。基于同震形變場(chǎng)采用Sobel算子推導(dǎo)具有地學(xué)意義的位移梯度、旋度、膨脹量、剪切應(yīng)變等算子。最后以2021年青?，敹郙W7.4地震為例,驗(yàn)證了基于Sentinel-2光學(xué)影像的同震形變場(chǎng)提取以及精細(xì)化后處理優(yōu)化流程,并推導(dǎo)出瑪多地震的旋度場(chǎng)、膨脹量場(chǎng)和剪切應(yīng)變場(chǎng),由此得出瑪多地震的發(fā)震斷層為左旋走滑,斷層兩盤存在傾滑(逆沖)運(yùn)動(dòng)。研究表明,對(duì)亞像素相關(guān)性匹配結(jié)果的后處理可有效消除軌道誤差、條帶誤差、場(chǎng)地因素和傳感器因素導(dǎo)致的白噪聲和加性噪聲,精化形變場(chǎng),更加客觀地描述同震形變的空間分布特征。