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

詹自敏 陳連旺 陸遠忠

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

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

詹自敏 陳連旺 陸遠忠

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

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

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

1 引言

大華北地區多次發生7級以上強震，給人民的生命財產安全造成了極大的威脅。

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

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

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

2 大華北地區地殼動力學參數

2.1 華北構造分區及主要的活動斷裂

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

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

2.2 地形數據與深部介質物性參數

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

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

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

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

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

3 華北地區三維有限元模型的建立

3.1 模型邊界與斷層的選取

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

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

3.2 網格劃分

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

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

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

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

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

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

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

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

3.3 地形起伏

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

3.4 物性參數

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

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

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

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

本模型與前人模型的重要區別在于將物性參數賦值到每個單元之上，實現了物性參數較為連續的變化，相對來說更接近實際。因此，本模型介質屬性的組數等于單元個數，共有254 156組介質屬性。由于目前波速結構結果分辨率的限制，導致一些單元的參數是相同的。

4 華北地區GPS位移場模擬與分析

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

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

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

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

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

5 結論

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

致謝 感謝李延興、黃忠賢、馬保起研究員在GPS觀測數據、深部波速結構數據及地形高程數據等方面提供的幫助與指導。

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

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

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

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

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

6 陳連旺，等.1966年邢臺地震引起的華北地區應力場動態演化過程的三維黏彈性模擬［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 劉峽.華北地區現今地殼運動及形變動力學數值模擬［D］.中國科學技術大學，2007.

9 胡勐乾.并行計算三維數值模擬在華北地區現今構造變形分析中的應用研究［D］.中國地震局地質研究所，2010.

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

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

12 黃忠賢，等.中國東部海域巖石圈結構面波層析成像［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 臧紹先，等.華北巖石圈三維流變結構的一種初步模型［J］.中國科學(D輯)，2002，32(7)：588-597.

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

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

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

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

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

20 張靜華，等.用GPS測量結果研究華北現今構造形變場［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

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

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

P315.8

A