詹自敏 陳連旺 陸遠忠

(中國地震局地殼應(yīng)力研究所，北京 100085)

大華北地殼動力學參數(shù)及三維有限元模型的建立*

詹自敏 陳連旺 陸遠忠

(中國地震局地殼應(yīng)力研究所，北京 100085)

利用華北地區(qū)的斷層數(shù)據(jù)、高程數(shù)據(jù)、波速結(jié)構(gòu)數(shù)據(jù)、彈性常數(shù)以及多種流變機制參數(shù)，建立大華北地區(qū)三維有限元數(shù)值模型，并通過給每個單元賦屬性的方法，實現(xiàn)物性參數(shù)的連續(xù)變化。以模型邊界附近實測GPS數(shù)據(jù)插值作為邊界位移條件，進行數(shù)值模擬。模擬結(jié)果表明，大華北區(qū)域位移場與實測位移場基本一致。

華北;有限元模型;數(shù)值模擬;GPS;地殼動力學

1 引言

大華北地區(qū)多次發(fā)生7級以上強震，給人民的生命財產(chǎn)安全造成了極大的威脅。

有限元的方法是在有限差分的基礎(chǔ)上發(fā)展起來的數(shù)值模擬方法，運用有限元方法建立的地塊實體模型，可以模擬地塊深部應(yīng)力場的變化和強震的孕育的過程。建立大華北地區(qū)三位實體模型，需要比較詳盡的參數(shù)，而華北地區(qū)大量的監(jiān)測數(shù)據(jù)，對于建立華北地區(qū)有限元模型提供了良好的條件。

早在1980年王仁等［1］便對華北地區(qū)進行了二維有限元數(shù)值模擬，采用彈性材料，并用初始應(yīng)力場和應(yīng)力釋放的方法模擬了邢臺、河間、渤海、海城、唐山5次地震的發(fā)生;2003年，文獻［2，3］用非連續(xù)變形和有限元相結(jié)合的方法對華北地區(qū)進行了二維數(shù)值模擬，同樣采用了彈性的本構(gòu)關(guān)系;1998年，劉昌銓等［4］建立了華北的三維有限元彈性模型;1999年及其后，陳連旺等［5，6］分別利用彈性、黏彈性本構(gòu)關(guān)系模擬了華北地區(qū)三維構(gòu)造應(yīng)力場以及強震應(yīng)力場隨時間的變化特征;2005年，劉勉等［7］建立了華北地區(qū)三維模型，分析了華北地區(qū)應(yīng)變及應(yīng)變能密度和震源機制;2007年，劉峽［8］在其博士論文中對華北地區(qū)進行了較詳細的三維有限元數(shù)值模擬，建立了均勻分層和實際分層兩個三維模型以及一個二維模型，使用了接觸單元模擬斷層，分析了華北地區(qū)的形變場和應(yīng)力場，并探討了其動力學機制;2010年，胡勐乾［9］也在其碩士學位論文中進行了華北地區(qū)的三維數(shù)值模擬，建立整個華北地區(qū)的分塊模型，選取了多條弱化斷裂，使用位移場作為邊界條件，模擬華北地區(qū)速度場和應(yīng)力場并進行了斷層性質(zhì)的分析。

然而以往的模型大多未能全面涵蓋本區(qū)的地形、物性、實際斷層幾何特性，有的模型邊界為簡單矩形，未考慮塊體構(gòu)造特性。本文將在建立較完善的華北地區(qū)地殼動力學參數(shù)庫的基礎(chǔ)上，建立華北地區(qū)三維黏彈性數(shù)值模型，并用此模型模擬華北地區(qū)的形變場，探索地震孕育規(guī)律。

2 大華北地區(qū)地殼動力學參數(shù)

2.1 華北構(gòu)造分區(qū)及主要的活動斷裂

根據(jù)中國區(qū)域構(gòu)造可以將華北地區(qū)劃分為鄂爾多斯活動地塊、華北平原活動地塊和魯東黃海活動地塊，其中華北平原活動地塊又可以分為3個三級活動地塊：太行山次級活動地塊、冀魯次級活動地塊和豫淮次級活動地塊［10］。

華北地塊的周圍被深大斷裂所圍限，其內(nèi)部也發(fā)育多條斷裂，主要走向有北北東、北東、北西和近東西向4組［11］，其主要斷裂如表1所示。

2.2 地形數(shù)據(jù)與深部介質(zhì)物性參數(shù)

