无限长圆柱体直流电场求解的最优化波数选取及应用方法与流程

本发明涉及环境保护与监测技术领域，尤其涉及一种无限长圆柱体直流电场求解的最优化波数选取及应用方法。

背景技术：

圆柱体截面上的电阻率成像广泛应用于树干木质评价、营养液传输监测等许多领域，其2.5d电位分布的正演计算是实现电阻率成像的基础。

目前，2.5d直流电阻率正演计算方法多是针对平坦地形或有一定地形起伏的地下二维半空间模型，而无限长圆柱体模型直流电场的求解域是一个封闭的圆，当前有许多研究者直接利用传统的2.5d直流电阻率正演方法计算无限长圆柱体模型上的电位分布，然而，这种传统的2.5d数值模拟方法计算得到的无限长圆柱体模型上电位分布具有较大误差，其主要原因就在于波数及对应步长选取的不合理。

技术实现要素：

本发明的目的是提供一种无限长圆柱体直流电场求解的最优化波数选取及应用方法，可以提高无限长圆柱体模型上电位分布的计算精度。

本发明的目的是通过以下技术方案实现的：

一种无限长圆柱体直流电场求解的最优化波数选取方法，包括：利用数值模拟计算方法得到无限长圆柱体模型横截面上电位的近似解析解；

将无限长圆柱体模型横截面上电位的近似解析解，代替平地条件下的解析解，然后利用最优化计算方法，来获得最优的波数序列及对应的步长序列。

一种无限长圆柱体直流电场求解的最优化波数选取及应用方法，其基于前述方法来获得最优的波数序列及对应的步长序列，并基于最优的波数序列及对应的步长序列，进行无限长圆柱体模型电阻率成像正演模拟或反演模拟，获得正演视电阻率图或反演电阻率成像结果图。

由上述本发明提供的技术方案可以看出，能够有效提高无限长圆柱体模型上直流电场求解的正演计算精度，继而进一步提高反演精度，促进电阻率成像在环境保护、环境监测等多个领域的应用效果。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例提供的一种无限长圆柱体模型直流电场求解的最优化波数选取方法的流程图；

图2为本发明实施例提供的无限长圆柱体模型的横截面示意图；

图3为本发明实施例提供的改进后的最优化波数方法计算结果与一般方法以及实验结果对比图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

本发明实施例提供一种无限长圆柱体模型直流电场求解的最优化波数选取方法，目的在于提高无限长圆柱体模型上电位分布的计算精度，该方法可通过常规方式编程，最终以计算机程序的方式来实现。

如图1所示，该方法主要包括如下步骤：

1、利用数值模拟计算方法得到无限长圆柱体模型横截面上电位的近似解析解。

在介绍无限长圆柱体模型横截面上电位的近似解析解之前，首先针对无限长圆柱体模型进行直流点电源电场求解原理进行介绍。

在对无限长圆柱体模型进行直流点电源电场求解时，由于模型为二度体，而源是三维的点电源，因此模型的场实际是三维分布的，称为二度半(2.5d)问题，通常变换到波数域进行二维求解。

本发明实施例中，将无限长圆柱体模型所满足的电位方程组进行傅立叶变换，其公式为：

其中，u(x,y,z)表示电位值，v(x,λ,z)是波数域中电位的二维分布，λ是波数，(x,y,z)是测量点的三维位置坐标，(x,z)是测量点的二维位置坐标。

本发明实施例中，利用有限单元法或积分方程法解傅立叶变换后的电位方程组所满足的拉普拉斯方程，获得不同波数情况下的v(x,λ,z)值；

再进行反傅立叶变换：

从而得到测量点(x,z)上的电位值。但是反傅立叶变换不是一个直接求解的过程，因为电位值是离散的，所以不能直接对其进行积分，为了求解电位值，将上式改写为：

其中，gj表示第j个波数λj对应的步长因子，是一个未知量，为了求出gj，必须要知道均匀无限长圆柱体上电位值的解析解。

本发明实施例中，所述无限长圆柱体模型横截面上电位的近似解析解为：

