一种针对RMT数据的双参数反演算法的制作方法

本发明涉及一种应用地球物理领域的高频电磁数据反演算法，特别涉及射频大地电磁法数据的反演。

背景技术：

目前，各类电磁资料的反演均是以寻求满足目标拟合差的地下介质电阻率分布为目的，忽略了地下介电常数差异对观测数据的影响。然而，实际中观测数据是由地下电阻率和介电常数共同决定的，因此采用传统的单参数电阻率反演来拟合观测数据是存在一定偏差的。尤其对于射频大地电磁法(radiomagnetotelluric,rmt)数据，由介电常数所引起的波动场在总场中的占比可达20％以上，这种情况下，仅通过电阻率参数的反演来拟合目标函数势必引起一定偏差，降低资料解释的准确性。

因此，有必要开发一种综合考虑地下电阻率和介电常数的双参数同步反演算法。

技术实现要素：

本发明所解决的技术问题是，针对现有技术的不足，提供一种双参数同步反演算法，本发明能够同时得到地下电阻率和介电常数分布，降低了数据拟合偏差，提高了反演准确性。

本发明的技术方案为：

一种针对rmt数据的双参数反演算法，包括以下步骤：

步骤1、读取测点的位置、测量频率以及各个测点的观测数据；各个测点的观测数据构成n维观测数据向量d＝(d1,d2,…,dn)^t，其中di为第i个测点的观测数据，n为测点的个数；

步骤2、为待反演的模型的两个参数m1(k)和m2(k)赋初值(一般取初始模型为均匀半空间模型，那么所有反演网格上的初始参数为同一常数，m1(k)和m2(k)的初值为常数向量)；k为迭代次数，其初始值为1；

步骤3、根据测点的位置和测量频率，计算由m1(k)和m2(k)共同作用得到的正演响应f(m1(k),m2(k))(计算正演响应为本领域公知技术)，f(m1(k),m2(k))为n维向量，f(m1(k),m2(k))＝(f1,f2,…,fn)^t，其中fi为第i个测点的正演响应值；

步骤4、计算数据拟合差rmsk，计算公式为

步骤5、判断rmsk是否大于目标拟合差(目标拟合差一般取为1)，若是，则进行第k次双参数反演迭代；否则退出迭代，反演结束；

进行第k次双参数反演迭代的步骤为：

(a)计算第k次迭代时的正演响应对两个模型参数的灵敏度矩阵j1(k)和j2(k)：

(b)求解如下双参数反演迭代方程组，得到第k次迭代的两个模型参数的修正量δm1(k)和δm2(k)；

其中，为观测数据向量d和正演响应f(m1(k),m2(k))的数据拟合差权重，其中diag(·)表示对角矩阵，即除了主对角线外的元素均为零的方阵，矩阵主对角线元素取值为1/εi，i＝1,2…,n，εi表示第i个测点的正演响应值fi与第i个测点的观测数据di之间的协方差，即和分别为与m1(k)和m2(k)相关的光滑度矩阵，其取值为相邻反演网格参数和的标准差，即其中x＝1,2；和的上标1/2表示矩阵的开方，即对矩阵中的所有元素分别开方；m10和m20分别为两个模型参数的先验信息，一般根据已知地质资料给出，若无，则通常令m10和m20为零；λ和γ为调节两个模型参数拟合差与数据拟合差之间权重的正则化因子；

(c)根据下式得到新的模型参数m1(k+1)和m2(k+1)：

m1(k+1)＝m1(k)+δm1(k)

m2(k+1)＝m2(k)+δm2(k)；

步骤6、令k＝k+1，返回步骤3。

进一步地，步骤1中，观测数据为视电阻率(ωm)或相位(°)；

进一步地，设置两个模型参数m1(k)和m2(k)分别为电阻率和介电常数。

进一步地，设置两个模型参数m1(k)和m2(k)分别为相对介电常数相对电导率；其中相对介电常数定义为εr＝ε/ε0，相对电导率定义为σr＝σ/σ0，ε和σ分别为均匀地下介质的介电常数和电导率，ε0为真空中的介电常数，σ0为参考场的电导率，σ0＝ω0ε0，ω0为参考频率。

进一步地，设置ω0的值，使j2(k)的数值略大于j1(k)，两者差值根据经验信息确定。

进一步地，正则化因子λ取j2(k)的最大特征值，γ取j1(k)的最大特征值。

本发明的原理为：

