CN112394420B - 一种基于复电阻率参数的矿体位置确定方法及系统 - Google Patents
一种基于复电阻率参数的矿体位置确定方法及系统 Download PDFInfo
- Publication number
- CN112394420B CN112394420B CN202011301082.5A CN202011301082A CN112394420B CN 112394420 B CN112394420 B CN 112394420B CN 202011301082 A CN202011301082 A CN 202011301082A CN 112394420 B CN112394420 B CN 112394420B
- Authority
- CN
- China
- Prior art keywords
- cole
- model
- parameterized
- parameter
- resistivity
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 39
- 230000035945 sensitivity Effects 0.000 claims abstract description 30
- 238000013016 damping Methods 0.000 claims abstract description 17
- 239000011159 matrix material Substances 0.000 claims description 35
- 230000002159 abnormal effect Effects 0.000 claims description 27
- 239000003607 modifier Substances 0.000 claims description 17
- 230000010287 polarization Effects 0.000 claims description 14
- 230000008859 change Effects 0.000 claims description 7
- 229910052500 inorganic mineral Inorganic materials 0.000 claims description 2
- 239000011707 mineral Substances 0.000 claims description 2
- 238000012545 processing Methods 0.000 abstract description 4
- 230000001419 dependent effect Effects 0.000 description 12
- 230000000694 effects Effects 0.000 description 11
- 230000014509 gene expression Effects 0.000 description 9
- 238000001228 spectrum Methods 0.000 description 7
- 238000004364 calculation method Methods 0.000 description 6
- 230000005284 excitation Effects 0.000 description 5
- 230000000875 corresponding effect Effects 0.000 description 4
- 230000007423 decrease Effects 0.000 description 3
- 239000006185 dispersion Substances 0.000 description 3
- 230000006870 function Effects 0.000 description 3
- 239000011435 rock Substances 0.000 description 3
- 230000007246 mechanism Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 241000120529 Chenuda virus Species 0.000 description 1
- 230000033558 biomineral tissue development Effects 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
- G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Measurement Of Resistance Or Impedance (AREA)
Abstract
本发明公开一种基于复电阻率参数的矿体位置确定方法及系统，涉及复电阻率数据处理技术领域，该方法包括：获取Cole‑Cole模型；将所述Cole‑Cole模型重参数化，得到重参数化的Cole‑Cole模型；根据所述重参数化的Cole‑Cole模型，利用所述重参数化的Cole‑Cole模型中的复电阻率值确定所述重参数化的Cole‑Cole模型的灵敏度；利用所述灵敏度，采用阻尼最小二乘法反演理论进行反演，得到所述重参数化的Cole‑Cole模型中的参数反演值；根据所述参数反演值确定矿体位置。本发明提供的方法及系统能够提高利用复电阻率模型进行矿体位置确定的准确率。
Description
技术领域
本发明涉及复电阻率数据处理技术领域，特别是涉及一种基于复电阻率参数的矿体位置确定方法及系统。
背景技术
国内外许多研究者利用不同的参数组合描述线性时不变系统，合理的表达各类岩(矿)石情况的导电机理，先后提出多种等效电路模型。近年，又提出了新的复电阻率模型，或者在已证实的模型基础上提出改进模型。其中，Dias模型较好地描述了各项极化机理的影响，但模型参数较多，具有很强的多解性，并且在实际复电阻率法反演中应用较少；Debye模型对测量频散范围超过两个数量级的电频散数据拟合较差，因此在实际数据处理中受到了很大的限制；当前，在复电阻率法反演中，Cole-Cole模型应用最为广泛，通过反演得到Cole-Cole模型四个参数效果并不特别理想，其中零频电阻率和极化率反演效果较好，频率相关系数和时间常数的反演效果较差。因此，在根据零频电阻率、极化率、频率相关系数和时间常数进行指导矿体位置确定时，得到的矿体位置可能并不准确。
发明内容
本发明的目的是提供一种基于复电阻率参数的矿体位置确定方法及系统，以提高利用复电阻率模型进行矿体位置确定的准确率。
为实现上述目的，本发明提供了如下方案：
一种基于复电阻率参数的矿体位置确定方法，包括：
获取Cole-Cole模型；
将所述Cole-Cole模型重参数化，得到重参数化的Cole-Cole模型；
根据所述重参数化的Cole-Cole模型，利用所述重参数化的Cole-Cole模型中的复电阻率值确定所述重参数化的Cole-Cole模型的灵敏度；
利用所述灵敏度，采用阻尼最小二乘法反演理论进行反演，得到所述重参数化的Cole-Cole模型中的参数反演值；
根据所述参数反演值确定矿体位置。
可选的，所述将所述Cole-Cole模型重参数化，得到重参数化的Cole-Cole模型，具体包括：
将所述Cole-Cole模型中的时间常数替换为设定时间常数，将所述Cole-Cole模型中的极化率替换为相位角，得到重参数化的Cole-Cole模型。
可选的，所述将所述Cole-Cole模型中的时间常数替换为设定时间常数，将所述Cole-Cole模型中的极化率替换为相位角，得到重参数化的Cole-Cole模型，具体包括：
根据如下公式替换所述Cole-Cole模型中的时间常数和所述Cole-Cole模型中的极化率：
其中，m0表示极化率，表示相位角，a(n)表示第一参数，b(n)表示第二参数，n表示模型参数的序号，τρ表示时间常数，ρ表示电阻率，/>表示设定时间常数，/>表示相位，Re()表示取复数的实部，Im()表示取复数的虚部，i表示复数的虚部，c表示频率相关系数。
可选的，所述根据所述重参数化的Cole-Cole模型，利用所述重参数化的Cole-Cole模型中的复电阻率值确定所述重参数化的Cole-Cole模型的灵敏度，具体包括：
改变所述重参数化的Cole-Cole模型中的任一参数，计算参数改变后的复电阻率；
根据所述参数改变后对应的复电阻率确定所述重参数化的Cole-Cole模型的灵敏度。
可选的，所述利用所述灵敏度，采用阻尼最小二乘法反演理论进行反演，得到所述重参数化的Cole-Cole模型中的参数反演值，具体包括：
根据所述灵敏度，确定所述重参数化的Cole-Cole模型参数反演的初值；
以单个Cole-Cole模型表示复电阻率正演模型，利用相对偏差确定所述重参数化的Cole-Cole模型参数反演的初值处的泰勒展开式；所述相对偏差为所述复电阻率正演模型理论正演场值和实测场值之间的拟合程度；
根据所述泰勒展开式确定拟合误差；
判断所述拟合误差是否小于设定拟合误差，得到第一判断结果；
若所述第一判断结果表示为是，则确定所述重参数化的Cole-Cole模型参数反演的初值为所述重参数化的Cole-Cole模型中的参数反演值；
若所述第一判断结果表示为否，则根据所述拟合误差确定异常体电阻率模型参数；
根据所述异常体电阻率模型参数和所述重参数化的Cole-Cole模型参数反演的初值之和更新所述重参数化的Cole-Cole模型参数反演的初值，并返回步骤“以单个Cole-Cole模型表示复电阻率正演模型，利用相对偏差确定所述重参数化的Cole-Cole模型参数反演的初值处的泰勒展开式”。
可选的，所述根据所述泰勒展开式确定拟合误差，具体包括：
根据所述泰勒展开式，利用如下公式确定拟合误差：
其中，表示拟合误差，/>表示泰勒展开式的一阶偏导数，pjk表示雅可比矩阵中的元素，Δxk表示模型修改量，xk表示第k个模型参数，/>表示异常复电阻率的参数数组，/>表示所述重参数化的Cole-Cole模型参数反演的初值，X表示模型参数。
可选的，所述根据所述拟合误差确定异常体电阻率模型参数，具体包括：
根据所述拟合误差和所述拟合误差取极小值的条件确定模型修改量的线性方程组的右端矢量；
根据所述模型修改量的线性方程组的右端矢量和雅可比矩阵确定异常体电阻率模型参数。
可选的，所述根据所述模型修改量的线性方程组的右端矢量和雅可比矩阵确定异常体电阻率模型参数，具体包括：
根据所述模型修改量的线性方程组的右端矢量和雅可比矩阵，利用如下公式确定异常体电阻率模型参数：
(PTP+λD)·ΔX＝S
其中，PT表示雅可比矩阵的转置，P表示雅可比矩阵，λ表示阻尼因子，D表示N×N阶对角矩阵，N表示对角矩阵的行列数，ΔX表示异常体电阻率模型参数，S为右端矢量。
一种基于复电阻率参数的矿体位置确定系统，包括：
获取模块，用于获取Cole-Cole模型；
重参数化模块，用于将所述Cole-Cole模型重参数化，得到重参数化的Cole-Cole模型；
灵敏度确定模块，用于根据所述重参数化的Cole-Cole模型，利用所述重参数化的Cole-Cole模型中的复电阻率值确定所述重参数化的Cole-Cole模型的灵敏度；
参数反演值确定模块，用于利用所述灵敏度，采用阻尼最小二乘法反演理论进行反演，得到所述重参数化的Cole-Cole模型中的参数反演值；
矿体确定模块，用于根据所述参数反演值确定矿体位置。
可选的，所述重参数化模块，具体包括:
重参数化单元，用于将所述Cole-Cole模型中的时间常数替换为设定时间常数，将所述Cole-Cole模型中的极化率替换为相位角，得到重参数化的Cole-Cole模型。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明提供的一种基于复电阻率参数的矿体位置确定方法及系统，通过对Cole-Cole模型重参数化，得到重参数化的Cole-Cole模型，将重参数化的Cole-Cole模型，采用阻尼最小二乘法反演理论进行反演，得到重参数化的Cole-Cole模型中的参数反演值，从而提高了复电阻率模型参数的分辨率，再根据所述参数反演值确定矿体位置，进而提高利用复电阻率模型进行矿体位置确定的准确率。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明基于复电阻率参数的矿体位置确定方法流程图；
图2为本发明重参数化的Cole-Cole模型不同时间常数复电阻率频谱图；
图3为本发明重参数化的Cole-Cole模型不同频率系数复电阻率频谱图；
图4为本发明重参数化的Cole-Cole模型不同相位角复电阻率频谱图；
图5为本发明ρ0＝100Ωm，C＝0.2时重参数化的Cole-Cole模型频谱图；
图6为本发明重参数化的Cole-Cole模型反演示意图；
图7为本发明基于复电阻率参数的矿体位置确定系统示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种基于复电阻率参数的矿体位置确定方法及系统，以提高利用复电阻率模型进行矿体位置确定的准确率。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
实施例一
如图1所示，本发明提供的一种基于复电阻率参数的矿体位置确定方法，包括：
步骤101：获取Cole-Cole模型。
步骤102：将所述Cole-Cole模型重参数化，得到重参数化的Cole-Cole模型。步骤102，具体包括：
将所述Cole-Cole模型中的时间常数替换为设定时间常数，将所述Cole-Cole模型中的极化率替换为相位角，得到重参数化的Cole-Cole模型。所述将所述Cole-Cole模型中的时间常数替换为设定时间常数，将所述Cole-Cole模型中的极化率替换为相位角，得到重参数化的Cole-Cole模型，具体包括：
根据如下公式替换所述Cole-Cole模型中的时间常数和所述Cole-Cole模型中的极化率：
其中，m0表示极化率，表示相位角，a(n)表示第一参数，b(n)表示第二参数，n表示模型参数的序号，τρ表示时间常数，ρ表示电阻率，/>表示设定时间常数，/>表示相位，Re()表示取复数的实部，Im()表示取复数的虚部，i表示复数的虚部，c表示频率相关系数。
步骤103：根据所述重参数化的Cole-Cole模型，利用所述重参数化的Cole-Cole模型中的复电阻率值确定所述重参数化的Cole-Cole模型的灵敏度。步骤103，具体包括：
改变所述重参数化的Cole-Cole模型中的任一参数，计算参数改变后的复电阻率。其中，分别改变重参数化的Cole-Cole模型四个参数中的每一参数，计算改变每一参数对应的复电阻率。采用控制变量法，每次只改变一个参数。
根据所述参数改变后对应的复电阻率确定所述重参数化的Cole-Cole模型的灵敏度。
步骤104：利用所述灵敏度，采用阻尼最小二乘法反演理论进行反演，得到所述重参数化的Cole-Cole模型中的参数反演值。步骤104，具体包括：
根据所述灵敏度，确定所述重参数化的Cole-Cole模型参数反演的初值。
以单个Cole-Cole模型表示复电阻率正演模型，利用相对偏差确定所述重参数化的Cole-Cole模型参数反演的初值处的泰勒展开式；所述相对偏差为所述复电阻率正演模型理论正演场值和实测场值之间的拟合程度。
根据所述泰勒展开式确定拟合误差；所述根据所述泰勒展开式确定拟合误差，具体包括：
根据所述泰勒展开式，利用如下公式确定拟合误差：
其中，表示拟合误差，/>表示泰勒展开式的一阶偏导数，pjk表示雅可比矩阵中的元素，Δxk表示模型修改量，xk表示第k个模型参数，/>表示异常复电阻率的参数数组，/>表示所述重参数化的Cole-Cole模型参数反演的初值，X表示模型参数。
判断所述拟合误差是否小于设定拟合误差，得到第一判断结果。若所述第一判断结果表示为是，则确定所述重参数化的Cole-Cole模型参数反演的初值为所述重参数化的Cole-Cole模型中的参数反演值。
若所述第一判断结果表示为否，则根据所述拟合误差确定异常体电阻率模型参数；所述根据所述拟合误差确定异常体电阻率模型参数，具体包括：
根据所述拟合误差和所述拟合误差取极小值的条件确定模型修改量的线性方程组的右端矢量。
根据所述模型修改量的线性方程组的右端矢量和雅可比矩阵确定异常体电阻率模型参数。所述根据所述模型修改量的线性方程组的右端矢量和雅可比矩阵确定异常体电阻率模型参数，具体包括：
根据所述模型修改量的线性方程组的右端矢量和雅可比矩阵，利用如下公式确定异常体电阻率模型参数：
(PTP+λD)·ΔX＝S
其中，PT表示雅可比矩阵的转置，P表示雅可比矩阵，λ表示阻尼因子，D表示N×N阶对角矩阵，N表示对角矩阵的行列数，ΔX表示异常体电阻率模型参数，S为右端矢量。
根据所述异常体电阻率模型参数和所述重参数化的Cole-Cole模型参数反演的初值之和更新所述重参数化的Cole-Cole模型参数反演的初值，并返回步骤“以单个Cole-Cole模型表示复电阻率正演模型，利用相对偏差确定所述重参数化的Cole-Cole模型参数反演的初值处的泰勒展开式”。
步骤105：根据所述参数反演值确定矿体位置。
实施例二
本发明还提供一种基于复电阻率参数的矿体位置确定方法的具体方式，对经典Cole-Cole模型进行重参数化，包括以下步骤：
(1)Pelton等人通过对大量的岩、矿石标本和露头测量，证明Cole-Cole模型确实可以近似描述激电效应，其复电阻率数学表达式为：
其中，ω为角频率，ρo称为零频电阻率，m0称为极化率，τ称为时间常数，c称为频率相关系数，ρ(ω)为复电阻率，i为复数的虚部。角频率、零频电阻率、极化率、时间常数和频率相关系数统称为Cole-Cole模型或复电阻率频谱参数。
将(1)式进行实虚部结合形式分解，得到
定义两个变量R和I，将其表示为
则复电阻率ρ(ω)实虚部结合形式的表达式为：
其中，ρ(iω)为复电阻率。
(2)对上述经典Cole-Cole模型重参数化，复电阻率相关参数 中用相位角/>设定时间常数/>来代替经典Cole-Cole模型中复电阻率相关参数{ρo，m0，c，τρ}中的极化率m0和时间常数τρ，设定时间常数/>和时间常数τρ之间关系如公式(6)。
相位角与极化率m0的关系表达式如下：
式中，Re()为取复数的实部，a(ω)表示第一参数。
定义两个参数第一参数a和第二参数b：
式中，Im()为取复数的虚部，b(ω)表示第二参数，i为复数的虚部。
根据经典Cole-Cole模型复电阻率表达式(1)可以改写为如下表达式：
ρ(ω)＝ρ0[1-m0(1-(a(ω)+ib(ω)))] (9)
其中，极化率m0采用如下迭代方法计算：
其中，Δm表示极化率的修改量，mo(n)表示极化率的第n个数值，mo(n-1)表示极化率的第n-1个数值。
其中，当m0(0)＝0，经迭代计算，τρ(n)和m0(n)是可以表达为：
则m0(n)如下表达：
将公式(12)代入经典Cole-Cole模型复电阻率表达式(1)，得到下列重参数化的复电阻率表达式：
其中，R，I分别为：
(3)重参数化的复电阻率频谱
在重参数化的Cole-Cole模型中，为了对复电阻率相关参数 影响程度进行分析，采用对四个参数的正常分布范围控制其中三个参数数值不变，研究某一参数变化时，计算其复电阻率值大小，统计分析从而得出参数变化对复电阻率值得影响程度，也就是灵敏度。参数的灵敏度高，在反演中能够恢复的越好，这也为后续反演认识打下基础。初始参数设置为：零频电阻率ρo＝100Ωm，相位角/>频率相关系数c＝0.2，设定时间常数/> 参数的变化范围：零频电阻率ρo：10-1～103Ωm，相位角8.34～71.20mrad，频率相关系数c：0.1～0.5，设定时间常数/>0.077～770s。
图2表示重参数化的Cole-Cole模型不同时间常数τ对复电阻率的影响，图2(a)为本发明重参数化的Cole-Cole模型不同设定时间常数复电阻率频谱图，其中，是以频率Hz为横坐标，以复电阻率ρ(iω)的虚部为纵坐标；图2(b)为本发明重参数化的Cole-Cole模型不同设定时间常数/>复电阻率频谱图，其中，是以频率Hz为横坐标，以复电阻率ρ(iω)的实部为纵坐标；图2(c)为本发明重参数化的Cole-Cole模型不同设定时间常数/>复电阻率频谱图，其中，是以频率Hz为横坐标，以相位φ(ω)为纵坐标；图2(d)为本发明重参数化的Cole-Cole模型不同设定时间常数/>复电阻率频谱图，其中，是以频率Hz为横坐标，以复电阻率ρ(iω)的振幅|ρ(iω)|为纵坐标；随着时间常数的增大，复电阻率的振幅和实分量衰减率逐渐减少，数值也逐渐减小；虚分量与相位在相应频段10-3～103Hz出现极小值，曲线呈现近似对称，频段范围原因未显示10-3Hz以下低频；随着时间常数的数值增大，虚分量极小值与对称区间向低频段移动。与经典Cole-Cole模型相似。
图3表示重参数化的Cole-Cole模型不同频率相关系数c对复电阻率的影响，图3(a)为本发明重参数化的Cole-Cole模型不同设定频率相关系数c复电阻率频谱图；其中，是以频率Hz为横坐标，以复电阻率ρ(iω)的虚部为纵坐标；图3(b)为本发明重参数化的Cole-Cole模型不同设定频率相关系数c复电阻率频谱图；其中，是以频率Hz为横坐标，以复电阻率ρ(iω)的实部为纵坐标；图3(c)为本发明重参数化的Cole-Cole模型不同设定频率相关系数c复电阻率频谱图，其中，是以频率Hz为横坐标，以相位φ(ω)为纵坐标；图3(d)为本发明重参数化的Cole-Cole模型不同设定频率相关系数c复电阻率频谱图，其中，是以频率Hz为横坐标，以复电阻率ρ(iω)的振幅|ρ(iω)|为纵坐标；随着频率相关系数c值增大，实分量与振幅的曲线趋缓，且频率相关系数c值越大，影响越小，变化范围更小。虚分量与相位在频段范围内出现极小值区间，且曲线存在轴对称现象，在频率100～101Hz之间存在几乎同一极小值。此外，频率相关系数c值越高，则复电阻率虚部以及相位变化范围越大。
图4表示重参数化的Cole-Cole不同相位角对复电阻率的影响，图4(a)为本发明重参数化的Cole-Cole模型不同设定相位角复电阻率频谱图；其中，是以频率Hz为横坐标，以复电阻率ρ(iω)的虚部为纵坐标；图3(b)为本发明重参数化的Cole-Cole模型不同设定相位角/>复电阻率频谱图；其中，是以频率Hz为横坐标，以复电阻率ρ(iω)的实部为纵坐标；图3(c)为本发明重参数化的Cole-Cole模型不同设定相位角/>复电阻率频谱图，其中，是以频率Hz为横坐标，以相位φ(ω)为纵坐标；图3(d)为本发明重参数化的Cole-Cole模型设定相位角/>复电阻率频谱图，其中，是以频率Hz为横坐标，以复电阻率ρ(iω)的振幅|ρ(iω)|为纵坐标；随着相位角的增大，复电阻率的振幅和实分量曲线的衰减强度逐渐加强；虚分量与相位在频段10-3～103Hz之间呈近似对称现象，且在该频段出现极小值。
综上所示图2至图4对比分析，复电阻率虚分量以及相位均恒为负值；复电阻率实分量与振幅随频率增大而减小。在零频电阻率变化时，影响复电阻率最大，振幅与实分量大小随零频电阻率大小而变，因此复电阻率影响程度依次为
相较于经典Cole-Cole模型，两者频率相关系数值影响结果有所不同。
重参数化的模型等效参数为：ρo＝100Ωm，C＝0.2绘制复电阻率频谱特性，如图5所示，图5(a)为本发明重参数化的Cole-Cole模型ρo＝100Ωm，/>C＝0.2时复电阻率频谱图，其中，是以频率为横坐标，以复电阻率ρ(iω)的振幅为纵坐标；图5(b)为本发明本发明重参数化的Cole-Cole模型ρo＝100Ωm，/> C＝0.2时复电阻率频谱图，其中，是以频率为横坐标，以复电导率相位为纵坐标；图5(c)为本发明重参数化的Cole-Cole模型ρo＝100Ωm，/>C＝0.2时复电阻率频谱图，其中，是以频率为横坐标，以复电导率虚分量为纵坐标；图5(d)为本发明本发明重参数化的Cole-Cole模型ρo＝100Ωm，/>C＝0.2时复电阻率频谱图，其中，是以频率为横坐标，以复电阻率虚分量为纵坐标；得到与经典Cole-Cole模型一致的曲线特征。
(4)反演理论及雅克比矩阵计算
①最小二乘反演理论
使用单个Cole–Cole模型表示复电阻率正演模型，选用阻尼最小二乘法(马奎特法)进行反演。用fsj表示实测场值，Fj(r)表示理论正演场值，表示异常体复电阻率的参数{ρo，m，τ，c}数组，理论和实测场值之间的拟合程度用相对偏差/>表示：
这样反演拟合误差为：
上式中下标j＝1、2、·····、m表示第j个工作频率点。
由于正演函数是非线性的，所以偏差函数和拟合误差也是非线性的，为了克服求解非线性方程组的难点，需要对偏差函数作线性化近似处理，先对复电阻率重参数化的-Cole-Cole模型四个参数赋予初值/>再对/>在/>处作泰勒展开，并忽略二阶以上的各阶偏导数，于是有：
式中下标k表示第k个模型参数，为模型修改量，/>表示相对偏差在/>处的值。再设/>由此可得到拟合误差表达式：
此时，拟合误差被表示成为复电导率模型修改量Δx1,Δx2,...,Δxn的多元函数，其取极小值的条件为：
由此推导得，
这样分别取j＝1,2,...,n,可以推导出如下求解模型修改量的线性方程组
(PTP+λD)·ΔX＝S (22)
式中，P为雅克比矩阵，其元素为pjk，ΔX＝(Δx1,Δx2,...,Δxn)T，右端矢量λ为阻尼因子，是一个大于零的常数，D为一个N×N阶对角矩阵，N表示对角矩阵的行列数，其对角线上元素为/>
根据(22)式求出异常体电阻率模型参数ΔX，并以作为新的模型参数初值，重新计算拟合误差。这样反复迭代以达到拟合误差小于设定拟合误差ε，此时的/>即为所求的反演结果，也即得到重参数化的Cole-Cole模型四个参数/>
②雅克比矩阵计算
最小二乘反演中需要涉及求取雅克比矩阵，首先需要求取复电阻率对重参数化的Cole-Cole模型四个参数的偏导数，从而得到反演过程中的雅克比矩阵，(22)式只涉及求取复电导率的雅克比矩阵，而复电阻率反演中需要求取重参数化的Cole-Cole模型四个参数/>的雅克比矩阵，其计算推导如下：
a.复电阻率ρ(iω)对相位角求偏导：
b.复电阻率ρ(iω)对频率相关系数c求偏导：
c.复电阻率ρ(iω)对时间常数求偏导：
d.复电阻率ρ(iω)对零频时电阻率ρo求偏导：
其中，推导如下：
/>
将求取得到的雅克比矩阵作为变量P代入到(22)式中，采用阻尼最小二乘法反演理论以及不同Cole-Cole模型参数雅可比矩阵，使用Fortran语言编程，得到正演理论值，再进行参数赋初值，然后用Fortran语言实现的反演程序进行反演计算，加入高斯误差，通过迭代计算，最终得到重参数化的Cole-Cole模型四个参数的反演值采用matlab程序语言对反演结果进行成图显示，将其重参数化的Cole-Cole模型反演结果与典型Cole-Cole模型参数反演结果进行对比，最后进行分析反演结果的分辨率。
本发明的效果在于：经典Cole-Cole模型经过重参数化得到重参数化的Cole-Cole模型，其等效参数有效解决了极化率m与频率相关系数c的强相关性。重参数化的Cole-Cole模型反演的各参数值更接近真值，反演效果更好，提高了参数分辨率，其对于时域激电和频域激电数据同样适用。
如图6所示，基于经典Cole-Cole模型，进行重参数化计算，得到重参数化的重参数化的Cole-Cole模型，对重参数化的-Cole-Cole模型中，复电阻率相关参数影响程度进行分析，采用对四个参数的正常分布范围控制变量研究某一参数的影响程度。然后，采用阻尼最小二乘法，结合推导得到的重参数化的Cole-Cole模型参数雅克比矩阵，通过Fortran语言编程实现的反演程序进行反演得到复电阻率模型参数，有效改善了复电阻率模型参数的分辨率，从而大大改善了复电阻率数据反演的精度，为复电阻率法野外实测数据反演提高技术支撑。对野外实际采集得到的复电阻率数据，通过采用重参数化后的复电阻率模型也即重参数化的Cole-Cole模型，由已得到的各模型参数的灵敏度，结合已知地质钻孔信息，设置反演的各模型参数初始模型，进行最小二乘反演，能够有效提高获得模型参数的分辨率，从而提供更可靠的地下介质的参数信息，更好地为解决地质问题服务。
实施例三
如图7所示，本发明提供的一种基于复电阻率参数的矿体位置确定系统，包括：
获取模块201，用于获取Cole-Cole模型。
重参数化模块202，用于将所述Cole-Cole模型重参数化，得到重参数化的Cole-Cole模型。
灵敏度确定模块203，用于根据所述重参数化的Cole-Cole模型，利用所述重参数化的Cole-Cole模型中的复电阻率值确定所述重参数化的Cole-Cole模型的灵敏度。
参数反演值确定模块204，用于利用所述灵敏度，采用阻尼最小二乘法反演理论进行反演，得到所述重参数化的Cole-Cole模型中的参数反演值。
矿体确定模块205，用于根据所述参数反演值确定矿体位置。
其中，所述重参数化模块202，具体包括:
重参数化单元，用于将所述Cole-Cole模型中的时间常数替换为设定时间常数，将所述Cole-Cole模型中的极化率替换为相位角，得到重参数化的Cole-Cole模型。
频率相关系数与时间常数在CR或SIP勘探中应作为重要的研究常数，在电法勘探中，零频电阻率可用于表征介质导电性强弱；极化率表示介质激电效应强弱；频率相关系数与时间常数在复电阻率法勘探中应作为重要的研究常数，时间常数可以直接区分极化体，并有可能在激电强度参数没有明显异常的情况下找到深部矿；频率相关系数的极化特征也能够从矿化围岩中划分出局部矿化体。为改善反演得到复电阻率模型模型参数的分辨率，提出先将典型Cole-Cole模型进行重参数化，再用复电阻率法数据进行反演，从而使得反演得到复电阻率模型参数的分辨率提高，从而根据重参数化的Cole-Cole模型的四个参数确定矿体位置。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。
Claims (6)
1.一种基于复电阻率参数的矿体位置确定方法，其特征在于，包括：
获取Cole-Cole模型；
将所述Cole-Cole模型重参数化，得到重参数化的Cole-Cole模型；
根据所述重参数化的Cole-Cole模型，利用所述重参数化的Cole-Cole模型中的复电阻率值确定所述重参数化的Cole-Cole模型的灵敏度；
利用所述灵敏度，采用阻尼最小二乘法反演理论进行反演，得到所述重参数化的Cole-Cole模型中的参数反演值；
根据所述参数反演值确定矿体位置；
所述将所述Cole-Cole模型重参数化，得到重参数化的Cole-Cole模型，具体包括：
将所述Cole-Cole模型中的时间常数替换为设定时间常数，将所述Cole-Cole模型中的极化率替换为相位角，得到重参数化的Cole-Cole模型；
所述将所述Cole-Cole模型中的时间常数替换为设定时间常数，将所述Cole-Cole模型中的极化率替换为相位角，得到重参数化的Cole-Cole模型，具体包括：
根据如下公式替换所述Cole-Cole模型中的时间常数和所述Cole-Cole模型中的极化率：
其中，m0表示极化率，mo(n)表示极化率的第n个数值，表示相位角，a(n)表示第一参数，b(n)表示第二参数，n表示模型参数的序号，τρ表示时间常数，ρ表示电阻率，/>表示设定时间常数，/>表示相位，Re()表示取复数的实部，Im()表示取复数的虚部，i表示复数的虚部，c表示频率相关系数；
所述根据所述重参数化的Cole-Cole模型，利用所述重参数化的Cole-Cole模型中的复电阻率值确定所述重参数化的Cole-Cole模型的灵敏度，具体包括：
改变所述重参数化的Cole-Cole模型中的任一参数，计算参数改变后的复电阻率；
根据所述参数改变后对应的复电阻率确定所述重参数化的Cole-Cole模型的灵敏度；
所述利用所述灵敏度，采用阻尼最小二乘法反演理论进行反演，得到所述重参数化的Cole-Cole模型中的参数反演值，具体包括：
根据所述灵敏度，确定所述重参数化的Cole-Cole模型参数反演的初值；
以单个Cole-Cole模型表示复电阻率正演模型，利用相对偏差确定所述重参数化的Cole-Cole模型参数反演的初值处的泰勒展开式；所述相对偏差为所述复电阻率正演模型理论正演场值和实测场值之间的拟合程度；
根据所述泰勒展开式确定拟合误差；
判断所述拟合误差是否小于设定拟合误差，得到第一判断结果；
若所述第一判断结果表示为是，则确定所述重参数化的Cole-Cole模型参数反演的初值为所述重参数化的Cole-Cole模型中的参数反演值；
若所述第一判断结果表示为否，则根据所述拟合误差确定异常体电阻率模型参数；
根据所述异常体电阻率模型参数和所述重参数化的Cole-Cole模型参数反演的初值之和更新所述重参数化的Cole-Cole模型参数反演的初值，并返回步骤“以单个Cole-Cole模型表示复电阻率正演模型，利用相对偏差确定所述重参数化的Cole-Cole模型参数反演的初值处的泰勒展开式”。
2.根据权利要求1所述的基于复电阻率参数的矿体位置确定方法，其特征在于，所述根据所述泰勒展开式确定拟合误差，具体包括：
根据所述泰勒展开式，利用如下公式确定拟合误差：
其中，表示拟合误差，/>表示泰勒展开式的一阶偏导数，pjk表示雅可比矩阵中的元素，Δxk表示模型修改量，xk表示第k个模型参数，/>表示异常复电阻率的参数数组，/>表示所述重参数化的Cole-Cole模型参数反演的初值，X表示模型参数。
3.根据权利要求2所述的基于复电阻率参数的矿体位置确定方法，其特征在于，所述根据所述拟合误差确定异常体电阻率模型参数，具体包括：
根据所述拟合误差和所述拟合误差取极小值的条件确定模型修改量的线性方程组的右端矢量；
根据所述模型修改量的线性方程组的右端矢量和雅可比矩阵确定异常体电阻率模型参数。
4.根据权利要求3所述的基于复电阻率参数的矿体位置确定方法，其特征在于，所述根据所述模型修改量的线性方程组的右端矢量和雅可比矩阵确定异常体电阻率模型参数，具体包括：
根据所述模型修改量的线性方程组的右端矢量和雅可比矩阵，利用如下公式确定异常体电阻率模型参数：
(PTP+λD)·ΔX＝S
其中，PT表示雅可比矩阵的转置，P表示雅可比矩阵，λ表示阻尼因子，D表示N×N阶对角矩阵，N表示对角矩阵的行列数，ΔX表示异常体电阻率模型参数，S为右端矢量。
5.一种基于复电阻率参数的矿体位置确定系统，其特征在于，所述基于复电阻率参数的矿体位置确定系统应用权利要求1所述的基于复电阻率参数的矿体位置确定方法，所述基于复电阻率参数的矿体位置确定系统包括：
获取模块，用于获取Cole-Cole模型；
重参数化模块，用于将所述Cole-Cole模型重参数化，得到重参数化的Cole-Cole模型；
灵敏度确定模块，用于根据所述重参数化的Cole-Cole模型，利用所述重参数化的Cole-Cole模型中的复电阻率值确定所述重参数化的Cole-Cole模型的灵敏度；
参数反演值确定模块，用于利用所述灵敏度，采用阻尼最小二乘法反演理论进行反演，得到所述重参数化的Cole-Cole模型中的参数反演值；
矿体确定模块，用于根据所述参数反演值确定矿体位置。
6.根据权利要求5所述的基于复电阻率参数的矿体位置确定系统，其特征在于，所述重参数化模块，具体包括:
重参数化单元，用于将所述Cole-Cole模型中的时间常数替换为设定时间常数，将所述Cole-Cole模型中的极化率替换为相位角，得到重参数化的Cole-Cole模型。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011301082.5A CN112394420B (zh) | 2020-11-19 | 2020-11-19 | 一种基于复电阻率参数的矿体位置确定方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011301082.5A CN112394420B (zh) | 2020-11-19 | 2020-11-19 | 一种基于复电阻率参数的矿体位置确定方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112394420A CN112394420A (zh) | 2021-02-23 |
| CN112394420B true CN112394420B (zh) | 2023-10-20 |
Family
ID=74607495
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011301082.5A Active CN112394420B (zh) | 2020-11-19 | 2020-11-19 | 一种基于复电阻率参数的矿体位置确定方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112394420B (zh) |
Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5121363A (en) * | 1990-12-26 | 1992-06-09 | Conoco Inc. | Fracture detection logging tool |
| US6485703B1 (en) * | 1998-07-31 | 2002-11-26 | The Texas A&M University System | Compositions and methods for analyte detection |
| RU2235347C1 (ru) * | 2003-08-11 | 2004-08-27 | РЫХЛИНСКИЙ Николай Иванович | Способ геоэлектроразведки (варианты) |
| JP2006177802A (ja) * | 2004-12-22 | 2006-07-06 | Central Res Inst Of Electric Power Ind | 地盤探査方法および地盤探査装置 |
| KR100806207B1 (ko) * | 2007-02-08 | 2008-02-22 | 한국지질자원연구원 | 다중 주파수 유도분극 역산 방법 |
| CN101520517A (zh) * | 2008-02-25 | 2009-09-02 | 中国石油集团东方地球物理勘探有限责任公司 | 一种能准确评价碎屑岩盆地含油气目标的方法 |
| CN101706587A (zh) * | 2009-11-24 | 2010-05-12 | 中南大学 | 一种电法勘探激电模型参数的提取方法 |
| CN101937105A (zh) * | 2010-04-28 | 2011-01-05 | 中国石油大学(北京) | 通过低频信号检测油气层的方法及装置 |
| CN103207413A (zh) * | 2011-11-16 | 2013-07-17 | 中国地质大学（北京） | 电法勘探装置及系统 |
| CN105044782A (zh) * | 2015-07-09 | 2015-11-11 | 成都理工大学 | 一种海洋地下介质总有机碳含量的获取方法 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7728596B2 (en) * | 2005-12-15 | 2010-06-01 | Evgenij Dmitrievich Lisitsyn | Method and device for sea electrical survey of oil-and-gas deposits |
| EP2102688B1 (en) * | 2007-01-03 | 2013-06-19 | Council of Scientific & Industrial Research | A process and device for measurement of spectral induced polarization response using pseudo random binary sequence (prbs) current source |
| US7863901B2 (en) * | 2007-05-25 | 2011-01-04 | Schlumberger Technology Corporation | Applications of wideband EM measurements for determining reservoir formation properties |
-
2020
- 2020-11-19 CN CN202011301082.5A patent/CN112394420B/zh active Active
Patent Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5121363A (en) * | 1990-12-26 | 1992-06-09 | Conoco Inc. | Fracture detection logging tool |
| US6485703B1 (en) * | 1998-07-31 | 2002-11-26 | The Texas A&M University System | Compositions and methods for analyte detection |
| RU2235347C1 (ru) * | 2003-08-11 | 2004-08-27 | РЫХЛИНСКИЙ Николай Иванович | Способ геоэлектроразведки (варианты) |
| JP2006177802A (ja) * | 2004-12-22 | 2006-07-06 | Central Res Inst Of Electric Power Ind | 地盤探査方法および地盤探査装置 |
| KR100806207B1 (ko) * | 2007-02-08 | 2008-02-22 | 한국지질자원연구원 | 다중 주파수 유도분극 역산 방법 |
| CN101520517A (zh) * | 2008-02-25 | 2009-09-02 | 中国石油集团东方地球物理勘探有限责任公司 | 一种能准确评价碎屑岩盆地含油气目标的方法 |
| CN101706587A (zh) * | 2009-11-24 | 2010-05-12 | 中南大学 | 一种电法勘探激电模型参数的提取方法 |
| CN101937105A (zh) * | 2010-04-28 | 2011-01-05 | 中国石油大学(北京) | 通过低频信号检测油气层的方法及装置 |
| CN103207413A (zh) * | 2011-11-16 | 2013-07-17 | 中国地质大学（北京） | 电法勘探装置及系统 |
| CN105044782A (zh) * | 2015-07-09 | 2015-11-11 | 成都理工大学 | 一种海洋地下介质总有机碳含量的获取方法 |
Non-Patent Citations (3)
| Title |
|---|
| 含甲烷水合物多孔介质的复电阻率模型对比;王彩程;陈强;邢兰昌;刘昌岭;郑金吾;;海洋地质前沿(05);全文 * |
| 复电阻率法二维正演并行计算研究;张志勇;谭捍东;王春阳;汪茂;何中波;梁鹏;;物探化探计算技术(03);全文 * |
| 频谱激电法的发展与展望;杨振威;郑伟;李晓斌;王华峰;;物探与化探(01);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112394420A (zh) | 2021-02-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106845010B (zh) | 基于改进SVD降噪和Prony的低频振荡主导模式辨识方法 | |
| Cohen et al. | A statistical study of temporally smoothed wavelet coherence | |
| US6366858B1 (en) | Method of and apparatus for independently determining the resistivity and/or dielectric constant of an earth formation | |
| Luo et al. | Interpolated DFT algorithms with zero padding for classic windows | |
| Yin et al. | Imaging the continuous conductivity profile within layered metal structures using inductance spectroscopy | |
| WO2017034924A1 (en) | Borehole acoustic logging receiver quality control and calibration | |
| US20110166842A1 (en) | Layer stripping method | |
| CN113591028A (zh) | 适用于直流偏磁风险评估的41点快速汉克尔变换的方法 | |
| Ferreira et al. | Comparison of simultaneous bender elements and resonant column tests on Porto residual soil | |
| Hui et al. | Some problems with the method of fundamental solution using radial basis functions | |
| CN114137616B (zh) | 一种利用量子力学原理进行储层含气性检测的方法 | |
| CN112394420B (zh) | 一种基于复电阻率参数的矿体位置确定方法及系统 | |
| CN113447992B (zh) | 一种运用时间域激发极化法进行矿体勘探的方法及系统 | |
| Watzenig et al. | Adaptive regularization parameter adjustment for reconstruction problems | |
| CN109541304A (zh) | 基于六项最小旁瓣窗插值的电网高次弱幅值谐波检测方法 | |
| Zorin | Spectral induced polarization of low and moderately polarizable buried objects | |
| WO2002003100A2 (en) | Determining earth formation resistivity by attenuation and phase shift measurements | |
| CN110320400B (zh) | 准同步采样和改进能量算子的电压闪变包络参数提取方法 | |
| CN116859457A (zh) | 时空-频率波数域同步挤压叠置致密河道储层表征方法 | |
| Stoynov | Structural spectral analysis of electrochemical impedance | |
| CN116127287A (zh) | 一种电阻率法勘探信号的降噪方法 | |
| Kelter et al. | Quantitative imaging of spectral electrical properties of variably saturated soil columns | |
| CN110750893B (zh) | 一种基于小波分析的时变电推力器建模方法 | |
| CN116203316A (zh) | 基于复频谱共轭插值的电网信号参数估计方法及系统 | |
| CN107917733A (zh) | 一种基于模型的导电结构厚度与电导率涡流检测方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |