具体实施方式

一种三维非对称结构土壤模型下基于场路耦合的直流偏磁计算方法，通过需要计及的直流偏磁影响范围确定土壤模型的尺寸，将模型离散为边长2m的小立方体；通过实际的测点数据，确定模型的区块划分；将反演的三维电阻率映射至各区块，非三维电阻率数据通过转换后映射；通过注流试验数据确定注流点坐标，施加激励，确定边界条件，划分网格，进行地表电位计算；确定观测路径，与注流试验结果进行对比，修正土壤模型；获取接地极周边电力系统接线图，坐标等参数，搭建直流电路模型；输入各节点电势，进行直流偏磁计算。本发明从初始模型出发，提高了直流偏磁的计算精度，在直流工程投运后减少了保护误动或拒动的可能，提高了电力系统运行的稳定性。

一种三维非对称结构土壤模型下基于场路耦合的直流偏磁计算方法，包括以下步骤：

步骤1：根据需要评估的直流偏磁影响范围，确定待建立的土壤模型尺寸，并对土壤模型进行离散化操作。

进行离散化操作的目的是将整个求解区域的土壤块离散为小的土壤块，以减少求解难度。

具体操作如图1所示，根据土壤模型尺寸及求解精度，首先在x轴一端建立一个小立方体(元)，然后通过复制操作在x轴上形成一列小立方体(线)，然后再通过复制操作将这列小立方体在xy平面上形成一面域内小立方体(面)，最后通过复制操作将这一面域内的小立方体在xyz直角坐标系中形成以小立方体为基的大立方体(体)。上述过程只是离散化的一种方法，但不局限于该方法。

步骤2：根据浅层大地电阻率和深层大地电阻率测点数据，对土壤模型进行区块划分；

步骤3：读取电阻率数据，判断是否为三维电阻率数据，若不是则将其转换为三维电阻率；

步骤4：将三维电阻率数据映射到相应的土壤模型区块中；

步骤5：读取注流试验相关设置参数，确定注流点坐标及注流试验测点路径；

步骤6：确定土壤模型边界条件，施加激励，将模型离散为四面体有限元网格；

步骤7：计算地表电位分布；

步骤8：取与注流试验一致的测点路径，并将其地电位与注流试验结果对比，判断|V_数-V_注|<ΔV；

步骤9：当判断结果为否时，对土壤电阻率进行修正，并再次进行计算；

步骤10：当判断结果为是时，读取周边电力系统接线数据，确定各变压器坐标及直流电路模型，然后获取相应坐标地电位分布；

步骤11：将获取的地电位分布映射至直流电路模型中，并进行直流偏磁计算；

步骤12：输出结果。

步骤1中，根据需要评估的直流偏磁影响范围，确定待建立的土壤模型尺寸，并对土壤模型进行离散化操作，土壤模型离散化程度可以根据计算精度需求自由控制。

步骤2中，根据浅层大地电阻率和深层大地电阻率测点数据，将浅层土壤中电阻率变化较小的区域合并为一个区块，减小深层土壤的离散程度，将其划为一个较大的区块。

步骤3中，读取电阻率数据，判断是否为三维电阻率数据，若不是则将其通过线性插值或二次函数插值的方式转换为三维电阻率。

步骤4中，三维电阻率数据为：由大地电磁法或大地音频电磁法经三维反演得到的电阻率数据。

步骤8中，当判断结果为否时，对土壤电阻率进行修正，当V_数<V_注时，增大所取路径上的土壤电阻率；当V_数>V_注时，减小所取路径上的土壤电阻率，然后再次进行计算。

本发明所述的三维非对称结构下基于场路耦合的直流偏磁计算方法，其计算流程图如图2所示，主要分为数据预处理部分，几何建模部分，地表电位计算部分，直流偏磁计算部分。

本发明所述计算方法由以下步骤组成：

1)根据需要评估的直流偏磁影响范围：根据DL/T5224-2014，应对极址周围地电位升大于3V的变电站，以及与其有直接电气联系的变电站进行仿真计算。现随着直流输电技术的发展，直流接地极的入地电流也逐步增大，其产生偏磁影响的范围也随之增大，所以建议将150km范围内的变电站、电厂等纳入计算网络。

