尹迪，董培育，石耀霖

(1 中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院 中國科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室, 北京 100049；2 中國地震局地震研究所 地震大地測量重點(diǎn)實(shí)驗(yàn)室， 武漢 430071)

川滇地區(qū)地處青藏高原東南緣，位于西部青藏高原和東部揚(yáng)子大陸地臺(tái)之間的過渡地帶。由于受到印度板塊的推擠作用，青藏高原內(nèi)物質(zhì)側(cè)向逃逸，其東緣受到華南板塊的阻擋作用，東流物質(zhì)向東南向—南向轉(zhuǎn)移，在川滇地區(qū)形成“順時(shí)針”旋轉(zhuǎn)的運(yùn)動(dòng)趨勢[1]。GPS觀測資料進(jìn)一步揭示了該區(qū)域的形變特征[2-3]，在塊體北部主要為東向運(yùn)動(dòng)，中部運(yùn)動(dòng)方向逐漸轉(zhuǎn)為東南向—南向，南部則呈現(xiàn)分流趨勢，即西南側(cè)向西南方向運(yùn)動(dòng)，東南側(cè)向東南方向運(yùn)動(dòng)。川滇地區(qū)復(fù)雜的構(gòu)造變形活動(dòng)，形成了川滇地區(qū)錯(cuò)綜復(fù)雜的斷層，主要發(fā)育有金沙江、怒江、瀾滄江斷裂帶、鮮水河、安寧河、則木河、小江斷裂帶以及紅河斷裂帶等具有深部構(gòu)造背景的大型活動(dòng)斷裂帶。它們組成了川滇塊體的邊界帶[4]。

川滇地區(qū)地震震源深度主要集中在10～20 km深度處的上地殼[5]，地震波速資料揭示該區(qū)域下地殼為速度較低的低速區(qū)[6]，大地?zé)崃鳒y量資料顯示川滇地區(qū)平均大地?zé)崃髦灯遊7]，這些證據(jù)均表明川滇地區(qū)上地殼為較薄的脆性層，而中下地殼為高溫的柔性層，且更容易發(fā)生流動(dòng)。前人對(duì)川滇地區(qū)數(shù)值模擬的研究成果[1, 8-10]均認(rèn)為青藏高原內(nèi)物質(zhì)被側(cè)向擠出過程中，川滇地區(qū)的柔性下地殼流動(dòng)速度較脆性上地殼快，導(dǎo)致脆性上地殼底部受到下地殼呈“喇叭狀”的拖曳力，該拖曳作用是造成繞喜馬拉雅東構(gòu)造結(jié)順時(shí)針轉(zhuǎn)動(dòng)的主要因素，形成川滇地區(qū)現(xiàn)今復(fù)雜的運(yùn)動(dòng)變形場。

川滇地區(qū)地表形變及深部結(jié)構(gòu)的復(fù)雜性，導(dǎo)致區(qū)域強(qiáng)震頻發(fā)，是中國地震最活躍、地震災(zāi)害最為嚴(yán)重的地區(qū)之一[4-5]。2008年汶川地震發(fā)生后，中國地震學(xué)界認(rèn)為地震預(yù)報(bào)不能僅僅依靠經(jīng)驗(yàn)預(yù)報(bào)，應(yīng)該開始探索數(shù)值地震預(yù)報(bào)的方法，石耀霖等[11]提出地震數(shù)值預(yù)報(bào)的路線圖，并具體提出利用二維有限單元法計(jì)算中長期時(shí)間相關(guān)的地震概率預(yù)報(bào)的方法。董培育等[12]利用該方法在青藏高原地區(qū)進(jìn)行初次嘗試，在計(jì)算過程中，將下地殼對(duì)上地殼的拖曳力等效為二維平面應(yīng)力問題中的體力[8]，然而在等效體力的估算中，主要依據(jù)經(jīng)驗(yàn)分區(qū)及試錯(cuò)法的設(shè)定。盡管這種方法可以在一定程度上較好地反映出研究區(qū)域所受的下地殼的拖曳力，但也存在一定的局限性，即可能會(huì)出現(xiàn)不同區(qū)域之間拖曳力的不連續(xù)，同時(shí)采用試錯(cuò)法進(jìn)行分區(qū)并尋找關(guān)鍵節(jié)點(diǎn)上的拖曳力過于依賴人為經(jīng)驗(yàn)。

