GPS約束下川滇地區下地殼拖曳作用及斷裂活動性有限元模擬

劉昌偉 常祖峰 李春光 曾平 馬琳

摘要：提取川滇地區GPS速度場邊界值作為約束條件，依據地質資料，利用ANSYS有限元模擬軟件建立川滇地區三維有限元模型，將川滇地區劃分為4個塊體，利用靜力學分析模塊模擬下地殼拖曳作用及2條主要活動斷裂的走滑速率。結果顯示：模型四周僅加載GPS速度場邊界值，川滇菱形塊體位移速度較實測速度偏差明顯，不滿足此地區地球動力學特征，增加下地殼拖曳作用荷載后速度場模擬結果得到優化，川滇菱形塊體位移速度差值減小;研究模擬下地殼拖曳作用產生的剪切力發現，拖曳作用在小金河斷裂附近由南東轉向正南時才能得到最優模擬結果，暗示其對川滇菱形塊體拖曳方向發生了偏轉。鮮水河、小江斷裂模擬走滑速率分別為85 mm/a，6 mm/a，與實際走滑速率具有較好一致性，優化后可對斷裂帶閉鎖研究提供參考。

關鍵詞：川滇地區;有限元模擬;拖曳作用;斷層走滑速率

中圖分類號：P315727?? 文獻標識碼：A?? 文章編號：1000-0666（2019）03-0385-08

0 引言

川滇地區位于青藏高原東南緣，處于羌塘、巴彥喀拉以及華南塊體交接部位，在印度板塊與歐亞板塊相互碰撞擠壓下其地質構造錯綜復雜。曾融生和孫為國（1992）研究青藏高原下地殼及地幔時發現，東鄰高原的青海、川滇地區上地幔存在低速帶，分析認為下地殼物質向東南向流動并堆積于此;相關數值模擬、重力異常以及橫波分裂等數據均顯示青藏高原東南部下地殼較軟，更易于流動，導致上地殼與地幔存在解耦（Royden et al，1997;熊熊等，2001;Flesch et al，2005）;受到NE向揚子板塊的阻擋，青藏高原物質向ES方向流動，盡管前人提出的下地殼通道流模型如 Couette流以及Poiseuille流等還具有爭議，但川滇地區下地殼作為青藏高原物質EN向流動的通道已基本不是爭論的焦點（Klemperer，2006，Beaumont et al，2004）;部分研究者通過給定地殼粘滯系數模擬得出了較為可信的川滇地區下地殼流動速度，認為其較上地殼運動快約10 mm/a（Wang，2007;曹建玲等，2013），朱守彪和石耀霖（2004）、王輝等（2007）等在川滇地區上、下地殼相互作用的研究過程中提出下地殼產生了拖曳作用。前人針對川滇地區數值模擬取得了豐富的研究成果，主要思路分為2類：一類將川滇地區四周GPS位移作為荷載，斷裂帶設置為軟弱帶，模擬給出應力分布、斷層相互作用等，并與發震規律進行對比分析（王凱英，馬瑾，2004;陳連旺等，2008），再根據模擬應力場分布劃分出了地震危險區（廖思佩等，2016）;另一類將川滇地區劃分為多個塊體，分析各個塊體內部應變狀態與相互作用（蔣鋒云等，2013），或將斷裂帶模擬走滑速率與實際走滑速率進行對比，分析斷裂帶活動性或閉鎖狀態（王閻昭等，2008;申重陽等，2002;趙靜等，2015;李長軍等，2018）。本文結合以上研究思路，利用ANSYS有限元模擬軟件與地質資料建立川滇地區三維模型，以近期GPS觀測數據作為邊界條件，模擬對比下地殼拖曳荷載對模型產生的影響;斷裂設置為軟弱帶，與各塊體間摩擦接觸，將模擬與實際斷裂帶走滑距離進行對比，探討摩擦接觸定義下斷裂帶有限元模擬與真實斷裂的可比性。

1 川滇地區構造特征及模型建立

11 主要斷裂與塊體劃分