确定待建立的土壤模型尺寸：为了保证直流偏磁的计算精度，所建立的土壤模型要大于需要计算偏磁电流的计算网络。所以，土壤模型的尺寸确定为400km，如图3所示。将接地极作为点元处理，其坐标设为(0,0,0)，根据模型尺寸建立电流场计算几何模型，即为原始的土壤模型，如图3所示。其中原点(0,0,0)为直流电流注入点，除地表外其余五面为零电势面。并将模型离散为边长为2m的小立方体，离散后的土壤模型如图4所示。

2)根据浅层大地电阻率测点数据，如表1所示，即为某接地极的浅层大地电阻率数据，所有测点位置分布如图5所示。现以此为例进行浅层大地电阻率的区块划分，但不限于此种划分方式。

a，首先确定在深度方向上的区块层深，依次为2m、3m、5m、5m、……、300m，共计13层；

b，确定在地表平面上的区块个数为49个，故设置每个区块的长、宽分别为150m、150m。

综上，浅层土壤区块划分完毕。其中部分区块无对应的电阻率数据，采用对应层深的视电阻率填充。获取分层坐标(x_i,y_i,z_i)，对接地极(0,0,0)邻近区域进行区块(x_i-1,x_i,y_i-1,y_i,z_i-1,z_i)划分。

表1某接地极浅层大地电阻率数据

3)据深层大地电阻率测点数据，获取分层坐标(x,y,z)，对划分区域进行区块(x_-1,x,y_-1,y,z_-1,z)划分。表2某地大地深层电阻率数据，进行深层大地电阻率的区块划分，划分。1)的400km，层区块的分400km 400km0.1km，2)进行浅层大地电阻率区块划分的区域接层区块的分400km400km0.9km，2)进行浅层大地电阻率区块划分的区域层区块的分400km400km2km，层的分400km400km310km。，深层区块划分。

表2某接地极深层大地电阻率数据

4)读取电阻率数据，判断是否为三维电阻率数据。如表3所示，为某局部区域的三维电阻率数据，表3中层厚即为该土壤层在深度方向的厚度，x、y、z为指定的直角坐标系对应的三个坐标轴。实际的大地因为土壤电阻率各向异性特性所体现出的电性在各方向上并不一致，如表3所示，各轴上电阻率的值及变化各不相同，且同一轴的正负半轴其值及变化也不相同，这体现了所建模型的非对称特性。若不是则将其通过线性插值或二次函数插值转换为三维电阻率。

表3某三维电阻率数据

上述表格为本发明所建议采用的三维电阻率数据格式。实际工程中所采用的勘测方法为一维大地电磁法及一维反演方法，使得所得到的电阻率数据不能通过插值等方式实现三维电阻率转换。但针对三维大地电磁法及三维反演所得到的电阻率数据则可以通过插值等方式实现三维转换。附录一中列举了一种插值方法，但不限于该方法。

5)将三维电阻率(ρ_xi,ρ_yi,ρ_zi)数据映射到相应的土壤模型区块(x_i-1,x_i,y_i-1,y_i,z_i-1,z_i)中。该部分的实质为对划分区块后的土壤块进行材料属性设置，该步骤可直接在相关软件中实现。

6)读取注流试验相关设置参数。注流试验是将一试验电极埋设在极址中心区域，在离开试验电极一定距离处埋设辅助电极，在二电极间用电流线将直流电源、电流表串联，并与大地构成试验回路。通过改变测量电极的位置测得不同位置处的地表电位，注流试验法模拟了接地极运行情况，获得的参数真实可靠。注流试验参数主要有：注入电流大小、注流深度、测量路径。确定注流点坐标(x_m,y_m,z_m)及注流试验测点路径(x_m-1,y_m-1,z_m-1；x_m,y_m,z_m)。

7)确定土壤模型边界条件：所谓边界条件即指在接地极真实运行的情况下，由直流电流散流所形成的恒定电流场，在土壤中产生的电位分布情况。因此，易知在无穷远边界处的电位为零。即如图3所示的，除地表外的其余五面为零电势面。

施加激励：所谓激励即为接地极运行时注入大地的直流电流，即在如图3所示原点(0,0,0)处，设置注入电流为6250A。

