利用光學(xué)影像相關(guān)獲取2019年Ridgecrest地震序列同震形變

王樂洋,鄒阿健,許光煜

(東華理工大學(xué) 測繪工程學(xué)院，江西 南昌 330013)

地震從古至今對人類的危害都十分巨大，而地震同震形變是地震最直觀的表現(xiàn)，能用其評估地震災(zāi)害等級，研究斷層構(gòu)造和震源機(jī)制。獲取同震形變的方法有很多，目前主要還是使用GPS技術(shù)和InSAR技術(shù)[1-4]。GPS技術(shù)獲取同震形變在遠(yuǎn)離板塊邊界和活躍斷裂帶地區(qū)GPS連續(xù)觀測站分布較為稀疏，空間分辨率較低，不能很好地體現(xiàn)地震產(chǎn)生的地表形變。在獲取有地表破裂的地震同震形變時，InSAR技術(shù)容易受相位失相干等因素的影響，從而造成靠近斷層區(qū)域形變數(shù)據(jù)的缺失。而光學(xué)影像在氣候良好的情況下，對于監(jiān)測較大地表形變能夠很好地彌補(bǔ)兩者的不足。

隨著科技的發(fā)展，高質(zhì)量的光學(xué)影像逐漸增多，它們可以用來監(jiān)測由于地質(zhì)過程、氣候變化或者人類活動造成的地表變化?？茖W(xué)家使用光學(xué)影像測量地表形變最早是在1991年Bindschadler等用Landsat TM影像觀察南極冰川流速[5]。2000年P(guān)uymbroeck等首次將光學(xué)影像相關(guān)應(yīng)用在地震形變監(jiān)測中，他們使用SPOT影像獲取了1994年Landers地震形變[6]。光學(xué)影像的分辨率和質(zhì)量，成像系統(tǒng)引起的幾何畸變，光學(xué)影像相關(guān)(Optical image Correlation，OIC)技術(shù)實(shí)行的困難等都是光學(xué)影像地表形變監(jiān)測所需要解決的問題。隨著科技的進(jìn)步，影像處理方法的改進(jìn)和光學(xué)影像的配準(zhǔn)和相關(guān)(Co-registration of Optically Sensed Images and Correlation,COSI-Corr)工具的開發(fā)[7]，這些問題已經(jīng)得到解決。文獻(xiàn)[8]介紹了從一對衛(wèi)星圖像中準(zhǔn)確提取地表位移的方法流程，COSI-Corr實(shí)現(xiàn)了這些流程。在此之后，使用光學(xué)影像監(jiān)測地表形變變得相對容易，且精度更高[9-10]。

2019年7月6日3時19分(UTC)，美國加利福尼亞Ridgecrest東北處發(fā)生Mw7.1地震(震中位于35.770°N，117.599°W)，震源深度8 km，該地震發(fā)生前兩天(7月4日)發(fā)生了一次Mw6.4地震，此次事件被稱為Ridgecrest地震序列。盡管加州地區(qū)經(jīng)常發(fā)生地震，但已經(jīng)有20年未發(fā)生過如此強(qiáng)烈的地震，根據(jù)USGS(United States Geological Survey)地震目錄記載，1999年1月1日至2019年7月31日。加利福尼亞東部發(fā)生了7次6級以上地震，其中2019年Ridgecrest地震序列是1999年Hector Mine Mw 7.1地震以來震級最大、最強(qiáng)烈的地震。本文首先分析哨兵2號影像相關(guān)結(jié)果中表現(xiàn)為明顯橫向條帶的姿態(tài)角誤差的處理，在已有方法上進(jìn)一步改進(jìn)；再使用哨兵2號影像相關(guān)獲取2019年Ridgecrest地震序列同震地表形變，分析地震的地表形變特征；最后對Landsat 8和Landsat 7影像分別進(jìn)行相關(guān)，并與哨兵2號結(jié)果進(jìn)行比較。

1 研究區(qū)域和數(shù)據(jù)處理

1.1 研究區(qū)域構(gòu)造背景

