聯合slipBERI與D-InSAR技術的青海省瑪多縣地震形變場提取與模擬

付玉蒂, 薛東劍, 張小軒

(成都理工大學地球科學學院, 成都 610059)

青海省瑪多縣于2021年5月22日2:04發生里氏7.4級地震,震源深度可達17 km,震中坐標為(34.59°N,98.34°E)。由于地震震源淺、震級大、地面震動強烈,使得震中地區房屋和牲畜棚遭到破壞,道路、橋梁和其他基礎設施出現不同程度的隆起或倒塌。對地震形變進行監測能有效地評估震區破壞分布以及滑動斷層的特征,用于地震形變監測的傳統技術包括GPS技術和傳統的水準監測技術等,但其對大面積形變場的監測始終存在空間分辨率不足等問題[1-2]。基于這一問題,應用差分干涉合成孔徑雷達(differential interferometric synthetic aperture radar,D-InSAR)進行大面積地表形變監測成為近年來的一大趨勢。它具有毫米級的視線精度,特別適用于監測快速而劇烈的地表形變,如地震、火山等。

Massonnet等[3]論證了采用D-InSAR提取地震同震位移場的有效性,并且成功提取了加州Landers地區形變場,讓InSAR技術開始受到廣泛關注;Budi等[4]采用D-InSAR技術提取龍博地區的地震形變場,分析得到其地震分布特征。余祥偉等[5]采用D-InSAR方法于2020年提取并模擬宜賓市地震形變場,驗證了觀測結果的可靠性;李成龍等[6]于2021年結合升降軌InSAR形變場,反演新疆伽師縣地震斷層參數以及滑動分布,并研究了其發震構造;韓炳權等[7]提取了2022年四川瀘定地震形變場,并且借鑒Okada彈性位錯模型思想,進行滑動斷層反演,以驗證觀測值的有效性;高二濤等[8]采用D-InSAR提取了九寨溝升降軌地震形變場,并對其位置、大小進行分析,為下一步研究斷層位置提供了條件。

綜上所述,D-InSAR在監測瞬時而劇烈的形變方面(如地震)具有精度高、結果準確等優勢,現將slipBERI方法與D-InSAR技術相結合,提取青海省瑪多縣地震同震形變場,采用slipBERI方法對地震滑動斷層進行反演,結合模型和形變參數,提高低相干區域的精度,為研究此次地震的斷層滑動、瑪多地區地震防災減災等方面提供一定依據。

1 研究區概況

1.1 地理位置

瑪多縣地處青海省果洛藏族自治州西北部,地理坐標為33°50′N～35°40′N,96°50′E～99°20′E,北鄰都蘭縣,東臨海南自治州興海縣、瑪沁縣,南臨達日縣、四川省石渠縣,西部為玉樹自治州曲麻萊縣,西南面則與巴音喀拉山、玉樹藏族自治州稱多縣接壤,如圖1所示。瑪多縣境內有扎陵湖和鄂陵湖兩大主要湖泊,中國第二大河黃河也流經其境內。

圖1 瑪多縣地理位置示意圖Fig.1 A map of Maduo County’s geographic location

1.2 地質構造及地震背景

瑪多縣地屬高平原區,大部分地區海拔在4 500 m以上,地勢較為平坦。大地構造帶屬巴顏喀拉褶皺帶,構造線自西北向東南延伸。近年來,巴顏喀拉地塊先后經歷了多次大地震,使其及周邊地塊成為地震高發區,如2001年青海昆侖山口Ms8.1地震(邊界帶中部)、2008年新疆和田Ms7.3地震(西段)、2010年玉樹Ms7.1地震(區塊南部邊界)、2013年蘆山Ms7.0地震(區塊東部邊界附近)等[9]。處于青藏高原東北部的瑪多地震,位于瑪多-甘德斷裂帶上,距巴顏喀拉塊體邊緣的東昆侖斷裂帶70 km[10],是發生于巴顏喀拉塊體的又一次大地震,與板塊周邊斷裂帶高度相關[11],瑪多地震及其周圍地質構造如圖2所示。

