CN110068873A - 一种基于球坐标系的大地电磁三维正演方法 - Google Patents

Abstract

本发明公开了一种基于球坐标系的大地电磁三维正演方法，属于地球物理勘探技术领域，包括以下步骤：a、建立大地电磁控制方程；b、在球坐标系中，划分成若干个小的倒立四棱柱网格单元；c、设置球坐标模型参数，构建磁场离散的球坐标交错网格单元；d、对单个频率进行循环，将大地电磁控制方程在球坐标交错网格单元进行数值离散；e、求解线性方程组；f、由球坐标张量阻抗公式计算阻抗，再带入卡尼亚视电阻率计算公式求取测点处的视电阻率和相位。本发明能够有效克服地球曲率对大地电磁三维深部探测的干扰，避免因地球曲率所带来的计算误差，适用于大地电磁三维正演数值模拟，能够与地球模型更好的匹配，精确度高。

Description

一种基于球坐标系的大地电磁三维正演方法

技术领域

本发明涉及到地球物理勘探技术领域，尤其涉及一种基于球坐标系的大地电磁三维正演方法。

背景技术

大地电磁测深是探测地球深部电性结构的主要方法，通过在地表同步观测电场、磁场分量，定性或定量的分析获取地球一定深度范围内的电性结构模型，在油气勘探、固体矿产资源勘察、深部地质结构探测、地热和地下水资源调查、地震预测和地质灾害防治等领域应用广泛。

正演模拟是地球物理方法的核心，大地电磁三维正演问题的本质是求复杂介质的三维电磁扩散方程的数值解。大地电磁三维正演主要计算方法有积分方程法、有限差分法和有限单元法等，其中有限差分法以实现相对简单、计算量较小、速度较快的优点，在三维正演模拟中使用最为广泛。(Mackie，1993，具体参见：Mackie R L.Three-dimensionalmagnetotelluric modeling using difference equations-Theory and comparisons tointegral equation solutions[J].Geophysics，1993，58(2)：215-226；谭捍东等，2003，具体参见：大地电磁法三维交错采样有限差分数值模拟[J].地球物理学报，2003，46(5)：705-711)。

这些传统的大地电磁三维正演模拟技术都是基于笛卡尔坐标系，将有曲率的地球局部区域模型拉伸为方形的六面体块状区域，忽略了地球曲率的影响。

随着仪器和计算机技术的全面发展，大地电磁在地球深部探测中发挥着重要的作用，如美国地球透镜计划、欧洲地球探测计划、加拿大岩石圈探测计划以及我国的深部探测技术与实验研究专项等，都包含了大量研究地球深部电性结构的大地电磁测深工作。其中的诸多大地电磁研究深度都深达上地幔，剖面长度可达数千公里，这种探测尺度情况下可能已无法忽略地球曲率的影响，继续采用现有的笛卡尔坐标系三维正演模拟将可能会造成较大误差。

基于球坐标系下有限差分的地磁测深三维正演；李建平；翁爱华；李世文；李大俊；李斯睿；杨悦；唐裕；张艳辉；吉林大学学报(地球科学版)；2018-03-26；该期刊文献公开了：“为了计算全球尺度电磁感应的响应，本文介绍地磁测深频率域三维正演。正演算法采用球坐标系下的交错网格有限差分方法，从Maxwell方程的积分形式出发，采用PARDISO对离散后的方程组求解，避免了迭代求解的散度校正。为了验证本文结果的正确性和精度，与前人的有限元和有限差分方法进行了对比，一维层状模型的三维交错网格有限差分数值结果和解析解相对误差小于5％，双半球模型的计算结果与前人的计算结果完全吻合。三维"棋盘模型"计算表明磁场分量对异常体的大小和位置具有很好的分辨能力”。

该期刊文献公开的基于球坐标系下有限差分的地磁测深三维正演，其中地磁测深是只观测和分析磁场信号来探测地球深部电性结构，其有效信号频率范围一般介于10⁵～10⁷秒，探测深度主要从约200千米的地球上地幔至中下地幔范围；因此，并不适用于大地电磁三维正演数值模拟，不能与地球模型更好的匹配，精确度低。

发明内容

