一种基于偏差原理的大地电磁正则化反演方法及处理终端与流程

1.本发明涉及地球科学中的地球物理探测技术领域，具体涉及一种基于偏差原理的大地电磁正则化反演方法及处理终端。


背景技术：

2.大地电磁反演问题属于不适定问题，通常引入正则化反演方法来进行求解，通过加入先验的约束条件来增强反演过程的稳定性质，减少反演结果的非唯一性。常用的正则化反演算法有occam反演、rebocc反演、nlcg反演等，occam反演能引入先验信息来去掉对已知构造边界处的光滑度或者加上对同类电性单元的电阻率差的限制，具有其他反演算法不具有的优点。但目前基于普通occam反演算法进行大地电磁反演，仍然存在以下不足，主要包括仍然需要依赖反演这的经验和未考虑收敛速度的问题。例如，吴小平等认为，可以在主要拟合观测数据的同时得到一个光滑而不是最光滑的地电模型，提出了正则化因子递减方案，具体做法是给一个初始正则化因子和步长进行每次迭代，虽避免了线性搜索正则化因子，减少了冗余计算，但初始正则化因子和步长的选取仍然存有一定的随意性。若想得到一个满意结果，仍然需要依赖反演者的经验，未能摆脱经验选择方式。张罗磊等提出，初始正则化因子每次迭代随模型粗糙度的变化，在总体粗糙度没用减少时，以粗糙度变化率减少正则化因子。这种方法运算速度较快，但容易陷入局部极小的缺点。朴英哲等在搜索正则化因子的过程中加上了观测数据拟合差与其期望值的比较，能够排除多余的正演计算，使得反演速度加快，但未考虑收敛速度问题。
3.相关的参考文献如下：
4.吴小平，徐果明.大地电磁数据的occam反演改进[j].地球物理学报,1998,41(4):547
‑
554；
[0005]
张罗磊，于鹏，王家林，等.光滑模型与尖锐边界结合的mt二维反演方法[j].地球物理学报,2009,52(6):1625
‑
1632；
[0006]
朴英哲,李桐林,刘永亮.在大地电磁二维occam反演中求取拉格朗日乘子方法改进[j].吉林大学学报:地球科学版,2014,44(2):660
‑
667。


技术实现要素：

