傅志豪 王 鑫 張景發(fā)

1 應(yīng)急管理部國家自然災(zāi)害防治研究院，北京安寧莊路1號，100085

2019-07-04美國加利福尼亞州發(fā)生MW6.4地震，34 h后該區(qū)域又發(fā)生MW7.1地震，震中位于小湖斷裂帶及Airport Lake斷裂帶交匯區(qū)域[1]，該地區(qū)構(gòu)造活動強(qiáng)烈，曾發(fā)生多次7級以上地震，包括1952年克恩縣MW7.5地震、1992年蘭德斯MW7.3地震及1999年MW7.1赫克托礦地震等[2]。地震發(fā)生后，許多學(xué)者利用InSAR、GPS數(shù)據(jù)提取了本次地震的形變場，但由于加州地震的震中區(qū)域地表破壞嚴(yán)重，D-InSAR結(jié)果在震中區(qū)域形變梯度過大，相位解纏失敗，無法提取準(zhǔn)確的斷層破裂軌跡[3]；此外，基于升降軌InSAR數(shù)據(jù)的offset-tracking方法雖可避免失相干影響提取方位向、距離向形變場，但在三維解算同震三維形變場時仍受D-InSAR結(jié)果控制；GPS數(shù)據(jù)由于點位稀疏，在地震等大區(qū)域形變監(jiān)測領(lǐng)域受到很大限制。相對而言，光學(xué)遙感數(shù)據(jù)具有覆蓋范圍全面、分辨率高、獲取方便等優(yōu)勢，其局限性在于只適合監(jiān)測水平向形變，對地表垂向沉降不敏感，同時形變測量精度受限于影像匹配精度及像元分辨率，但在監(jiān)測以水平位移為主、地表破裂明顯的形變領(lǐng)域，高分辨率光學(xué)影像同樣可以提供高精度地表形變場。目前，該方法廣泛應(yīng)用于冰川、滑坡及同震位移形變場監(jiān)測領(lǐng)域[4-7]。

本文利用地震前后的Landsat-8和Sentinel-2光學(xué)影像數(shù)據(jù)提取加州地震同震形變場，通過光學(xué)形變后處理算法去除系統(tǒng)誤差項，并利用形變誤差的方差-協(xié)方差函數(shù)及融合offset-tracking技術(shù)的InSAR同震三維形變場綜合驗證光學(xué)形變測量的精度，最終根據(jù)形變信息提取地震的發(fā)震斷層幾何展布及水平位移信息，包括斷層長度、最大水平位移量等參數(shù)，為該地震的震源參數(shù)反演及發(fā)震背景研究提供依據(jù)。

1 數(shù)據(jù)處理與方法

1.1 光學(xué)影像處理

表1 光學(xué)數(shù)據(jù)參數(shù)

圖1 加利福尼亞地震遙感影像概況Fig.1 Remote sensing image of the 2019 MW7.1 and MW6.4 earthquakes in California

分別選用地震前后的兩景Landsat-8和Sentinel-2光學(xué)影像第8波段數(shù)據(jù)，數(shù)據(jù)來源于USGS網(wǎng)站，參數(shù)如表1所示。原始影像如圖1所示，選取的兩類影像含云量均小于3%，時間基線分別為16 d及10 d，可有效避免地物變化造成的時間失相干及云霧所引起的噪聲。Sentinel-2光學(xué)數(shù)據(jù)經(jīng)大氣校正及輻射定標(biāo)處理，Landsat-8數(shù)據(jù)屬于1級產(chǎn)品，由于兩者數(shù)據(jù)覆蓋范圍不一致，為減少運(yùn)算量，將Landsat-8影像裁剪至相同范圍。最后，以COSI-Corr為數(shù)據(jù)處理平臺，基于頻率域互相關(guān)算法，通過逆傅里葉變換計算兩景影像的相對位移，經(jīng)亞像元匹配可達(dá)到1/20亞像元精度，即同震位移形變場精度可提高20倍[8](Landsat-8為0.75 m，Sentinel-2為0.5 m)。假設(shè)i1、i2分別為震前、震后兩景影像，其相對位移量為Δx、Δy，對其進(jìn)行傅里葉變換后存在以下關(guān)系式：

i2(x,y)=i1(x-Δx,y-Δy)

(1)

I2(wx,wy)=I1(wx,wy)e-j(wxΔx+wyΔy)

(2)

式中，I1、I2為兩幅影像的傅里葉變換，wx、wy分別為行和列的頻率變化。兩幅影像的歸一化互譜Ci1,i2為：

ej(wxΔx+wyΔy)

(3)

式中，符號*表示復(fù)數(shù)共軛相乘?；谠摵瘮?shù)進(jìn)行逆傅里葉變化可得到與位移量有關(guān)的互相關(guān)函數(shù)F-1，之后取二維脈沖函數(shù)的峰值點即可得到兩景影像相對位移及信噪比信息：