(一)构建反演目标函数时，同时考虑地下电阻率和介电常数两个参数，构建双参数目标函数如下所示：

其中，目标函数φ(m1,m2)是与介电常数m1和电阻率m2均相关的函数；m1和m2为m^inv维的待反演的模型的两个参数；d为n维观测数据向量，d＝(d1,d2,…,dn)^t，其中di为第i个测点的观测数据，n为测点的个数；观测数据一般为视电阻率(ωm)和相位(°)；f(m1,m2)是由m1和m2共同作用得到的正演响应，f(m1,m2)＝(f1,f1,…,fn)^t；为观测数据向量d和正演响应f(m1,m2)的数据拟合差权重，其中diag(·)表示对角矩阵，即除了主对角线外的元素均为零的方阵，矩阵主对角线元素取值为1/εi，i＝1,2…,n，εi表示第i个正演响应值fi与第i个观测数据di之间的协方差，即和分别为与m1和m2相关的光滑度矩阵，其取值为相邻反演网格参数mxi和mxj的标准差，即其中x＝1,2；和的上标1/2表示矩阵的开方，即对矩阵中的所有元素分别开方；m10和m20分别为两个模型参数的先验信息，λ和γ为调节两个模型参数拟合差与数据拟合差之间权重的正则化因子；

f(m1(k),m2(k))表示第k次时两个模型参数的m1(k)和m2(k)共同作用得到的正演响应；对正演响应f(m1,m2)进行二元taylor展开，表达式如下：

f(m1(k+1),m2(k+1))＝f(m1(k),m2(k))+j1(k)(m1(k+1)-m1(k))+j2(k)(m2(k+1)-m2(k)),(2)

其中，mx(k)和mx(k+1)分别表示第k次和第k+1次迭代获得的模型参数，其中x＝1,2；j1(k)和j2(k)分别表示第k次迭代时的正演响应对两个模型参数的灵敏度矩阵，即正演响应对两个模型参数一阶偏导数：本发明中采用互换原理进行求解(siripunvaraporn,2012)；

将双参数目标函数φ(m1,m2)分别对m1和m2求导数，公式如下：

令(3a)和(3b)分别为零，得到如下等式，

联立(4a)和(4b)后可得如下双参数反演迭代方程组，

其中，j1(k)和j2(k)分别表示第k次迭代时的正演响应对两个模型参数的灵敏度矩阵，式中和表示j1(k)和j2(k)的转置矩阵；

通过求解方程(5)即可得到第k次迭代的两个模型参数的修正量δm1(k)和δm2(k)，根据下式得到新的模型参数m1(k+1)和m2(k+1)：

m1(k+1)＝m1(k)+δm1(k),(6)

m2(k+1)＝m2(k)+δm2(k).(7)。

进行反演迭代得到两个模型参数的分布结果。

(二)在进行多参数同步反演前，首先要研究不同地电参数对正演响应的灵敏度，然后根据参数灵敏度的不同对电导率σ和介电常数ε进行合理的尺度变换以保证反演的稳定性。

这里首先以均匀半空间模型为例来讨论电导率和介电常数的改变对电场场值的影响。假设均匀半空间电导率为σ^b，介电常数为ε^b，那么地下某处的电场ex可表示为：

ex＝ex0e^iωte^-kz,(8)

其中，ex0为地表电磁场值；e为数学中的欧拉数；ω为角频率；t为时间；k为波数，其计算式为k²＝ω²με-iωμσ，式中μ、σ、ε分别为均匀地下介质的磁导率、电导率和介电常数；z为电磁场的传播距离。

将ex分别对电导率σ和介电常数ε求导可得：

其中，

根据差分公式，电导率改变量δσ^b和介电常数改变量δε^b引起的电场场值改变量和分别为：

将(9)(10)式代入(11)式，并将(11a)除以(11b)得，

通常，令电导率和介电常数改变量均为原来的n倍，即δσ^b＝σ^b·n、δε^b＝ε^b·n，那么(12)式可转换为：

从(13)式可看出，在mt的勘探频率下σ^b>>ωε^b，即也就是说此时观测数据对电导率更为灵敏；类似地，在gpr的勘探频率下σ^b<<ωε^b，即这就意味着gpr观测数据对介电常数更为灵敏(尤其是在高阻地区)；而在rmt勘探频段，σ^b≈ωε^b，此时，rmt观测数据对电导率和介电常数的灵敏度相当。然而，考虑到实际中电导率的变化范围更广(10^-4～0.1s/m)，而介电常数的变化范围为1～81，这使得实际反演中，电导率更容易引起观测数据的改变，而削弱了介电常数的影响。此外，从(12)式中我们可看出，观测数据对电导率的灵敏度与频率的一次方成正比，而观测数据对介电常数的灵敏度与频率的二次方成正比，因此，随着频率的升高，介电常数对观测数据的影响更大。为此，我们期望寻求一种电导率变换式，其数据尺度缩小，同时保证频率对电导率灵敏度的影响与对介电常数灵敏度的影响相当。

