《大地电磁学_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 二维介质的张量阻抗旋转 二维介质张量阻抗的旋转不变量（由前述的坐标 旋转时张量阻