F-1{ej(ωxΔx+ωyΔy)}=δ(x+Δx,y+Δy)

(4)

圖2為Sentinel-2、Landsat-8融合EW向、NS向形變場及信噪比的三波段融合影像，其中area1、area2、area3為失相干區(qū)、軌道區(qū)及條帶誤差區(qū)，area4為本文選定的精度評定區(qū)域。

1.2 誤差后處理

1)光學(xué)影像特征追蹤法通過計算同名點像元值的偏移量提取形變信息，在植被、建筑、雨雪覆蓋區(qū)域影像DN值隨時間變化劇烈，容易丟失相關(guān)性。如圖2(a)和2(b)所示，area1位于Ridgecrest小鎮(zhèn)湖泊及農(nóng)田附近，受時間失相干影響，在遠(yuǎn)離震中區(qū)域仍存在大量異常形變值，在形變圖上表現(xiàn)為信噪比均低于0.9。本文通過掩膜處理去除整景形變場中存在的低相干區(qū)域，并對該區(qū)域進(jìn)行插值平滑處理。

圖2 2019年加利福尼亞MW7.1及MW6.4地震同震形變場Fig.2 Co-seismic deformation field of the 2019 MW7.1 and MW6.4 earthquakes in California

2)衛(wèi)星成像過程中由于拍攝位置的變化，形變場全局存在一個固有偏移量，即軌道誤差。如圖2(a)所示，受軌道誤差影響，形變場存在一個明顯的線性趨勢，常采用多項式曲面擬合模型方法去除軌道誤差[9]。在非形變區(qū)域建立參考點，根據(jù)最小二乘原理解算待求參數(shù)，用原始形變場減去模擬的整幅圖像軌道誤差項，即可得到地表真實形變場：

φorbit=a0+a1x+a2y+a3xy

(5)

式中，φorbit為軌道趨勢項誤差，x為水平向坐標(biāo)，y為垂直向坐標(biāo)，a0、a1、a2、a3為待求參數(shù)。

3)推掃式衛(wèi)星成像過程中傳感器保持固定，CCD(電荷耦合器件)線性陣列探測器所導(dǎo)致的非線性失真小于1 mm，在形變監(jiān)測的精度范圍之內(nèi)，因此誤差主要來源于飛行平臺抖動導(dǎo)致的影像幾何變形。如圖2(b)所示，由于Landsat-8、Sentinel-2均為處理后的一級產(chǎn)品，無法獲取傳感器軌道參數(shù)及姿態(tài)參數(shù)進(jìn)行校正，沿形變場交叉軌道方向存在大量平行排列的帶狀條紋。常規(guī)的均值相減法容易引入地形、植被造成的異常形變值，因此本文采用改進(jìn)均值相減法[10]，將旋轉(zhuǎn)后的形變圖均分成12等份，使條帶誤差呈豎直分布，再以每一等份的均值代表該行像素所受到的姿態(tài)誤差影響水平，最終減去該誤差即可去除形變場中帶狀條紋。

處理后的影像如圖3所示，整體形變場信噪比均高于0.96，Sentinle-2的EW向形變場中異常低值及Landsat-8影像中殘存的條帶誤差均被有效去除。由于地震發(fā)生于巴拿明山谷及死亡谷附近，受地形誤差影響，在Sentinel-2形變場左上區(qū)域仍存在部分異常低值；Landsat-8影像由于只經(jīng)過粗略的輻射校正，且數(shù)據(jù)質(zhì)量、分辨率較低，同時衛(wèi)星姿態(tài)誤差并未完全去除，遠(yuǎn)場區(qū)域存在大量異常形變值。

圖3 2019年加利福尼亞MW7.1及MW6.4地震同震形變場Fig.3 Co-seismic deformation field of the 2019 MW7.1 and MW6.4 earthquakes in California

2 同震形變場分析與評定

2.1 地表形變場分析

通過數(shù)據(jù)處理流程，可獲取整個研究區(qū)EW向及NS向的同震形變場。如圖3所示，加州地震由AB、CD兩段地表形變場中的斷層破裂跡線組成，其中Landsat-8因受分辨率及數(shù)據(jù)質(zhì)量制約，遠(yuǎn)場區(qū)域存在異常形變值，但所揭露的斷層軌跡與Sentinel-2一致。在地表形變場中，AB段發(fā)震斷層跡線整體呈NNW向展布，斷層南端距加洛克斷層僅3 km，在南部三分之一處走向轉(zhuǎn)變?yōu)镹W，由南向北逐漸延伸。該段發(fā)震斷層最大水平滑移量為2.82 m，長度達(dá)55 km，且上盤形變明顯大于下盤，在遠(yuǎn)離斷層方向形變值逐漸減小，斷層呈右旋走滑運(yùn)動特征，推測MW7.1主震是由NNW向的斷層運(yùn)動所引起的。距NNW破裂跡線南部18 km處存在一個與之相交的NE向破裂跡線，該斷層跡線整體走向NE，最大滑移量為1.05 m，地表破裂軌跡15 km，運(yùn)動學(xué)特征以左旋走滑為主，推測MW6.4前震是由NE向斷層運(yùn)動所致。由于本次地震震源機(jī)制復(fù)雜，地震序列中包含多達(dá)20個2 km以上的斷層破裂，受到MW7.1主震及余震影響，NS向形變場中CD段斷層上、下兩盤均表現(xiàn)為正N向運(yùn)動，但斷層上盤形變值遠(yuǎn)低于下盤；EW向形變場中斷層上下兩盤仍呈相對EW向運(yùn)動趨勢，說明實際斷層呈現(xiàn)走滑特征，滑動方向為左旋。

