本发明属于工程、水文及环境地球物理勘探技术领域，尤其涉及一种接触式激发极化法有限单元数值模拟方法。
背景技术：
目前，充电法是利用天然的或者人工揭露的良导体露头、地下水出露点，直接接上供电电极(一般为正极)，同时将另外一个供电电极置于满足“无穷远”要求的位置，通过两个测量电极观测充电电场变化及分布规律，推测良导体的物性及空间分布。由于充电介质范围与观测电位异常范围形状较为相似，通常可以根据观测电位异常推断实际充电介质范围，常应用于金属勘探详查和勘探阶段、水文地质工程地质调查中地下流体追索等勘探领域。然而，传统充电法原理基于稳定电流场中理想导体(电阻率为零)的等位体前提，如金属矿体或高矿化度地下水，相对周围岩石电阻率很低，可近似看成是理想导体，通常也将非理想导体(不等位体或有损耗介质)近似为理想导体。但实际应用中充电介质的电阻率通常不为零，即实际应用均涉及有损耗介质目标体，而目前尚缺乏一种有效针对有损耗介质接触式激发极化法有限单元数值模拟的方法。
在实践勘探应用中，若充电体的规模较小或其中心埋深较大时，其充电电场与位于中心的点电源电场较为相近，因而可以使用点电源场电位曲线分布进行中心埋深推断。当充电体规模较大或者其中心埋深较浅，地表观测的充电电场畸变将明显不同于地下点电源电场分布。通常，针对上述情况，传统方法是根据电位曲线的极大值与充电点在地面投影位置不重合关系，判断为非理想导体(不等位体)。而针对不等位体的范围推断，则是建议不同充电点位置供电，配合其他物探方法综合判断不等位体的范围。因而，发展针对任意形状有损耗介质有效的接触式激发极化法有限单元数值模拟方法，为实际勘探应用提供复杂形状有损耗介质充电激发极化法理论支撑，具有非常重要的现实意义。
当前工程、水文及环境地球物理勘探技术领域还未见有针对任意形状有损耗介质接触式激发极化法有限单元数值模拟的研究报道，因此迫切需要开发一种创新性针对任意形状有损耗介质接触式激发极化法有限单元数值模拟方法，用以对工程、水文及环境地球物理勘探涉及的复杂有损耗介质电源条件下目标体分布探测及研究工作。
通过上述分析，现有技术存在的问题及缺陷为：
现有接触式激发极化法实际勘探涉及的复杂有损耗介质缺乏有效数值模拟方法，导致方法理论适用性差。
技术实现要素：
针对现有技术存在的问题，本发明提供了一种接触式激发极化法有限单元数值模拟方法。
本发明是这样实现的，一种接触式激发极化法有限单元数值模拟方法包括：
步骤一，任意形状有损耗介质区域剖分：针对任意形状三维有损耗介质区域，采用四面体交错网格进行剖分，将接触式充电三维区域离散为多个四面体网格；
步骤二，有限单元法数值方程求解：建立接触式激发极化法电位满足的方程组，求解方程组获取剖分网格总电位分布；
步骤三，接触式极化率计算：以每个剖分网格等效电阻率换算相对应的极化率值，再次进行所有剖分网格单位接触式激发极化有限单元法求解计算，获取考虑接触式激发极化的情况下总电位分布；
步骤四，接触式激发极化法视极化率计算：根据视极化率求取公式计算各个地面测点的视极化率值。
进一步，所述步骤一中，所述任意形状有损耗介质区域剖分具体包括：
采用笛卡尔直角坐标系，取坐标系z轴垂直向上，根据右手螺旋法则确定x轴和y轴方向；确定包含任意形状三维有损耗介质和围岩介质计算区域在坐标轴x、y、z方向规则空间分布范围，采用规则长方体网格将上述规则计算区域在x方向上按一定间距剖分为m个网格节点，其y方向上按一定间距剖分为有n个网格节点，及其z方向上按一定间距剖分为有q个网格节点，再针对每个规则长方体网格进行交错四面体网格剖分。
进一步，坐标系三轴方向剖分网格节点长度允许各不相同，其间距剖分长度d应满足d≤l/10，其中，l为任意形状三维有损耗介质在各坐标轴方向的展布长度，如此，包含任意形状三维有损耗介质和围岩介质计算区域被剖分为(m-1)×(n-1)×(q-1)×5个网格单元；
按任意形状三维有损耗介质在坐标系内的空间分布，确定其在(m-1)×(n-1)×(q-1)×5个网格单元中对应自身空间分布的nm个网格单元，nm个任意形状三维有损耗介质剖分网格单元对应nm个单元点电源。
进一步，任意形状三维有损耗介质在坐标轴x方向上分布节点坐标为xi，y方向上分布节点坐标为yj，z方向上分布节点坐标为zk，其中，i＝1,…,mx；j＝1,…,ny；z＝1,…,qz，mx、ny及qz分别为x、y及z方向分布节点数，任意形状三维有损耗介质剖分单元网格中心坐标为(xi dx,i/2,yj dy,j/2,zk dz,k/2)，dx,i、dy,j及dz,k为坐标轴x、y及z方向第i、j及k个剖分网格间距；
充电点位置s坐标为(xs,ys,zs)，供电点供电电流强度为itotal，任意形状三维有损耗介质区域剖分电阻率可允许随剖分单元网格位置变化而变化，记为ρq，q＝1,…,nm；ρb为围岩介质电阻率；pi＝3.1415926；剖分单位网格视为点电源，电流强度记为iq，q＝1,…,nm，地面观测点为坐标为(xa,ya,za)，a＝1,…,a；a为观测点数量；u为剖分单位网格节点接触式供电电源条件下引起的总电位值，即待求解任意形状三维有损耗介质在计算区域所有网格节点的总电位值。
进一步，所述步骤二中，所述有限单元法数值方程求解具体包括：
根据公式(1)计算在供电电流强度itotal充电时，任意形状三维有损耗介质区域剖分单元网格点电源电流强度系数iq，q＝1,…,nm：
其中，rs,q为充电点到第q个剖分单元网格点电源中心的距离；
其中ωm为区域ω对电流源im点所张的立体角，γs和γ∞分别为区域ω的地面边界和地下无穷远边界，n为边界的法向向量方向，δ为狄拉克函数，r1和r2为计算点至供电点a/b的空间距离，
通过采用有限单元法的线性插值、单元积分、边界积分，最终总体合成、再扩展成由全体节点组成的矩阵，进而全部单元相加，得
其中，u是全部节点的u组成的列向量；ke＝1/ρ(k1e k2e)，
令(3)式的变分为0，得线性方程组：
ku＝p(4)
解方程组，得到各节点的电位。
进一步，所述步骤三中，所述接触式极化率计算具体包括：
根据步骤二所计算得到的任意形状三维有损耗介质区域剖分单元网格点总电位，按等效电阻率公式(5)，将考虑极化率模型的激化效应换算为各剖分单元电阻率值：
其中，η为模型极化率值；
进一步，所述步骤四，所述接触式激发极化法视极化率计算具体包括：
根据步骤三得到考虑极化率模型的所有剖分单元网格点总电位u2，根据极限极化率求取公式(6)进行换算接触式激发极化法地面观测的视极化率值：
其中，δu2和δu分别为地面相邻两个网格节点针对考虑极化率模型、不考虑极化率模型时有限单元法求解的总电位，ηs对应为上述两个相邻网格节点中间点的视极化率值。
本发明的另一目的在于提供一种接触式激发极化法有限单元数值模拟系统，所述接触式激发极化法有限单元数值模拟系统包括：
有损耗介质区域剖分模块，用于对任意形状三维有损耗介质区域，采用四面体交错网格进行剖分，将接触式充电三维区域离散为多个四面体网格；
有限单元法数值方程求解模块，用用于建立接触式激发极化法电位满足的方程组，求解方程组获取剖分网格总电位分布；
接触式极化率计算模块，用于以每个剖分网格等效电阻率换算相对应的极化率值，再次进行所有剖分网格单位接触式激发极化有限单元法求解计算，获取考虑接触式激发极化的情况下总电位分布；
接触式激发极化法视极化率计算模块，用于根据视极化率求取公式计算各个地面测点的视极化率值。
本发明的另一目的在于提供一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如下步骤：
步骤一，任意形状有损耗介质区域剖分：针对任意形状三维有损耗介质区域，采用四面体交错网格进行剖分，将接触式充电三维区域离散为多个四面体网格；
步骤二，有限单元法数值方程求解：建立接触式激发极化法电位满足的方程组，求解方程组获取剖分网格总电位分布；
步骤三，接触式极化率计算：以每个剖分网格等效电阻率换算相对应的极化率值，再次进行所有剖分网格单位接触式激发极化有限单元法求解计算，获取考虑接触式激发极化的情况下总电位分布；
步骤四，接触式激发极化法视极化率计算：根据视极化率求取公式计算各个地面测点的视极化率值。
本发明的另一目的在于提供一种信息数据处理终端，所述信息数据处理终端包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行所述的接触式激发极化法有限单元数值模拟方法。
结合上述的所有技术方案，本发明所具备的优点及积极效果为：
与现有技术相比，针对地球物理勘探领域复杂有损耗介质开发接触式激发极化法数值模拟方法是本专利发明点之一，它为涉及复杂有损耗介质电源地球物理勘探领域提供了一种有效的数值模拟方法。
本发明是一个新型数值模拟算法，根据步骤二所设计的一种特殊的充电电流强度满足的有损耗介质剖分网格有限单元法边值问题，实现了针对地球物理勘探领域复杂有损耗介质开发接触式激发极化法数值模拟在可行性及精度上的统一。
本发明解决了传统任意形状有损耗介质接触式激发极化法有限单元数值模拟计算方法的难题，算法结构简单，实现过程简便，设计合理，计算结果稳定可靠，可广泛应用于工程、水文及环境地球物理勘探领域涉及的任意形状有损耗介质充电法勘探极化率数值模拟计算，为基于复杂有损耗介质电源的工程、水文及环境地球物理勘探数值模拟、反演及解释等提供了方法支撑。
附图说明
为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图做简单的介绍，显而易见地，下面所描述的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下还可以根据这些附图获得其他的附图。
图1是本发明实施例提供的接触式激发极化法有限单元数值模拟方法的计算流程示意图。
图2是本发明实施例提供的设计模型三维示意图。
图3是本发明实施例提供的设计模型三维网格剖分图。
图4是本发明实施例提供的设计模型及观测系统示意图。
图5是本发明实施例提供的设计模型接触式激发极化法数值模拟计算结果(numerical_new_ground)与物理模拟结果(physical_new_ground)、传统激发极化法(physical_old_ground)对比图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
针对现有技术存在的问题，本发明提供了一种接触式激发极化法有限单元数值模拟方法，下面结合附图对本发明作详细的描述。
如图1所示，本发明实施例提供的接触式激发极化法有限单元数值模拟方法，步骤为：
步骤一：任意形状有损耗介质区域剖分：
参照图2，采用笛卡尔直角坐标系，取坐标系z轴垂直向上，则根据右手螺旋法则确定x轴和y轴方向。参照图3，确定包含任意形状三维有损耗介质和围岩介质计算区域在坐标轴x、y、z方向规则空间分布范围为0～100cm、0～120cm、0～30cm，采用规则长方体网格将上述规则计算区域在x方向上按一定间距剖分为m＝58个网格节点，其y方向上按一定间距剖分为有n＝58个网格节点，及其z方向上按一定间距剖分为有q＝17个网格节点，再针对每个规则长方体网格进行交错四面体网格剖分。坐标系三轴方向剖分网格节点长度允许各不相同，其间距剖分长度d应满足d≤l/10，其中，l为任意形状三维有损耗介质在各坐标轴方向的展布长度，如此，包含任意形状三维有损耗介质和围岩介质计算区域被剖分为(m-1)×(n-1)×(q-1)×5个网格单元。
按任意形状三维有损耗介质在坐标系内的空间分布，确定其在(m-1)×(n-1)×(q-1)×5个网格单元中对应自身空间分布的nm＝57188个网格单元，nm个任意形状三维有损耗介质剖分网格单元对应nm个单元点电源。
任意形状三维有损耗介质在坐标轴x方向上分布节点坐标为xi，y方向上分布节点坐标为yj，z方向上分布节点坐标为zk，其中，i＝1,…,mx；j＝1,…,ny；z＝1,…,qz，mx、ny及qz分别为x、y及z方向分布节点数。任意形状三维有损耗介质剖分单元网格中心坐标为(xi dx,i/2,yj dy,j/2,zk dz,k/2)，dx,i、dy,j及dz,k为坐标轴x、y及z方向第i、j及k个剖分网格间距。
充电点位置s坐标为(xs,ys,zs)，供电点供电电流强度为itotal＝1安培，参照图4，任意形状三维有损耗介质区域剖分电阻率可允许随剖分单元网格位置变化而变化，记为ρq＝15ω.m，q＝1,…,nm；ρb＝50ω.m为围岩介质电阻率；pi＝3.1415926；剖分单位网格视为点电源，电流强度记为iq，q＝1,…,nm，地面观测点为坐标为(xa,ya,za)，a＝1,…,a；a＝30为观测点数量；u为剖分单位网格节点接触式供电电源条件下引起的总电位值，即步骤三待求解任意形状三维有损耗介质在计算区域所有网格节点的总电位值。
步骤二：有限单元法数值方程求解：
根据公式(1)计算在供电电流强度itotal充电时，任意形状三维有损耗介质区域剖分单元网格点电源电流强度系数iq，q＝1,…,nm：
其中，rs,q为充电点到第q个剖分单元网格点电源中心的距离。
其中ωm为区域ω对电流源im点所张的立体角，γs和γ∞分别为区域ω的地面边界和地下无穷远边界，n为边界的法向向量方向，δ为狄拉克函数，r1和r2为计算点至供电点a/b的空间距离，
通过采用有限单元法的线性插值、单元积分、边界积分，最终总体合成、再扩展成由全体节点组成的矩阵，进而全部单元相加，得
其中，u是全部节点的u组成的列向量；ke＝1/ρ(k1e k2e)，
令(3)式的变分为0，得线性方程组：
ku＝p(4)
解方程组，得到各节点的电位。
步骤三：接触式极化率计算
根据步骤二所计算得到的任意形状三维有损耗介质区域剖分单元网格点总电位，按等效电阻率公式(5)，将考虑极化率模型的激化效应换算为各剖分单元电阻率值：
其中，η为模型极化率值为6.8％，背景介质极化率值为1％。再根据公式(5)求解得到考虑极化率模型的所有剖分单元网格点总电位u2，根据极限极化率求取公式(6)进行换算接触式激发极化法地面观测的视极化率值：
其中，δu2和δu分别为地面相邻两个网格节点针对考虑极化率模型、不考虑极化率模型时有限单元法求解的总电位，ηs对应为上述两个相邻网格节点中间点的视极化率值。
本发明还提供一种接触式激发极化法有限单元数值模拟系统，所述接触式激发极化法有限单元数值模拟系统包括：
有损耗介质区域剖分模块，用于对任意形状三维有损耗介质区域，采用四面体交错网格进行剖分，将接触式充电三维区域离散为多个四面体网格；
有限单元法数值方程求解模块，用用于建立接触式激发极化法电位满足的方程组，求解方程组获取剖分网格总电位分布；
接触式极化率计算模块，用于以每个剖分网格等效电阻率换算相对应的极化率值，再次进行所有剖分网格单位接触式激发极化有限单元法求解计算，获取考虑接触式激发极化的情况下总电位分布；
接触式激发极化法视极化率计算模块，用于根据视极化率求取公式计算各个地面测点的视极化率值。
将本发明提供的针对土槽有损耗介质接触式激发极化法数值模拟计算的视极化率曲线及相同模型(图2、图3、图4)土槽物理模拟观测结果。参照图5所示为本发明算例设计模型接触式激发极化法数值模拟计算结果(numerical_new_ground)与物理模拟结果(physical_new_ground)、传统激发极化法(physical_old_ground)对比图。图中，cm表示厘米，chargeability表示极化率。对比可见，从形态上看，本发明计算结果与物理模拟观测结果吻合较好，仅在模型边界附近存在微弱差异，这是由于物理模拟环境与数值模拟环境不能完全一致所导致。以两种数据的拟合相对误差衡量计算结果吻合程度，图5所示电位曲线拟合相对误差为2％，同时，接触式激发极化法观测结果相比传统激发极化法异常幅值高出许多，验证了本发明实施算例的准确性及可行性。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。
1.一种接触式激发极化法有限单元数值模拟方法，其特征在于，所述接触式激发极化法有限单元数值模拟方法包括：
步骤一，任意形状有损耗介质区域剖分：针对任意形状三维有损耗介质区域，采用四面体交错网格进行剖分，将接触式充电三维区域离散为多个四面体网格；
步骤二，有限单元法数值方程求解：建立接触式激发极化法电位满足的方程组，求解方程组获取剖分网格总电位分布；
步骤三，接触式极化率计算：以每个剖分网格等效电阻率换算相对应的极化率值，再次进行所有剖分网格单位接触式激发极化有限单元法求解计算，获取考虑接触式激发极化的情况下总电位分布；
步骤四，接触式激发极化法视极化率计算：根据视极化率求取公式计算各个地面测点的视极化率值。
2.如权利要求1所述的接触式激发极化法有限单元数值模拟方法，其特征在于，所述步骤一中，所述任意形状有损耗介质区域剖分具体包括：
采用笛卡尔直角坐标系，取坐标系z轴垂直向上，根据右手螺旋法则确定x轴和y轴方向；确定包含任意形状三维有损耗介质和围岩介质计算区域在坐标轴x、y、z方向规则空间分布范围，采用规则长方体网格将上述规则计算区域在x方向上按一定间距剖分为m个网格节点，其y方向上按一定间距剖分为有n个网格节点，及其z方向上按一定间距剖分为有q个网格节点，再针对每个规则长方体网格进行交错四面体网格剖分。
3.如权利要求2所述的接触式激发极化法有限单元数值模拟方法，其特征在于，坐标系三轴方向剖分网格节点长度允许各不相同，其间距剖分长度d应满足d≤l/10，其中，l为任意形状三维有损耗介质在各坐标轴方向的展布长度，如此，包含任意形状三维有损耗介质和围岩介质计算区域被剖分为(m-1)×(n-1)×(q-1)×5个网格单元；
按任意形状三维有损耗介质在坐标系内的空间分布，确定其在(m-1)×(n-1)×(q-1)×5个网格单元中对应自身空间分布的nm个网格单元，nm个任意形状三维有损耗介质剖分网格单元对应nm个单元点电源。
4.如权利要求3所述的接触式激发极化法有限单元数值模拟方法，其特征在于，任意形状三维有损耗介质在坐标轴x方向上分布节点坐标为xi，y方向上分布节点坐标为yj，z方向上分布节点坐标为zk，其中，i＝1,…,mx；j＝1,…,ny；z＝1,…,qz，mx、ny及qz分别为x、y及z方向分布节点数，任意形状三维有损耗介质剖分单元网格中心坐标为(xi dx,i/2,yj dy,j/2,zk dz,k/2)，dx,i、dy,j及dz,k为坐标轴x、y及z方向第i、j及k个剖分网格间距；
充电点位置s坐标为(xs,ys,zs)，供电点供电电流强度为itotal，任意形状三维有损耗介质区域剖分电阻率可允许随剖分单元网格位置变化而变化，记为ρq，q＝1,…,nm；ρb为围岩介质电阻率；pi＝3.1415926；剖分单位网格视为点电源，电流强度记为iq，q＝1,…,nm，地面观测点为坐标为(xa,ya,za)，a＝1,…,a；a为观测点数量；u为剖分单位网格节点接触式供电电源条件下引起的总电位值，即待求解任意形状三维有损耗介质在计算区域所有网格节点的总电位值。
5.如权利要求1所述的接触式激发极化法有限单元数值模拟方法，其特征在于，所述步骤二中，所述有限单元法数值方程求解具体包括：
根据公式(1)计算在供电电流强度itotal充电时，任意形状三维有损耗介质区域剖分单元网格点电源电流强度系数iq，q＝1,…,nm：
其中，rs,q为充电点到第q个剖分单元网格点电源中心的距离；
其中ωm为区域ω对电流源im点所张的立体角，γs和γ∞分别为区域ω的地面边界和地下无穷远边界，n为边界的法向向量方向，δ为狄拉克函数，r1和r2为计算点至供电点a、b的空间距离，
通过采用有限单元法的线性插值、单元积分、边界积分，最终总体合成、再扩展成由全体节点组成的矩阵，进而全部单元相加，得
其中，u是全部节点的u组成的列向量；ke＝1/ρ(k1e k2e)，
令(3)式的变分为0，得线性方程组：
ku＝p(4)
解方程组，得到各节点的电位。
6.如权利要求1所述的接触式激发极化法有限单元数值模拟方法，其特征在于，所述步骤三中，所述接触式极化率计算具体包括：
根据步骤二所计算得到的任意形状三维有损耗介质区域剖分单元网格点总电位，按等效电阻率公式(5)，将考虑极化率模型的激化效应换算为各剖分单元电阻率值：
其中，η为模型极化率值。
7.如权利要求1所述的接触式激发极化法有限单元数值模拟方法，其特征在于，所述步骤四，所述接触式激发极化法视极化率计算具体包括：
根据步骤三得到考虑极化率模型的所有剖分单元网格点总电位u2，根据极限极化率求取公式(6)进行换算接触式激发极化法地面观测的视极化率值：
其中，δu2和δu分别为地面相邻两个网格节点针对考虑极化率模型、不考虑极化率模型时有限单元法求解的总电位，ηs对应为上述两个相邻网格节点中间点的视极化率值。
8.一种接触式激发极化法有限单元数值模拟系统，其特征在于，所述接触式激发极化法有限单元数值模拟系统包括：
有损耗介质区域剖分模块，用于对任意形状三维有损耗介质区域，采用四面体交错网格进行剖分，将接触式充电三维区域离散为多个四面体网格；
有限单元法数值方程求解模块，用用于建立接触式激发极化法电位满足的方程组，求解方程组获取剖分网格总电位分布；
接触式极化率计算模块，用于以每个剖分网格等效电阻率换算相对应的极化率值，再次进行所有剖分网格单位接触式激发极化有限单元法求解计算，获取考虑接触式激发极化的情况下总电位分布；
接触式激发极化法视极化率计算模块，用于根据视极化率求取公式计算各个地面测点的视极化率值。
9.一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如下步骤：
步骤一，任意形状有损耗介质区域剖分：针对任意形状三维有损耗介质区域，采用四面体交错网格进行剖分，将接触式充电三维区域离散为多个四面体网格；
步骤二，有限单元法数值方程求解：建立接触式激发极化法电位满足的方程组，求解方程组获取剖分网格总电位分布；
步骤三，接触式极化率计算：以每个剖分网格等效电阻率换算相对应的极化率值，再次进行所有剖分网格单位接触式激发极化有限单元法求解计算，获取考虑接触式激发极化的情况下总电位分布；
步骤四，接触式激发极化法视极化率计算：根据视极化率求取公式计算各个地面测点的视极化率值。
10.一种信息数据处理终端，其特征在于，所述信息数据处理终端包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行权利要求1～7任意一项所述的接触式激发极化法有限单元数值模拟方法。
技术总结