本文提出一種新的計(jì)算思路，根據(jù)川滇地區(qū)的活動(dòng)地塊構(gòu)造，建立研究區(qū)域的二維彈性有限元模型，以GPS觀測資料為約束，計(jì)算不考慮拖曳力條件下模擬的位移場與實(shí)際GPS觀測的位移場之間的差值，將該差值與坐標(biāo)值進(jìn)行多項(xiàng)式擬合，該擬合關(guān)系式可以在一定程度上反映研究區(qū)域所需要的拖曳力作用，以此確定不同地點(diǎn)所受拖曳力大小的合理取值。該方法建立了拖曳力分布與空間坐標(biāo)的數(shù)學(xué)關(guān)系，保證了拖曳力的連續(xù)性且降低了拖曳力反演試錯(cuò)法的主觀性，為探討川滇地區(qū)地殼運(yùn)動(dòng)變形的動(dòng)力學(xué)機(jī)制和地震數(shù)值預(yù)報(bào)計(jì)算提供了良好基礎(chǔ)。

1 模型及邊界條件

1.1 區(qū)域模型

參考研究區(qū)域地質(zhì)資料以及數(shù)值模擬的成果[13-14]，同時(shí)考慮模型的邊界影響，確定研究區(qū)域范圍為96°～106°E, 22°～34°N，其中由圖1可知，研究區(qū)域西南角為構(gòu)造環(huán)境極其復(fù)雜的阿薩姆構(gòu)造結(jié)區(qū)域，該地區(qū)有限的GPS觀測數(shù)據(jù)顯示塊體運(yùn)動(dòng)方向?yàn)镹NE向，與其東側(cè)(中國境內(nèi))塊體運(yùn)動(dòng)方向(S或SW向)相反。為了減小計(jì)算結(jié)果的影響，本文模型去掉阿薩姆構(gòu)造結(jié)區(qū)域，在西南角大致以國境線為邊界，建立有限元數(shù)值計(jì)算模型展開研究。介質(zhì)分區(qū)參照活動(dòng)地塊劃分的結(jié)果[4, 15-17]，活動(dòng)塊體主要包括：羌塘、巴顏喀拉、川滇、滇西、滇南及華南地塊等Ⅱ級(jí)塊體，其中川滇菱形塊體由于受北東向麗江—小金河斷裂的切割，可進(jìn)一步劃分為川西北和滇中兩個(gè)次級(jí)塊體(圖1)。

GPS觀測數(shù)據(jù)引自文獻(xiàn)[3]。其中，B1：巴顏喀拉塊體，B2：華南塊體，B3：羌塘塊體，B4：川西北塊體，B5：滇中塊體，B6：滇西塊體，B7：滇南塊體。圖1 川滇地區(qū)活動(dòng)地塊及主要斷裂帶Fig.1 Active blocks and faults in the Sichuan-Yunnan region

模型中包括的斷裂帶主要有：甘孜—玉樹斷裂帶、鮮水河斷裂帶、安寧河斷裂帶、則木河斷裂帶、小江斷裂帶、金沙江斷裂帶、紅河斷裂帶、瀾滄江斷裂帶、怒江斷裂帶、南汀河斷裂帶、理塘斷裂帶、楚雄—建水?dāng)嗔褞?、岷江斷裂帶及龍門山斷裂帶。在計(jì)算中，將斷裂帶設(shè)置為具有一定寬度的軟弱帶，其中龍門山斷裂帶是由3條斷層組成，故將其寬度設(shè)為10 km，其余斷層的寬度設(shè)定為5 km。模型網(wǎng)格剖分采用三角形網(wǎng)格，在斷裂帶及邊界處進(jìn)行加密處理，最終得到的有限元模型共包含88 498個(gè)節(jié)點(diǎn)，174 280個(gè)單元。

1.2 材料參數(shù)