本发明为了克服上述现有技术的缺陷，提供一种基于球坐标系的大地电磁三维正演方法，本发明能够有效克服地球曲率对大地电磁三维深部探测的干扰，避免因地球曲率所带来的计算误差，适用于大地电磁三维正演数值模拟，能够与地球模型更好的匹配，精确度高。

本发明通过下述技术方案实现：

一种基于球坐标系的大地电磁三维正演方法，其特征在于，包括以下步骤：

a、建立大地电磁控制方程；

b、在球坐标系中，沿r、θ、三个坐标轴方向，分别用若干平行的球面以不同的间距划分成若干个小的倒立四棱柱网格单元；

c、设置球坐标模型参数，包括网格节点坐标、单元、节点编号和单元电阻率，构建磁场离散的球坐标交错网格单元；

d、对单个频率进行循环，将步骤a中的大地电磁控制方程在步骤c中的球坐标交错网格单元进行数值离散，得到系数矩阵和右端项，组装到线性方程组中；

e、求解线性方程组，迭代过程中开展散度校正，计算得到各个节点的场值；

f、根据大地电磁场源，由球坐标张量阻抗公式计算阻抗，再带入卡尼亚视电阻率计算公式求取测点处的视电阻率和相位。

所述步骤a中，大地电磁控制方程采用忽略位移电流后的麦克斯韦尔积分方程形式。

所述步骤b中，划分成若干个小的倒立四棱柱网格单元是指沿θ轴方向剖分成N_θ段，每段的编号i沿θ轴方向序号递增，i＝1,2,…,N_θ，网格弧度为Δθ(i)(1,...,N_θ)；沿轴方向被剖分成段，网格弧度为沿r轴方向被剖分成N_r段，网格间距为Δr(k)(1,...,N_r)。

所述步骤f中，大地电磁场源，分解为两个正交源场等效作用的结果，表征为Sθ和两个极化模式。

本发明的基本原理如下：

采用麦克斯韦尔积分方程形式建立大地电磁控制方程，基于球坐标系交错网格剖分，在球坐标系中采用倒立四棱柱网格单元进行网格剖分，构建磁场离散的球坐标交错网格。采用预条件双共轭梯度迭代法解线性方程组求取每个采样点上的磁场分量值，代入球坐标张量阻抗公式求得球坐标系大地电磁阻抗张量，再由卡尼亚视电阻率计算公式求取地面测点处的视电阻率和相位。

本发明的有益效果主要表现在以下方面：

1、本发明，“a、建立大地电磁控制方程；b、在球坐标系中，沿r、θ、三个坐标轴方向，分别用若干平行的球面以不同的间距划分成若干个小的倒立四棱柱网格单元；c、设置球坐标模型参数，包括网格节点坐标、单元、节点编号和单元电阻率，构建磁场离散的球坐标交错网格单元；d、对单个频率进行循环，将步骤a中的大地电磁控制方程在步骤c中的球坐标交错网格单元进行数值离散，得到系数矩阵和右端项，组装到线性方程组中；e、求解线性方程组，迭代过程中开展散度校正，计算得到各个节点的场值；f、根据大地电磁场源，由球坐标张量阻抗公式计算阻抗，再带入卡尼亚视电阻率计算公式求取测点处的视电阻率和相位”，相较于背景技术中“基于球坐标系下有限差分的地磁测深三维正演”这篇期刊文献而言，尽管正演模拟都基于球坐标系，但地磁测深法与大地电磁测深法属于完全不同的地球物理方法。其中地磁测深是只观测和分析磁场信号来探测地球深部电性结构的方法，其有效信号频率范围一般介于10⁵～10⁷秒，探测深度主要从约200千米的地球上地幔至中下地幔范围；而对于本发明所涉及的大地电磁测深法，是需要同步观测和分析电场与磁场来探测地球电性结构的方法，其有效信号频率范围一般介于10^-4～10⁵秒，探测深度范围从地表至不超过约200千米深度的地球上地幔范围。两种方法所研究的深度范围和地球物理场类型均不相同，因此，地磁测深和大地电磁测深本质上属于不同的地球物理方法，地磁测深不适用于大地电磁三维正演数值模拟。地磁测深因其研究区域范围较大，一直以来都是基于球坐标系，而大地电磁测深法过去认为研究尺度较小，为简化正演模拟，在笛卡尔坐标系中将有曲率的地球局部区域模型拉伸为方形的六面体块状区域，从而忽略了地球曲率的影响；本发明针对大地电磁测深法探测尺度较大情况下存在的地球曲率影响问题，同时观测和分析磁场及电场分量，能够有效克服地球曲率对大地电磁三维深部探测的干扰，避免因地球曲率所带来的计算误差，适用于大地电磁三维正演数值模拟，能够与地球模型更好的匹配，精确度高。