美國加州位于太平洋板塊和北美洲板塊之間的交匯處，活躍的板塊運(yùn)動導(dǎo)致了加州成為美國地震最活躍的地區(qū)。持續(xù)的大地測量觀測數(shù)據(jù)表明太平洋板塊相對于北美洲板塊以約48 mm/a的速度向西北移動[11]。2019年Ridgecrest地震序列的構(gòu)造背景見圖1。

圖1 2019年Ridgecrest地震序列構(gòu)造背景圖

圖中黑線為板塊運(yùn)動碰撞的邊界線，白色方框是哨兵2號影像覆蓋區(qū)域，黃色五角星表示的是2019年7月4日發(fā)生的Mw6.4地震，藍(lán)色五角星表示的是2019年7月5日發(fā)生的Mw7.1地震(震源機(jī)制球來自USGS)，白色五角星表示的是1999年10月16日發(fā)生的Hector Mine Mw7.1地震。紅色圓點(diǎn)表示地震前后三周發(fā)生的大于4級的地震。黃色正方形表示Ridgecrest所在位置。

1.2 實(shí)驗(yàn)數(shù)據(jù)處理

本文研究使用歐洲空間局(European Space Agency,ESA)的哨兵2號第八波段影像數(shù)據(jù)，以及美國國家航空航天局(National Aeronautics and Space Administration,NASA)的Landsat 7和Landsat 8衛(wèi)星第八波段影像數(shù)據(jù)，影像對的相關(guān)信息見表1。

表1 文中采用的光學(xué)影像信息表

使用基于ENVI平臺的COSI-Corr軟件包，通過光學(xué)影像相關(guān)技術(shù)獲取同震地表形變場。由于歐空局的哨兵2號數(shù)據(jù)是L2A級數(shù)據(jù)，USGS的Landsat 7和Landsat 8數(shù)據(jù)是L1C級數(shù)據(jù)，都是經(jīng)過粗略的幾何校正和正射校正的產(chǎn)品，預(yù)處理只需保證影像包含整個震區(qū)以及影像對大小一致，相關(guān)參數(shù)設(shè)置如下: 配準(zhǔn)運(yùn)算窗口為32像素×32像素，偏移量步長為6像素×6像素(最終的分辨率結(jié)果為60 m)；掩膜閾值設(shè)為0.9；迭代次數(shù)為2次。哨兵2號影像對(2019-06-13—2019-07-03)相關(guān)的結(jié)果見圖2。其中，圖2(a)中正值代表東向形變，負(fù)值代表西向形變；圖2(b)中正值代表北向形變，負(fù)值代表南向形變。

圖2 哨兵2號影像對(2019-6-13—2019-7-03)得到的未處理誤差的東西向和南北向形變場

如圖2所示，處理數(shù)據(jù)得到的水平形變場中存在各種誤差，其中包括失相關(guān)噪聲、軌道誤差、條帶誤差、衛(wèi)星姿態(tài)角誤差、地形陰影誤差。由于地形陰影誤差較為復(fù)雜，本文未進(jìn)行考慮。為減弱或消除以上誤差帶來的影響，一般情況下采用以下方法：對于失相關(guān)噪聲，在軟件中掩膜掉形變場中信噪比(Signal Noise Ratio，SNR)低于經(jīng)驗(yàn)閾值0.9的區(qū)域，將過大和過小的值替代為空值(Nan值)。由于下載的影像沒有進(jìn)行嚴(yán)格的校正，得到的形變場中依舊會存在軌道誤差，對于軌道誤差，用一次多項(xiàng)式曲面擬合來消除。在衛(wèi)星成像過程中，傳感器的電荷耦合器件(Charge Coupled Device，CCD)陣列的錯位引起衛(wèi)星條帶偽影即條帶誤差[12]，對于條帶誤差，旋轉(zhuǎn)影像使條帶呈豎直分布，對列采用均值相減法[10]。在結(jié)果中由于衛(wèi)星成像過程中姿態(tài)角的變化，導(dǎo)致形變結(jié)果中出現(xiàn)衛(wèi)星姿態(tài)角誤差，對于衛(wèi)星姿態(tài)角誤差，一般表現(xiàn)為微弱的橫向帶狀信號，與條帶誤差相互垂直，旋轉(zhuǎn)相同角度使得條帶呈水平分布，再對行采用均值相減法。最后再用3像素×3像素的中值濾波窗口進(jìn)行進(jìn)一步降噪。

