蘆山MS7.0地震同震形變的地形效應

王金馳 段虎榮,2 張成浩 梁文康 劉 鵬

1 西安科技大學測繪科學與技術學院,西安市雁塔中路58號,710054 2 中國科學院精密測量科學與技術創新研究院大地測量與地球動力學國家重點實驗室,武漢市徐東大街340號,430077

前人根據地震同震形變反演斷層參數時,為了簡化模型,常常會忽略真實地形因素。在計算同震形變時地表真實地形會產生多大影響至今仍然不明確,部分學者認為地形效應對同震形變影響較大[1-3],但也有學者持相反意見[4-6]。

2013-04-20四川蘆山發生MS7.0地震,震中位于四川盆地西緣。四川盆地平均海拔為250～750 m,而青藏高原平均海拔為3 000～5 000 m,地形梯度變化較大,因此此次地震是研究地形對地震同震形變影響的較好案例。譜元法(spectral element method,SEM)可兼顧有限元法網格的靈活性和偽譜法的高精度、快速收斂特性,在復雜地震動數值模擬方面得到廣泛應用[7-9]。本文利用譜元法模擬2013年蘆山MS7.0地震同震形變,探討地形對地表同震形變的影響。

1 譜元法計算同震形變

1.1 控制方程

圖1 準半空間中斷層示意圖

區域Ω的控制方程可表示為:

?·T+f=0

(1)

式中,T為應力張量,f為外力。

在自由表面Γ0上滿足邊界條件:

(2)

同震形變由彈性本構關系T=c:ε控制,c為4階彈性張量,ε為應變張量。

1.2 分裂節點法

地震震源由斷層面Σ上的滑動矢量Δs定義。如圖1所示,當從斷層Σ的正側向負側移動時,滑動矢量捕獲到位移場s的不連續性[10]:

(3)

本文采用分裂節點法對模型域進行網格劃分,使元素邊界滿足斷層面[11](圖2)。在分裂節點法中,通過在斷層面上規定一個不連續的位移場來引入載荷,因此控制方程(1)中f為0。

圖中使用譜元素(淺灰色單元)離散化研究區域Ω,在區域外添加一層無限元素(深灰色單元);黑色粗實線表示斷層

2 數值模擬

通過模擬基于彈性半空間各向同性介質的垂直斷層走滑和傾滑,驗證譜元法在地震同震形變模擬中的可靠性。采用圖3所示模型,尺寸為100 km×80 km×30 km,泊松比為0.25,楊氏彈性模量為5.68 GPa。圖中深灰色矩形表示斷層,其頂部距離地表4 km,底部距離地表14 km,中心距離模型左右兩側邊界50 km。使用近似為2 km×2 km×2 km的六面體單元對模型進行網格劃分,共包含30 000個單元和253 611個節點,節點自由度為725 370。模型頂面施加應力自由邊界條件,側邊和底面施加Dirichlet邊界條件。

表1 不同傾角下計算的同震位移

圖3 垂直斷層模型

2.1 走滑斷層引起的同震形變

設定斷層右旋純走滑量為4 m,并使用譜元法計算地表三維同震形變,將數值模擬結果與Okada解析法結果進行對比(圖4)。從圖4(a)和4(b)可以看出,采用2種方法計算得到的水平位移場沿斷層兩側呈反對稱分布,靠近斷層位置的水平位移較大且與斷層走向平行。從圖4(c)和4(d)可以看出,采用2種方法計算得到的垂直位移場呈四象限對稱分布。水平和垂直位移場的分布均符合右旋走滑斷層特征,并且譜元法與Okada解析法結果在數值上具有很好的一致性。

粗黑直線代表斷層所在位置,紅色箭頭代表斷層走向

為詳細研究譜元法與Okada解析法的差異,選擇位移場剖面X=-10 km位置,將2種方法的模擬結果進行對比。從圖5可以看出,對于水平位移x分量,除邊界區域外,兩者的振幅和相位基本一致。對于水平位移y分量,在斷層附近,左側Okada解析法結果大于譜元法,右側譜元法結果大于Okada解析法;而在超過7 km之后,Okada解析法結果大于譜元法。對于垂直位移z分量,在斷層附近,譜元法結果大于Okada解析法。總體而言,兩者差異最大不超過10%。