李玶和汪良謀（1975），闞榮舉等（1977）首先提出了川滇菱形塊體的概念，以麗江—小金河斷裂為界將川滇菱形塊體劃分為川西北與滇中2個次級塊體（王鎧元等，1987;李鐵明等，2003）。巴彥喀拉塊體以東昆侖斷裂、龍門山斷裂和甘孜—玉樹斷裂為界，早期也被稱為川青塊體（韓渭賓，蔣國芳，2003），一般根據龍日壩斷裂帶和虎牙斷裂帶將其劃分為阿壩次級塊體、馬爾康次級塊體和龍門山次級塊體（徐錫偉等，2008，陳長云等，2013）。滇東塊體活動性斷裂分布較為有限，以大涼山斷裂為界劃分出大涼山次級塊體和川中次級塊體（宋方敏等，2002）。滇西南塊體劃分觀點較多，主流的劃分方式是以瀾滄江—紅河斷裂為界劃分出印支板塊和滇緬泰板塊（蘇有錦，秦嘉政，2001）。綜上，本文將川滇地區劃分成川滇菱形塊體、巴顏喀拉塊體、滇東塊體和滇西南塊體4個一級塊體（圖1）。將連接一級塊體的斷裂定義為一級斷裂，有甘孜—玉樹—鮮水河斷裂、龍門山斷裂、安寧河—則木河—小江斷裂、怒江斷裂及紅河斷裂。次級塊體間接觸的為次級斷裂，有東昆侖斷裂、里塘斷裂、金沙江斷裂、麗江小金河斷裂、南華—楚雄—建水斷裂、南汀河斷裂以及怒江斷裂南段。

12 介質參數與網格劃分

GPS 觀測結果顯示，川滇地區上地殼形變滿足彈性形變的特征（Shen et al，2005）。本文ANSYS模擬采用線彈性靜力學分析模塊，取川滇地區平均上地殼厚度為30 km，參考前人對川滇地區地殼泊松比研究，將上述4個一級塊體賦予不同的密度、彈性模量以及泊松比，斷裂設置為軟弱帶，其力學性質參考可塑狀態粘土（含水），各塊體介質參數如表1所示。川滇地區斷裂以走滑為主，龍門山推覆構造帶產狀具有淺部陡深部緩的特點（張竹琪等，2010），因此本文將龍門山斷裂傾角簡化設置為65°，其余斷裂帶傾角設置為直立，斷裂帶傾角與塊體為摩擦接觸，當滑動速率<1 mm/s且走滑距離小于1 mm時質密石英巖摩擦系數為06（Dieterich，1978;Dieterich，Kilgore，1996），在不影響計算結果的前提下斷裂帶寬度均簡化為2 km，模型中共計12條斷裂帶、15個摩擦面。本文模型均采用四面體網格劃分方法，塊體網格尺寸為40 km，斷裂帶網格尺寸為5 km，斷裂帶與塊體接觸位置進行局部網格加密，共劃分出3 458個網格，28 567個節點（圖2）。

13 邊界條件

本文GPS數據來自中國地震局GNSS服務平臺提供的基礎數據產品 http：//wwwcgpsaccn。中國大陸構造環境監測網絡共建設有260個基準站和2 000多個區域站，其中66個基準站位于研究區（96°～107°E，21°～35°N）內。本文采用ANSYS靜力學分析模塊進行求解，時間在此求解模塊中無實際物理意義，故荷載均設置為1年產生的位移量。建模過程中對地形起伏與部分斷裂構造進行了簡化處理，斷裂模型均采用了非線性摩擦接觸分析，由于具有大量非線性分析，荷載步設置為5，每次加載20%，使其荷載平滑加載以避免難以收斂的問題。為得到最理想的模擬結果，共設計3個模型。模型Ⅰ在不考慮位移隨深度發生變化的基礎上，將川滇地區66個GPS基準站以及142個流動站（共計208個）實測的年均速度值（2013—2015年、2015—2017年2期）作為控制點，利用Matlab插值功能求出模型四周邊界處的位移速度，依次將位移荷載加載在模型四周（圖1b），將模型Ⅰ作為初始參考模型。模型Ⅱ在模型Ⅰ的基礎上使川滇菱形塊體底部受到下地殼拖曳產生的剪切力作用（圖3a），位移方向與川滇菱形塊體長軸方向一致（約140°方向），位移大小通過多次模擬試驗確定理想值。模型Ⅲ在模型Ⅱ的基礎上考慮川滇菱形塊體底部受到的剪切力在小金河斷裂附近發生由東向南偏轉，?偏轉角度與地表實測GPS偏轉角平均值一致約為30°（圖3b），位移值大小同樣通過多次模擬實驗確定。

2 模擬結果