2 數據概況及研究方法

2.1 數據概況

(1)合成孔徑雷達(synthetic aperture radar,SAR)數據。本次研究獲取覆蓋瑪多縣2景地震前、后Sentinel-1A衛星降軌SAR數據,詳細參數如表1所示。

表1 SAR數據參數Table 1 SAR data parameters

(2)精密軌道數據及數字高程模型(digital elevation model,DEM)數據。精密定軌星歷數據(precise orbit ephemerides,POD)是sentinel-1衛星最精確的軌道數據,在全球導航衛星系統(global navigation satellite system,GNSS)下行21 d后才能使用。每天生成一個文件,可覆蓋26 h,定位精度優于5 cm。POD精密定軌星歷數據用于參考校正,能有效地減少形變相位誤差,提高監測精度。同時,處理過程中采用DEM數據進行地形校正和地理編碼,基于研究區域的地理位置,選取30 m分辨率的航天飛機雷達地形測繪使命-數字高程模型(shuttle radar topography mission-digital elevation model,SRTM-DEM)數據作為參考DEM數據,如圖3所示。

圖3 瑪多縣DEM圖Fig.3 DEM of Maduo County

2.2 研究方法

2.2.1 雙軌差分測量原理

1.2.2 計量分析法 主要使用CiteSpaceV分析系統繪制知識圖譜，它是基于JAVA平臺的可視化應用軟件，通過該軟件對文獻數據分析處理繪制出直觀易懂的科學圖譜，分析探測相關研究的重點、熱點、前言、變化趨勢等[2]，也是一款著眼于分析科學分析中蘊含的潛在知識，并在科學計量學、數據和信息可視化背景下逐漸發展起來的一款引文可視化分析軟件[3]。借助知識圖譜可視化技術，探究我國體育教學評價領域的研究現狀、熱點及發展趨勢，把握該領域發展動態及演化規律，從而為我國體育教學評價的進一步研究提供理論支撐。

針對地震這種具有大尺度、瞬時等特點的形變,采用雙軌差分干涉測量方法對其形變場進行提取。雙軌法需要地震前后各一景SAR影像,基本思想是對影像做干涉處理,然后,以外部DEM為參考,對地形相位進行反演,在干涉相位圖中將其去除,提取出震區的形變信息。基本原理如圖4所示。

A1、A2為傳感器在地震前后兩次觀測P點的位置;R1、R2為傳感器到P點之間的距離;A1、A2之間的線性距離是基線L;α為基線L與水平面的夾角;γ為入射角;β⊥為垂直基線;β∥為平行基線;h為觀測點P的高程;H為衛星飛行高度;δ為R1、R2之間的距離差圖4 D-InSAR基本原理圖Fig.4 D-InSAR basic schematic diagram

獲取形變信息的關鍵是從干涉圖中去除地形信息,以數字高程模型數據為參考對地形相位進行反演,并在干涉圖中進行二次差分,以消除地形相位,獲得震區的形變相位。采用雙軌差分提取地震形變場,具有數據易獲取、運算速率快、精度高等優點。記Δφ為兩幅SAR圖像的相位差,在SAR影像的獲取過程中,大氣等各種因素也會影響電磁波的傳輸,即Δφ表示為

Δφ=φflat+φtopography+φmove+φatomosphere+φnoise

(1)

式(1)中:φflat為干涉相位受平地效應的影響;φtopography為受地形相位影響的部分;φmove為遙感平臺移動時引起的干涉相位誤差;φatomosphere為大氣相位延遲;φnoise為噪聲引起的相位誤差。其具體流程如圖5所示。

圖5 D-InSAR技術流程圖Fig.5 D-InSAR technical flow chart

2.2.2 數據處理關鍵步驟

1)濾波