對于哨兵2號光學(xué)影像相關(guān)的結(jié)果，姿態(tài)角誤差表現(xiàn)為明顯的橫向條帶，使用傳統(tǒng)的均值相減法并不能去除，文獻(xiàn)[10]提出了改進(jìn)的均值相減法，將研究區(qū)域分成12等份再進(jìn)行均值相減，很好地去除了凱庫拉地震區(qū)域哨兵2號影像相關(guān)結(jié)果中的橫向條帶，但使用此方法對本次加州地震區(qū)域哨兵2號影像相關(guān)結(jié)果中橫向條帶的去除效果不是很理想(見圖3(a))，本文對其進(jìn)行進(jìn)一步的改進(jìn)。由于地形起伏，植被覆蓋，存在云和河流湖泊等因素，形變場結(jié)果中各區(qū)域相干性有差別，導(dǎo)致傳統(tǒng)的均值相減并不能完全去除橫向條帶，所以把結(jié)果分成若干區(qū)域，增強(qiáng)均值相減區(qū)域的相干性，再對各區(qū)域進(jìn)行均值相減來去除條帶。從圖3(a)中可以看出，雖然將橫向條帶分成12份使得各自區(qū)域內(nèi)相干性變強(qiáng)，但處理后還存在一些橫向的條帶，仔細(xì)觀察發(fā)現(xiàn)這些未去除的條帶存在于各豎向條帶的兩旁，可能是由于各豎向條帶兩旁區(qū)域相干性較差，但兩豎向條帶之間區(qū)域的橫向條帶相干性較好。基于這點(diǎn)，將橫向條帶按照豎向條帶所處位置進(jìn)行劃分，然后對劃分后的各區(qū)域進(jìn)行均值相減，并按照去除情況對豎向條帶內(nèi)的區(qū)域進(jìn)行一定數(shù)量的均分。

圖3 哨兵2號影像對(2019-6-13—2019-7-03)得到的處理完誤差的東西向和南北向形變場

由于結(jié)果中各處相干性有差別，各處的值是不均勻的，很難依靠旋轉(zhuǎn)后的列均值比較等方法來分辨圖中各豎向條帶所處的位置。本文根據(jù)旋轉(zhuǎn)前兩豎向條帶頂端經(jīng)差來確定兩豎向條帶之間區(qū)域所占比例，確定旋轉(zhuǎn)前形變場豎向條帶頂端各點(diǎn)在形變場所代表矩陣中的位置，然后得到旋轉(zhuǎn)后形變場中豎向條帶的大致位置。相關(guān)過程可表示為：

(1)

式中:L1，L2表示旋轉(zhuǎn)前第一個和第二個豎向條帶頂點(diǎn)的經(jīng)度；l表示旋轉(zhuǎn)前第一個和第二個豎向條帶頂點(diǎn)的經(jīng)度差；k表示旋轉(zhuǎn)前形變場所代表矩陣的列數(shù)；L表示整個結(jié)果圖中最右端與最左端的經(jīng)度差；k1表示旋轉(zhuǎn)前兩個豎向條帶頂點(diǎn)處的列數(shù)差；round(·)表示按照四舍五入原則將結(jié)果轉(zhuǎn)化成整數(shù)；按旋轉(zhuǎn)前形變場設(shè)立坐標(biāo)系，令左下角坐標(biāo)為(0,0)，(x0,y0)表示形變場的中心位置坐標(biāo)；(x1,y1)表示旋轉(zhuǎn)前第一個豎向條帶頂點(diǎn)的坐標(biāo)，(x2,y2)表示旋轉(zhuǎn)后該點(diǎn)的坐標(biāo)；b表示逆時針旋轉(zhuǎn)使得橫向條帶呈水平分布的角度。按求得k1相同的方法可得到形變圖中豎向條帶所劃分的各區(qū)域的列數(shù)差，從而得到旋轉(zhuǎn)前各豎向條帶頂點(diǎn)在所設(shè)坐標(biāo)系中的坐標(biāo)，根據(jù)求得(x2,y2)相同的方法可以得到旋轉(zhuǎn)后各豎向條帶頂點(diǎn)坐標(biāo)。