定义复介电常数ε'为：

(14)式的实部代表位移电流的贡献，它在电磁波的传播过程中不会产生能量损耗；而虚部代表传导电流的贡献，由于传导电流的传播过程中会产生焦耳热，从而使能量有所损耗。根据物理学中的定义，介质损耗角正切tanδ0是表征电磁波在介质中传播时将电磁能量转化为热能所消耗的能量，其数学表达式为：

令(15)式中的tanδ0＝1，此时，位移电流与传导电流对场的贡献相等，定义tanδ0＝1时的地电参数为参考场。σ0为参考场的电导率，ω0为参考频率，ε0为真空中的介电常数。由此可得参考场的电导率σ0＝ω0ε0。与相对介电常数εr＝ε/ε0的定义类似，定义一个相对电导率参数σr＝σ/σ0，然后在双参数反演中反演相对电导率和相对介电常数。经过这一变换后，观测数据对σr和εr的灵敏度的频率依赖性相当，并且缩小了电导率参数的变化范围，提高了反演稳定性。

从相对电导率的定义可看出，相对电导率的大小随着参考频率的变化而改变，也就是说，参考频率直接影响观测数据对相对电导率的灵敏度的大小。为此，有必要研究参考频率的选取对反演过程的影响。为简单起见，在研究参考频率时令λ和γ均为零，即忽略模型拟合项，同时令等于单位矩阵i。经过上述简化后，(5)式左端的系数矩阵可表示为：

其中，系数矩阵的维数为2m^inv×2m^inv，m^inv为反演网格数，j1(k)和j2(k)分别为第k次迭代时的正演响应对相对介电常数εr和相对电导率σr的灵敏度矩阵。由于参考频率的大小只影响j2(k)，因此(16)式中j1(k)与参考频率无关。理论而言，合理的参考频率ω0应保证j2(k)的数值与j1(k)相当，即使观测数据对相对电导率σr和相对介电常数εr的灵敏度相当，而考虑到实际中相对电导率参数对rmt观测数据影响比相对介电常数更显著，因此，在rmt实测数据反演中应使j2(k)的数值略大于j1(k)。具体分析见图4。

有益效果：

本发明提供了一种新的双参数同步反演算法。建立了关于双参数的二元目标函数，并进行二元极小化得到由两个参数耦合的反演方程组。

本发明提供的双参数反演算法能够同时得到地下电阻率和介电常数分布，丰富了反演解释资料，提高了现有的资料处理与解释水平。本算法主要应用于rmt数据反演及资料解释。

与传统的电阻率单参数反演相比，本发明中正演计算和反演迭代过程考虑了地下介电常数的影响，降低了数据拟合偏差，从理论上更符合实际情况，提高了反演准确性；在不增加数据采集工作量的情况下，仅通过反演计算即可得到地下电阻率和介电常数两个电性参数分布，提供了更丰富的地电信息；采用本算法进行rmt数据反演，可有效压制传统准静态近似下反演所带来的浅地表高祖冗余构造，不仅增强了rmt数据反演理论水平，而且提高了rmt实际资料处理解释准确性。

使用本发明，通过反演计算得到地下电阻率和介电常数的电性分布，可以查明地下地电特征、地质构造及矿产分布或者解决其它浅地表工程、水文及环境地质问题。

附图说明

图1为本发明原理图。

图2为均匀半空间中赋存一矩形异常体的模型示意图。其中背景电阻率为1000ωm，相对介电常数σr0为5，矩形异常体电阻率为1000ωm，相对介电常数σr1为1。

图3为图1中模型的te模式下视电阻率及相位灵敏度。其中观测数据点位于(0,0)处，计算频率为250khz。图3(a)～图3(c)分别为视电阻率对对数电阻率、相对电导率及相对介电常数的灵敏度，图3(d)～图3(f)分别为相位对三个参数的灵敏度。