圖5 走滑斷層在X=-10 km處地表位移剖面

2.2 傾滑斷層引起的同震形變

設定斷層純傾滑量為4 m,同樣采用上節中所構建的模型區域和介質屬性。對比譜元法結果與Okada解析法結果(圖6)可知,水平和垂直位移場的分布均符合斷層傾滑特征,并且譜元法與Okada解析法結果在數值上具有很好的一致性。

粗黑實線代表斷層所在位置,紅色圓圈和叉代表斷層傾向

選取X=10 km地表位移剖面,對比譜元法數值模擬結果與Okada解析法結果。由圖7可知,x和z分量呈中心對稱分布,y分量呈軸對稱分布。靠近斷層時,2種方法所有分量的振幅和相位吻合較好;遠離斷層時,2種方法所有分量均存在一定差異,這可能是由建模誤差或數值方法的近似性所致。在下一步研究中,將考慮添加無限元素以更好地探究這些差異的原因。

圖7 傾滑斷層在X=10 km處地表位移剖面

綜合走滑運動和傾滑運動數值模擬結果可知,譜元法與Okada解析法計算結果一致性良好。因此,認為使用譜元法模擬地震地表同震形變具有可行性。

3 地形因素對2013年蘆山MS7.0地震同震形變的影響

3.1 數 據

本文使用的GPS同震位移資料來自于33個GPS連續站,GPS觀測數據最大水平誤差小于3 mm,最大垂直誤差小于12.4 mm[12]。圖8為本文研究區,震源中心位于30.3°N、103.0°E。使用單一平面斷層生成格林函數,斷層參數參照Duan等[13]的研究成果。

圖中沙灘球為CENC CMT震源機制解;紅色五角星為震中位置;U-V為譜元法計算時所采用的坐標系,U為沿斷層走向方向,V為垂直斷層走向方向

將蘆山地區真實地表-地殼結構納入模型中,地形網格文件使用NOAA發布的ETOPO1數據[14],分辨率為1′,包含陸地地形和海洋水深數據。模型介質參數參考談洪波等[15]的研究成果。

3.2 模型構建

為研究地形因素對同震形變的影響,構建圖9中2種模型,其中A為地形模型(圖9(a)),B為無地形模型(圖9(b))。本文模型基于均勻彈性介質,泊松比為0.25,楊氏彈性模量為79.8 GPa,水平尺寸為116 km×133 km,深度為40 km(不包括地形)。離散模型使用六面體單元,所有模型平均網格間距設置為2 km,模型A共包含79 750個單元和675 608個節點,模型B共包含81 200個單元和686 868個節點。

粉紅色矩形為斷層位置;U為沿斷層走向方向,V為沿垂直斷層走向方向

3.3 基于分裂節點法計算地形對蘆山地震同震形變的影響

2013年蘆山MS7.0地震是發生在龍門山斷裂帶西南段的一次淺源逆沖型地震,破裂長度20 km[16],地震引起明顯的地表水平位移及垂直形變,可為研究地形效應對同震形變的影響提供理想的條件。

本文分別基于模型A和模型B,利用譜元法計算地表同震位移,并探討地形對地表同震位移的影響。圖10為2種模型生成的地表同震位移場,可以看出,模型A計算的沿斷層走向方向位移為-71～0 mm,垂直斷層走向方向位移為-18～28 mm,垂直方向位移為-42～67 mm;模型B則分別為-80～2.5 mm、-19～31 mm和-47～75 mm。由圖9可知,同震形變場分布形態符合逆沖型地震特征,距離震中越近,同震形變越大。地形的存在不僅會改變地表同震形變場分布形態,也會改變同震位移的大小。從圖10(a)和10(d)可以看出,在震中上方區域,模型A形變場分布不如模型B平滑,且模型A的最大值比模型B少9 mm,對比圖10(b)與10(e)、10(c)與10(f),也有類似的情況。

第1、2行分別為基于模型A、B計算的同震位移,第3行為模型A與模型B之差;第1、2、3列分別代表沿斷層走向的水平位移、垂直斷層走向的水平位移、垂直位移