使用GPS位移作為邊界條件加載在模型Ⅰ（圖4a）四周，為了進一步討論其是否滿足川滇地區地球動力學特征，需與GPS實測點位移進行對比驗證。圖4a為GPS實測值（基準站點）與同坐標模擬值對比結果，圖中箭頭指向表示位移矢量方向。從圖中可以看出，巴顏喀拉塊體、滇東塊體以及滇西南塊體模擬值能夠反映出西北部青藏高原的擠壓以及東部四川盆地阻擋的構造格局，與GPS實測值較為接近，而川滇菱形塊體模擬值卻與實測值存在較大偏差，其模擬—實測速度差分布顯示，滇中次級塊體（菱塊南部）差值尤為明顯。

模型Ⅱ（圖4b）在川滇菱形塊體底部添加一個與剪切力等效的位移荷載，加載方向為川滇菱形塊體長軸方向（140°），荷載大小以川滇菱形塊體GPS實測位移的平均值148 mm/a為基數進行調整，經過多次模擬實驗，確定最佳位移量為9 mm/a，但滇中次級塊體殘差并未發生明顯改善。模型Ⅲ（圖4c）保持底部最佳位移量不變，在小金河斷裂附近對荷載進行偏轉調整后，川滇菱形塊體位移速度差分布如圖4左側橙色點所示，較GPS實測東西向、南北向位移速度差在±5 mm/a以內，無系統性誤差并呈現呈正態分布，因此模型Ⅲ為最佳擬合結果。

川滇地區平均震源深度為15 km，截取模型Ⅲ在該深度的等效應力并去除邊界效應影響（圖5）。由圖5可見，巴顏喀拉塊體應力主要集中于北部以及東昆侖斷裂兩端;川滇菱形塊體應力分布由WN向ES遞減，其中小金河以及鮮水河斷裂交匯處出現局部應力集中，滇中次級塊體內基本無應力分布;滇西南塊體應力主要集中于怒江斷裂、南汀河斷裂兩端，塊體南緣速度場分布復雜，局部地區出現應力集中。斷裂帶走滑速率如圖6所示，鮮水河斷裂、則木河斷裂、小江斷裂、龍門山斷裂以及紅河斷裂走滑速率明顯。其中鮮水河斷裂走滑速率由WN至ES表現出先增大后減小的特點，走滑速率變化范圍較大，為11～104 mm/a;則木河斷裂、小江斷裂走滑速率分布較為均勻，分別為26～56 mm/a，52～65 mm/a;龍門山斷裂走滑速度整體偏小并集中于中段，為11～25 mm/a;紅河斷裂走滑速率為32～55 mm/a，集中于與小江斷裂交匯處，其余分段未發生明顯走滑。

3 討論

31 川滇菱形塊體下地殼拖曳作用

模型Ⅰ在排除了加載方式與介質參數設置不合理的情況下，偏差依然存在，筆者認為是川滇菱形塊體的模擬中缺少位移約束所致。下地殼流動對川滇菱形塊體產生的“拖曳作用”可能是缺少的主要約束之一。巴顏喀拉塊體、滇東塊體以及滇西南塊體的模擬值與實測值吻合度較高，“拖曳作用”可能并未在這些塊體產生顯著影響，其作用范圍僅在川滇菱形塊體內，支持川滇菱形塊體下地殼作為青藏高原東南方向流動通道的觀點。“拖曳作用”產生的剪切力方向發生了偏轉，如圖4c所示，川滇菱形塊體的2個次級塊體（川西北塊體、滇中塊體）GPS實測位移方向由SE向偏轉至近S向，偏轉處大致沿小金河斷裂分布。前人研究川滇地區地殼結構時發現鮮水河斷裂與小金河斷裂的交叉處（深紫色區域）地殼泊松比大于03（圖7），認為下地殼存在低速層甚至局部熔融（Xu et al，2007;Wang et al，2010）。在此區域出現2種異常并非巧合，筆者推測軟弱的下地殼物質由青藏高原沿ES向擠入川滇菱形塊體（紅色箭頭），受到四川盆地剛性塊體的“阻撓”，局部因高溫擠壓形成低速層或熔融，剩余物質轉向繼續向南流動（紫色箭頭），進而下地殼牽引作用對川滇菱形塊體產生的剪切力也隨之發生改變（圖7）。與模型Ⅰ對比顯示，模型Ⅲ在川滇菱形塊體底部添加剪切力后，其模擬速度場得到優化，EW向與NS向位移差明顯改善，十分接近GPS實測值。其次，模型Ⅲ川滇菱塊底部剪切力方向由NE向正S偏轉30°后才與GPS實測值擬合最佳，拖曳作用隨下地殼流向改變而發生偏轉。

