付玉蒂, 薛東劍, 張小軒

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

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

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

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

1 研究區(qū)概況

1.1 地理位置

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

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

1.2 地質(zhì)構(gòu)造及地震背景

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

2 數(shù)據(jù)概況及研究方法

2.1 數(shù)據(jù)概況

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

表1 SAR數(shù)據(jù)參數(shù)Table 1 SAR data parameters

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

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

2.2 研究方法

2.2.1 雙軌差分測量原理

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

針對地震這種具有大尺度、瞬時(shí)等特點(diǎn)的形變,采用雙軌差分干涉測量方法對其形變場進(jìn)行提取。雙軌法需要地震前后各一景SAR影像,基本思想是對影像做干涉處理,然后,以外部DEM為參考,對地形相位進(jìn)行反演,在干涉相位圖中將其去除,提取出震區(qū)的形變信息?；驹砣鐖D4所示。

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

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

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

(1)

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

圖5 D-InSAR技術(shù)流程圖Fig.5 D-InSAR technical flow chart

2.2.2 數(shù)據(jù)處理關(guān)鍵步驟

1)濾波

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

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

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

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

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

3 形變場分析

采用雙軌差分干涉方法對瑪多震前震后兩景SAR數(shù)據(jù)進(jìn)行處理后,通過干涉、濾波、相位解纏等各項(xiàng)優(yōu)化,最終得到瑪多地震形變場,如圖9所示,其形狀近似于一個(gè)呈現(xiàn)西北-東南向分布的不規(guī)則橢圓,這與瑪多縣的構(gòu)造線(西北至東南)走向一致。以AA′為參考界線,AA′下側(cè)為隆升區(qū),剖面上側(cè)為沉降區(qū)。從中部往南北方向,形變增大,且北部形變小于南部形變,最大視線向(line of sight,LOS)向上升形變量為0.81 m,最大下沉量為0.65 m。AA′形變剖面(圖10)沿發(fā)震的斷層展布,其形變值在0左右,變化極小;BB′形變剖面(圖11)穿過隆升和沉降區(qū),最大沉降量為0.55 m左右,從B開始穿過發(fā)震斷層,形變量從負(fù)值快速增加到該剖面最大隆升值約0.45 m。這與華俊等[12]的InSAR形變場結(jié)果(最大形變約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]的研究一致,破裂帶北盤沿視線方向呈拉伸運(yùn)動(dòng)(變形約0.65 m),南盤沿視線向呈縮短運(yùn)動(dòng)(變形約0.81 m),反映了瑪多地震的主要變形場為東西向水平運(yùn)動(dòng),具有明顯的左旋走滑特征。斷層兩側(cè)的相對視線向位移可達(dá)1.50 m,同震形變場的長軸方向(斷層走向)整體上為北-西-西(NWW),表明地震引起的地表破裂錯(cuò)動(dòng)位移明顯。斷層滑動(dòng)特征與楊君妍等[14]和王守文等[15]研究結(jié)果一致。

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

4 瑪多地震斷層滑動(dòng)分布反演

4.1 反演方法

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

4.2 斷層滑動(dòng)分布分析

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

awid=1.860+0.34alen

(2)

adip=-0.390+0.43awid

(3)

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

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

表2 滑動(dòng)斷層參數(shù)Table 2 Slip fault parameters

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

(4)

(5)

4.3 反演結(jié)果分析

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

圖13 slipBERI模擬結(jié)果Fig.13 slipBERI simulation results

5 結(jié)論

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

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

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