本文計(jì)算程序采用的是FEPG自動(dòng)生成的二維有限元彈性平面計(jì)算開源程序包，主要用到的材料參數(shù)為楊氏模量E和泊松比υ。根據(jù)區(qū)域?qū)游龀上裱芯砍晒@取的地殼速度結(jié)構(gòu)信息[6, 18]，可計(jì)算區(qū)域各個(gè)塊體的楊氏模量E和泊松比υ，同時(shí)參考前人的研究成果[9, 19]，將斷層處理為具有一定寬度的軟弱帶，其楊氏模量為其所在塊體的1/3左右，泊松比稍大，各塊體及斷裂帶的具體參數(shù)選取見表1。

表1 介質(zhì)參數(shù)表Table 1 Material parameters

1.3 邊界條件

已有研究表明，GPS觀測給出的現(xiàn)今川滇地區(qū)的地殼運(yùn)動(dòng)圖像與晚第四紀(jì)構(gòu)造變形的圖像一致[20]。假設(shè)板塊運(yùn)動(dòng)在百年尺度范圍內(nèi)保持恒定，在計(jì)算中選用歐亞參考框架下川滇地區(qū)1999—2015年長期GPS觀測平均速度場數(shù)據(jù)作為穩(wěn)定的速度場[3]，將觀測臺(tái)站數(shù)據(jù)插值到模型邊界上，作為邊界約束，具體邊界條件如圖2所示。

圖2 模型網(wǎng)格剖分及邊界條件圖Fig.2 The numerical model with meshs and boundary condition

2 下地殼拖曳力的計(jì)算

2.1 未考慮拖曳力的初步結(jié)果

僅加載GPS位移邊界條件，不考慮柔性下地殼的拖曳作用，計(jì)算結(jié)果如圖3所示。在巴顏喀拉塊體和四川盆地內(nèi)觀測值與模擬值吻合較好，差異較小；川滇菱形塊體及滇南塊體內(nèi)觀測值與模擬值存在明顯差異：1)菱形塊體北區(qū)，觀測值表現(xiàn)為明顯的東南向運(yùn)動(dòng)特征，而模擬值南向運(yùn)動(dòng)分量較小，整體表現(xiàn)為東向運(yùn)動(dòng)為主。2)菱形塊體中南區(qū)，模擬值雖與觀測值方向基本一致，但模擬值的運(yùn)動(dòng)量較觀測值小。3)菱形塊體南區(qū)，觀測值基本為南向運(yùn)動(dòng)，而模擬值仍偏東南向，差異明顯。4)滇南塊體，主要差異表現(xiàn)為模擬值運(yùn)動(dòng)分量較小于觀測值。另外值得注意的是，在滇南地區(qū)觀測值有明顯的東西向分流的運(yùn)動(dòng)趨勢，在西南側(cè)呈向西南方向運(yùn)動(dòng)，在東南側(cè)呈向東南方向運(yùn)動(dòng)。尤其在東喜馬拉雅構(gòu)造結(jié)附近，表現(xiàn)為明顯的順時(shí)針旋轉(zhuǎn)的運(yùn)動(dòng)形態(tài)，而計(jì)算結(jié)果并沒有反映出川滇菱形塊體順時(shí)針旋轉(zhuǎn)的趨勢。綜合觀測結(jié)果與初步計(jì)算結(jié)果表明：在構(gòu)造復(fù)雜的川滇地區(qū)僅考慮脆性上地殼彈性變形是不夠的，因此，需要在彈性模型中考慮柔性下地殼流動(dòng)對(duì)脆性上地殼產(chǎn)生的附加拖曳作用[1, 8]。前人在川滇地區(qū)數(shù)值模擬過程中采用不同的方法考慮了拖曳力的影響，朱守彪和石耀霖[8]采用遺傳算法反演，利用逐漸逼近的方法，進(jìn)行大量迭代計(jì)算搜索得到最佳匹配模型；王輝等[9]利用在滇南地區(qū)邊界上加大其速度場，模擬區(qū)域形變能夠與實(shí)際吻合，上述方法取得了一定成果但缺乏數(shù)學(xué)物理規(guī)律且需要大量試錯(cuò)并不斷重復(fù)計(jì)算。

2.2 拖曳力計(jì)算方法及結(jié)果