32 斷層走滑速率對比分析

大地構造應力場與模擬得出的理想應力場主要區別在于斷裂帶局部形成的“閉鎖”，閉鎖位置、程度等均不可測，而閉鎖勢必產生局部應力集中，直接利用模擬得出的應力云圖去討論似乎并不合理。本文以鮮水河斷裂、小江斷裂為例，對斷裂帶模型摩擦滑動速率與實際活動斷層走滑速率進行對比研究。鮮水河斷裂是川滇地區現今強烈活動的一條大規模左旋走滑斷裂帶，附近地震活動頻率高，以惠遠寺為界分為NW，SE兩段。前人對鮮水河斷裂進行了大量研究（表2），其整體走滑速率呈中間快兩端慢的特點，為888～973 mm/a;本文模擬結果如圖6所示，鮮水河斷裂SE段與龍門山斷裂交匯處走滑速率為104 mm/a，SE段均值為85 mm/a。小江斷裂是云南境內地震活動最為頻繁的斷裂帶之一，主要表現為左旋走滑斷裂，毛玉平和韓新民（2003）、施發奇等（2012）以及何宏林等（2002）計算出小江斷裂走滑速率約7 mm/a，本文模擬走滑速率均值為6 mm/a，與前人結果對比較為接近。值得注意的是，在對龍門山斷裂及紅河斷裂等多條斷裂帶的對比中發現大部分模擬走滑速率高于實際走滑速率（未發表數據），筆者認為原因除模型簡化造成的誤差之外，斷裂帶局部閉鎖也是造成實際與模擬斷裂走滑速率不一致的重要因素：若能精細化斷裂帶模型的介質參數與邊界條件，有限元法可以得出更真實有效的模擬結果，為研究斷裂帶閉鎖提供更多手段。

33 地質體仿真模擬的局限

盡管本文利用ANSYS模擬得出了部分結果，但地質體仿真模擬存在的局限主要有3點：①深部地應力數據匱乏。地殼具有漫長的構造演化史，后期的改造均疊加在前期地殼固有的地應力之上，但我們對這部分應力分布及大小知之甚少。巖土、采礦領域在地下1 km內取得了部分地應力成果，但其研究過程已十分困難（趙德安等，2007;汪斌等，2012）。②斷裂帶延伸情況或相互穿插造成的影響，真實地質環境中斷裂往往表現為雁列式、共軛式等不規則分布，局部次生斷裂的加密很可能對模型整體力學性質產生影響，大量隱伏斷裂的存在使得斷裂的延伸情況變得更加復雜，單一摩擦模型很難還原真實斷裂構造過程。③由于運算速度的限制，若在模型內考慮下地殼粘性、塑形、蠕變等非線性分析以及地幔對流、殼幔作用等熱力學條件，運算時間會成倍增加，在不影響地球動力學特征的前提下恰當優化模型變得必要。

4 結論

本文以川滇地區GPS速度場邊界值為約束條件，依據地質資料建立川滇地區三維有限元模型，使用ANSYS靜力學分析模塊進行模擬，主要結論如下：

（1）川滇地區下地殼物質流動對川滇菱形塊體產生“拖曳作用”，其產生的剪切力荷載隨下地殼流動發生了偏轉，仿真模擬考慮“拖曳作用”的影響后能在一定程度上優化川滇地區位移速度場模擬結果。

（2）鮮水河、小江斷裂模擬走滑速率分別為85 mm/a和6 mm/a，與實測走滑速率具有較好一致性，對斷裂優化后可為“閉鎖”研究提供更多手段。

（3）地質體仿真模擬具有局限性，應該從殼內地應力、斷裂帶精細分布以及模型效率方面進行改進。

參考文獻：

曹建玲，王輝，張晶2013青藏高原下地殼流動方式的數值模擬研究[J].地震，33（4）：55-63

陳長云，任金衛，孟國杰，等2013巴顏喀拉塊體東部活動塊體的劃分、形變特征及構造意義[J].地球物理學報，56（12）：4125-4141

陳連旺，李紅，陸遠忠2008汶川MS80地震同震效應的三維非線性數值模擬分析[J].國際地震動態，（11）：10

韓渭賓，蔣國芳2003川青塊體及其向南東方向運動的新證據[J].地震工程學報，25（2）：175-178