去平后的干涉圖(圖6)中還有大量的噪聲斑點,這難免會對形變信息的提取產生一定影響,導致干涉圖的相位值產生偏差或不連續,這不僅會影響后續解纏質量,還會影響最終結果,影響監測精度。因此,應通過濾波去除噪聲斑點。3種常用濾波方法包括Adaptive、Goldstein以及Boxcar,其中Adaptive濾波用于處理高分辨率影像;Boxcar濾波采用局部濾波原理,可以有效地保留微小條紋;Goldstein濾波方法使用率最高,它之所以能有效抑制噪聲斑點、改善條紋可見度,是因為其具有可變的濾波器。在處理過程中,考慮到圖中條紋的清晰度和細條紋的存在,選擇了Goldstein方法,將窗口值擴大至4,改善濾波效果,濾波后的差分干涉圖(圖7)噪聲顯著減少,干涉條紋清晰,有利于后續相位解纏處理。

圖6 瑪多地震區差分干涉圖Fig.6 Differential interferogram of Maduo earthquake area

圖7 Goldstein濾波后的差分干涉圖Fig.7 The difference interferogram after Goldstein filtering

相位解纏在整個D-InSAR處理過程中十分關鍵,相位解纏的結果與形變信息提取的效果密切相關。在未解纏時,干涉圖的相位介于[-π,π],此時相位與實際情況相差2nπ,對它解纏是將此時的相位值轉換為真實相位值的過程,解纏效果也會直接影響形變信息的質量。根據研究區域的特點,本文研究中采用最小費用流(minimum cost flow,MCF)法進行相位解纏,得到解纏相位(圖8)。

圖8 相位解纏圖Fig.8 Phase unwrapping

3 形變場分析

采用雙軌差分干涉方法對瑪多震前震后兩景SAR數據進行處理后,通過干涉、濾波、相位解纏等各項優化,最終得到瑪多地震形變場,如圖9所示,其形狀近似于一個呈現西北-東南向分布的不規則橢圓,這與瑪多縣的構造線(西北至東南)走向一致。以AA′為參考界線,AA′下側為隆升區,剖面上側為沉降區。從中部往南北方向,形變增大,且北部形變小于南部形變,最大視線向(line of sight,LOS)向上升形變量為0.81 m,最大下沉量為0.65 m。AA′形變剖面(圖10)沿發震的斷層展布,其形變值在0左右,變化極小;BB′形變剖面(圖11)穿過隆升和沉降區,最大沉降量為0.55 m左右,從B開始穿過發震斷層,形變量從負值快速增加到該剖面最大隆升值約0.45 m。這與華俊等[12]的InSAR形變場結果(最大形變約0.9 m)較一致。

圖9 青海省瑪多縣地震形變圖Fig.9 Earthquake deformation map of Maduo County, Qinghai Province

圖10 AA′形變剖面圖Fig.10 AA′ deformation profile

圖11 BB′形變剖面圖Fig.11 BB′ deformation profile

瑪多地震有明顯的地表破裂痕跡,這與李智敏等[13]的研究一致,破裂帶北盤沿視線方向呈拉伸運動(變形約0.65 m),南盤沿視線向呈縮短運動(變形約0.81 m),反映了瑪多地震的主要變形場為東西向水平運動,具有明顯的左旋走滑特征。斷層兩側的相對視線向位移可達1.50 m,同震形變場的長軸方向(斷層走向)整體上為北-西-西(NWW),表明地震引起的地表破裂錯動位移明顯。斷層滑動特征與楊君妍等[14]和王守文等[15]研究結果一致。

根據此次地震的形變信息、形變特征以及斷層走向,可以判斷該破裂帶位于巴顏喀拉塊體,為昆侖山口-江口斷裂。此次地震導致巴顏喀拉地塊向東強烈擠壓,從而導致北部走滑斷層在構造作用下左旋運動,與之前發生在巴顏喀拉地塊外圍邊界的中強地震不同,瑪多7.4級地震是近年來該區塊內發生的唯一一次強震。

4 瑪多地震斷層滑動分布反演