图4为不同参考频率的选取对反演系数矩阵数值分布的影响。各子图的横纵坐标取系数矩阵的角标，图中符号s1、s2、s3、s4分别表示系数矩阵的左上部分j1k、右上部分各个j2k、左下部分j1k和右下部分j2k；图4(a)～图4(c)分别为参考频率ω0＝10khz、50.1khz和100khz时反演系数矩阵数值分布；

图5为不同参考频率下σr和εr同步反演结果。三个参考频率分别为10khz、50.1khz和100khz；图5(a)～图5(c)分别为三个参考频率下反演得到的电导率图像，图5(d)为100khz下同步反演得到的相对介电常数分布图。图中下三角为测点位置，虚线框为异常体所在位置。

图6为只考虑介电常数模型约束(即令γ为零)的反演迭代收敛曲线。图中λ分别取5，0.5，0.05和0.005。

图7为λ＝0.5,γ＝0.01时双参数εr和σr同步反演结果，图7(a)和图7(b)分别为εr和σr同步反演结果。

具体实施方式

以下结合附图和具体实施方式对本发明作进一步的说明。

本发明的双参数反演算法的实施包括以下步骤：

1、读取观测数据文件(视电阻率和相位)、测量频率文件、地形文件及测点位置信息等；其中观测数据d是由仪器测得，用于反演计算；测点是仪器的布设位置，频率是从仪器中所设置，测点和频率用于计算正演响应f(m1k,m2k)；地形是观测区域的地表信息，用于进行反演网格离散；

2、根据测点位置信息、目标深度以及观测区域的地表信息确定反演模型的规模；基于开源代码triangle，对待反演的模型进行非结构网格(三角形网格)剖分，生成反演网格，然后通过限定每个反演单元的最大面积对反演网格进行加密生成正演网格(shewchuk,1996)；

3、为待反演的模型的两个参数m1k和m2k赋初值(一般取初始模型为均匀半空间模型，那么所有反演网格上的初始参数为同一常数)；为了计算正演响应，需要将反演网格参数映射到正演网格上才能进行正演计算，因此采用遍历法将反演网格参数映射到正演网格上，得到正演网格参数；

采用遍历法将反演网格参数映射到正演网格上的方法为：设正演网格单元数为m^fwd，反演网格单元数为m^inv，针对每一个正演网格对所有反演网格进行搜索，当第i个正演网格的三角形中心点与某一个反演网格的三个节点相连得到的三个角度之和为360°，那么这个正演网格就位于该反演网格内，反之则位于该反演网格之外；若确定第i个正演网格位于第j个反演网格内，则将第j个反演网格的电性参数(包括电阻率及介电常数)赋值给第i个正演网格(正演网格参数和反演网格参数相同，均为电阻率和介电常数)针对所有正演网格进行上述搜索和赋值后，完成网格参数映射；

得到正演网格参数后，计算正演响应；然后计算初始数据拟合差rms，计算公式为

4、进行第k次双参数反演迭代。

(a)计算灵敏度矩阵j1k和j2k；

(b)组建反演方程组；

(c)求解方程组，得修正后的模型m1k+1、m2k+1；

(d)计算rms，若rms大于目标拟合差(目标拟合差一般取为1)，则继续执行下一次反演迭代，否则退出迭代，反演结束；

5、双参数地电模型可视化及资料解释。

图2为均匀半空间中赋存一矩形异常体的模型示意图。该模型用于研究视电阻率和相位对电导率及介电常数参数的灵敏度。同时作为理论模型测试本发明所提出的双参数同步反演算法的有效性。

图3(a)～图3(f)为te模式下计算得到的视电阻率及相位对不同反演参数的灵敏度。观测数据点位于(0,0)处，计算频率为250khz，计算耗时12.76s。图中，图3(a)～图3(c)为视电阻率对以10为底的对数电阻率log10(ρ)、相对电导率σr及相对介电常数εr的灵敏度；图3(d)～图3(f)为相位对以10为底的对数电阻率log10(ρ)、相对电导率σr及相对介电常数εr的灵敏度。从图3(a)和图3(d)两图中可看出，观测数据对log10(ρ)的灵敏度在10的-2～-3次方数量级上，而图3(c)和图3(f)图中，观测数据对εr的灵敏度较小，为10的-4～-5次方数量级上，在这种情况下，由于观测数据对log10(ρ)的灵敏度远大于εr，也就是说εr对观测数据的影响微乎其微，导致反演迭代中对相对介电常数的修正几乎为零。图3(b)和图3(d)为本发明中所定义的相对电导率σr对电阻率参数进行变换后的灵敏度，参数变换中参考频率ω0取100khz。从图中可看出，经过尺度变换后，视电阻率和相位对σr的灵敏度也处于10的-4～-5次方数量级上，从而保证σr和εr对观测数据有相当的贡献，为双参数的反演奠定基础。

