戴前偉，張富強(qiáng)，楊坤平，李蕓蕓

（中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083）

電導(dǎo)率分塊線性變化的二維高頻大地電磁法有限元正演模擬

戴前偉，張富強(qiáng)，楊坤平，李蕓蕓

（中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083）

為了適應(yīng)電性參數(shù)在水平方向和垂直方向是連續(xù)變化的實(shí)際需要，這里采用矩形剖分，線性插值的有限元法研究了電導(dǎo)率分塊線性變化的高頻大地電磁高精度，快速正演模擬。首先給出了二維高頻大地電磁測深的變分問題以及電導(dǎo)率分塊線性變化時(shí)的有限元數(shù)值解法，編制了相應(yīng)的有限元程序。利用該程序?qū)訝钅Ｐ瓦M(jìn)行了計(jì)算，并與解析解的結(jié)果做了對比分析，證明了程序的正確性和有效性，然后對典型電導(dǎo)率隨深度線性變化的模型和典型地塹模型進(jìn)行了正演模擬，結(jié)果表明能夠有效地解決電導(dǎo)率在水平和垂直方向分塊線性變化的高頻大地電磁正演問題。

高頻大地電磁法；電導(dǎo)率分塊線性變化；正演模擬

0 前言

隨著中深度地球物理勘探的迫切需要，以美國EH－4電導(dǎo)率成像系統(tǒng)為代表的高頻大地電磁法，在我國地球物理勘探行業(yè)應(yīng)用越來越廣泛［1、2］，其勘探深度在地下1 000m 以內(nèi)，采集的信號頻率范圍在10Hz～64KHz之間，開展高頻大地電磁場的數(shù)值模擬有著重要的理論意義與實(shí)際意義［3］。國內(nèi)、外眾多學(xué)者，都對大地電磁法的正演模擬進(jìn)行了研究：Wannamaker［4］等采用矩形對角線二次劃分的剖分方式研究了二維帶地形問題的大地電磁的有限元正演模擬；陳樂壽［5］等采用有限元法對二維大地電磁正演計(jì)算并提出了一些改進(jìn)的方法；徐世浙［6、7］等采用矩形剖分的方式對電導(dǎo)率分塊線性變化的大地電磁法進(jìn)行了正演計(jì)算；陳小斌［8］研究了大地電磁正演計(jì)算中地形的影響；劉云［9］等采用自適應(yīng)地形四邊形網(wǎng)格剖分進(jìn)行了大地電磁場的模擬；湯井田［10］等采用矩形剖分的有限單元法對高頻大地電磁進(jìn)行了有限元正演計(jì)算；阮百堯等［11、12］采用矩形剖分方式的有限單元法，對電導(dǎo)率分塊線性變化地電斷面電阻率測深進(jìn)行了有限元數(shù)值模擬；歐東新［13］等對二維電導(dǎo)率線性變化測井進(jìn)行了正演計(jì)算，取得了良好的效果。

大地電磁側(cè)重研究區(qū)域地質(zhì)構(gòu)造劃分，而高頻大地電磁（EH－4）側(cè)重于研究局部區(qū)域的電性變化，相比于大地電磁正演模擬的區(qū)域大小，高頻大地電磁的計(jì)算區(qū)域相比要小，網(wǎng)格的大小要比較小的反應(yīng)局部地質(zhì)體電性的有微小變化。此外，高頻大地電磁正演模擬所采用的頻點(diǎn)是固定的，而大地電磁采用的頻率是天然電磁場的隨機(jī)頻點(diǎn)，受野外電磁場干擾比較大。

前人在研究大地電磁或高頻大地電磁有限元模擬時(shí)，大都假設(shè)電導(dǎo)率分塊均勻，但實(shí)際工作中的巖礦石，由于組成、溫度、濕度及孔隙度等的變化，其電性參數(shù)在水平方向和垂直方向是連續(xù)變化的。因此，研究采用矩形剖分，線性插值的有限元法進(jìn)行高頻大地電磁的數(shù)值模擬，有著重要的理論和實(shí)際意義。