将模型离散为四面体有限元网格：与步骤1)中离散化类似，当求解整个模型时使十分困难的，为了简化求解的难度，通过将整个模型剖分为若干个小的单元，分别求解这样单元是比较简单的。因此，完成网格剖分后的模型如图7所示，其中图7-a为单个四面体网格的几何示意图。

8)对以下电流场方程进行数值计算，求得地表电位分布；

第一类边界条件：

第二类边界条件：

在上式中，(x,y,z)为单位元的坐标，V为该单位元的电势，ρ为该单位元的电阻率，f(x,y,z)为该单位元的电流函数。

9)取与注流试验一致的测点路径(x_m-1,y_m-1,z_m-1；x_m,y_m,z_m)，并将其地电位V_数与注流试验结果V_注对比，判断|V_数-V_注|<ΔV。

10)当判断结果为否时，确定电位偏移区域，对该区域土壤电阻率进行修正，例如：注流试验的测点路径为(x_m-1,0,0；x_m,0,0)，若V_数-V_注<0，则增大该区域(z＝0,x_m-1<x<x_m,)的土壤电阻率ρ_x；若V_数-V_注>0，则减小该区域(z＝0,x_m-1<x<x_m,)的土壤电阻率ρ_x，然后重复步骤(8)进行计算。

11)当判断结果为是时，读取周边电力系统接线数据：当接地极运行时，强大的直流电流改变了地表电位的分布，使得不同位置的地表电位各不相同。因此使得周边通过中性点接地的电力设施间存在电位差，从而使得在交流电网产生一定的直流分量，而该直流分量受到交流系统参数的影响。所以主要包括周边厂、站的地理接线图即相对于接地极的坐标，厂、站变压器接地电阻，线路分裂数、回数、距离。确定各变压器坐标(x_n,y_n,z_n)及直流电路模型：如图6所示，为某接地极周边两变电站构成的直流电路模型。同上所述，因为A、B两站间存在电位差，使得两站通过输电线路及大地构成了直流通路，由此产生了直流分量，即偏磁电流。偏磁电流的大小与回路中的电阻参数直接相关。

然后获取相应坐标的地电位V_n。

12)将获取的地电位分布V_n映射至直流电路模型中，按如下电路公式进行直流偏磁计算；

其中a为接地点，b为非接地点。

13)、输出计算结果：如图8所示，为采用本发明方法计算得到的某接地极广域范围内的直流偏磁影响结果。通过采用本发明所提出的计算方法，可以直接求得相关厂、站的偏磁电流大小，可以直接为评估偏磁影响提供基础数据。如图8中第一行所示，2.75A即为从站点1流入大地的偏磁电流，而在相关设计中，该站允许通过的最大中性点电流为2.6A。通过计算，可知在接地极正常运行的情况下，该站的偏磁电流已经超标6％，这将严重影响该站变压器的正常运行，并且威胁电力系统的正常运行。所以，通过本发明所提出方法，可以在接地极选址阶段对周边的直流偏磁影响进行有效的评估，避免由接地极投运后偏磁电流超标产生的经济损失，并对接地极的优化设计提供有效的参考。

附录一：

插值拟合:

将该问题用数学方法描述，即已知某一区域(R³是三维欧式空间)内n个点(x_i,y_i,z_i)的测量值ρ_i(i＝1,2,…,n)，对需要求出该点的值ρ。可以采用四维数据插值方法求解这一数学问题，进行土壤电阻率数据的拟合。

(1)二次样条函数边界条件的确定

记I＝[a,b]，△:a＝x₀＜x₁＜…＜x_n＝b为I的一个分划。若给定型值序列及m₀，m_n，则存在唯一的半节点二次插值样条S(x)，它适合

S(x_i)＝y_i(i＝0,1,…,n)

S′(x₀)＝m₀ S′(x_n)＝m_n

若记m_i＝S′(x_i),M_i＝S″(x_i)(i＝0,1,…,n)，则在上S(x)可以表示成

由于

所以二次插值样条S(x)的二阶导数在各内半节点处有跳跃。如果利用最小二乘法，使得内半节点处二阶导数的变化最小，即

这相当于保证了曲率的变化最小。

S(x)的M-关系式为

__iM_i-1+3M_i+λ_iM_i+1＝d_i(1≤i≤n-1) (3)

其中，