图4(a)～4(c)为三个不同参考频率下计算得到的系数矩阵。系数矩阵的维数为2m^inv×2m^inv，子图的横纵坐标取系数矩阵的角标，图中符号s1、s2、s3、s4分别表示系数矩阵的左上部分j1(k)、右上部分各个j2(k)、左下部分j1(k)和右下部分j2(k)。图4(a)中参考频率ω0＝10khz，由于参考频率较小，使得系数矩阵中s4的数值小于s1，也就是说εr对观测数据的贡献更大；图4(b)中参考频率ω0＝50.1khz，在这个参考频率下，系数矩阵中s4的数值和s1相当，也就是说此时εr和σr对观测数据的贡献差不多；而图4(c)中，当参考频率ω0＝100khz时，系数矩阵中s4的数值要大于s1，这意味着此时σr对观测数据的贡献更大。合理的参考频率应保证j2(k)的数值与j1(k)相当，即使双参数εr和σr对观测数据的贡献相当。具体分析见图5。

图5(a)～图5(d)为不同参考频率下的双参数同步反演结果。三个参考频率分别取10khz、50.1khz和100khz，反演迭代次数均为15次。反演的初始模型取背景电阻率10000ωm和相对介电常数5。图5(a)～图5(c)是参考频率分别取10khz、50.1khz和100khz时反演得到的电导率分布图，从图中可看出，当参考频率ω0为10khz时，反演结果无法得到异常电导率分布信息(图5(a))，这是因为参考频率较小时会降低相对电导率的灵敏度，使得每次反演迭代中电导率的改变非常微弱，进而无法得到异常信息；当ω0增加到50.1khz时，电导率异常有所显示，但是异常仍不太明显(图5(b))；而当参考频率增至100khz时，电导率异常分布位置和大小均能准确的反演出来(图5(c))，这是因为较大的参考频率增加了相对电导率的灵敏度。图5(d)为参考频率为100khz时反演得到的相对介电常数分布，这里仅给出100khz下的相对介电常数反演结果是因为介电常数灵敏度大小不受参考频率的影响，不同参考频率反演得到的介电常数分布类似。

综上分析，可得到以下结论：在双参数同步反演时，ω0的选择至少要保证σr对观测数据的贡献大于或等于εr才能得到正确的反演结果。同时，考虑到rmt数据对电导率参数的依赖程度更高，因此在实际中取略大的ω0，使σr的贡献略大于εr能够得到更为合理的反演结果。

图6和图7用于研究权重因子λ和γ的选择对反演结果的影响。从图5(d)可看出，当令λ和γ均为零时，得到的相对介电常数模型左右边界分辨率较差，为此，我们对目标函数增加介电常数约束。图6为只考虑介电常数模型约束的反演迭代收敛曲线。从图中可看出，相对介电常数模型约束的权重因子λ的取值既不能过大，使得收敛速度非常缓慢，在合理的迭代步中无法突显出异常体，同时，λ的取值也不能太小，虽然当λ取值很小时会得到更小的数据拟合差，但是，过小的λ使得反演结果不可靠，甚至是错误的。图7(a)和图7(b)为λ取0.5，γ取0.01时εr和σr同步反演的结果。反演迭代次数均为15次。从图中可看出，γ取0.01时电导率模型和介电常数模型能够准确反映出异常体的位置，同时数据拟合差也逐步收敛。此外，我们也测试了γ取0.1，0.001和0.0001的反演结果，当γ取0.1时，电导率模型在异常体上方出现了高阻异常，且数据拟合差无法稳定收敛；而当γ取0.001和0.0001时得到的模型分布与图7类似。

根据前面的对比分析及测试，在进行εr和σr同时反演时，一般情况下选取λ为j2的最大特征值，γ为j1的最大特征值能够得到较为合理的结果。

参考文献：

[1]siripunvaraporn,w.three-dimensionalmagnetotelluricinversion:anintroductoryguidefordevelopersandusers[j].surveysingeophysics,2012,33(1):5-27.

[2]shewchukj.triangle:engineeringa2dqualitymeshgeneratoranddelaunaytriangulator[j].lecturenotesincomputerscience,1996,1148:203-222.