按照上述方法得到兩個豎向條帶旋轉(zhuǎn)后的橫坐標(biāo)之差可以得到兩豎向條帶的間距，根據(jù)此間距可以對旋轉(zhuǎn)后的形變場進(jìn)行劃分，再對每個區(qū)域進(jìn)行均值相減。本次事件覆蓋區(qū)域哨兵2號影像相關(guān)結(jié)果旋轉(zhuǎn)后最左端至第一個豎向條帶間距與兩豎向條帶間距相同，將其與其他4個兩豎向條帶之間區(qū)域形變場提取出來，按照下述過程去除橫向條帶。

(2)

式中:φ表示橫向條帶誤差水平；δ表示原形變場所代表矩陣；β表示減去橫向條帶誤差的形變場所代表矩陣；m,n表示提取出的形變場所代表矩陣的行數(shù)和列數(shù)；a,b表示分成的每個列等份的初始像素點(diǎn)和最后的像素點(diǎn)在矩陣中的位置；Dk表示該像素點(diǎn)的值；ceil(·)表示大于某一整數(shù)的最小值;c表示將形變場分成5個列等份，可以按照去除情況對各區(qū)域內(nèi)部進(jìn)一步均分，即可以將其分成5的倍數(shù)等份，但不宜過多。具體流程見圖4，結(jié)果見圖3(c)、圖3(d)。

圖4 本文方法去除橫向條帶數(shù)據(jù)處理流程

在圖3(c)、圖3(d)中，沒有看到明顯的橫向條帶，從視覺上看去除效果較為良好。選取豎向條帶旁區(qū)域作剖面線(見圖3(a)黑色線條)，其各點(diǎn)誤差值的情況見圖5。可以看出東西向形變場中直接均分的均值相減法得到的值明顯較為分散，而按本文方法則基本集中在-0.2～0.2 m之間，兩者精度分別為0.215 m和0.133 m，本文方法改善較大。而南北向形變場中可能由于其豎向條帶兩旁相干性差距不大，直接均分的均值相減法得到的結(jié)果中也沒有看到明顯的橫向條帶，兩種方法精度分別為0.148 m和0.138 m，從精度上看，本文方法略好于直接均分的均值相減法。

圖5 直接均分的均值相減法與本文方法得到的形變場剖面

2 同震地表形變場特征分析

在2019年Ridgecrest地震序列研究中，本文采用基于ENVI的COSI-Corr軟件包，通過光學(xué)影像相關(guān)技術(shù)，使用震前震后的哨兵2號影像對(2019-06-28—2019-07-18)，使用前節(jié)所述處理方法獲取該次地震的同震地表形變場如圖6(c)、圖6(d)所示，具體破裂區(qū)域如圖7所示。從圖中可以看出Ridgecrest Mw7.1地震的破裂跡線從東南到西北呈近直線形狀，Mw6.4地震的地表形變主要表現(xiàn)在Mw7.1地震破裂的東南方向，破裂跡從東北到西南呈近直線形狀，兩次破裂相互交叉。在東西形變場中，Mw7.1地震的破裂跡以北形變值為正，向東位移，破裂跡以南為負(fù)，向西位移，兩側(cè)形變較為平均，呈現(xiàn)拉伸趨勢，最大形變值在1.7 m左右；Mw6.4地震的破裂跡東西兩側(cè)均為負(fù)。南北形變場中，Mw7.1地震破裂線北側(cè)形變?yōu)樨?fù)，向南位移，而南側(cè)形變?yōu)檎虮蔽灰?，南?cè)形變主要集中在中部，北側(cè)形變比較平均，北側(cè)形變明顯高于南側(cè)，呈擠壓趨勢，最大形變值在2.3 m左右；Mw6.4地震的破裂跡東西兩側(cè)均為正。Mw6.4破裂跡兩側(cè)形變情況可能受在其之后發(fā)生的Mw7.1地震的影響，導(dǎo)致其破裂跡兩側(cè)位移方向相同，但還能看到明顯破裂跡線，說明雖然位移方向相同，但兩側(cè)大小還是有一定程度的落差。在結(jié)果圖中還可以看到在主破裂西方向的北部和南部均有一個小破裂。形變結(jié)果不僅能提供形變數(shù)據(jù)，還可以提供地震的斷層幾何參數(shù)。從本次地震的形變結(jié)果來看，可以得到該地震序列的主破裂的走向角大約為320°，破裂長度約為45 km。次破裂的走向角約為230°，破裂長度約為15 km。相關(guān)文獻(xiàn)中使用PlanetScope影像[13]，worldview影像[14]得到的Ridgecrest地震序列同震形變場結(jié)果，在形變大小的量級，以及斷層跡線細(xì)節(jié)上均與本文相似，幾者間存在著良好的一致性，可見本文結(jié)果較為合理。