其中，i表示电流强度，单位为安培(a)，ρ表示电阻率值，单位为欧姆·米(ω·m)，r表示无限长圆柱体模型的半径，单位为米(m)，fθab(θ)是一系列通过正演模拟计算得到的经验值，θ表示测量点与供电电极之间的夹角。

如图2所示，为无限长圆柱体模型的横截面示意图。图2中的n与m为测量点，a与b为供电电极，θab为供电电极a与供电电极b之间的夹角，θm、θn分别为测量点n、测量点m与供电电极a之间的夹角，可统称为夹角θ。取较多的波数进行正演计算得到不同波数λ对应的v(λ,x,z)，对这计算到的这些v(x,λ,z)进行拟合，将得到的函数表达式进行积分即可得到空间域的电位值u0(θ)，就等于u0(θ)r/(iρ)，它与夹角θab以及夹角θ有关，对其进行一次计算后可存储起来，之后的计算只需要通过查表即可。

2、将无限长圆柱体模型横截面上电位的近似解析解，代替平地条件下的解析解，然后利用最优化计算方法，来获得最优的波数序列及对应的步长序列。

综合反傅立叶变换结果与近似解析解得到：

其中，x＝r·cos(θ)，z＝r·sin(θ)，对于任意的测量点(xi,zi)或者(r·cos(θi),r·sin(θi))都有：

上式的矩阵形式为：

ag≈i

其中，a是由组成的n×m阶矩阵，n为参加计算的测量点总个数，m为所使用的波数个数，g是由gj组成的m维列向量，i为n维单位列向量。

本发明实施例中，最优化计算方法中，首先计算波数序列对应的最优步长序列，再判断相应波数序列是否为最优的波数序列，如果是，则计算结束，获得最终结果；否则，需要对波数序列进行优化；主要过程如下：

1)计算波数序列对应的最优步长序列，并判断波数序列是否为最优的波数序列。

对于一组波数序列λ^k，能够计算出相应的a；k表示迭代次数，k＝0，1,2,..，在初次计算时，k＝0，λ⁰表示给定的一组随机的波数序列；

引入限定条件：i-ag的二范数最小，得到：

dg^t不为零，则有：

a^tag＝a^ti

即：

g＝(a^ta)^-1(a^ti)

通过求解上述方程，对于一组波数序列λ^k，能够得到一组最优步长序列g^k。

再利用来计算二范数的大小，如果二范数足够小(即二范数满足要求)，则将波数序列λ^k作为最优的波数序列，g^k作为最优步长序列；否则，转入步骤2)。

2)对波数序列进行优化。

本发明实施例中，需要对波数序列进行优化，使得二范数达到全局最小，记v＝ag，将v进行关于变量λj的泰勒级数展开：

上式中，δ为数学中的变分符号；v0通过一组波数序列λ^k及对应的最优步长序列g^k计算得到；

通过上式得到：

二范数是关于δλ的函数，要使二范数取极值，则有：

mδλ＝h

上式中，

通过求解上式得到δλ的值，将求得的δλ值代入下式，得到一组新的波数序列λ^k+1，并返回步骤1)：

λ^k+1＝λ^k+δλ，k＝0，1,2,...。

另一方面，本发明实施例还提供一种基于上述方案所得到的最优的波数序列及对应的步长序列进行相关应用的方法。该方法中，基于最优的波数序列及对应的步长序列，进行无限长圆柱体模型电阻率成像正演模拟或反演模拟，从而获得正演视电阻率图或反演电阻率成像结果图。该方案可以提高无限长圆柱体模型的正演计算精度，最终提高反演计算精度，得到更加准确的反演电阻率成像。例如，可以应用于树干木质评价、营养液传输监测等方面。

需要说明的是，与平地条件下不同，无限长圆柱体的最优化波数序列并不是固定不变的，而是随着供电电极的夹角θab和圆柱体的半径的变化而变化的，如果采取不变的波数序列对不同的供电点角度的所有情况进行计算，那么误差还会进一步增大。另外，为了避免半径对波数选取的影响，可以先计算半径为1米的圆柱体模型上的电位值，然后利用电位值与半径成反比这一规律将其变换到目标半径值对应的电位分布，即在进行数值模拟时将模型缩放到1米半径后代入最优化波数序列进行正演计算，这样就可以直接利用半径为1米时对应的最优化波数序列，而不需要单独对某一半径条件下对应的最优化波数序列进行计算，大大节省了计算时间。