地形數(shù)據(jù)來源于 http：//srtm.csi.cgiar.org/。地形數(shù)據(jù)包括經(jīng)度、緯度和高程值，高程的基準點為海平面。

地殼上地幔S波速度成果為文獻［12，13］通過面波層析成像的方法得到的結(jié)果。依據(jù)理想各向同性介質(zhì)中彈性常數(shù)剪切模量、泊松比和P波、S波速度之間的定量關(guān)系，計算得到的楊氏模量在三維立體空間內(nèi)畫散點圖如圖1所示。

巖石圈隨深度不同具有不同的溫壓條件，巖石在高溫高壓下會出現(xiàn)一定的流變性質(zhì)。2002年，臧紹先等［14］建立了一種華北巖石圈流變結(jié)構(gòu)的初步模型，給出了3種不同巖石強度的計算公式(表2)。 s-1［15，16］。對于大陸巖石圈應(yīng)變率一般取10-14～10-15s-1，然而深部巖石應(yīng)變率的變化比較復(fù)雜，想要確定深部巖石的應(yīng)變率比較困難，所以本文模型暫不考慮應(yīng)變率隨深度的變化，計算時使用全球平均結(jié)果10-15

圖1 華北地區(qū)楊氏模量分布散點圖Fig.1 Scattered diagram of Young modulus in North China

本文根據(jù)華北巖石組成的一般模型［17，18］，結(jié)合地下巖石的流變特性等性質(zhì)［19］，再綜合利用華北巖石圈三維流變結(jié)構(gòu)的初步模型［14］，給出華北地下巖石分層以及巖石特性(表2)。

3 華北地區(qū)三維有限元模型的建立

3.1 模型邊界與斷層的選取

本研究區(qū)約為北緯30～43°、東經(jīng)105～125°，北邊界依燕山構(gòu)造帶和陰山一帶而行，直至遼東;南邊界以秦嶺-大別山為界，直至黃海中部;西邊界以賀蘭山、六盤山為界;東邊界位于124.5°E附近的海域中。模型邊界全長5 320 km，東西向長度最大為1 660 km，南北向長度最大為1 350 km，深度為200 km。

斷層在模型中的位置如圖2所示。根據(jù)表1中各條斷層的走向、傾向和傾角等數(shù)據(jù)，實現(xiàn)斷層的傾斜，每一條斷層用一個小體表示。

3.2 網(wǎng)格劃分

網(wǎng)格為上層細網(wǎng)格、下層粗網(wǎng)格，模型主要采用六面體網(wǎng)格劃分的方法，這種劃分有兩個難題：一是每一層網(wǎng)格要求的尺寸大小不一致，如何實現(xiàn)層與層之間六面體網(wǎng)格的過渡和轉(zhuǎn)換;二是斷層的傾斜，使得斷層附近不可能建立完全直立的六面體單元，如何即保證網(wǎng)格是六面體，又實現(xiàn)斷層的傾斜。

對于第一個問題，不同層之間不同尺寸的單元網(wǎng)格的連接，采用金字塔過渡單元來解決(圖3 (a))。模型中，結(jié)合表3的分層，0～30 km為第一大層，此層又可分為4個小層，每層約為8 km;40～80 km為第二大層，100～200 km為第三大層，30～40 km和80～100 km為過渡層。

表1 華北地區(qū)主要斷裂一覽表Tab.1 Major faults in North China

表2 模型選用巖石分層及巖石性質(zhì)參數(shù)列表Tab.2 Parameters of rocks in different layers of this model

關(guān)于第二個問題，在斷層附近采取掃略方法實現(xiàn)六面體網(wǎng)格的劃分，而在距離斷層較遠的地方，使用正常的六面體網(wǎng)格劃分(圖3(b)、(c))。

網(wǎng)格劃分的整體效果如圖3(a)所示。模型共有254 156個單元，234 056個節(jié)點。

圖2 模型中斷層的位置及名稱示意圖Fig.2 Location and name of selected faults in the model

圖3 華北模型網(wǎng)格劃分及傾斜斷層示意圖Fig.3 Meshing and the inclined fault of the model

3.3 地形起伏

為了在模型頂面實現(xiàn)地形的起伏，并考慮到地表起伏的高度相對模型的范圍來說較小，一般不會出現(xiàn)夾角過大的單元，所以采用直接改變地表節(jié)點Z坐標值的方法，即將模型地表的每一個節(jié)點的坐標值提取出來，在地形參數(shù)庫中運用插值的方法得到相應(yīng)的高程值，然后改變地表節(jié)點的Z坐標，實現(xiàn)模型的地表高程起伏(圖4，將模型地表的高程值放大50倍的效果圖)。