圖6 哨兵2號影像獲取的Ridgecrest地震序列的同震形變場

圖7 哨兵2號數(shù)據(jù)獲取的Ridgecrest地震序列的同震形變場具體區(qū)域

本文通過選取遠(yuǎn)場部分區(qū)域(見圖6(c)黑色方框)，計(jì)算其標(biāo)準(zhǔn)差作為精度評定的指標(biāo)。由于結(jié)果中的大部分誤差都已經(jīng)去除，理論上選定區(qū)域的形變值為0。但誤差的處理會受到軟件亞像素相關(guān)過程的精度的限制，光學(xué)影像用COSI-Corr進(jìn)行相關(guān)性匹配監(jiān)測形變精度至少為1/20像素[8]，對于哨兵2號結(jié)果來說形變監(jiān)測精度為50 cm。按圖6(c)方框分別選取范圍內(nèi)的東西向和南北向形變，得到它們東西向和南北向的均方根誤差分別為0.287 m和0.3 m，精度水平比較合理[9-10]。

加州地區(qū)有數(shù)以千計(jì)的GNSS站，但是盡管有這么多精度較高的GNSS站，但其空間覆蓋率仍較低。目前Ridgecrest 地震序列周圍地區(qū)GNSS站點(diǎn)之間的平均間距為20～30 km[15]，與加州西部板塊邊界區(qū)域相比，本次地震周邊的GPS站點(diǎn)更少，站點(diǎn)之間的間距更遠(yuǎn)，與這次地震的斷裂長度相當(dāng)，不足以用來研究本次地震形變的細(xì)節(jié)。如果相鄰像素間的應(yīng)變梯度超過一個相位周期，InSAR測量通常會在破裂跡線附近產(chǎn)生失相關(guān)，相關(guān)文獻(xiàn)中哨兵1號與ALOS-2 PLAYSAR-2數(shù)據(jù)得到的Ridgecrest地震序列地表形變，在地表破裂跡線附近區(qū)域失相干情況較為嚴(yán)重，且不能清晰的看到前震破裂和主震破裂之外的破裂跡線[16]。相較于GPS和InSAR結(jié)果，哨兵2號影像相關(guān)的結(jié)果雖然損失了一定精度，但它既具有較高的空間分辨率，又能夠彌補(bǔ)InSAR結(jié)果中失相干的情況，更能體現(xiàn)該次地震的同震形變情況。

3 Ridgecrest地震序列的不同光學(xué)影像相關(guān)結(jié)果

為了探究使用其他開源光學(xué)影像獲取的本次地震的形變場情況，本文選取Landsat 7和Landsat 8兩種光學(xué)影像分別進(jìn)行相關(guān)，其結(jié)果見圖8和圖9。

圖9 Ridgecrest地震序列的Landsat 7影像相關(guān)結(jié)果

從圖6和圖8中可以看出，Landsat 8與哨兵2號相關(guān)結(jié)果中均能明顯看出地表破裂跡線，均能得出大致斷層走向和長度。但Landsat 8相關(guān)結(jié)果中破裂跡線沒有那么清晰，整個形變場中異常值更多，存在的噪聲誤差也更加明顯。選取圖中黑色方框區(qū)域，東西向和南北向均方根誤差為0.345 m和0.334 m，顯示效果和精度均稍差于哨兵2號結(jié)果。