2、本发明，在球坐标系中，沿r、θ、三个坐标轴方向，分别用若干平行的球面以不同的间距划分成若干个小的倒立四棱柱网格单元；设置球坐标模型参数，包括网格节点坐标、单元、节点编号和单元电阻率，构建磁场离散的球坐标交错网格单元；对单个频率进行循环，将大地电磁控制方程在球坐标交错网格单元进行数值离散，得到系数矩阵和右端项，组装到线性方程组中；针对传统笛卡尔坐标系六面体剖分与实际大地电磁模型存在的拉伸效应，采用球坐标系中交错倒立四棱柱网格单元剖分方式，能够更好的匹配大地电磁实际模型，避免坐标拉伸投影所引入的误差，从而能够克服传统大地电磁正演方法存在的计算偏差，具有符合地球形态的三维大地电磁正演模拟的特定效果。

3、本发明，较现有技术而言，针对传统笛卡尔坐标系的阻抗张量公式存在的误差问题，构建了球坐标张量阻抗公式，使得求得的阻抗更加精准；对工程应用尺度的音频大地电磁、地壳尺度的大地电磁以及上地幔尺度的长周期大地电磁均具有良好的适用效果。

附图说明

下面将结合说明书附图和具体实施方式对本发明作进一步的具体说明：

图1是本发明基于球坐标系的大地电磁三维正演模拟流程图；

图2是本发明的球坐标交错网格倒立四棱柱网格单元示意图；

图3是本发明具体实例1对国际标准测试DTM1模型正演视电阻率与现有笛卡尔坐标对比图；

图4是本发明具体实例1对国际标准测试DTM1模型正演相位与现有笛卡尔坐标对比图；

图5是本发明具体实例2对放大后的国际标准测试DTM1模型正演视电阻率与现有笛卡尔坐标对比图；

图6是本发明具体实例2对放大后的国际标准测试DTM1模型正演相位与现有笛卡尔坐标对比图。

具体实施方式

实施例1

一种基于球坐标系的大地电磁三维正演方法，包括以下步骤：

a、建立大地电磁控制方程；

b、在球坐标系中，沿r、θ、三个坐标轴方向，分别用若干平行的球面以不同的间距划分成若干个小的倒立四棱柱网格单元；

c、设置球坐标模型参数，包括网格节点坐标、单元、节点编号和单元电阻率，构建磁场离散的球坐标交错网格单元；

d、对单个频率进行循环，将步骤a中的大地电磁控制方程在步骤c中的球坐标交错网格单元进行数值离散，得到系数矩阵和右端项，组装到线性方程组中；

e、求解线性方程组，迭代过程中开展散度校正，计算得到各个节点的场值；

f、根据大地电磁场源，由球坐标张量阻抗公式计算阻抗，再带入卡尼亚视电阻率计算公式求取测点处的视电阻率和相位。