3.4 物性參數(shù)

使用PRONY模型模擬黏彈性，PRNOY模型的核方程為：

圖4 華北模型地表起伏示意圖Fig.4 Surface topographic undulation of the model

其中G為剪切模量，K為體變模量，G∞和K∞為時間無窮大時剪切模量和體變模量的取值，而Gi和Ki則分別為第i個PRONY元件的剪切和體變模量。nG和nK為PRONY元件的個數(shù)(簡單情況下取1)。t為時間和分別為剪切模量和體變模量的松弛時間。

模擬中用弱化帶模擬斷層，將斷層的物性參數(shù)賦為周邊塊體的十分之一。

本模型與前人模型的重要區(qū)別在于將物性參數(shù)賦值到每個單元之上，實現(xiàn)了物性參數(shù)較為連續(xù)的變化，相對來說更接近實際。因此，本模型介質(zhì)屬性的組數(shù)等于單元個數(shù)，共有254 156組介質(zhì)屬性。由于目前波速結(jié)構(gòu)結(jié)果分辨率的限制，導(dǎo)致一些單元的參數(shù)是相同的。

4 華北地區(qū)GPS位移場模擬與分析

為了驗證模型內(nèi)部材料屬性是否能夠反應(yīng)真實情況，本文將中國地震局第一形變監(jiān)測中心獲得的GPS年平均位移數(shù)據(jù)［20］插值得到三維模型邊界處的位移數(shù)值作為三維模型6個側(cè)面的邊界條件。由于深部運動速率與GPS速率之間的差異尚無明確結(jié)論，作為一種近似，本文所施加的側(cè)面邊界條件未考慮沿深度的變化。模型的底面在垂直方向固定不動，這近似于一種重力均衡的表現(xiàn)。在上述邊界條件作用及考慮重力載荷下，對模型進行靜態(tài)計算，并將計算得到的地表位移場與GPS實測位移場進行比較(圖5，紅色箭頭表示華北地區(qū)GPS年位移場，每一個數(shù)據(jù)值均為一個觀測點的數(shù)值，沒有進行插值;藍色箭頭為數(shù)值模擬結(jié)果)。

為進一步分析華北塊體內(nèi)部的模擬情況，將華北地區(qū)分為西北、中北、東北、西南、中南、東南6塊，分別對每個塊區(qū)的模擬值與GPS觀測值的區(qū)域平均結(jié)果進行比較(表3)。

圖5 華北地區(qū)GPS實測位移場與模型模擬結(jié)果對照圖Fig.5 Comparison between GPS observed displacement values and the simulated values

表3 GPS觀測值與模型模擬值分區(qū)比較Tab.3 Comparsion between GPS observed values and the simulated values in different districts

從圖5及表3可以看出，除河北中東部分地區(qū)差別較大以及山西和河南境內(nèi)部分觀測點方向不一致外，大部分地區(qū)位移場的大小和方向的分布與觀測GPS位移場基本一致，其總體走勢也基本相同，模型模擬出來的位移場與真實的GPS位移場分布具有很高的相似性，除東北區(qū)的角度相差較大外，絕大部分區(qū)域的數(shù)值大小和角度方向均相差不大。

5 結(jié)論

本文建立的大華北地區(qū)三維精細結(jié)構(gòu)模型所模擬出來的華北地區(qū)位移場與GPS觀測位移場符合較好，說明本文所建立的大華北地區(qū)參數(shù)數(shù)據(jù)來源可靠，具有真實的地理信息依據(jù)，能夠真實地反應(yīng)華北地區(qū)深部構(gòu)造狀況。

致謝 感謝李延興、黃忠賢、馬保起研究員在GPS觀測數(shù)據(jù)、深部波速結(jié)構(gòu)數(shù)據(jù)及地形高程數(shù)據(jù)等方面提供的幫助與指導(dǎo)。

1 王仁，等.華北地區(qū)地震遷移規(guī)律的數(shù)學模擬［J］.地震學報，1980，2(1)：32-42.

2 白武明，林邦慧，陳祖安.1976年唐山大震發(fā)生對華北地區(qū)各地塊運動與變形影響的數(shù)值模擬研究［J］.中國科學(D輯)，2003，33(增刊)：99-107.