4.1 反演方法

將slipBERI[16-17](slip from bayesian regularized inversion)滑動斷層反演方法用于瑪多地震滑動斷層的反演,該方法與常用斷層反演方法如梯度下降法、Okada彈性位錯模型不同,其使用Okada模型計算表面變形的解析解,再引入馮·卡曼自相關檢驗與地震的分型特性相結合,利用貝葉斯方法對輸入的地震斷層數據(InSAR)的滑動進行求解,在滑動模型的建立中,能對光滑度以及反演模型進行選擇,并且其自帶模擬退火算法用于選擇貝葉斯方法起始值,整體操作簡便,適合用作斷層滑動分布反演。

4.2 斷層滑動分布分析

反演原始數據使用D-InSAR處理得到的瑪多地震形變數據,在反演過程中采用馮·卡曼正則化結合貝葉斯方法,通過馬爾科夫鏈蒙特卡羅(Markov Chain Monte Carlo,MCMC)鏈進行數字近似求解,光滑方式選擇馮·卡曼正則化配套光滑方式進行反演。斷層相關參數如寬度及深度使用Mai等[18]的公式,計算公式為

awid=1.860+0.34alen

(2)

adip=-0.390+0.43awid

(3)

式中:alen為斷層長度;awid為斷層寬度;adip為斷層深度。

斷層滑動分布反演結果如圖12所示,從圖12中可見,左旋走滑是瑪多地震同震滑動主要走向,與D-InSAR得到形變場結果相同,沿斷層走向破裂為85 km,接近地表滑動量大多在0～1 m,最大滑動量約為5 m處于斷裂帶東部,位于地下10 km左右,東側滑移量明顯大于西側滑移量,斷層主體破裂帶位于地下5～15 km,斷層參數與全球質心矩張量目錄(global centroid-moment-tensor project,GCMT)給出的斷層參數基本一致,反演得出地震矩約為1.108 8×1020NM,根據式(4)和式(5)得出震級約為Mw7.36,與GCMT所給出震級Mw7.4基本相當,如表2所示。

表2 滑動斷層參數Table 2 Slip fault parameters

圖12 反演得到的瑪多地震滑動斷層分布圖Fig.12 Earthquake slide fault distribution map obtained by inversion

(4)

(5)

4.3 反演結果分析

InSAR數據滑動斷層分布反演擬合情況如圖13所示。圖13(a)為地震原始形變圖,圖像中紅線為此次斷層軌跡,沿著這一軌跡附近反演效果良好;圖13(b)為slipBERI方法結合InSAR數據模擬得到的形變圖,斷層上部為沉降區,下部為隆升區,斷層附近形變量在-0.2 m和0.2 m左右,與原始數據基本吻合;圖13(c)為前兩幅圖像擬合得到的殘差圖,大部分區域殘差小于10 cm,但圖像右下角區域擬合效果不佳,對照地震形變數據,該區域形變復雜,斷層可能還存在另一條分支斷層,這與余鵬飛等[19]的研究結果一致。

圖13 slipBERI模擬結果Fig.13 slipBERI simulation results

5 結論

綜合D-InSAR技術監測和slipBERI模擬得到的瑪多地震形變場結果,可以得出以下結論。

(1)瑪多地震形變場形狀近似于一個不規則橢圓,受西北-東南走向的左旋走滑斷層控制,斷層下側為隆升區,最大視線向隆升量達0.81 m;斷層上側為沉降區,其最大視線向沉降量為0.65 m。此次形變場斷層整體呈北-西-西(NWW)走向,根據形變信息、形變特征和斷層走向,可以判斷出該斷裂帶位于巴顏喀拉塊體內,屬昆侖山口-江口斷裂。

(2)slipBERI模擬得到的形變場與觀測結果相似,大部分地區的殘差小于10 cm,這不僅驗證了觀測結果的可靠性,而且反演結果連續完整,形變場邊界清晰,有助于恢復大氣造成的觀測結果缺失等問題。