1 變分問題

假定地電結(jié)構(gòu)是二維的，選取右旋直角坐標(biāo)系的原點(diǎn)在地面上，y軸正方向垂直向上，z軸方向平行于走向方向（見圖1），即介質(zhì)的電性參數(shù)僅是x、y兩個(gè)坐標(biāo)的函數(shù)。當(dāng)平面電磁波以任意角度向地面入射時(shí)，地下介質(zhì)中的電磁波總以平面波形式，幾乎垂直地向下傳播。這就使電磁場的各分量沿z方向都沒有變化，即z方向的偏導(dǎo)數(shù)為零。根據(jù)Maxwell方程，并考慮到?／?z＝0，得到兩個(gè)獨(dú)立的方程組，并以z分量為準(zhǔn)，分別稱為TE模式和TM模式。

圖1 大地電磁二維正演研究區(qū)域Fig.1 Studied area of 2Dforward modeling in magnetotelluric

（1）TE模式：

式（3）和式（4）可以統(tǒng)一表示為：

對于TE模式：

對于TM模式：

相應(yīng)的邊值問題歸結(jié)為：

方程組（6）與下列變分問題等價(jià)［14］：

2 高頻大地電磁法有限元求解

用有限單元法解變分問題（7）的步驟如下：

（1）網(wǎng)格剖分。首先將求解的二維區(qū)域剖分成矩形單元，如圖2所示（在圖2中，1、2、…代表節(jié)點(diǎn)；①、②、…代表單元）。

圖2 網(wǎng)格剖分及節(jié)點(diǎn)編號示意圖Fig.2 Sketch map of mesh division and node number

（2）線性插值。圖3（a）是母單元示意圖，圖3（b）是子單元示意圖，兩個(gè)單元間的坐標(biāo)變換關(guān)系為：

其中 x0、y0是子單元中點(diǎn)的坐標(biāo)；a、b是子單元的兩個(gè)邊長，微分關(guān)系為：

圖3 母單元、子單元示意圖Fig.3 Sketch map of parent element and sub－element

雙線性插值的形函數(shù)為：

在單元e上對電磁場u和參數(shù)τ、λ、k進(jìn)行雙線性插值。用ui、τi、λi（i＝1、2、3、4）表示各節(jié)點(diǎn)的u、τ、λ，單元中的u、τ、λ可表示為式（11）。

（3）單元分析。根據(jù)有限單元法的原理，對于區(qū)域的積分，可以劃分為各個(gè)矩形單元的積分之和。

式（7）第一項(xiàng)：

其中 K1e＝（k1ij）；k1ij＝k1ji；ue＝ （ui）T；i、j＝1、2、3、4

式（7）第二項(xiàng)：

其中 K2e＝（k2ij）；k2ij＝k2ji；ue＝1、2、3、4

式（7）第三項(xiàng)：

式（2）右側(cè)最后一項(xiàng)線積分只對CD邊界進(jìn)行，但單元的23邊落在CD邊界上時(shí)，線積分為式（14）。

其中 K3e＝（k3ij）；k3ij＝k3ji；ue＝ （ui）T；i、j＝1、2、3、4

（4）總體合成。將以上各單元積分相加，地單元 （i，j）的Fe（u），將各單元的Fe（u）相加，便得到總體系數(shù)矩陣的F（u）為式（15）。

（5）視電阻率和阻抗相位的計(jì)算。令式（15）為零，得線性方程組，加入上邊界條件u｜AB＝1，解線性代數(shù)方程組即可得各節(jié)點(diǎn)的u值，它代表各節(jié)點(diǎn)的（對TE模式）或（對TM模式）。在求得地表各節(jié)點(diǎn)的u后，再利用數(shù)值方法求出場值沿地表的偏導(dǎo)數(shù)?u／?y，代入式（16）、式（17），便可計(jì)算視電阻率和相位。

TE模式：

3 模型的計(jì)算