本文利用林曉光等[3]的方法計算地形效應的影響。結果表明,沿斷層走向方向、垂直斷層走向方向和垂直方向,地形對同震形變的影響最大可達約16.8%、19.2%和18.3%。

4 討 論

本文基于平面斷層模型研究地形對同震形變的影響,斷層參數在一定程度上會影響地表同震位移,這里主要討論斷層傾角和深度對地表同震位移的影響。

在Duan等[13]的研究結果基礎上,構建6種模型,采用控制變量法研究斷層傾角和深度對地表同震位移的影響,所有模型使用與§3.2相同的介質參數。數值模擬結果如表1、2所示。

根據所選研究區域,使用8個GPS站點觀測值,為了與譜元法結果進行對比,將這些GPS觀測值轉換到U-V坐標系下(譜元法計算時所采用的坐標系),沿斷層走向方向U分量變化范圍為-63.8～2.4 mm,沿垂直斷層走向方向V分量變化范圍為-17.6～22.1 mm,垂直方向Z分量變化范圍為-4.9～83.6 mm。

通過表1可以看出,隨著斷層傾角增大,沿斷層走向方向分量、垂直斷層走向方向分量和垂直方向分量位移均減小。

由表2可知,隨著斷層深度變淺,地表同震位移逐漸變大。這是因為淺層斷層距離地表更近,斷層運動能更直接地影響地表,從而使地表同震形變更為明顯。

表2 不同深度下計算的同震位移

在原始斷層幾何參數基礎上,將斷層走滑分量和傾滑分量分別設置為0.15 m和0.70 m,計算得到地形效應對沿斷層走向方向同震形變的影響最大可達約10.3%,對垂直斷層走向方向同震形變的影響最大可達約12.4%。這與Lin等[2]計算汶川地震地形效應對同震位移的影響(9.05%)不同,可能是由于其使用的是二維有限元模型,而本文采用的是三維譜元法模型。Langer等[1]計算智利地震地形效應對同震位移的影響為30%,遠大于本文計算值,這可能是由于秘魯-智利海溝地形急劇變化,相對高差過大,而且斷層兩側介質不同,一側為海域,而另一側為陸地。

值得一提的是,GPS結果是某一確定點的形變數據,而譜元法結果是整個網格區域的形變情況。為了方便二者對比,譜元法結果使用與GPS觀測點鄰近區域的位移值。本文選取同震區域位移較大的LS05、LS07、SCTQ測站進行對比。LS05測站東、北分量觀測值分別為-9.9 mm、-66.8 mm,譜元法模擬值分別為-10.4 mm、-53.3 mm;LS07測站東、北分量觀測值分別為24.1 mm、-17.8 mm,譜元法模擬值分別為17.1 mm、-17.0 mm;SCTQ測站東、北分量觀測值分別為-9.8 mm、-18.7 mm,譜元法模擬值分別為-7.0 mm、-20.1 mm。除LS05測站北分量差異較大外(13.5 mm),其他結果二者較為接近。

5 結 語

在不考慮地形條件下,采用譜元法與Okada解析法分別計算垂直斷層走滑、傾滑運動引起的地表同震位移。結果表明,無論對于走滑斷層還是傾滑斷層,2種方法結果的一致性均較好,表明譜元法模擬同震形變具有可靠性。

針對2013-04-20蘆山MS7.0地震,本文分別構建蘆山區域地形模型和無地形模型,探討地形效應對同震形變的影響。結果表明,地形因素不僅會改變地表位移場分布形態,而且量值在不同方向上也有所差異。對于水平同震形變分量而言,沿斷層走向和垂直斷層走向的影響最大分別約為10.3%、12.4%;對于垂直同震形變分量而言,最大影響可達約11.9%。反演此次地震斷層滑動速率時,若不考慮地形因素的影響,計算結果會低估斷層走滑分量和高估傾滑分量。

本文基于平面斷層幾何模型探究地形對地震同震形變的影響,然而蘆山地區斷層幾何復雜,為了更全面地理解同震形變的影響,后續將進一步收集同震DInSAR數據,并考慮更精細的斷層幾何模型,以便更準確地研究地形對同震形變的影響。另外,對于復雜的斷層幾何模型,不同構建方法所得結果的差異也值得進一步探索。

致謝：感謝Gharti教授提供SPECFEM-X軟件。