“a、建立大地电磁控制方程；b、在球坐标系中，沿r、θ、三个坐标轴方向，分别用若干平行的球面以不同的间距划分成若干个小的倒立四棱柱网格单元；c、设置球坐标模型参数，包括网格节点坐标、单元、节点编号和单元电阻率，构建磁场离散的球坐标交错网格单元；d、对单个频率进行循环，将步骤a中的大地电磁控制方程在步骤c中的球坐标交错网格单元进行数值离散，得到系数矩阵和右端项，组装到线性方程组中；e、求解线性方程组，迭代过程中开展散度校正，计算得到各个节点的场值；f、根据大地电磁场源，由球坐标张量阻抗公式计算阻抗，再带入卡尼亚视电阻率计算公式求取测点处的视电阻率和相位”，相较于背景技术中“基于球坐标系下有限差分的地磁测深三维正演”这篇期刊文献而言，尽管正演模拟都基于球坐标系，但地磁测深法与大地电磁测深法属于完全不同的地球物理方法。其中地磁测深是只观测和分析磁场信号来探测地球深部电性结构的方法，其有效信号频率范围一般介于10⁵～10⁷秒，探测深度主要从约200千米的地球上地幔至中下地幔范围；而对于本发明所涉及的大地电磁测深法，是需要同步观测和分析电场与磁场来探测地球电性结构的方法，其有效信号频率范围一般介于10^-4～10⁵秒，探测深度范围从地表至不超过约200千米深度的地球上地幔范围。两种方法所研究的深度范围和地球物理场类型均不相同，因此，地磁测深和大地电磁测深本质上属于不同的地球物理方法，地磁测深不适用于大地电磁三维正演数值模拟。地磁测深因其研究区域范围较大，一直以来都是基于球坐标系，而大地电磁测深法过去认为研究尺度较小，为简化正演模拟，在笛卡尔坐标系中将有曲率的地球局部区域模型拉伸为方形的六面体块状区域，从而忽略了地球曲率的影响；本发明针对大地电磁测深法探测尺度较大情况下存在的地球曲率影响问题，同时观测和分析磁场及电场分量，能够有效克服地球曲率对大地电磁三维深部探测的干扰，避免因地球曲率所带来的计算误差，适用于大地电磁三维正演数值模拟，能够与地球模型更好的匹配，精确度高。

实施例2

一种基于球坐标系的大地电磁三维正演方法，包括以下步骤：

a、建立大地电磁控制方程；

b、在球坐标系中，沿r、θ、三个坐标轴方向，分别用若干平行的球面以不同的间距划分成若干个小的倒立四棱柱网格单元；

c、设置球坐标模型参数，包括网格节点坐标、单元、节点编号和单元电阻率，构建磁场离散的球坐标交错网格单元；

d、对单个频率进行循环，将步骤a中的大地电磁控制方程在步骤c中的球坐标交错网格单元进行数值离散，得到系数矩阵和右端项，组装到线性方程组中；

e、求解线性方程组，迭代过程中开展散度校正，计算得到各个节点的场值；

f、根据大地电磁场源，由球坐标张量阻抗公式计算阻抗，再带入卡尼亚视电阻率计算公式求取测点处的视电阻率和相位。

所述步骤a中，大地电磁控制方程采用忽略位移电流后的麦克斯韦尔积分方程形式。

实施例3

一种基于球坐标系的大地电磁三维正演方法，包括以下步骤：

a、建立大地电磁控制方程；

b、在球坐标系中，沿r、θ、三个坐标轴方向，分别用若干平行的球面以不同的间距划分成若干个小的倒立四棱柱网格单元；

c、设置球坐标模型参数，包括网格节点坐标、单元、节点编号和单元电阻率，构建磁场离散的球坐标交错网格单元；

d、对单个频率进行循环，将步骤a中的大地电磁控制方程在步骤c中的球坐标交错网格单元进行数值离散，得到系数矩阵和右端项，组装到线性方程组中；

e、求解线性方程组，迭代过程中开展散度校正，计算得到各个节点的场值；

f、根据大地电磁场源，由球坐标张量阻抗公式计算阻抗，再带入卡尼亚视电阻率计算公式求取测点处的视电阻率和相位。

所述步骤a中，大地电磁控制方程采用忽略位移电流后的麦克斯韦尔积分方程形式。

所述步骤b中，划分成若干个小的倒立四棱柱网格单元是指沿θ轴方向剖分成N_θ段，每段的编号i沿θ轴方向序号递增，i＝1,2,…,N_θ，网格弧度为Δθ(i)(1,...,N_θ)；沿轴方向被剖分成段，网格弧度为沿r轴方向被剖分成N_r段，网格间距为Δr(k)(1,...,N_r)。