λ_i＝h_i+1/(h_i+h_i+1),__i＝1-λ_i,h_i＝x_i-x_i-1,d_i＝8(E_i+1-E_i)/(h_i+h_i+1),E_i＝(y_i-y_i-1)/h_i

由_₁M₀+3M₁+λ₁M₂＝d₁

得M₁+a₁M₂＝b₁+c₁M₀ (4)

其中，a₁＝λ₁/3,b₁＝d₁/3,c1＝-_₁/3。

由_₂M₁+3M₂+λ₂M₃＝d₂

再将式(4)代入，得M₂+a₂M₃＝b₂+c₂M₀ (5)

其中，a₂＝λ₂/(3-_₂a₁),b₂＝(d₂-_₂b₁)/(3-_₂a₁),c₂＝-_₂c₁/(3-_₂a₁)。

反复运用上述方法推下去，可得一般式

M_i+a_iM_i+1＝b_i+c_iM₀(1≤i≤n-1) (6)

其中，a_i＝λ_i/(3-__ia_i-1),b_i＝(d_i-_₃b_i-1)/(3-__ia_i-1),c_i＝-__ic_i-1/(3-__ia_i-1)。

若令

e_i＝(3-__ia_i-1)^-1,a₀＝b₀＝0,c₀＝1

则a_i＝λ_ie_i,b_i＝(d_i-__ib_i-1)e_i,c_i＝-__ic_i-1e_i，特别地，由

M_n-1+a_n-1M_n＝b_n-1+c_n-1M₀

得M_n-1＝F_n-1M₀+G_n-1M_n+H_n-1 (7)

其中，F_n-1＝c_n-1,G_n-1＝-a_n-1,H_n-1＝b_n-1。

由

M_n-2+a_n-2M_n-1＝b_n-2+c_n-2M₀

及式(7)，得

M_n-2＝F_n-2M₀+G_n-2M_n+H_n-2

其中，F_n-2＝c_n-2-a_n-2F_n-1,G_n-2＝-a_n-2G_n-1,H_n-2＝b_n-2-a_n-2H_n-1，反复运用这种方法，可得一般式

M_i＝F_iM₀+G_iM_n+H_i(0≤i≤n) (8)

其中，F_i＝c_i-a_iF_i+1,G_i＝-a_iG_i+1,H_i＝b_i-a_iH_i+1，并规定：F_n＝H_n＝0,G_n＝1，由式(8)，式(2)可写成

若令可得关于M₀，M_n的线性代数方程组

其中，

由柯西不等式，即得

AD-B²≥0

此不等式等号仅在{F_i-F_i-1}与{G_i-G_i-1}(1≤i≤n)线性相关时，等号成立。

当AD-B²≠0时，从方程组(9)可以解出

再由式(8)就可以解得M₁,…,M_n-1。

由二次样条的I型边界条件，可得

即

式(12)即为边界条件。

(2)三元二次样条函数

定义：设在三维直角坐标系o-xyz中的立方体区域上给定一个立方体网格分割其中

△_x:a＝x₀＜x₁＜…＜x_L＝b,

△_y:c＝y₀＜y₁＜…＜y_M＝d,

△_z:e＝z₀＜z₁＜…＜z_N＝f

在R上满足下述两个条件的函数T(x,y,z)称为三元二次样条函数，简称三二次样条。

在每个子立方体

上关于x,y,z都是二次多项式函数，即

其中半节点并约定

在整个R³上，连续，简记作T(x,y,z)∈C^I,I,I(R³)。

此外，给定一组数据{T_ijk}(0≤i≤L；0≤j≤M；0≤k≤N)，若T(x_i,y_j,z_k)＝T_ijk(0≤i≤L；0≤j≤M；0≤k≤N)则称T(x,y,z)为三二次插值样条函数。

在x轴上关于△_x的二次样条函数的全体组成L+3维的线性空间S₂(x；△_x)，它的基样条{h_r(x)}(r＝0,1,…,L+2)满足条件

在y轴上关于分割△_y的二次样条函数的全体组成M+3维线性空间S₂(y；△_y)；在z轴上关于分割△_z的二次样条的全体组成N+3维的线性空间S₂(z；△_z)，它们的基样条{j_s(y)}与{e_t(z)}也像{h_r(x)}一样，在节点处适合类似的条件。