本文提出一種新的思路來計(jì)算川滇地區(qū)下地殼流對(duì)上地殼的拖曳作用。將未考慮拖曳力的模擬結(jié)果與實(shí)際GPS觀測資料結(jié)果之間的位移場差值(如圖3(b)所示)記為ΔUx和ΔUy，其與坐標(biāo)x(東西向)，y(南北向)之間的關(guān)系在一定程度上可以反映研究區(qū)域需要施加的拖曳力與坐標(biāo)的關(guān)系。因此嘗試將ΔUx，ΔUy與x，y進(jìn)行多項(xiàng)式擬合并將該擬合關(guān)系式乘以適當(dāng)?shù)谋壤蜃觼矸从逞芯繀^(qū)域所需要的適當(dāng)大小的拖曳力。通過試錯(cuò)法，找到模型中各個(gè)點(diǎn)位上的拖曳力，令最終計(jì)算結(jié)果可與實(shí)測值最佳吻合，具體計(jì)算如下：

本文共選用遠(yuǎn)離邊界區(qū)域的357個(gè)均勻節(jié)點(diǎn)的值(ΔUx，ΔUy，x，y)，進(jìn)行多項(xiàng)式擬合，經(jīng)過嘗試，認(rèn)為3次擬合精度即可滿足計(jì)算需求。x方向上任意點(diǎn)位的擬合公式如式(1)所示，多點(diǎn)位的擬合公式(矩陣形式)如式(2)所示，y方向上的擬合公式同x方向擬合公式類似，如式(3)所示，其中式(2)和式(3)中的n代表點(diǎn)位號(hào)(n=1,2，…，357)。

(1)

(2)

圖3 僅加載GPS邊界條件的計(jì)算結(jié)果Fig.3 The results only with GPS boundary conditions

(3)

經(jīng)計(jì)算得到的x方向和y方向的擬合系數(shù)見表2。

表2 多項(xiàng)式擬合系數(shù)Table 2 Polynomial fitting coefficients

(4)

圖4 加載至模型中的拖曳力Fig.4 Drag force in this model

3 考慮下地殼拖曳作用的計(jì)算結(jié)果

考慮脆性上地殼變形和柔性下地殼拖曳共同作用的計(jì)算結(jié)果與實(shí)際GPS觀測結(jié)果的比較如圖5所示，模擬得到的GPS速度值與觀測值之間可比性較高，先前速度場差異較大的川滇菱形塊體內(nèi)部在加載拖曳力后得到了很大改善，僅在個(gè)別地區(qū)存在誤差(東喜馬拉雅構(gòu)造結(jié)附近的觀測點(diǎn)不列入考慮)。根據(jù)式(5)計(jì)算二者之間的均方根誤差，計(jì)算結(jié)果如表3所示，加載拖曳力后的模擬值均方根誤差明顯小于未加載拖曳力的情況。

圖5 施加拖曳力的模型計(jì)算結(jié)果與觀測值對(duì)比圖Fig.5 Comparison between the observed and the simulated value with drag force in the model

(5)

式中:n是所有加載拖曳力的節(jié)點(diǎn)的個(gè)數(shù)，Δi是每個(gè)節(jié)點(diǎn)的計(jì)算值和觀測值之差，RMES是均方根誤差。

