《大地电磁学_chp4水平非均匀介质》由会员分享,可在线阅读,更多相关《大地电磁学_chp4水平非均匀介质(61页珍藏版)》请在金锄头文库上搜索。
1、大地电磁学 Geo-electromagnetism Magnetotellurics(MT) 地球物理专业用成都理工大学 Mao Lifeng 2011年3月21日Chp4 水平不均匀介质大地电磁测深理论4.1 二维介质大地电磁场4.1.1 二维介质大地电磁场基本公式及边界条件 4.1.2 二维介质大地电磁场有限元正演 4.2 水平非均匀介质的阻抗张量4.2.1 二维介质的张量阻抗定义4.2.2 二维介质的坐标旋转下的张量阻抗4.2.3 三维介质的张量阻抗 4.3 一维、二维、三维介质的判别4.1 二维介质大地电磁场 二维介质是指介质的电学性质除沿深度方向上变 化外,还沿一个水平方向上有变化
2、,而沿另一个 方向保持不变,这个方向称为走向。 二维介质情况下,大地电磁的解析求解就变得十 分困难,除极少数情况外,一般不能给出解析解 ,只能借助微分方程的数值计算方法求出近似解 。即数值解法,常用的数值解法有:有限元法、 有限差分法、积分方程法、有限体积法和边界元 法等。 计算二维介质的大地电磁场-阻抗-视电阻率和 相位的过程成为二维大地电磁正演。4.1.1 二维介质大地电磁场基本公式及边界条件 二维介质大地电磁场基本公式: 假定y轴指向构造走向,z轴指向地下,介质电性 参数仅是x、z两个坐标的函数。忽略位移电流, 大地电磁场满足的麦克斯韦方程组为:4.1.1 二维介质大地电磁场基本公式及边
3、界条件二维介质极化模式:4.1.1 二维介质大地电磁场基本公式及边界条件 二维介质大地电磁场的亥姆霍茨方程: 当电导率和磁导率是均匀的区域时,两式可简化 为:4.1.1 二维介质大地电磁场基本公式及边界条件 TM模式下的边界条件: 空间介质电导率为0S/m,等 效阻抗为0,所以可以将Hy 分量视为常数,考虑到MT计 算的电磁场分量之比,所以 可将TM上边界z=zmin取为地 空界面,且Hy=1. TE模式下的边界条件: 上边界取高空z=zmin,满足 地下异常体的电磁感应可忽 略不计,其上电磁场水平分 量为常数,可取为1。且有: |zmin|=ymax-ymin。两种模式下,下边界要在地下一定
4、 距离的地方,保证下部介质已属一 维的或均匀半空间的地方。由于其 阻抗容易得到,故可用第三类边界 条件;两侧边界条件应离横向不均 匀异常体足够远的地方,使得该处 电磁场可用一维介质正演得到,所 以可用第二类边界条件。4.1.1 二维介质大地电磁场基本公式及边界条件4.1.2 有限元法正演有限元法的基本思想:将大地电磁场的开域问题离散化 为有限个单元并通过有限个节点相连接的等效集合体。 有限元法利用在每一个单元内假设的近似函数来分片地 表示全求解域上待求的未知场函数,单元内的近似函数 由未知函数在单元各个节点的数值和其插值函数来表达 。这样未知函数在各个节点上的数值成为未知量,使一 个连续的无限
5、自由度的电磁场计算问题变成离散的有限 自由度问题。求出这些未知量,就可计算出整个求解域 上的近似电磁场值。 显然,随着单元数目的增加,也即单元尺寸的缩小,或 者随着单元自由度的增加及插值函数精度的提高,解的 近似程度将不断改进,如果单元是满足收敛要求的,近 似解最后将收敛于精确解。4.1.2 有限元法正演有限元法的步骤: 1. 计算域离散化。对单元编号,包含局部编号 和全局编号。 2. 采用合适的插值函数,单元节点电磁场插值 。插值函数要力求简单,如采用多项式方案 ,但要满足一定的收敛要求。 3. 计算单元刚度矩阵、右端项。 4. 组装整体刚度矩阵和右端项,形成大型带状 、对称线性方程组。 5
6、. 求解上述线性方程组,即得到全部节点场值 。得到计算域空间电磁场。4.1.2 有限元法正演 国内外大地电磁有限元法正演的情况 1971 年Coggon首先将有限单元法应用在电磁法正演模拟中 ,他从电磁场能量最小原理出发,实现了二维地电断面有 限单元法正演计算;1976 年William 等发展了矩形网格剖 分,以解决二维大地电磁测深正演问题,使有限单元法向 前发展了一步;1987 年Wannamaker 等采用矩形单元中的 三角形剖分,形成了二维大地电磁测深正演模拟的经典程 序PW2D。1977 年Reddy首先探讨了有限单元法的三维大 地电磁正演问题; Pridmore D F等人(198
7、1)研究了利用 二次场进行三维模拟时代散度校正;Zyserman 等和 Mitsuhata 等用不同形式的有限单元法在三维大地电磁的正 演模拟中取得了成功。Badea E A等(2001)研究了库仑规 范下的可控源音频大地电磁法的有限元模拟。Yuji Mitsuhata和Toshihiro Uchida(2004)研究了基于T-的三维 大地电磁有限元正演模拟;Kerry Key和Chester Weiss( 2006)研究了二维大地电磁的非结构化网格的自适应有限 单元法并对两种后验误差进行了比较;韩国学者Myung、 Jin Nam, Hee Joon kim(2007)等人实现了矢量有限元模
8、拟 三维大地电磁文章。4.1.2 有限元法正演 国内学者徐世浙对有限单元法进行了深入研究, 发表了一系列有关有限单元法在大地电磁测深正 演计算的论文或专著,网格剖分也由原来的简单 剖分发展到三角单元和三角-矩形综合剖分,提高 了计算精度和速度;刘小军等用有限单元法模拟 了基于二次场的二维大地电磁。对于三维大地电 磁正演,国内采用有限单元法模拟的较少,王若 等用有限单元法正演模拟了三维CSAMT。闫述( 2003)和师学明(2004)应用矢量有限元法模拟 大地电磁响应。4.1.2 有限元法正演 大地电磁有限元法数值模拟的边值问题: (3)、(4)可统一表示成:4.1.2 有限元法正演 边值问题(
9、8)与下列变分问题等价: 有限单元计算计算域非均匀网格剖分4.1.2 有限元法正演单元计算 这里使用双二次插值,采用矩形单元,即在每个 单元上取8个点,4个顶点和4条边的中点,编号及 坐标见图a(母单元)和图b(子单元)。4.1.2 有限元法正演 两个坐标系的关系为:4.1.2 有限元法正演 单元中的可表示成:式中, 是单元8个节点的待定值,对 应的插值函数分别是U,x,z的双二次函数。单元 积分其中, (i , j =1,2,8)4.1.2 有限元法正演 其中: 单元积分 其中,4.1.2 有限元法正演其中,l为边长。将K1e,K2e 和K3e 扩展成全体节点 组成的矩阵,由全部单元相加得:
10、求变分得: 故4.1.2 有限元法正演加上相应的边界条件,解线性方程组,即可得到计算域各 节点场值u。再用数值方法计算垂向导数du/dz,代人极化 模式下的视电阻率关系式,即可得到视电阻率和相位:4.1.2 有限元法正演线性方程组求解(不完全LU 分解预处理 BICGSTAB 算法) 系数矩阵K是对称正定的稀疏矩阵。 算法流程: 1.给定系数矩阵A,右端项b,初始值x0,最大迭代次数 kmax,精度要求tol0,及预优矩阵M,计算r0=b-Ax0 ,并令k=1,r0=r0; 2.若k=tol0,则转到3,否则,停止,输 出xk。 3.ck-1=rT rk-1。若ck-1=0或者wk-1=0,则
11、算法终止,输出 失败信息,否则转4.4.1.2 有限元法正演4. 若k=1,则p1=r0,否则,计算Sk-1=ck-1ak-1/(ck-2wk-1),pk=rk+1+Sk-1(pk-1-wk-1Vk-1). 求解方程Mp=p,计算Vk=Ap,ak=pk-1 /(rTVk),gk=rk-1-akVk。tol=|gk|,若toltol0,则xk=xk-1+akp,否 则,停止输出xk. 4. 求解Mg=g,t=Ag,wk=tTs/(tTt),xk=xk-1 +akpk+wkg,rk=g-wkt,tol=|rk|,k+,转 到2.4.1.2 有限元法正演不完全LU 分解预处理 预处理矩阵M 采用不完
12、全 LU 分解。若系数矩阵A 的顺序主子式矩阵都是非 奇异的,则矩阵A 一定能 进行LU 分解。若A 是一 个大型稀疏矩阵,则它的 因子L 和U 的下、上三角 部分一般来说都是满的矩 阵。若取预处理矩阵 M=LU,则理论上最好(不 需任何迭代),但实际上这 完全行不通。 使用无填充的不完全LU 分解,记为ILU(0)。右图 是一个五点矩阵ILU(0)分 解示意图。4.1.2 有限元法正演ILU(0)算法: For i=2, ,n, do: For k=1, ,i-1 and for A(i,k)是A的非0值,Do: 计算a(i,j)=a(I,k)/a(k,k); For j=k+1, ,n a
13、nd for A(i,k)是A的非0值,Do: 计算a(i,j):=a(i,j)-a(i,k)*a(k,j); End do; End do; End do4.1.2 有限元法正演模拟结果: 例1:2层G型介质,电阻率100/600欧姆米, 厚度1.8km,用双线性插值、双二次插值有限 元法正演与解析解法的结果对比如下图:(a)、(b)为TE模式下的视电阻率和相位,(c)、 (d)为 TM模式下的视电阻率和相位。4.1.2 有限元法正演 二维模型正演结果 1为TE模式,2为 TM模式,图(a)为 视电阻率,图(b) 为相位。(a)(b)小结 二维介质大地电磁正演结果示例理论模型TE模式视电阻率
14、TE模式相位TM模式视电阻率TM模式相位4.2.1 二维介质的张量阻抗定义 4.2.1 二维介质的张量阻抗定义 二维介质的波阻抗、视电阻率、相位 一维、二维波阻抗的区别:4.2.1 二维介质的张量阻抗定义 以上公式是坐标轴的y轴沿构造走向、x轴垂直构造走向时 的结果,当测量坐标轴为地面任意正交方向时,电磁场各 分量是对应上述特殊坐标轴时对应TE、TM分量的两种模式 下的场的矢量和,所以在一般坐标系下,二维介质波阻抗 可表示为:x-y为地面任意正交测量轴,波阻抗也相应地用张量描述。张量阻抗元素可通过观测各分量电磁场数据解线性方程组而得到4.2.2 二维介质的张量阻抗旋转 任意坐标轴方向下的张量阻
15、抗与在y轴选为构造走向的坐标 轴方向下的张量阻抗之间的关系可通过坐标轴旋转时对应的 张量阻抗矩阵旋转变换得到:设测量轴为x-y,构造走向为y轴, 倾向为x轴,倾向与x轴夹角为xita 角,则有:4.2.2 二维介质的张量阻抗旋转4.2.2 二维介质的张量阻抗旋转对两组测量坐标系x-y和x-y,假定x轴与x轴夹角为xita(x轴逆时 针旋转xita角到x轴),则:4.2.2 二维介质的张量阻抗旋转给定任意正交测量轴的张量阻抗结果后,可通 过上述的坐标变换求得二维介质的主轴及主轴 阻抗张量元素。主轴的求法: 1.试探法 将测量坐标逐次旋转一个角度(如5度),计算相应 的阻抗张量元素,直到满足: ,即为主轴相对于测量轴的方 位角。实际工作中存在测量误差、噪声或地下结构 是三维的, 一般不为0,但在非三维 结构下它们是极小值。4.2.2 二维介质的张量阻抗旋转2.解析法 阻抗元素在主轴上有以下特点:4.2.2 二维介质的张量阻抗旋转据此求得方位角后,将该角值代入旋转后的张量阻抗 式中便得到对应两偏振波分量的波阻抗。该方法得到的解共8个,彼此夹角45度,无法得到那个 是走向,那个是倾向。:4.2.2 二维介质的张量阻抗旋转 二维介质张量阻抗的旋转不变量(由前述的坐标 旋转时张量阻