1.本发明涉及时频电磁勘探数据反演技术领域，尤其是一种时频电磁频率域数据自约束稳健电阻率反演方法。
背景技术：
2.时频电磁方法(tfem)是一种地球物理勘探的新方法，广泛用于油气勘探、矿产勘探以及环境监测等。基于电磁感应原理，该方法将频率域测深与时间域测深统一在一个系统中。由于时频电磁收发距较近，接地电性源(一般在10公里左右)无法等效为电偶极子，激发的电磁场为非平面波。因此，在定义电阻率时，需要考虑近场效应和远场效应，这是一个繁琐的过程。全区视电阻率的计算是目前较普遍的解决办法，但是全区视电阻率对地下异常介质分辨率依然较低，参与反演时计算速度较慢。虽然原始电磁场数据包含丰富的地下异常信息，但是由于电磁场衰减较快，数量级变化较大，如果使用电磁场幅值直接成像，其对地下异常介质的分辨率更低。反演是一种逆向求解问题，通过拟合寻优的数学方法，可以逼近地下真实模型。
3.例如“专利公开号为：cn104407393a、名称为：基于时频电磁的自适应遗传模拟退火反演方法以及系统”的中国发明专利，其包括：根据时频电磁的发射频率范围以及分量类型获取反演数据；获取预先设定的反演层数、初始温度、初始步长、最大循环次数；根据所述的反演层数、反演数据确定多个个体模型；对所述的个体模型进行遗传算法运算生成个体新模型；根据所述的初始步长对所述的个体新模型强制进行退火运算生成预测模型；确定所述预测模型与所述个体新模型之间的拟合误差；根据所述的拟合误差判定所述的预测模型是否符合退火标准。
4.但需要指出的是，上述专利公开的方法是通过事先生成多个个体模型，然后由模拟退火算法，通过实际数据对比大量的个体模型生成新的个体模型。需要通过正演事先生成大量的个体模型，同时需要设定初始温度，模拟退火，计算较繁琐，通过预先生成的模型与反演数据得到的拟合差，对初始模型依赖较大。与此同时，由于需要事先生成个体模型，对于复杂模型，需要大量甚至海量的模型，这导致运算时间较长，如果事先生成的个体模型数量较少或较简单，可能会导致反演结果分辨率低，甚至无法逼近真实模型。
5.因此，急需要提出一种原理简单、工艺简单、准确可靠的时频电磁频率域数据自约束稳健电阻率反演方法，其既能保留原始数据丰富的异常信息，又能较准确提取电阻率。
技术实现要素：
6.针对上述问题，本发明的目的在于提供一种时频电磁频率域数据自约束稳健电阻率反演方法，本发明采用的技术方案如下：
7.一种时频电磁频率域数据自约束稳健电阻率反演方法，其包括以下步骤：
8.采集获得地下频率域数据；
9.利用地下频率域数据中的电磁场幅值进行迭代反演；构建模型粗糙度矩阵；
10.迭代反演过程中，预设雅可比矩阵权重系数，采用变步长拉格朗日乘子反演方法进行迭代计算，并更新反演模型；
11.基于柱坐标系计算时频电磁频率域电磁场响应，并对正演算法的验证。
12.进一步地，所述雅可比矩阵权重系数的目标函数u的表达式为：
[0013][0014]
其中，其表示每个频率下的场值与反演模型之间的权重；m表示反演模型参数；m表示反演数据对应的频率个数；λ表示拉格朗日乘子；dj表示实测数据；ε
*
表示拟合差阈值；fj[m]表示正演函数；表示模型粗糙度矩阵。
[0015]
进一步地，所述采用变步长拉格朗日乘子反演方法进行迭代计算，包括以下步骤：
[0016]
预设一拉格朗日变步长反演迭代的阈值、最大迭代次数、反演最大拟合差和初始反演迭代步长；
[0017]
结合初始反演迭代步长，并进行反演迭代；
[0018]
若反演迭代过程中拉格朗日乘子λ小于预设的拉格朗日变步长反演迭代的阈值，则将拉格朗日乘子λ除以k1，并代入反演迭代过程中；所述k1取值为大于1且小于1.5；
[0019]
若反演迭代过程中拉格朗日乘子λ大于预设的拉格朗日变步长反演迭代的阈值，则将拉格朗日乘子λ除以k2，并代入反演迭代过程中；所述k2的取值为大于等于1.6且小于3。
[0020]
进一步地，所述更新反演模型中：
[0021]
求得反演结果的拟合差；若拟合差小于反演最大拟合差或迭代次数大于预设的最大迭代次数，则停止反演。
[0022]
进一步地，所述反演模型的更新表达式为：
[0023][0024]
其中，m
k 1
(λ)表示第k 1次迭代的反演模型参数；jk表示第k次迭代的雅可比矩阵；dk表示第k次迭代的实测数据和反演拟合数据的取自然对数后的绝对误差矩阵；mk表示第k次迭代的反演模型参数；t表示矩阵的转置。
[0025]
进一步地，所述基于柱坐标系计算时频电磁频率域电磁场响应，包括：
[0026]
推导柱坐标系下电偶极子激发的电磁场方程；
[0027]
基于电偶极子激发的电磁场方程，并采用高斯积分求得接地电性源激发的电磁场方程；
[0028]
利用汉克尔积分求解接地电性源电磁场分量。
[0029]
进一步地，所述柱坐标系下电偶极子激发的电磁场方程和其直角坐标系下的水平电场e
x
的关系式为：
[0030]
[0031][0032][0033]ex
＝ercosθ-e
θ
sinθ
[0034][0035]
其中，i表示电流强度；θ表示测点与电偶极子的夹角；μ0表示真空磁导率；r表示测点与源之间的距离；j1、j0表示1阶和0阶的贝塞尔函数；m'表示汉克尔积分项；i表示虚数；w表示圆频率；dx表示偶极子长度；dm'表示汉克尔积分项的微分；ρ1表示第一层电阻率；r表示与地层厚度和地层电阻率相关的传递函数；er表示柱坐标下的电场r方向分量；e
θ
表示柱坐标下的电场θ方向分量；hz表示垂直磁场分量；r*表示与地层厚度相关的传递函数。
[0036]
与现有技术相比，本发明具有以下有益效果：
[0037]
(1)本发明通过预设雅可比矩阵权重系数，采用变步长拉格朗日乘子反演方法进行迭代计算，并更新反演模型，其好处在于，反演迭代开始阶段，由于误差过大，使步长较快减小，使迭代过程加速收敛，在反演后期拟合差较小，使用小步长，在保证反演稳定的情况下，持续逼近真实模型。与现有反演方法相比，具有简单，高效且稳健的特点。
[0038]
(2)本发明通过基于柱坐标系计算时频电磁频率域电磁场响应，其好处在于，使用柱坐标系推导电性源激发的电磁场表达式，能够很好的利用场的轴对称特点，可以更高效的计算电磁场，最后由柱坐标系下的解，通过坐标变换得到直角坐标系下的电磁场解。
[0039]
(3)本发明在正演过程中，通过推导柱坐标系下电偶极子激发的电磁场方程；基于电偶极子激发的电磁场方程，并采用高斯积分求得接地电性源激发的电磁场方程；利用汉克尔积分求解接地电性源电磁场分量，其好处在于，通过汉克尔变化数值逼近，解决层状模型无法求得解析解的难题，同时由高斯积分求解接地线源激发的电磁场，其通过电偶极子表达式，转换成沿着线源方向的积分，再通过高斯积分公式计算，具有简单高效且精确度高的特点。
[0040]
综上所述，本发明具有原理简单、工艺简单、准确可靠等优点，在时频电磁勘探数据反演技术领域具有很高的实用价值和推广价值。
附图说明
[0041]
为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需使用的附图作简单介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对保护范围的限定，对于本领域技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
[0042]
图1本发明中的反演逻辑流程图。
[0043]
图2本发明中的正演逻辑流程图。
[0044]
图3本发明中的正演结果ex与dipole1d对比曲线。
[0045]
图4本发明中的正演结果hz与dipole1d对比曲线。
[0046]
图5本发明的反演模型示例。
具体实施方式
[0047]
为使本技术的目的、技术方案和优点更为清楚，下面结合附图和实施例对本发明作进一步说明，本发明的实施方式包括但不限于下列实施例。基于本技术中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本技术保护的范围。
[0048]
本实施例中，术语“和/或”，仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，a和/或b，可以表示：单独存在a，同时存在a和b，单独存在b这三种情况。
[0049]
本实施例的说明书和权利要求书中的术语“第一”和“第二”等是用于区别不同的对象，而不是用于描述对象的特定顺序。例如，第一目标对象和第二目标对象等是用于区别不同的目标对象，而不是用于描述目标对象的特定顺序。
[0050]
在本技术实施例中，“示例性的”或者“例如”等词用于表示作例子、例证或说明。本技术实施例中被描述为“示例性的”或者“例如”的任何实施例或设计方案不应被解释为比其它实施例或设计方案更优选或更具优势。确切而言，使用“示例性的”或者“例如”等词旨在以具体方式呈现相关概念。
[0051]
在本技术实施例的描述中，除非另有说明，“多个”的含义是指两个或两个以上。例如，多个处理单元是指两个或两个以上的处理单元；多个系统是指两个或两个以上的系统。
[0052]
如图1至图5所示，本实施例提供了一种时频电磁频率域数据自约束稳健电阻率反演方法，其包括以下步骤：
[0053]
第一步，准备时频电磁频率域数据，确定接地线源长度以及电流强度。
[0054]
第二步，给定反演层数，给定拉格朗日变步长反演迭代的阈值，给定反演最大的迭代次数，给定反演的最大拟合差。
[0055]
第三步，给定初始反演迭代步长(拉格朗日初值，给10左右，无量纲)，给定模型粗造度矩阵，给定初始模型(由于计算中取自然对视，所以初值要大于1)，开始反演。
[0056]
第四步，如果迭代中拉格朗日乘子小于给定的阈值，将其除1.1，如果大于给定的阈值，将其值除2，保证随着反演迭代，步长逐渐减小。
[0057]
第五步，计算雅可比矩阵，其中在雅可比计算中，加入权重系数，平衡电磁场场值的量级变化较大问题。
[0058]
第六步，计算新的模型，同时在模型更新中保证迭代前后的连续性，计算反演结果的拟合差(均方根误差，rms)，如果拟合差小于给定的阈值，或者迭代次数大于给定的最大迭代次数，则反演结束；否则，继续反演迭代，直到达到满足条件。
[0059]
最后，保存反演结果。
[0060]
本实施例在反演过程中，将时间域响应数据转化为频率域数据，利用不同频率上的电磁响应信息之差，来估计地下的电导率分布。初始时，通过滤波、去噪等处理后，为反演做准备。
[0061]
其中，基于对初始模型依赖较小的occam反演理论，构建目标函数，其表达式为：
[0062][0063]
其中，其表示每个频率下的场值与反演模型之间的权重；m表示反演模型参数；m表示反演数据对应的频率个数；λ表示拉格朗日乘子；dj表示实测数据；ε
*
表示拟合差阈值；fj[m]表示正演函数；表示模型粗糙度矩阵。
[0064]
本实施例中，采用为每个频率下的场值与反演模型之间的权重，减小由于场值的量级变化过大而引起的收敛过慢或者不收敛问题。
[0065]
在本实施例中，对目标函数求梯度并令其等于0，得到反演模型的更新表达式，即可求得反演结果。
[0066]
通过求取目标函数的梯度，并令其等于0，求其最小值，将问题转化成线性问题，得到模型的更新表达式：
[0067][0068]
其中，m
k 1
(λ)表示第k 1次迭代的反演模型参数；jk表示第k次迭代的雅可比矩阵；dk表示第k次迭代的实测数据和反演拟合数据的取自然对数后的绝对误差矩阵；mk表示第k次迭代的反演模型参数；t表示矩阵的转置。
[0069]
另外，在正演过程中：首先，推导柱坐标系下电偶极子激发的电磁场表达式；其次，基于电偶极子激发的电磁场表达式，通过高斯积分计算接地电性源激发的电磁场表达式；最后，通过汉克尔积分变化求解接地电性源电磁场分量。
[0070]
下面列举以实际案例：
[0071]
首先建立地质模型，地质模型为四层层状地质模型，每层地层具有电阻率和层厚两个参数。这里设定电流强度为1a，接地电性源长度为4km，地层厚度自地表向下分别为400m、200m、2000m、第四层为均匀半空间；电阻率分别为100ω
·
m、10ω
·
m、500ω
·
m、100ω
·
m；接收器平行接地电性源在8km处。
[0072]
本实施案例参照图1，准备好时频电磁勘探的频率域数据，时频电磁勘探频率域数据自约束稳健电阻率反演包含如下步骤：
[0073]
第一步，给定反演层数(本实施案例反演层数为30层)，设定地层厚度按比例增加，比例系数为1.1，给定首层厚度(本实施案例首层厚度为50m)；设置反演跌打步长阈值和初值(本实施案例拉格朗日乘子初始值为10，阈值为0.0005，大于该阈值，步长除2，小于该阈值步长除1.1，保证迭代步长随迭代次数而降低)；给定最大迭代次数；给定最大拟合误差。在本实施例中，k1和k2为经验取值，可以根据拉格朗日初值以及实际反演效果进行调整，其中，k1建议取值大于1且小于1.5，以保证反演稳定。另外，k2建议取大于等于1.6且小于3，以保证前期快速收敛，后期反演稳定。
[0074]
第二步，设置模型粗糙度矩阵，给定反演迭代步长拉格朗日乘子初值(本实施案例为10)，给定模型初始值(迭代运算，需要给予初始值才能运行，本实施案例所有层的电阻率初始值均为10ω
·
m)。
[0075]
第三步，计算雅可比矩阵，设置雅可比矩阵权重系数。
[0076]
第四步，计算模型更新表达式，将迭代前后的模型进行相关联，保证结果的连续性。
[0077]
第五步，计算拟合误差rms，如果rms小于最大拟合误差，或者迭代次数大于最大迭代次数，则保存反演数据，反演终止；否则返回第三步，继续反演，直到达到满足的条件。
[0078]
在本实施例中，理论计算步骤如图2所示，理论模型通过正演，得到理论数据。同时与业界内广泛认可的dipole1d正演程序进行对比，如图3和图4所示，结果表明，本发明的正演算法与其结果基本一直，相对误差在(10-11-10-1
)％以内，说明了本发明正演算法的正确性。
[0079]
如图5所示(rms曲线表明反演稳定；拟合曲线表明反演结果与原始数据完全拟合；反演模型对比表明反演模型与真实模型吻合)。其中，图5(a)表示反演迭代的收敛过程；图5(b)表示反演模型正演电场分量与实测数据电场分量结果的拟合情况；图5(c)表示反演模型与真实模型的对比图。本发明的时频电磁勘探频率域数据自约束稳健电阻率反演结果基本恢复了真实模型的电阻率，可直接使用反演结果进行定量分析。
[0080]
综上所述，本实施案例提供的一种时频电磁勘探频率域数据自约束稳健电阻率反演方法，反演结果与初始模型相互独立，对初始模型依赖度较低，基本原理简单，且反演算法稳定性较高。由场值数据直接反演电阻率，保留了地下异常介质足够多的响应信息，避免了对近区和远区分析的问题，进一步提高电阻率值对地下导电介质的分辨率。
[0081]
上述实施例仅为本发明的优选实施例，并非对本发明保护范围的限制，但凡采用本发明的设计原理，以及在此基础上进行非创造性劳动而作出的变化，均应属于本发明的保护范围之内。