3.1 模型1－水平層狀模型

水平均勻?qū)訝钅Ｐ头炙膶樱旱谝粚与娮杪蕿?00Ω·m，厚 度 為 50m；第 二 層 電 阻 率 為1 000Ω·m，厚度為100m；第三層電阻率為500Ω·m，厚度為400m；基底電阻率為10Ω·m。正演所采用的頻率為10Hz～10 000Hz。

用基于Fortran語言程序編寫的電導(dǎo)率分塊線性變化的高頻大地電磁進(jìn)行模擬計(jì)算，得出視電阻率和相位值并繪制出圖4。如圖4所示，在所研究的頻點(diǎn)上，程序計(jì)算的結(jié)果與解析解都很接近，誤差在3%以內(nèi)。這說明程序在TE、TM二種模式下，計(jì)算均勻?qū)訝钅Ｐ头矫娴恼_性和有效性。在TM模式下，視電阻率值稍大于解析解，相位值稍小于解析解；在TE模式下，電阻率稍小于解析解，相位值稍大于解析解?？v觀兩種模式下的誤差，高頻段相對明顯，高頻段的誤差比低頻段的誤差要稍大些，如果網(wǎng)格間距能夠更小，更易于接近解析解。

圖4 模型1曲線圖Fig.4 Curves of model 1

3.2 模型2－電導(dǎo)率線性變化的層狀模型

圖5 為三層地電模型示意圖，第一層的電阻率ρ1＝100Ω·m，厚度為h1＝30m；第二層的厚度h2＝150m；第三層的電阻率ρ3＝20Ω·m。其中第二層電阻率變化過渡層，其電阻率是深度線性關(guān)系為：

圖6是用所編制的程序計(jì)算的TE及TM模式的視電阻率和相位擬斷面圖。從圖6上可以看出：電阻率值隨著頻率的減小，從約100Ω·m逐漸線性減小到20Ω·m，與計(jì)算模型的電阻率變化相一致。此外，無論是TE模式還是TM模式，都能清楚地反映出中間層電阻率線性減小的變化和地塹構(gòu)造形成的異常。尤其在TE模式的視電阻率和相位的擬斷面圖上，更能反映中間層電阻率線性變化的趨勢。在TM模式下的視電阻率斷面圖上，中間層能明顯地反應(yīng)出來。但是在深度方向上，形態(tài)被明顯向下拉伸，比較而言，TM模式的相位斷面圖對模型的反應(yīng)不如TE模式的效果好。

圖5 模型2Fig.5 The parameters of model 2

圖6 模型2視電阻和相位擬斷面圖Fig.6 The apparent resistivity and phase pseudo－section map of model 2

3.3 模型3－地塹模型

圖7 為地塹模型，區(qū)域大小為1 200m×500m。第一層的電阻率為10Ω·m，厚度為50m；在厚度為200m的中間存在一個(gè)地塹，長為150m，寬為150m；其中第二層電阻率變化過渡層，其電阻率是深度線性關(guān)系為：

具體參數(shù)詳見圖7。

圖8是用所編制的程序計(jì)算的TE，TM模式的視電阻率和相位擬斷面圖。從圖8上可以看出：電阻率值隨著頻率的減小從約10Ω·m逐漸線性減小到100Ω·m。在碰到地塹構(gòu)造時(shí)，電阻率曲線明顯發(fā)生變化，形成一個(gè)凹陷區(qū)域，與計(jì)算模型的電阻率變化相一致。無論是TE模式還是TM模式，都能清楚地反映出第二層電阻率線性減小的變化和地塹結(jié)構(gòu)形成的異常。相比之下TE和TM模式的視電阻率擬斷面圖，都能夠清楚地反映出地塹構(gòu)造模型的電阻率實(shí)際分布情況。在TM模式下的相位斷面圖上，第二層能明顯地反應(yīng)出來，在地塹位置形成了向深部延伸的低阻異常。比較而言，TE模式的相位斷面圖對模型的反應(yīng)不如TM模式的效果好。