何宏林，池田安隆，宋方敏，等2002小江斷裂帶第四紀晚期左旋走滑速率及其構造意義[J].地震地質，24（1）：14-26

蔣鋒云，朱良玉，王雙緒2013川滇地區地殼塊體運動特征研究[J].地震研究，36（3）：263-268

闞榮舉，張四昌，晏鳳桐，等1977我國西南地區現代構造應力場與現代構造活動特征的探討[J].地球物理學報，20（2）：96-109

李長軍，甘衛軍，秦姍蘭，等2018滇中主要活斷層現今活動性研究[J]地震研究，41（3）：381-389

李玶，汪良謀1975云南川西地區地震地質基本特征的探討[J].地質科學，10（4）：308-326

李鐵明，鄧志輝，呂弋培2003川滇地區現今地殼形變及其與強震時空分布的相關性研究[J].中國地震，19（2）：132-147

李延興，楊國華，李智，等2003中國大陸活動地塊的運動與應變狀態[J].中國科學：地球科學，33（增刊1）：65-81

廖思佩，侯強，杜永超2016基于GPS形變資料的川滇地區應力場數值模擬研究[J].大地測量與地球動力學，36（7）：645-649

毛玉平，韓新民2003云南地區強震（M≥6）研究[M].昆明：云南科技出版社

申重陽，王琪，吳云，等2002川滇菱形塊體主要邊界運動模型的GPS數據反演分析[J].地球物理學報，45（3）：352-361

施發奇，尤偉，付云文2012GPS資料揭示的小江斷裂近期運動特征[J].地震研究，35（2）：207-212

宋方敏，李如成，徐錫偉2002四川大涼山斷裂帶古地震研究初步結果[J].地震地質，24（1）：27-34

蘇有錦，秦嘉政2001川滇地區強地震活動與區域新構造運動的關系[J].中國地震，17（1）：24-34

唐文清，陳智梁，劉宇平，等2005青藏高原東緣鮮水河斷裂與龍門山斷裂交會區現今的構造活動[J].地質通報，24（12）：1169-1172

唐文清，劉宇平，陳智梁，等2007基于GPS技術的活動斷裂監測——以鮮水河、龍門山斷裂為例[J].山地學報，25（1）：103-107

汪斌，李新平，周桂龍2012我國大陸實測深部地應力分布規律研究[J].巖石力學與工程學報，31（1）：2875-2880

王輝，曹建玲，張懷，等2007川滇地區下地殼流動對上地殼運動變形影響的數值模擬[J].地震學報，29（6）：581-591

王凱英，馬瑾2004川滇地區斷層相互作用的地震活動證據及有限元模擬[J].地震地質，26（2）：259-272

王鎧元，孫克祥，薛嘯鋒1987川滇交界金河—洱海深斷裂帶的形成機理和構造屬性探討[J].地震研究，10（3）：92-105

王敏，沈正康，牛之俊，等2004現今中國大陸地殼運動與活動塊體模型[J].國際地震動態，（增刊1）：21-32

王閻昭，王恩寧，沈正康2008基于GPS資料約束反演川滇地區主要斷裂現今活動速率[J].中國科學：地球科學，38（5）：582-597

熊探宇，姚鑫，張永雙2010鮮水河斷裂帶全新世活動性研究進展綜述[J].地質力學學報，16（2）：176-188

熊熊，許厚澤，滕吉文2001青藏高原物質東流的巖石層力學背景探討[J].大地測量與地球動力學，21（2）：1-6

徐錫偉，聞學澤，陳桂華，等2008巴顏喀拉地塊東部龍日壩斷裂帶的發現及其大地構造意義[J].中國科學：地球科學，38（5）：529-542

曾融生，孫為國1992青藏高原及其鄰區的地震活動性和震源機制以及高原物質東流的討論[J].地震學報，（增刊1）：534-564

張竹琪，張培震，王慶良2010龍門山高傾角逆斷層結構與孕震機制[J].地球物理學報，53（9）：2068-2082

趙德安，陳志敏，蔡小林，等2007中國地應力場分布規律統計分析[J].巖石力學與工程學報，26（6）：1265-1271

趙靜，江在森，牛安福，等2015川滇菱形塊體東邊界斷層閉鎖程度與滑動虧損動態特征研究[J].地球物理學報，58（3）：872-885

朱守彪，石耀霖2004用遺傳有限單元法反演川滇下地殼流動對上地殼的拖曳作用[J].地球物理學報，47（2）：232-239