圖8 Landsat 8獲取的Ridgecrest地震序列處理完誤差的同震形變場

從圖8和圖9可以看出，Landsat 7結(jié)果中兩圖均未發(fā)現(xiàn)Landsat 8結(jié)果中出現(xiàn)的明顯破裂跡線。在Landsat 8結(jié)果中，東西向形變場和南北向形變場中異常值和空值主要集中在瑟爾斯湖以及一些植被茂密、地形起伏的區(qū)域；在Landsat 7結(jié)果中，只有東西向形變場中充滿空值的瑟爾斯湖區(qū)域與Landsat 8結(jié)果同樣明顯，但其他區(qū)域充滿了異常值，在南北向形變場中則是充滿了空值。由此可見在本次地震中，Landsat 8結(jié)果與Landsat 7結(jié)果存在著較大的差距，且用Landsat 7影像相關(guān)無法獲取本次地震的同震形變。從結(jié)果可以看出，光學(xué)影像的分辨率對結(jié)果的顯示效果和精度有一定的影響，影像分辨率越高，相關(guān)結(jié)果的顯示效果和精度越好。而即使是同種分辨率的光學(xué)影像，結(jié)果也會有差別。Landsat 7和Landsat 8同樣為15 m空間分辨率，但兩者的傳感器不同，Landsat 7攜帶增強(qiáng)型專題制圖儀(Enhanced Thematic Mapper, ETM+)傳感器，而Landsat 8則攜帶陸地成像儀(Operational Land Imager, OLI)和熱紅外傳感器(Thermal Infrared Sensor, TIRS)，兩者發(fā)射時間相差14年，各項(xiàng)技術(shù)Landsat 8均較為成熟。與Landsat 7第八波段相比，Landsat 8的第八波段范圍較窄，可以更好地區(qū)分植被和無植被特征，數(shù)據(jù)較好。而且2003年之后Landsat 7衛(wèi)星受到損壞，導(dǎo)致之后的數(shù)據(jù)產(chǎn)生缺失，雖然已有方法對其數(shù)據(jù)進(jìn)行補(bǔ)充，但質(zhì)量依舊與損壞前有一定的差距。由于Landsat 7影像數(shù)據(jù)質(zhì)量相對較差，橫向的條帶誤差大小間隔過大，地形起伏、河流湖泊、植被覆蓋等因素造成的失相關(guān)情況更加嚴(yán)重，導(dǎo)致誤差不能很好地削弱。2021年3月將發(fā)射的Landsat 9衛(wèi)星將為光學(xué)影像地表監(jiān)測提供更加有效且豐富的數(shù)據(jù)。

4 結(jié) 論

本文使用哨兵2號光學(xué)影像數(shù)據(jù)獲取了2019年Ridgecrest地震序列的同震地表形變場，并對表現(xiàn)為明顯橫向條帶的姿態(tài)角誤差去除方法進(jìn)行改進(jìn)。同時對Landsat 7和Landsat 8兩種光學(xué)影像分別進(jìn)行相關(guān)，并與哨兵2號影像相關(guān)結(jié)果進(jìn)行比較。得到了以下的主要結(jié)論：

1) 在哨兵2號影像相關(guān)的結(jié)果中，豎向條帶會影響形變場相干性的連續(xù)，相較于傳統(tǒng)的均值相減法和直接均分的均值相減法，本文按照豎向條帶劃分的均值相減法對橫向條帶的去除效果更加良好。

2) Ridgecrest地震序列主要存在一個主破裂、一個次破裂和主破裂上的兩個小破裂，東西向呈拉伸趨勢，南北向呈擠壓趨勢，兩方向形變特征表明Ridgecrest地震序列主震是一個右旋走滑地震。

3) 在使用光學(xué)影像相關(guān)的方法進(jìn)行地表形變監(jiān)測時，光學(xué)影像的空間分辨率越高，顯示效果和精度越好；同分辨率的光學(xué)影像，由于影像質(zhì)量的不同，得到的結(jié)果依舊會有差別。