对R³上固定分割由定义的三二次样条的全体构成一个(L+3)(M+3)(N+3)维的线性空间S₂(x,y,z；△)，即

并且三个一元二次基样条的直积

{h_r(x)j_s(y)e_t(z)}(0≤r≤L+2；0≤s≤M+2；0≤t≤N+2)

恰好构成它的一组基底，称之为三二次基样条。于是S₂(x,y,z；△_x)中任一个三二次样条函数可以表示为

上式共有(L+3)(M+3)(N+3)个待定系数{a_rst}，插值条件{T_ijk}有(L+1)(M+1)(N+1)个，所以T(x,y,z)的表达式中尚剩余2(LM+MN+NL)+8(L+M+N)+26个自由度，而这些可由边界条件来确定。由于上述一元二次基样条是针对I型边界给出的，所以这里很自然地应用下面形式的边界条件

立方体R³的六个边界平面上的所有节点处的一阶法向偏导数

其中，T＝0,L；U＝0,M；V＝0,N；0≤i≤L；0≤j≤M；0≤k≤N。

立方体R³的12条边界棱边上的所有节点处的二阶混合偏导数

T,U,V,i,j,k的取值范围与式(14)中相同。

立方体R³的8个顶点处的三阶混合偏导数

S_TUV＝T′″_xyz(x_T,y_U,z_V) (16)

T＝0,L；U＝0,M；V＝0,N。

在R上关于给定分割△的任意一个三二次样条T(x,y,z)，总可以表示为式(13)，直接代入插值条件(3)，并使之满足边界条件(14)、(15)和(16)，则从基样条的性质就可以得到

上面最后四式中I,J,K分别规定如下

上述八式的右端无重复地出现了式(13)中全部系数

a_rst(0≤r≤L+2；0≤s≤M+2；0≤t≤N+2)，并且函数T(x,y,z)的偏导数的连续性可由式(13)中每个一元基样条是C^I连续而保证。因此设给定了直角坐标系o-xyz中的立方体区域R³及其立方体网格分割△，并且给出(L+3)(M+3)(N+3)个常数T_ijk,P_Tjk,q_iUk,r_ijV,PP_iUV,qq_TjV,rr_TUk,S_TUV，

(0≤i≤L；0≤j≤M；0≤k≤N；T＝0,L；U＝0,M；V＝0,N)则存在唯一一个三二次样条函数T(x,y,z)，它以条件(3)为插值条件，以式(14-16)为边界条件。

(3)三二次样条函数边界条件的确定

按定义进行三二次插值需具体给定边界条件(14-16)。利用插值条件(3)来确定边界条件式(14-16)，方法如下：

设T(x,y,z)是这样的三二次插值样条，固定y＝y_j,z＝z_k(0≤j≤M；0≤k≤N)，则T(x,y_j,z_k)是关于x的二次样条函数，它在各节点x_i上的函数值T(x_i,y_j,z_k)(0≤i≤L)可以从插值条件(3)中获得，它在两端点x₀,x_L处的一阶导数值

T′_x＝(x₀,y_j,z_k)＝P_0jk

T′_x＝(x_L,y_j,z_k)＝P_Ljk

可按下式求出

其中M_ijk＝T″_xx(x_i,y_j,z_k)，

而F_i,G_i,H_i满足下列递推公式

且F_L＝H_L＝0,G_L＝1,a₀＝b₀＝0,c₀＝1，M_1jk＝F₁M_0jk+G₁M_Ljk+H₁,M_L-1,jk＝F_L-1M_0jk+G_L-1M_Ljk+H_L-1。这样由式(18)确定了边界条件式(14)中的第一个条件，同样方法可以确定另外两个条件。

固定x＝x_i,z＝z_V(0≤i≤L；V＝0,N)，T′₂(x_i,y,z_V)是关于y的二次样条函数，它在节点y_j处的函数值T′₂(x_i,y_j,z_V)可由刚确定的边界条件式(14)的第三个条件获得。它在两端点y₀,y_M处的一阶导数值

T″_y2(x_i,y₀,z_V)＝PP_i0V

T″_y2(x_i,y_M,z_V)＝PP_iMV

可由前面求边界条件式(14)的方法类似地求出，这样便可求得式(15)中的第一个条件，同理可求出另外两个条件。