[0007]
针对现有技术的不足，本发明的目的之一提供一种基于偏差原理的大地电磁正则化反演方法，其能够解决大地电磁反演的问题；
[0008]
本发明的目的之二提供一种处理终端，其能够解决大地电磁反演的问题。
[0009]
实现本发明的目的之一的技术方案为：一种基于偏差原理的大地电磁正则化反演方法，包括如下步骤：
[0010]
步骤1：获得大地电磁观测数据，并构建反演算子方程；
[0011]
步骤2：利用偏差原理求解出当前的反演算子方程中的正则化因子，得到新的反演算子方程；
[0012]
步骤3：求解出大地电磁观测数据中当前的灵敏度，
[0013]
步骤4：利用所述当前的灵敏度和新的反演算子方程对大地电磁观测数据进行反演，得到反演结果；
[0014]
步骤5：重复步骤2
‑
4，直至最后一次反演，最后一次反演的结果作为最终的反演结果。
[0015]
进一步地，所述正则化因子按公式
①
得到：
[0016][0017]
其中，α
k
表示第k次反演的正则化因子。
[0018]
进一步地，所述步骤2中，采用morozov偏差原理求解出正则化因子。
[0019]
实现本发明的目的之二的技术方案为：一种处理终端，其包括：
[0020]
存储器，用于存储程序指令；
[0021]
处理器，用于运行所述程序指令，以执行所述基于偏差原理的大地电磁正则化反演方法的步骤。
[0022]
本发明的有益效果为：本发明基于偏差原理来计算出正则化因子，进而得到反演结果，使得不依赖于经验值，反演效果更好，并且通过收敛速度更快的偏差原理高阶收敛算法计算得到灵敏度，使得收敛速度得到快速提高，也即提高了反演速度。
附图说明
[0023]
图1为较佳实施例的流程示意图；
[0024]
图2是在低阻模型的示意图；
[0025]
图3是在低阻模型下按现有的tm模型进行反演得到正则化因子的示意图；
[0026]
图4是在低阻模型下按本发明进行反演得到正则化因子的示意图；
[0027]
图5是在低阻模型下按tm模型和本发明进行反演得到均方根误差曲线图；
[0028]
图6是在低阻模型的示意图；
[0029]
图7是在低阻模型下按现有的te模型进行反演得到正则化因子的示意图；
[0030]
图8是在低阻模型下按本发明进行反演得到正则化因子的示意图；
[0031]
图9是在低阻模型下按te模型和本发明进行反演得到均方根误差曲线图；
[0032]
图10是在高低阻模型的示意图；
[0033]
图11是在高低阻模型下按现有的tm模型进行反演得到正则化因子的示意图；
[0034]
图12是在高低阻模型下按本发明进行反演得到正则化因子的示意图；
[0035]
图13是在高低阻模型下按tm模型和本发明进行反演得到均方根误差曲线图；
[0036]
图14是在高低阻模型的示意图；
[0037]
图15是在高低阻模型下按现有的te模型进行反演得到正则化因子的示意图；
[0038]
图16是在高低阻模型下按本发明进行反演得到正则化因子的示意图；
[0039]
图17是在高低阻模型下按te模型和本发明进行反演得到均方根误差曲线图；
[0040]
图18为处理终端的示意图。
具体实施方式
[0041]
下面，结合附图以及具体实施方案，对本发明做进一步描述。
[0042]
如图1
‑
图17所示，一种基于偏差原理的大地电磁正则化反演方法，包括如下步骤：
[0043]
步骤1：获得大地电磁观测数据，并构建反演算子方程，以对大地电磁观测数据进行反演。
[0044]
大地电磁反演算子方程如公式
①
：
[0045]
f(m)＝d
‑‑‑‑‑‑①
[0046]
式中，d表示大地电测观测数据，为向量形式，向量中的值为一系列大地电磁实际观测得到的数据。m表示电阻率模型参数，为向量形式，f表示大地电磁观测数据正演响应算子。通常，野外测量得到的大地电磁观测数据带有一定的误差，即有d
δ
＝d+δd，d
δ
表示带有误差后的大地电磁观测数据，δ表示误差系数，通常为常数。因此，公式
①
可以变成公式
②
：
[0047]
d
δ
＝f
α
(m)
‑‑‑‑‑‑②
[0048]
式中，f
α
(m)表示对大地电磁观测数据d
δ
正演响应算子。只需要对公式
②
进行求解，即可完成对大地电磁数据的反演。公式
②
的解为m
α
，
[0049]
为了求解公式
②
的m
α
，现有方法中通常会通过正则化进行求解，其中occam反演是一种较常用的正则化反演方法，可以用于求解公式
②
。在occam反演中，occam反演的展平泛函φ(m)如公式
③
：
[0050][0051]
其中，表示拟合差的期望值，w表示数据归一化对角矩阵，r(m)表示地下电阻率模型粗糙度，α表示正则化因子，也可称之为拉格朗日乘子。公式
③
为现有技术，在此不作过多的介绍。
[0052]
因此，只需要求解出公式
③
中的正则化因子α，即可求解公式
①
，从而得到大地电磁的反演结果。目前传统的方法求解公式
③
中的正则化因子α，大多数是凭借经验直接预设一个固定值。这导致存在一个明显的缺陷和不足：需要依赖反演者的经验，未能摆脱经验选择方式。一旦反演者的经验不足，导致反演结果出现很大的偏差，反演结果质量差。
[0053]
步骤2：利用偏差原理求解出正则化因子α，本实施例优选morozov偏差原理求解出正则化因子。偏差原理要求正则化因子α的选取与原始的观测数据误差匹配，也即满足如公式
④
的目标函数
[0054][0055]
式中，d
δ
表示大地电磁观测数据，z
α
表示正则化解，a表示大地电磁正演算子。公式
④
可通过具有二阶收敛速度的newton求解，即对公式
④
进行一阶泰勒展开，得到公式
⑤
：
[0056][0057]
令l(α)＝0，便得到newton迭代公式
⑥
：
[0058][0059]
为了进一步提高收敛效率，对上述目标函数在α＝α
k
展开并在(α
‑
α
k
)3截断，则
得到公式
⑦
：
[0060][0061]
式中，ξ
k
表示介于α和α
k
之间的常数。公式
⑦
也即是二阶公式，联合公式
⑥
和
⑦
得到三阶收敛的迭代公式，也即得到公式
⑧
：
[0062][0063]
其中，时，否则，否则，已由公式
⑤
给出，要求解出公式
⑧
，还需要求解出大地电磁反演目标函数的一阶导数和二阶导数其中，根据公式
④
，可得出，
[0064][0065][0066]
其中，是大地电磁反演过程中第k次迭代的正则化解(即反演得到的模型解)，和则分别是正则化解对应的一阶和二阶导数，它们在第k步解求方程组
⑨
得到：
[0067][0068]
其中ω
′
是稳定泛函，起着恢复解的稳定性的作用；f是公式
②
中的正演算子，f
*
是其伴随算子。方程组中的三个方程的系数矩阵是一样的，因此，只需要对系数矩阵做一次分解和回代后，便可以得到解和回代后，便可以得到和进而计算得到和从而可以求解出正则化因子α。
[0069]
当然，以上只是本实施例给出求解基于偏差原理下的正则化因子，具体求解过程也可以通过其他数学手段进行求解得到，本质上属于可通过现有数学方法进行求解得到，因此，也不进行赘述。
[0070]
所述α
k
表示第k次反演的正则化因子，α
k+1
表示第k+1次反演的正则化因子，也即每一次反演对应有一个正则化因子。因此，通过公式
⑧
，只需要给出初始值(也即给出α1的值)，即可求解出任意一次反演的正则化因子，得到每次反演均更新得到的正则化因子，而不依赖于人为给出的经验值。
[0071]
步骤3：在大地电磁观测数据反演中，灵敏度是一个重要的参数，灵敏度的计算最终可归结为电磁场对电导率的求导，在反演中，灵敏度有多个，并构成一个灵敏度矩阵，因此，灵敏度的计算也即是灵敏度矩阵的计算过程。通过基于偏差原理的高阶收敛算法，可以提高大地电磁观测数据反演过程的收敛速度，最终得到反演结果。灵敏度具体计算过程属于现有技术，在此不进行具体求解过程的介绍。
[0072]
本实施例中，在每一次将正则化因子代入大地电磁观测数据构建的反演算子方程
中进行反演过程中，均重新计算灵敏度，然后将灵敏度和对应的正则化因子代入到反演算子方程中进行反演，直至最后一次反演(也即最后一个k值)，并将最后一次反演的结果作为最终的反演结果，也即得到大地的电性结构，至于反演次数可以通过人为预设，例如预设反演30次，也即α
k
中的k的最大值可以预设得到。也即，每一次反演，均会计算一次正则化因子和灵敏度并代入反演中，得到反演结果。
[0073]
本发明基于偏差原理来计算出正则化因子，进而得到反演结果，使得不依赖于经验值，反演效果更好，并且通过收敛速度更快的偏差原理高阶收敛算法计算得到灵敏度，使得收敛速度得到快速提高，也即提高了反演速度。
[0074]
将背景技术中由吴小平提出的“大地电磁数据的occam反演改进”进行计算得到的正则化因子记为按方案一得到的正则化因子，按本发明基于偏差原理进行高阶收敛算法计算得到的正则化因子记为按方案二得到的正则化因子。
[0075]
图2
‑
图9是在低阻模型下分别按现有的te和tm模型进行反演得到的反演结果图。图2和图6部分是低阻模型示意图，图3和图7部分是按方案一得到的正则化因子，图4和图8部分是按方案二得到的正则化因子，图5和图9部分是均方根误差曲线图。其中，研究区域中为电阻率100ω
·
m的均匀半空间，存在电阻率为10ω
·
m的低阻异常体，长宽为2.4km x 1km。沿测线以300m的间距布置30个测点，正演频率以对数等间隔在0.01hz
‑
100hz的频段中取25个频率，反演的初始模型取均匀半空间。将两种正则化因子的选取方案应用到te、tm两种反演模式中。其中方案一的正则化因子的初始值为50，步长为2。
[0076]
从图中可以看出，tm和te反演结果均能显示异常体的顶部、底部埋深基本与理论模型一致。反演精度都较高。单独的模式下，将二种搜索正则化因子的方案加以对比，反演结果显示，方案二的反演效果最佳，精度最高。
[0077]
从均方根误差曲线来看，两种反演模式下，方案一迭代初期下降较慢，方案二迭代初期下降迅速，且方案二比方案一收敛更快。tm模式下初始均方根误差均为0.4716，迭代8次后，方案一的均方根误差为0.0295,方案二的均方根误差为0.0278。te模式下初始均方根误差均为0.4387，迭代8次后，方案一的均方根误差为0.1817，方案的二均方根误差为0.1417。这表明了方案二收敛迅速的优点，体现了算法的有效性和稳健性。
[0078]
图2
‑
图9均是低阻模型下的比较，下面给出高低阻模型的比较，如图10
‑
图17。图10和图14高低阻模型示意图，图11和图15是按方案一得到的正则化因子，图12和图16部分是按方案二得到的正则化因子，图13和图17部分是均方根误差曲线图。其中，电阻率100ω
·
m的均匀半空间中，存在高低组两个异常体，左低阻体长宽为1.5km x 1km，电阻率为10ω
·
m。右高阻体长宽为1.5km1km，电阻率为1000ω
·
m。频率、测点数目和距离和低阻模型的一致。反演初始模型为均匀半空间，将二种正则化因子的选取方案应用到te、tm两种反演模式中。其中方案一的初始正则化因子取50，步长为2。
[0079]
tm和te两种模式的反演结果分别如图4和图5所示。从图中可以看出te和tm两种反演模式基本都能反应出高低阻体的位置，但两种模式又有区别，te模式对高低阻体的反映在纵向分辨率上比tm模式的低。将二种搜索正则化因子的方案加以对比，反演结果显示，方案二的反演效果最佳，精度最高。方案一迭代初期下降较慢，方案二迭代初期下降迅速，且方案二比方案一收敛更快。tm模式下初始均方根误差均为0.4940，迭代10次后，方案一均方根误差为0.0213,方案二均方根误差为0.0188。te模式下初始均方根误差均为0.4784，迭代
10次后，方案一的均方根误差为0.1226，方案二的均方根误差为0.1001。表明了方案二收敛迅速的优点，体现了算法的有效性和稳健性。
[0080]
通过模型试算，均体现了本发明(也即方案二)的高效收敛。若将此算法进一步推广到三维大地电磁反演中，因此算法快速收敛的优点，有望进一步提高反演效率和缩短反演耗时。
[0081]
如图18所示，本发明还涉及一种处理终端100，其包括：
[0082]
存储器101，用于存储程序指令；
[0083]
处理器102，用于运行所述程序指令，以执行所述基于偏差原理的大地电磁正则化反演方法的步骤。
[0084]
本说明书所公开的实施例只是对本发明单方面特征的一个例证，本发明的保护范围不限于此实施例，其他任何功能等效的实施例均落入本发明的保护范围内。对于本领域的技术人员来说，可根据以上描述的技术方案以及构思，做出其它各种相应的改变以及变形，而所有的这些改变以及变形都应该属于本发明权利要求的保护范围之内。