与传统方法相比，本发明所提出的最优化波数选取方案更适用于圆柱体模型的2.5d电阻率成像正演模拟，能够有效提高圆柱体模型上电位值的正演计算精度，为了验证本发明上述方案的性能，还进行一组实验，对所测的数据进行两种方法的正演模拟，在其他条件相同的情况下，其中一组正演模拟采用本发明所提出的最优化波数选取方法进行波数计算，而另外一组的正演模拟采用传统的最优化波数方法进行波数计算。

主要实验过程如下：

(1)模型制作：首先，用胶水将pvc管与底部的塑料板黏结起来；其次，测量pvc管的长度、直径及周长，在pvc管纵向中心处的平面上等距离均匀标记上24个点作为电极位置，用电钻在这些标记的位置上打孔，电钻最小直径2mm。

(2)微电极制作：将铜导线剪短成长度约为2cm的细铜丝，每两根合并为一个微电极，插入到pvc管上的小孔中，用速干胶水将孔洞封闭，剪去管内过长的微电极使之等价于点电极，在微电极附近标上序号。

(3)制备导体水：水的电阻率与nacl浓度的关系可以通过查表获得，为了配备电阻率为50欧姆米的导体水溶液，需向30l水中加入2.85gnacl晶体。

(4)测量前的准备工作：将制备好的nacl水溶液倒入到模型中，用四根导线分别将供电电极和测量电极接入duk-2b高密度电法测量系统，将电池于主机连接起来，打开仪器设备。

(5)对三种不同模型分别进行温纳装置和偶极装置的测量，测量是逐点进行的，每测量一次都需要人工进行移动电极及记录数据，接入的电压为120v。

(6)测量完毕之后需要对水进行电阻率的测量，因为在测量过程中水会蒸发，并且所购买的蒸馏水中本身就存在一定量的杂质，使得水的电阻率要比配备的理论值偏低。对电阻率的测量是通过将水溶液注入到圆柱体模型中，圆柱体模型上下均粘上铜片接入电路，在保证上下底面与水溶液接触良好的情况下对圆柱体两端的电压和电流值进行测量，然后利用欧姆定律计算圆柱体的电阻，进而利用圆柱体的半径和长度计算水溶液的电阻率。

实验结果如图3所示。左上部分为实验的实测结果(实测电阻数据)，右上部分为平地条件下的最优化波数序列计算得到的电阻率值，右下部分为改进后的波数选取方法(即本发明方案)计算得到的电阻值，各部分色标相同，不难发现，平地条件下的最优化波数选取方法所计算得到的结果与实验数据相差较大，而改进后的方法计算得到的电阻率值与实测数据更为接近，不难看出，改进后的计算方法能够大大提高计算的精度。注意到改进后的方法计算得到的电阻值比实测数据稍大，这是由于圆周上的微电极的影响使得所测数据偏低。值得注意的是，在计算过程中，有一很有意义的发现，因为改进后的最优化波数序列是在r＝1m时计算得到的，当直接将模型半径设置为99.75mm时，利用改进后的波数方法计算得到的电阻值如图3左下部分所示，计算的误差比较大，这是因为，对于圆柱体而言，其最优化波数序列不仅仅跟供电电极的夹角有关，还与圆柱体的半径有关，而在野外条件下，树木的半径各不相同，在应用过程中需要不断重新对最优化波数进行计算，这样非常耗费计算时间，为了解决这一问题的方案是：在进行正演计算时，首先将模型半径统一设置为1m，这样需要使用的最优化波数序列就只与供电夹角有关而与半径无关，在计算得到半径为1m条件下的电位值分布之后，利用电位值与半径成反比的关系，将半径等距离缩放，同时也将对应的电位值等距离缩放，这样就可以计算任一半径条件下的电位值分布而避免了半径对最优化波数值选取的影响。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例可以通过软件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是cd-rom，u盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。