在球坐标系中，沿r、θ、三个坐标轴方向，分别用若干平行的球面以不同的间距划分成若干个小的倒立四棱柱网格单元；设置球坐标模型参数，包括网格节点坐标、单元、节点编号和单元电阻率，构建磁场离散的球坐标交错网格单元；对单个频率进行循环，将大地电磁控制方程在球坐标交错网格单元进行数值离散，得到系数矩阵和右端项，组装到线性方程组中；针对传统笛卡尔坐标系六面体剖分与实际大地电磁模型存在的拉伸效应，采用球坐标系中交错倒立四棱柱网格单元剖分方式，能够更好的匹配大地电磁实际模型，避免坐标拉伸投影所引入的误差，从而能够克服传统大地电磁正演方法存在的计算偏差，具有符合地球形态的三维大地电磁正演模拟的特定效果。

实施例4

一种基于球坐标系的大地电磁三维正演方法，包括以下步骤：

a、建立大地电磁控制方程；

b、在球坐标系中，沿r、θ、三个坐标轴方向，分别用若干平行的球面以不同的间距划分成若干个小的倒立四棱柱网格单元；

c、设置球坐标模型参数，包括网格节点坐标、单元、节点编号和单元电阻率，构建磁场离散的球坐标交错网格单元；

d、对单个频率进行循环，将步骤a中的大地电磁控制方程在步骤c中的球坐标交错网格单元进行数值离散，得到系数矩阵和右端项，组装到线性方程组中；

e、求解线性方程组，迭代过程中开展散度校正，计算得到各个节点的场值；

f、根据大地电磁场源，由球坐标张量阻抗公式计算阻抗，再带入卡尼亚视电阻率计算公式求取测点处的视电阻率和相位。

所述步骤a中，大地电磁控制方程采用忽略位移电流后的麦克斯韦尔积分方程形式。

所述步骤b中，划分成若干个小的倒立四棱柱网格单元是指沿θ轴方向剖分成N_θ段，每段的编号i沿θ轴方向序号递增，i＝1,2,…,N_θ，网格弧度为Δθ(i)(1,...,N_θ)；沿轴方向被剖分成段，网格弧度为沿r轴方向被剖分成N_r段，网格间距为Δr(k)(1,...,N_r)。

所述步骤f中，大地电磁场源，分解为两个正交源场等效作用的结果，表征为Sθ和两个极化模式。

较现有技术而言，针对传统笛卡尔坐标系的阻抗张量公式存在的误差问题，构建了球坐标张量阻抗公式，使得求得的阻抗更加精准；对工程应用尺度的音频大地电磁、地壳尺度的大地电磁以及上地幔尺度的长周期大地电磁均具有良好的适用效果。

本发明基于球坐标系的大地电磁三维正演方法，包括以下步骤：

1)在大地电磁研究的频率范围内，约10⁴～10^-4赫兹，由于频率相对较低，满足似稳电磁场，可以忽略位移电流；取地下介质的磁导率为真空磁导率，时谐因子为e^-iωt，此时大地电磁控制方程即为频率域电场强度矢量和磁场强度矢量满足的Maxwell方程积分形式：

其中，E为电场强度矢量，H为磁场强度矢量，μ₀为真空磁导率，μ₀＝4π×10^-7H/m，σ为介质电导率，ω为圆频率；

2)在球坐标系下沿r、θ、三个坐标轴方向，分别用若干平行的球面以不同的间距划分成若干个小的网格单元；在球坐标系中沿θ轴方向剖分成N_θ段，每段的编号i沿θ轴方向序号递增，i＝1,2,…,N_θ，网格弧度为Δθ(i)(1,...,N_θ)；沿轴方向被剖分成段，网格弧度为沿r轴方向被剖分成N_r段，网格间距为Δr(k)(1,...,N_r)；

3)对于某倒立四棱柱网格单元(i,j,k)，rΔθ(i)、rΔθ(i+1)、 Δr(k)和Δr(k+1)分别为此单元的各边长；ρ(i,j,k)为此单元的电阻率值，磁场H_θ(i,j,k)、H_r(i,j,k)分布于棱边中心，电场E_θ(i,j,k)、E_r(i,j,k)分布于侧面中心，电场和磁场分量形成交错采样分布以遵守能量守恒定律；

4)对单元某倒立四棱柱网格单元(i,j,k)，步骤1中的方程1按电流密度J的分量形式在该单元中离散化表示为：

同理，方程2在该单元中离散化表示为：

另外，电场E可用电阻率和电流密度J的离散表示为：

E_θ(i,j,k)＝ρ_θ(i,j,k)·J_θ(i,j,k) (9)