3 陳祖安，等.1966年以來華北地區(qū)一系列七級大震破裂過程的數(shù)值模擬［J］.地球物理學報，2003，46(3)：373-381.

4 劉昌銓，劉明軍，嘉世旭.利用華北北部深部地球物理資料數(shù)值模擬地殼應(yīng)力場［J］.地震學報，1998，20(3)：240-249.

5 陳連旺，等.華北地區(qū)斷層運動與三維構(gòu)造應(yīng)力場的演化［J］.地震學報，2001，23(4)：349-361.

6 陳連旺，等.1966年邢臺地震引起的華北地區(qū)應(yīng)力場動態(tài)演化過程的三維黏彈性模擬［J］.地震學報，2001，23 (5)：480-491.

7 Liu Mian and Yang Youqing.Contrasting seismicity between the north China and south China blocks：Kinematics and geodynamics［J］.Geophysics Research Letters，2005，32，L12310，doi：1029/2005GL023048.

8 劉峽.華北地區(qū)現(xiàn)今地殼運動及形變動力學數(shù)值模擬［D］.中國科學技術(shù)大學，2007.

9 胡勐乾.并行計算三維數(shù)值模擬在華北地區(qū)現(xiàn)今構(gòu)造變形分析中的應(yīng)用研究［D］.中國地震局地質(zhì)研究所，2010.

10 韓竹軍，等.華北地區(qū)活動地塊與強震活動［J］.中國科學(D輯)，2003，33(增刊)：108-118.

11 李鐵朋，等.華北地區(qū)Ms≥6.5級地震震源斷層參數(shù)的研究［J］.地球物理學進展，2007，22(1)：95-103.

12 黃忠賢，等.中國東部海域巖石圈結(jié)構(gòu)面波層析成像［J］.地球物理學報，2009，52(3)：653-662.

13 Zhongxian Huang，et al.The lithosphere of North China Craton from surface wave tomography［J］.Earth and Planetary Science Letters，2009，288：164-173)

14 臧紹先，等.華北巖石圈三維流變結(jié)構(gòu)的一種初步模型［J］.中國科學(D輯)，2002，32(7)：588-597.

15 魏榮強，臧紹先.大陸巖石圈流變結(jié)構(gòu)研究進展及存在的問題［J］.地球物理學進展，2007，22(2)：359-364.

16 魏榮強，臧紹先.巖石圈流變結(jié)構(gòu)的一種新的應(yīng)變率約束［J］.地球物理學報，2004，47(6)：1 029-1 034.

17 吳宗絮，等.華北大陸地殼——上地幔巖石學結(jié)構(gòu)與演化［J］.巖石礦物學雜志，1994，13(2)：106-115.

18 邢集善，劉建華，趙晉泉.華北板內(nèi)深部構(gòu)造［J］.山西地震，2002，4：3-12.

19 周永勝，何昌榮.地殼主要巖石流變參數(shù)及華北地殼流變性質(zhì)研究［J］.地震地質(zhì)，2003，25(1)：109-122.

20 張靜華，等.用GPS測量結(jié)果研究華北現(xiàn)今構(gòu)造形變場［J］.大地測量與地球動力學，2004，(3)：40-46.

CRUSTAL DYNAMIC PARAMETERS AND 3D FINITE ELEMENT MODEL OF NORTH CHINA

Zhan Zimin，Chen Lianwang and Lu Yuanzhong
(Institute of Crustal Dynamics，CEA，Beijing 100085)

Using the data of fault and terrain，seismic velocity structure and rarious rheological parameters of North China，to build a crustal dynamic parameters library and build a 3D finite element numerical simulation model of this region based on the library.Different from other models，an attribute is set for every element，so that the variation of attribute is continuous.By using the interpolated GPS data as displacement boundary conditions，the calculated displacement field is similar to the observational displacement field.This model has great significance for crustal dynamics study，for example，the research on the stress and strain fields of spatio-temporal evolution in North China，the characteristics of strong earthquakes activities，and the analysis of stress field after strong earthquakes and so on.

North China;finite element model;numerical simulation;GPS;crustal dynamics

1671-5942(2011)Supp.-0028-05

2011-01-14

中央級公益性科研院所基本科研業(yè)務(wù)專項重點項目(ZDJ2009-06);中央級公益性科研院所基本科研業(yè)務(wù)專項重大項目(ZDJ2007-01)

詹自敏，女，1987年生，碩士研究生，主要研究方向為地殼動力學數(shù)值模擬.E-mail：zhanzimin@gmail.com

P315.8

A