固定y＝y_U,z＝z_V(U＝0,M；V＝0,N)，则T″_y2(x,y_U,z_V)是关于x的二次样条函数，它在各节点x_i处的函数值T″_y2(x_i,y_U,z_V)可由刚确定的边界条件式(15)中第一个条件获得。

它在两端点x₀,x_L处的一阶导数值

T′″_xyz(x₀,y_U,z_V)＝S_0UV

T′″_xyz(x_L,y_U,z_V)＝S_LUV

可用求边界条件式(14)的方法类似地求出，这样就确定了边界条件式(16)，至此确定了全部边界条件。

(4)三元二次样条计算方法

对三元二次样条函数可以表示成

其中系数a_rst可由插值条件(3)及边界条件(2)-(4)完全确定。给定一点(x^*,y^*,z^*)∈R³，x^*∈[x_l-1,x_l],y^*∈[y_m-1,y_m],z^*∈[z_n-1,z_n]，(1≤l≤L,1≤m≤M,1≤n≤N)，则

其中，分别表示基样条ψ_s(y),σ_t(z)依次在x_l,y_m,z_n处的一阶导数值，可由二次样条的m-连续性方程组求出。而F₀,F₁,G₀,G₁为混合函数，代入式(20)可以求出T(x^*,y^*,z^*)。

(5)四维超曲面图形表示

三元函数T＝T(x,y,z)∈C(R³)对应高程值W(常数)的等值面图形绘制方法的基本思想如下：

1)将R³作网格分割与三二次样条定义域分划类似；

2)在每个平面z＝z_k(k＝0,1,…,N)上，等值面方程T(x,y,z_k)＝W实质上是关于变量x,y的等值线方程。求出平面z＝z_k上矩形网格上的全部等值点，再作正等轴测投影变换，得到它们在xoy平面上的投影点，然后按一定次序逐点连接这些投影点，绘出等值线的正等轴测投影图形。

3)按k＝1,2,…,N的次序，绘出所有平面z＝z_k上矩形区域内的等值线图形，就形成了一张三元函数T＝T(x,y,z)对应于高程值W的等值面图形。

其中等值线的绘制方法如下：

4)矩形域的网格划分

设其中对[a,b],[c,d]分别作均匀网格分划：

△_x:a＝x₀＜x₁＜…＜x_L＝b

△_y:c＝y₀＜y₁＜…＜y_M＝d

得到的分划若记子矩形的横边与纵边分别为inc,ind，则

x_i＝a+i·inc,(0≤i≤L)

y_j＝c+j·ind,(0≤j≤M)

记f_i,j＝f(x_i,y_j)，网格点(x_i,y_j)记为P_i,j(0≤i≤L；0≤j≤M)。

5)等值点确定

等值线与网格横边或纵边的交点称为等值点。为绘出等值线，必须先求出所有等值点。假设ind,inc充分小，等值点可用线性插值求出。检查每一网格横边或纵边是否有等值点的具体方法如下：

[1]当(W-f_i-1,j)(W-f_i,j)＜0时，在网格横边上有一等值点，其投影坐标为

[2]当(W-f_i,j-1)(W-f_i,j)＜0时，在网格横边上有一等值点，其投影坐标为

记

规定当(W-f_i-1,j)(W-f_i,j)>0时，xx(i,j)＝-1；当(W-f_i,j-1)(W-f_i,j)>0时，yy(i,j)＝-1。xx(i,j),yy(i,j)分别称为相对横、纵坐标。它们分别是横、纵边上有等值点的标志。

[3]等值点的追踪

建立规则，有次序地将等值点逐点连成等值线。

[4]等值线起点和终点的搜索

二元函数对应于高程值的等值线可能有若干条分支，其中有的是开曲线，有的是闭曲线。对每条支曲线，如果能找到其起点，按照等值点的追踪方法，求出其上所有等值点，并逐段连接相邻两等值点，直到其终点，则完成该支曲线的绘制。对于开曲线，其线头在R_xy的边界上，线尾也必在边界上。而对于闭曲线，其上任何一个等值点均可作为起点，该点也同时为终点。

按上述方法用Matlab编制程序可实现从大地电磁法三维反演模型到三维非对称结构土壤电阻率模型的数据拟合。