E_r(i,j,k)＝ρ_r(i,j,k)·J_r(i,j,k) (11)

再根据在各电场分量采样点处的电阻率值为相邻网格单元电阻率的网格间距加权平均值，将方程(3)～(11)联立消除电场分量后，对于研究区内部的每个H_θ(i,j,k)、和H_r(i,j,k)分量，都可以分别与其周围的12个磁场分量建立线性关系，即可满足(12)、(13)、(14)的线性方程：

其中方程(12)中各系数的表达式分别为：

a₇＝-a₆，a₉＝-a₈，

a₁₁＝-a₁₀，a₁₃＝-a₁₂；

其中方程(13)中各系数的表达式分别为：

b₃＝-b₁，b₄＝-b₂，

b₁₁＝-b₁₀，b₁₃＝-b₁₂；

其中方程(14)中各系数的表达式分别为：

c₃＝-c₁，c₄＝-c₂，c₇＝-c₅，c₈＝-c₆，

其中方程(12)中的a₀、(13)中的b₀和(14)中的c₀在模型内部单元离散时值为0，在边界单元由边界条件确立；

5)步骤4中的方程(12)～(14)可等价认为H_θ(i,j,k)可以根据周围的四个H_θ分量、四个分量及四个H_r分量求出，将所有单元的方程组装后，方程组形式为:

其中，H_θ、和H_r为求解向量，分别是离散磁场编号存储的一维向量的子向量，b_θ、和b_r分别是各子向量对应的右端项子向量； 分别为方程(12)～(14)不同单元相关的a、b、c系数组成；根据边界条件类型，将其在线性方程组中的贡献加入方程组的左边或右端项即形成定解问题；采用预条件双共轭梯度迭代法解线性方程组即可求取每个采样点上的H_r、H_θ、值；另外，在低频情况下往往无论如何迭代，残差都无法收敛到理想水平，还需要进行散度校正来校正预条件双共轭梯度迭代的中间解；

6)通过步骤5)求得各个节点的磁场后，再根据麦克斯韦方程电场和磁场的偏微分关系，即可求得电场值；对于大地电磁场源，可分解为两个正交源场等效作用的结果，在笛卡尔坐标中表征为Sx和Sy两个极化模式，在本发明中则表征为Sθ和两个极化模式，因此，球坐标张量阻抗公式为：

再由卡尼亚视电阻率计算公式，代入场值可求取地面测点处的视电阻率和相位。

具体实例1

为验证本发明算法的正确性，与一维理论模型和三维国际标准测试DTM1模型进行了对比测试。地球平均半径6371千米，因此在模型尺寸较小时，球坐标理论上应该与笛卡尔坐标正演结果一致。

首先，设计层状模型第一层电阻率100欧姆米，厚1000米，第二层电阻率1欧姆米，厚1000米，第三层电阻率100欧姆米。对于一维模型，利用理论公式计算出视电阻率和相位的解析解。其次，计算过程中球坐标和笛卡尔坐标系正演算法使用相同的网格剖分，网格数为30×30×40，背景模型是100欧姆米均匀半空间。计算了从1000赫兹到0.0001赫兹共8个频率的视电阻率和相位响应。表1是两种坐标系数值解与解析解的对比结果；球坐标解和笛卡尔坐标解完全一致，均与解析解吻合良好，视电阻率平均相对误差1.4％，相位平均误差0.23°。

表1

对于简单的一维层状模型，可以看出本发明计算结果与现有的笛卡尔坐标计算结果完全一致，但对不考虑地球曲率的复杂三维模型，还需采用复杂的三维模型进行测试。用于测试的DTM1模型为2008年国际三维大地电磁反演讨论会上提出的标准测试模型，该模型在100欧姆米的均匀半空间中设置了三个电阻率异常块，电阻率分别为10欧姆米，1欧姆米，10000欧姆米，三个异常体在空间上交错分布，电阻率差异明显，复杂性较强。对国际标准测试DTM1模型进行三维正演计算，图3示出了本发明与现有笛卡尔坐标正演结果xy和yx视电阻率对比曲线，图4示出了本发明与现有笛卡尔坐标正演结果角xy和yx相位对比曲线，计算频率从10～0.0001赫兹共21个频点，网格采用40×40×60。球坐标和笛卡尔坐标系正演结果重合度较高，均方误差小于0.1％，曲线基本重合。因此，结合一维简单层状模型和三维复杂标准测试模型计算结果分析，本发明基于球坐标系大地电磁正演算法是正确的，当模型尺度在可忽略地球曲率的前提下，球坐标系正演结果与现有笛卡尔坐标系基本相同。

