本发明涉及地球物理勘探领域，具体地来讲为一种基于粒子群优化算法的电导-极化率多参数成像方法。
背景技术：
：：在电法勘探领域，电阻率解释是一个非常重要的内容，通常将其用来估算地下介质的分布。受激发极化效应的影响，岩(矿)石电阻率是一个与频率有关的复数，因此在向地下供入稳定电流时，不仅能产生以二次涡流场为主的感应场，也能产生极化场。极化场的存在，使得电磁响应产生“反号”现象，该“反号”包含大量的极化信息。因此，研究极化现象，准确提取地下介质中的极化参数，进行极化率、电导率-深度成像是十分必要的。目前已经应用到实际电法勘探极化参数提取的方法中，最为广泛的是constable(1987)提出的奥克姆反演方法，但在存储反演时面临着灵敏度矩阵需要多次正演的问题。闫国祥(2015)应用最小二乘法对极化参数进行了分步提取，但结果过于依赖第一步反演结果，而且仅是理论模型研究，没有应用于实际例子。cn109376446a公开了一种基于粒子群优化的极化曲线拟合方法，特别是涉及在电化学领域腐蚀动力学参数的提取，待拟合的动力学电流曲线呈e指数非线性变化，但是其使用的是基本粒子群算法对整条曲线进行拟合，并没有对特征区间进行划分，不足以描述极化特征，而且随着种群数目、待求参数量的增大，时间成本也在累计增加，为此将粒子群算法进行优化是十分必要的。cn101706587a公开了一种电法勘探激电模型参数的提取方法，特别是涉及在柯尔-柯尔模型中的未知参数。该方法结合随机统计算法与最小二乘法，对单个柯尔-柯尔模型相对相位谱和振幅联合提取极化效应参数，既避免了随机统计算法中过于依赖初值的问题，又保留了最小二乘法的提取速度较快的特点。但以上参数提取方法针对的仅是复电导率，在实际野外测量中，不能直接测得复电导率值，更多的是得到电磁响应数据，因此严重影响参数提取在实际野外中的应用，所以研究实际极化地质的极化参数提取方法具有重要意义。技术实现要素：本发明的目的在于针对现有极化参数提取方法的不足，根据实际地质中存在的极化问题，建立柯尔-柯尔模型，提供一种基于粒子群优化算法的电导-极化率多参数成像方法。本发明是这样实现的，一种基于粒子群优化算法的电导-极化率多参数成像方法，该方法包括：1)基于地下极化介质的柯尔-柯尔复电导率模型，获得梯形波激励下的时域感应-极化磁场响应；2)采用超导量子传感器时域电磁探测系统进行野外磁场数据观测，并对实测数据进行叠加取样、滤波、以及一次场剔除处理；3)利用步骤2)数据的早期数据计算零频电导率值、晚期数据e指数多项式拟合得到时间常数，将获得的零频电导率值和时间常数对粒子群搜索范围进行约束；4)根据实测磁场数据与步骤1的时域感应-极化磁场响应构建目标函数，利用粒子群优化算法，根据步骤3约束的范围，提取极化参数，包括零频电导率σ0、时间常数τ、极化率η、频散系数c、深度d；5)利用步骤4提取到实测数据的极化参数，带入广义趋肤深度计算公式，进行极化介质的电导率、极化率-深度成像。进一步地，步骤3中，发射电流关断后，早期响应以感应场为主，利用早期数据计算零频电导率值σ1；晚期数据以极化场为主，对其晚期数据进行e指数拟合，得到时间常数值τ1；零频电导率值σ1和时间常数值τ1的±20％分别作为零频电导率和时间常数值的粒子群搜索范围上下限。进一步地，步骤4中具体包含以下步骤：4ⅰ、输入相关参数，随机初始化粒子的位置、速度和个体最优值pbest、种群最优值gbest；4ⅱ、构建目标函数：minφ(x)＝φ1(x)+λφ2(x)，x＝σ0、η、τ、c，其中φ2(x)＝ft(x)，ft(x)为步骤1的时域感应-极化磁场响应，bt为实测磁场数据，λ为收敛因子；4ⅲ、计算每个粒子当前适应度，与pbest、gbest进行比较，若当前适应度值较小，则令其替代pbest、gbest，反之继续保留pbest、gbest值；4ⅳ、实测曲线依据斜率划分为两段，据此分别按照式(1)和式(2)更新粒子的速度和位置：xij(k+1)＝xij(k)+vij(k+1)·t(2)式中，t0为曲线斜率由负变为正的时刻，pk为最大迭代次数，iter为当前迭代次数；迭代过程中，当粒子的速度或位置超出搜索范围，产生变异时，根据二者移动方向分为两种情况，若二者变异方向相反，则不做任何处理，反之，将由下式对位置和速度重新定义，当超出上边界时：vij(k+1)＝vmax-2·r3·vmax(3)xij(k+1)＝xmax-0.5·r4·xmax(4)当粒子的位置或速度超出下边界时：vij(k+1)＝vmin-2·r3·vmin(5)xij(k+1)＝xmin-0.5·r4·xmin(6)4ⅴ、判断是否达到最大迭代次数，若是，获得pbest、gbest的位置信息，输出σ0、τ、η、c、d；否则返回步骤4ⅲ，继续执行。进一步地，步骤5中包含以下步骤：5ⅰ、应用matlab编程，输入步骤4提取出的零频电导率σ0、时间常数τ、极化率η、频散系数c到程序中，根据公式，获得复电导率值；5ⅱ、根据广义趋肤深度公式，计算该复电导率σ(ω)的趋肤深度；5ⅲ、利用步骤5ⅰ中的零频电导率、极化率和步骤5ⅱ获得的趋肤深度，分别进行电导率-深度、极化率-深度成像。本发明与现有技术相比，有益效果在于：本发明的极化介质电导、极化率参数成像方法，相比传统电阻率成像方法，更符合实际地下极化介质的传播扩散规律，从而提高了电导率探测深度的解释精度。本方法为开展电磁探测寻找矿产资源提供新的技术保障，有利于瞬变电磁勘探方法的实用化。附图说明图1是基于地下极化介质模型的电导率、极化率多参数成像方法流程图；图2是粒子群优化算法提取极化参数的流程图；图3是根据地质资料建立的三层层状大地模型；图4是本发明一个实施例的层状模型提取结果成像对比图。具体实施方式为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。本发明在电流关断后100ms内，采用超导量子传感器时域电磁探测系统观测电磁数据，晚期数据采用e指数拟合确定时间常数τ1；根据发射电流计算一次磁场，从实测数据中减去一次磁场后，计算零频电导率σ1；二者大小的±20％作为搜索范围输入粒子群中，根据理论模型的磁场响应与实测磁场响应构建目标函数，最后采用粒子群优化算法进行零频电导率σ0、时间常数τ、极化率η、频散系数c、深度d的提取。本发明是这样实现的，一种基于粒子群优化算法的电导-极化率多参数成像方法包括：1)基于地下极化介质的柯尔-柯尔复电导率模型，获得梯形波激励下的时域感应-极化磁场响应表达式；先获取层状模型下负阶跃波的磁场响应，将梯形波发射电流一阶导数与磁场响应卷积，可得到梯形波激励下的感应磁场响应；再带入柯尔-柯尔模型频域表达式其中σ0为零频电导率、τ为时间常数、η为极化率、c为频散系数，最后采用基于贝塞尔函数的快速hankel变换、余弦变换算法，将结果从频域转换到时域，从而获得层状大地模型下时域感应-极化双场电磁响应，即理论模型磁场响应。2)使用超导量子传感器时域电磁探测系统接收磁场，对电磁实测数据进行处理，包括叠加取样、抽道、滤波、一次场剔除等处理；3)处理后的数据，利用早期数据(一般为100微秒～1毫秒)计算零频电导率值，利用晚期数据(一般二十个毫秒以后)e指数拟合得到时间常数，二者对粒子群搜索范围进行约束；电流关断后，晚期数据以极化场为主，对其进行e指数拟合，得到时间常数值τ1；早期响应以感应场为主，根据发射电流计算一次磁场，从实测数据中减去一次磁场后，计算零频电导率σ1。4)根据实测磁场数据与理论模型磁场响应构建目标函数，利用粒子群优化算法进行参数提取，包括零频电导率σ0、时间常数τ、极化率η、频散系数c、深度d；基于粒子群优化算法的电导-极化率多参数提取，具体步骤如下：4ⅰ、输入粒子群优化算法相关参数d、pk、bk、xlimit、vlimit，其中d为问题维数，bk为种群规模，pk为最大迭代次数，为加速度因子，xlimit为粒子的位置限制，vlimit为速度限制；随机初始化粒子的位置、速度、个体最优值pbest和种群最优值gbest；4ⅱ、构建目标函数：minφ(x)＝φ1(x)+λφ2(x)，x＝σ0、η、τ、c，其中φ2(x)＝ft(x)，bt为实测磁场数据，ft(x)为理论模型磁场响应，λ为收敛因子；4ⅲ、计算每个粒子当前适应度，与pbest、gbest进行比较，若当前适应度值较小，则令其替代pbest、gbest，反之继续保留pbest、gbest值；4ⅳ、实测曲线依据斜率划分为两段，据此更新粒子的速度(式(1))和位置(式(2))：xij(k+1)＝xij(k)+vij(k+1)·t(2)式中，t0为斜率由负变为正的时刻，pk为最大迭代次数，iter为当前迭代次数；下标i表示第i个粒子(i＝1,2,3......m)，下标j表示粒子的第j维(j＝1,2,3......d)，i表示当前的迭代次数，假设粒子群搜索空间为d维，粒子群有m个粒子(通常将粒子数设置为问题维数的5～10倍)，r1、r2是(0，1)之间的随机数。为了保证粒子在搜索空间内搜索，通常将粒子的速度限制在[-vmaxvmax]区间中，将步骤3得到的搜索范围作为粒子的位置限制[-xminxmax]，则粒子的最大速度可以取xmax的10％～20％。由于粒子的速度和位置均有方向性，迭代过程中，当其到达边界，即产生变异时，大致可分为两种情况，若二者变异方向相反，则不做任何处理，反之，将由下式对位置和速度重新定义，当粒子的位置或速度超出上边界时：vij(k+1)＝vmax-2·r3·vmax(3)xij(k+1)＝xmax-0.5·r4·xmax(4)当粒子的位置或速度超出下边界时：vij(k+1)＝vmin-2·r3·vmin(5)xij(k+1)＝xmin-0.5·r4·xmin(6)式中r3、r4是(0，1)之间的随机数。ⅴ、判断是否达到最大迭代次数，若是，获得pbest、gbest的位置信息，输出σ0、η、τ、c、d；否则返回步骤ⅲ，继续执行；5)利用提取到的实测数据的极化参数，带入深度计算公式中，绘制极化介质的电导率、极化率-深度图。包括5ⅰ、应用matlab编程，输入步骤4提取出的零频电导率、时间常数、充电率、频散系数到程序中，带入中获得复电导率值；5ⅱ、根据广义趋肤深度公式，计算该复电导率σ(ω)的趋肤深度；5ⅲ、利用步骤5ⅰ中的零频电导率、极化率和步骤5ⅱ获得的趋肤深度，分别进行电导率-深度、极化率-深度成像。实施例结合图1所示，一种基于粒子群优化算法的电导-极化率多参数成像方法，包括：1)基于地下极化介质的柯尔-柯尔复电导率模型，获得了梯形波激励下的时域感应-极化磁场响应表达式；如图1所示的步骤中，先根据图4(a)中的地质资料，建立一个三层大地模型，初始化各层电性参数，获取该层状大地模型的负阶跃波磁场响应，将梯形波发射电流与磁场响应卷积，得到梯形波下的感应磁场响应；再带入柯尔-柯尔模型表达式最后采用基于贝塞尔函数的快速hankel变换、余弦变换算法，将结果从频域转换到时域，从而获得层状大地模型下时域感应-极化双场电磁响应。2)使用超导量子传感器时域电磁探测系统接收磁场，对电磁实测数据进行处理，包括叠加取样、抽道、滤波、一次场剔除等处理；进一步地，超导量子传感器时域电磁探测系统接收到的原始数据需要进行转换及归一化，使原始数据变为有意义的电磁响应量。由于电磁勘探系统在实际测量过程中，有很多因素使得数据产生噪声，为了提高数据质量，要对噪声特点分析，并对实测数据进行噪声滤波，采用多次叠加消除漂移量。根据发射电流计算一次场其中r表达收发距，实测数据减去一次场。3)利用早期数据计算零频电导率值，晚期数据e指数多项式拟合得到时间常数，二者对粒子群搜索范围进行约束；进一步地，晚期以极化放电为主，因此利用晚期数据采用e指数拟合确定时间常数τ1＝0.0055ms，利用早期数据计算零频电导率σ1＝0.058s/m。4)根据实测磁场数据与理论模型磁场响应构建目标函数，利用粒子群优化算法进行极化参数提取，包括零频电导率、时间常数、充电率、频散系数、深度五个参数；进一步地，基于粒子群优化算法的柯尔-柯尔模型参数提取，如图2所示，具体步骤如下：ⅰ、输入粒子群优化算法相关参数d、pk、bk、xlimit、vlimit，其中d为5，bk为20，pk为300，均相等，大小均设定为2，xlimit为[0.04640.0696；0.00440.0066；01；01；-50-300]；vlimit为[-0.001,0.001；-0.001,0.001；-0.0001,0.0001；-0.002,0.002；-0.002,0.002]，依次为零频电导率、时间常数、充电率、频散系数、深度的搜索范围；随机初始化粒子的位置、速度和个体最优值pbest、种群最优值gbest；ⅱ、构建目标函数：minφ(x)＝φ1(x)+λφ2(x)，x＝σ0、η、τ、c，其中φ2(x)＝ft(x)，bt为实测数据，ft(x)为理论模型磁场响应，λ为收敛因子；ⅲ、计算每个粒子的当前适应度，并与pbest、gbest进行比较，若当前适应度值较小，则令当前适应值替代pbest、gbest；ⅳ、实测曲线依据斜率划分为两段，据此更新粒子的速度和位置：xij(k+1)＝xij(k)+vij(k+1)·t(2)式中，t0为斜率由负变为正的时刻，pk为最大迭代次数，iter为当前迭代次数，下标i表示第i个粒子(i＝1,2,3......m)，下标j表示粒子的第j维(j＝1,2,3......d)，假设粒子群搜索空间为d维，粒子群有m个粒子(通常粒子数设置为问题维数的5～10倍)，r1、r2是(0，1)之间的随机数。为了保证粒子在搜索空间内搜索，通常将粒子的速度限制在[-vmaxvmax]区间中，将步骤3得到的搜索范围做为粒子的位置限制[-xminxmax]，则粒子的最大速度可以取xmax的10％～20％。由于粒子的速度和位置均有方向性，迭代过程中，当其到达边界，即产生变异时，大致可分为两种情况，若二者变异方向相反，则不做任何处理，反之，将由下式对位置和速度重新定义，当粒子的位置或速度超出上边界时：vij(k+1)＝vmax-2·r3·vmax(3)xij(k+1)＝xmax-0.5·r4·xmax(4)当粒子的位置或速度超出下边界时：vij(k+1)＝vmin-2·r3·vmin(5)xij(k+1)＝xmin-0.5·r4·xmin(6)式中r3、r4是(0，1)之间的随机数。ⅴ、判断是否达到最大迭代次数，若是，获得pbest、gbest的位置信息，输出σ0＝0.051s/m，τ＝0.0050ms，η＝0.791，c＝1，d＝-107m；否则返回步骤ⅲ，继续执行；5)进行极化介质的电导率、极化率-深度成像。参见图1所示，利用提取到的极化参数，带入深度计算公式，绘制极化介质的电导率、极化率-深度图像。ⅰ、应用matlab编程，输入步骤4提取出的零频电导率、时间常数、充电率、频散系数到程序中，根据获得复电导率值；ⅱ、将复电导率带入广义趋肤深度公式中，可得令ic＝a+ib，可得其中ε为介电常数，ω为角频率，σ∞为高频交流电导率，μ0是磁导率，计算该复电导率σ(ω)的趋肤深度；ⅲ、利用步骤ⅰ中的零频电导率、极化率和步骤ⅱ获得的趋肤深度，分别进行电导率-深度、极化率-深度成像。图2是粒子群优化算法提取参数的流程图；图3是根据地质资料建立的三层层状大地模型；图4为采用图1所示的本发明一个实施例的层状模型成像对比图，其中a是提取结果的零频电导率-深度成像图，b是提取结果的极化率-深度成像图。结果表明，本方法的电导-极化率深度的成像结果，与实际大地一致，相比传统电阻率计算方法，更符合实际地下极化介质的传播扩散规律，提高了电导率探测深度的解释精度，为电磁勘探方法野外高精度解释提供了新的思路和方法。以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。当前第1页12当前第1页12