线源频率电磁测深二维正演(二)*

LINE SOURCE FREQUENCY ELECTROMAGNETIC SOUNDING TWO-DIMENSIONAL FORWARD MODELING

Yan Shu
(Xi'an Jiaotong University,Microwave and Optical Communication Institute)
Chen Mingseng
(Xi'an Branch,CCRI)

　　Abstract　Following the development of frequency electromagnetic sounding (including controlled source audio frequency magnetotelluric sounding)if is essential to study three-dimensional case.Line source 2-D modeling is a base of point source 3-D modeling.
　　Keywords　line source,frequency electromagnetic sounding,two-dimensional geoelectric model,finite element

2.3　人工边界条件
　　地球物理场问题是一个分布在无限空间的开放域问题，当用有限元这种场域无法解开放域问题时，需将无限空间截断为有限空间，在截断处设定人工边界，形成外边界。在求频率测深二维正演问题的数值解时，可将外边界取为矩形，如图4所示，在人工边界上确定边界条件，进而求解。

图3　几何关系

图4　坐标系和矩形区域

　　a.　下边界EF：在相当深的地下，认为电磁场已衰减为零。因此将下边界取为第一类边界条件

E_x｜_EF=0。

(31)

　　b.　地中侧边界CE、DF：对于地中侧边界，在离源足够远的地方，也就是远区，可认为是近于均匀的“平面波”^［1］，这时地中侧边界取为第二类边界条件

(32)

　　c.　空中边界CABD：在空中离地下二维地质异常体足够远处，可将大地的影响做为水平分层大地(一维情况)来考虑，此时在空中边界上可采用一维吸收边界条件。设在边界上

BE_x｜_CABD=0，

(33)

式中B是边界算子；设在空中边界上的场为强度1的向外传播的波与经大地反射的波的迭加：

E_x=e^ik₀r＋Re^ik₀r。

(34)

R是反射系数。将上式求导后有

(35)

用乘(34)式后与(35)式相加得：

θ是矢径和边界外法线方向之间的夹角。将上式与(33)式比较后，可知边界算子

故在空中边界上取第三类边界条件

(36)

边界条件(36)式满足一维情况下的Sommerfeld辐射条件这就是吸收边界条件的物理意义。这个边界条件与文献［10］中采用的空中边界条件
当式中修正贝塞尔函数中的宗量比较大时是一致的。
2.4　有限元方程
　　综上所述，与二维地电模型线源频率电磁测深正演问题对应的泛函极值问题为：

(37)

在上式中，关于δ函数的积分为：

(38)

　　用有限无法求解泛函极值问题时，要将研究的场域剖分成互不重迭的有限个数的单元，在每一基本单元内，由结点处的场值(即要求的未知量)表示出的插值函数来逼近各单元内场的分布，由此将求泛函极值问题变成了求多元函数极值的问题。这样可把(37)式中泛函的面积分表示成每个单元上面积分的总和(如M个单元)，将线积分表示成边界单元上边界线积分的总和(如N条边界)，即：

(39)

考虑到三角形单元的灵活性和有源情况下对称性的要求，可采用交叉对称网格(图5)。其中的基本单元可取为三角形6节点单元，即将三角形顶点及各边中点做为节点，排列顺序是i、j、k、I、J、K。

图5　交叉对称网格与三角形6节点单元

　　在三角形6节点单元上取E_x的二次函数

E_x=a₁＋a₂z＋a₃y＋a₄z²＋a₅zy＋a₆y²，

(40)

经过计算可求出式(39)中的面积分为：

(41)

　　为了求式(39)中的线积分，设在一段边界l上电场E_x是边长的二次函数(图6)：

图6　边界段示意图

E_x=b₀＋b₁S＋b₂S²，

经计算后

(42)

　　为了求得所讨论泛函极值问题的解，还应求出式(41)与式(42)之和相对节点上未知场值E^e_xmE^e_xn的偏导数，并令其等于零，得到下列有限元方程：

(源所在节点上)，

(43)

式中FG、FC、W、Δ_m的具体形式见参考文献［8］。为执行式中的求和，要将各单元矩阵扩充为总体矩阵之后，再把相同位置上的元素迭加。
　　在求解最终形成的线性方程组过程中，考虑到下一步三维问题巨大的计算量，我们采用了crount^［9］分解算法。这样可将系数矩阵分块求解，只要计算机容量满足公共三角块的要求，方程组的求解就可顺利进行。

3　计算结果

3.1　E_x单分量视电阻率
　　与以往线源正演计算中使用阻抗视电阻率^［10］不同，我们使用了由单分量电场定义的视电阻率。
　　已知均匀半空间电场E_x分量的表达式为^［11］：

(44)

式中　r——场点至源点的距离。
　　由于k₀<<k₁，K₁(－ik₁r)<<K₁(－ik₀r),故有则视电阻率

(45)

3.2　计算精度
　　表1是假设大地为电阻率ρ=100 Ω^.m的均匀媒质，在频率f=64 Hz时，有限元算法与精确式(44)计算结果的比较。由于实际勘查中习惯使用视电阻率，这两种计算结果已由(45)式转换成了视电阻率。其中有限元法单元的边长约为1／8波长，下边界约是3倍的趋肤深度，上边界是下边界的1.5倍。表1中所列15个测点的总平均相对误差为1.04%。从总体上讲，有限元算法的计算精度低于边界元^［12］。但用边界元法处理多种媒质是困难的，而有限元法则非常适合于解决多种介质的计算问题。

表1　精度比较