具体实例2

通过具体实例1验证了本发明算法的正确性，接下来就是通过计算较大尺寸模型即地球曲率不可忽略情况下，本发明的计算结果与现有笛卡尔坐标计算结果的差异。仍然以较复杂的国际标准测试DTM1模型为基础，将国际标准测试DTM1模型进行放大扩展，模型相对位置示意图与国际标准测试DTM1模型一致，但异常体尺寸变为表2所示，表2为国际标准测试DTM1模型放大尺寸。两种坐标系使用相同的网格剖分，使用网格是40×40×60，表层单元网格横向宽度均为100千米宽的等距网格，背景模型是100欧姆米均匀半空间，计算频率从10～0.00001赫兹共25个频点。

异常体编号	x范围/千米	y范围/千米	深度范围/千米
				异常体1	-1500～1500	-100～100	5～30
异常体2	-1500～0	-100～1500	30～90
				异常体3	0～1500	-1500～100	90～180

表2

对放大后国际标准测试DTM1模型进行三维正演计算，图5示出了本发明与现有笛卡尔坐标在模型测线X＝0处、计算周期为10000秒的视电阻率对比曲线，图6示出了本发明与现有笛卡尔坐标在模型测线X＝0处、计算周期为10000秒的相位对比曲线。两种坐标系视电阻率和相位计算结果差异明显，差异较大位置集中在模型电阻率突变区域。其中xy视电阻率、xy相位、yx视电阻率、yx相位最大误差分别为5.8％、5.5％、9.8％、2.7％。计算结果表明，对于较大尺度的大地电磁应用，地球曲率影响不能忽略，现有的笛卡尔坐标正演模拟在低频段会导致较大误差。

说明采用本发明，在大地电磁三维正演模拟中，能够有效克服传统基于笛卡尔坐标系正演模拟存在的计算偏差，具有符合地球形态的三维大地电磁正演模拟的特定效果；能够与地球模型更好的匹配，精确度高。

Claims (4)

1.一种基于球坐标系的大地电磁三维正演方法，其特征在于，包括以下步骤：

a、建立大地电磁控制方程；

b、在球坐标系中，沿r、θ、三个坐标轴方向，分别用若干平行的球面以不同的间距划分成若干个小的倒立四棱柱网格单元；

c、设置球坐标模型参数，包括网格节点坐标、单元、节点编号和单元电阻率，构建磁场离散的球坐标交错网格单元；

d、对单个频率进行循环，将步骤a中的大地电磁控制方程在步骤c中的球坐标交错网格单元进行数值离散，得到系数矩阵和右端项，组装到线性方程组中；

e、求解线性方程组，迭代过程中开展散度校正，计算得到各个节点的场值；

f、根据大地电磁场源，由球坐标张量阻抗公式计算阻抗，再带入卡尼亚视电阻率计算公式求取测点处的视电阻率和相位。

2.根据权利要求1所述的一种基于球坐标系的大地电磁三维正演方法，其特征在于：所述步骤a中，大地电磁控制方程采用忽略位移电流后的麦克斯韦尔积分方程形式。

3.根据权利要求1所述的一种基于球坐标系的大地电磁三维正演方法，其特征在于：所述步骤b中，划分成若干个小的倒立四棱柱网格单元是指沿θ轴方向剖分成N_θ段，每段的编号i沿θ轴方向序号递增，i＝1,2,…,N_θ，网格弧度为Δθ(i)(1,...,N_θ)；沿轴方向被剖分成段，网格弧度为沿r轴方向被剖分成N_r段，网格间距为Δr(k)(1,...,N_r)。

4.根据权利要求1所述的一种基于球坐标系的大地电磁三维正演方法，其特征在于：所述步骤f中，大地电磁场源，分解为两个正交源场等效作用的结果，表征为Sθ和两个极化模式。