圖7 模型3Fig.7 The parameters of model 3

圖8 模型3視電阻率和相位擬斷面圖Fig.8 The apparent resistivity and phase pseudo－section map of model 3

4 結(jié)論

通過對電導(dǎo)率分塊線性變化二維高頻大地電磁的正演計(jì)算，可以得出以下幾點(diǎn)認(rèn)識：

（1）通過模型1的計(jì)算和誤差分析表明，計(jì)算程序可靠，準(zhǔn)確，能夠有效地進(jìn)行高頻大地電磁正演模擬，能夠有效地解決電導(dǎo)率在水平和垂直方向分塊線性變化的高頻大地電磁正演問題。

（2）兩種模式都能清楚地反映模型中異常體形成的異常。相比較而言，TE模式和TM模式的視電阻率斷面圖都能對模型的電阻率分布進(jìn)行正確的反應(yīng)，但是TM模式的相位擬斷面圖比TE模式更能反應(yīng)模型的電阻率分布情況。

［1］ 鄧居智，李紅星，楊海燕，等.高頻大地電磁法在長大深埋隧道勘察中的應(yīng)用研究［J］.工程勘察，2010（9）：85.

［2］ 林昀.基于高頻大地電磁法對三峽某隧道不良地質(zhì)體的勘察［J］.工程地球物理學(xué)報(bào)，2010，7（4）：456.

［3］ 王燁.基于矢量有限元的高頻大地電磁法三維數(shù)值模擬［D］.長沙：中南大學(xué)，2008.

［4］ WANNAMAKER P E，STODT J A，RIJO L.Two－dimensional topographic responses in magnetotelluric model using finite elements［J］，Geophysics，1986（51）：2131.

［5］ 陳樂壽.有限元法在大地電磁測深正演計(jì)算中的應(yīng)用與改進(jìn)［J］.石油勘探，1981，20（3）：83.

［6］ 徐世浙，于濤，李予國，等.電導(dǎo)率分塊連續(xù)變化的二維 MT有限元模擬［J］.高校地質(zhì)學(xué)報(bào)，1995，1（2）：65.

［7］ 徐世浙.電導(dǎo)率分段線性變化的水平層的點(diǎn)電源電場的數(shù)值解［J］.地球物理學(xué)報(bào)，1986，21（1）：84.

［8］ 陳小斌.MT二維正演計(jì)算中地形影響的研究［J］.石油物探，2000，39（3）：112.

［9］ 劉云，王緒本。大地電磁自適應(yīng)地形有限元正演模擬［J］.地震地質(zhì)，2010，32（3）：382.

［10］湯井田，王燁，杜華坤，等.高頻大地電磁法有限元數(shù)值模擬［J］.物探化探計(jì)算技術(shù)，2009，31（4）：297.

［11］阮百堯，徐世浙.電導(dǎo)率分塊線性變化二維地電斷面電阻率測深有限元數(shù)值模擬［J］.中國地質(zhì)大學(xué)學(xué)報(bào)，1998，23（3）：303.

［12］阮百堯，熊彬.電導(dǎo)率連續(xù)變化的三維電阻率測深有限元模擬［J］.地球物理學(xué)報(bào)，2002，45（1）：131.

［13］歐東新，阮百堯，王家林.電導(dǎo)率線性變化測井二維有限單元法數(shù)值模擬［J］.同濟(jì)大學(xué)學(xué)報(bào)，2005，33（5）：692.

［14］徐世浙.地球物理中的有限單元法［M］.北京：科學(xué)出版社，1994.

P 631.3＋25

A

10.3969／j.issn.1001－1749.2012.05.09

1001—1749（2012）05—0552—07

中南大學(xué)碩士生學(xué)位論文創(chuàng)新資助項(xiàng)目（2011ssxt055）

2012－01－09 改回日期：2012－06－16

戴前偉（1968－），男，湖南漣源人，博士，教授，從事電磁法方法及理論、工程地球物理勘探等方向的研究。