Fielding等[3]利用升降軌ALOS-2和Sentinel-1影像，通過聯(lián)合offset-tracking及升降軌D-InSAR結(jié)果解算加州地震三維形變場，最后基于GNSS結(jié)果驗證同震形變場精度，提取出的形變場EW向形變范圍為-1.5～1.5 m，NS向形變范圍為-2～2 m。該結(jié)果與本文利用光學(xué)影像所提取的同震形變場具有很好的一致性，同時由于D-InSAR提取的LOS向形變在震中區(qū)域出現(xiàn)相位纏繞，沿發(fā)震斷層附近區(qū)域存在大量空值區(qū)，本文所提取的斷層軌跡可為后續(xù)震源參數(shù)分析及滑動模型反演提供約束條件。

2.2 精度評定

為估算形變場中存在的誤差，本文選取遠(yuǎn)離形變區(qū)分布的area4開展精度評定，該區(qū)域遠(yuǎn)離震中分布，去除誤差項后實際形變值應(yīng)為0。根據(jù)地統(tǒng)計學(xué)中的變異函數(shù)理論，分別采用方差、協(xié)方差函數(shù)對該區(qū)域形變場取樣，以推算整體形變場的誤差水平[11]。其中，變差函數(shù)及協(xié)方差函數(shù)可表示為：

(6)

(7)

式中，xi為影像中隨機(jī)采樣點的位置，f(xi)為在影像中隨機(jī)采樣點xi上的形變值，hi為距離間隔，N為給定距離h時采樣點對的數(shù)目。

從圖4可以看到，EW向形變協(xié)方差在15 km處逐漸趨于穩(wěn)定，NS向協(xié)方差在20 km才開始穩(wěn)定，表明形變場在大于20 km處的兩點空間上不存在相關(guān)性。由于Sentinel-2光學(xué)影像未借助衛(wèi)星輔助姿態(tài)文件進(jìn)行校正，同震形變場始終存在因地形因素造成的誤差形變值，在協(xié)方差曲線上反映為影像存在微弱的相關(guān)性。根據(jù)變差函數(shù)顯示，EW向及NS向的方差分別為35.17 cm和74.96 cm，說明NS向形變場受地形影響更為明顯。

圖4 同震形變場誤差評定Fig.4 Error assessment diagram of co-seismic deformation

3 結(jié) 語

本文以2019年加州MW7.1和MW6.4地震為例，詳細(xì)論述Landsat-8、Sentinel-2光學(xué)影像同震形變數(shù)據(jù)處理流程及誤差后處理方法，開展了地表同震形變場分析，并結(jié)合InSAR同震三維形變及形變誤差的方差-協(xié)方差函數(shù)綜合驗證光學(xué)形變精度。結(jié)果顯示，兩次地震的發(fā)震斷層分別為NNW向的左旋走滑斷裂及NE向的右旋走滑斷裂，形成了共軛關(guān)系，其中MW6.4前震以左旋走滑為主，最大滑移量約為1.05 m，斷層長度為15 km；MW7.1主震以右旋走滑為主，地表最大滑移量達(dá)到2.82 m，斷層長度約55 km。由于D-InSAR結(jié)果在缺少斷層近場數(shù)據(jù)的情況下常導(dǎo)致錯估滑動分布[12]，利用GPS融合InSAR的方法常導(dǎo)致斷層反演的最大滑動量小于實際值[13]?；诟叻直媛使鈱W(xué)影像獲取的同震形變場可有效避免D-InSAR形變梯度過大時出現(xiàn)的失相干現(xiàn)象，準(zhǔn)確揭示該地震的地表形變、發(fā)震斷層幾何展布和水平位移信息，從而提高加利福尼亞地震序列的斷層滑動分布和震源參數(shù)反演的準(zhǔn)確性，并為本次地震的發(fā)震斷層及運(yùn)動學(xué)背景分析提供依據(jù)，同時為國內(nèi)開展基于光學(xué)遙感影像的地表同震形變場研究奠定基礎(chǔ)。

致謝：本文所采用的Landsat-8、Sentinel-2光學(xué)影像來源于美國USGS，在此表示感謝。