在此基礎(chǔ)上計(jì)算得到川滇地區(qū)的水平主應(yīng)變率場如圖6所示，該計(jì)算結(jié)果與前人根據(jù)GPS觀測速度場數(shù)據(jù)利用最小二乘配置應(yīng)變率計(jì)算方法計(jì)算得到的研究區(qū)域主應(yīng)變場，在大小及方向上均具有較好的一致性[21-22]。進(jìn)一步分析研究區(qū)域應(yīng)變率，大致估計(jì)研究區(qū)域斷層錯(cuò)動(dòng)方向，由圖6可以看出在龍門山斷裂帶上擠壓應(yīng)變明顯高于拉伸應(yīng)變，處于逆沖型環(huán)境。而在川滇菱形塊體邊界帶上大部分區(qū)域表現(xiàn)為拉伸和擠壓應(yīng)變大小相當(dāng)，以走滑型為主。沿著甘孜—玉樹，鮮水河斷裂帶向南延伸，水平主壓應(yīng)變方向由近EW向逐漸變?yōu)镹WW向，以左旋走滑為主；南段安寧河—?jiǎng)t木河段主壓應(yīng)變方向持續(xù)變化為NW向，在小江斷裂帶上變化為NNW向，在其南端，接近NS向，均以左旋走滑為主。金沙江斷裂北段受到近EW向的擠壓和近NS向的拉張，其南段主壓應(yīng)變方向轉(zhuǎn)為NNW向，整體表現(xiàn)為右旋走滑；麗江—小金河斷裂主壓應(yīng)變方向?yàn)镹W向，張應(yīng)變方向?yàn)镹E向，主要表現(xiàn)為逆沖左旋走滑；紅河斷裂北段主壓應(yīng)變方向?yàn)镹W向，隨著斷裂帶向南延伸，主壓應(yīng)變轉(zhuǎn)為近水平向，以右旋走滑為主；瀾滄江斷裂帶北段主壓應(yīng)變方向?yàn)榻麰W向，南段轉(zhuǎn)為NE向，南汀河斷裂主要壓應(yīng)變方向?yàn)镹E向，成左旋走滑。模擬得到的主要斷層的滑動(dòng)性質(zhì)與實(shí)際斷層滑動(dòng)性質(zhì)均一致(具體見表4)，一定程度上驗(yàn)證了計(jì)算的正確性。

表3 計(jì)算結(jié)果與觀測值均方根誤差對(duì)比Table 3 Comparison of root mean square error between the calculated and observed values

藍(lán)色箭頭為壓縮分量，紅色箭頭為拉張分量。圖6 模擬計(jì)算的應(yīng)變率場Fig.6 Simulated strain rate field

表4 川滇地區(qū)主要活動(dòng)斷裂性質(zhì)Table 4 The characteristics of the main faults in the Sichuan-Yunnan region

4 結(jié)論與討論

本文以GPS觀測數(shù)據(jù)為約束條件，依據(jù)川滇地區(qū)構(gòu)造活動(dòng)特征建立二維有限元彈性平面模型，考慮柔性下地殼差異性快速流動(dòng)對(duì)脆性上地殼的拖曳作用，并以等效體力的形式加載至模型中。拖曳力的反演計(jì)算通過擬合彈性模型中計(jì)算結(jié)果與實(shí)測值之間的位移差和坐標(biāo)點(diǎn)之間的關(guān)系，并通過試錯(cuò)法計(jì)算得到，模擬的結(jié)果可以較好地反映川滇地區(qū)的實(shí)際運(yùn)動(dòng)變形和受力情況。本文計(jì)算結(jié)果表明，研究區(qū)域拖曳力可能主要集中在川滇菱形塊體的中北部地區(qū)，且主要為南向；在菱形塊體的南部則可能存在東南向的拖曳力；在其西側(cè)即滇西地區(qū)可能存在西南向的拖曳力。與前人[8，12]反演的川滇地區(qū)下地殼拖曳力方向、大小等結(jié)果基本保持一致。但本研究方法是基于位移殘差與坐標(biāo)的數(shù)學(xué)關(guān)系擬合得到，據(jù)此數(shù)學(xué)關(guān)系，僅需調(diào)節(jié)合理的比例因子，提高了人為試錯(cuò)分區(qū)計(jì)算的效率。同時(shí)該方法可以得到研究區(qū)域內(nèi)連續(xù)漸變的拖曳力分布，相比于人為試錯(cuò)尋找關(guān)鍵節(jié)點(diǎn)上拖曳力的做法，該方法在物理上更加合理，方法上更加客觀，可以為川滇地區(qū)的長期地表形變動(dòng)力學(xué)機(jī)制研究提供一種更簡潔的新思路。

本文是基于二維彈性模型基礎(chǔ)上的一種簡潔推演計(jì)算，僅是對(duì)真實(shí)地質(zhì)模型的一種近似最優(yōu)解。在相關(guān)問題的深入研究中，考慮真實(shí)的柔性下地殼與脆性上地殼之間的拖曳關(guān)系，仍需要通過建立更加接近實(shí)際地殼分層結(jié)構(gòu)的三維模型，采用更加精細(xì)的介質(zhì)和斷層參數(shù)，才能夠得到更加精確的結(jié)果。