基于场路耦合的偏磁电流计算方法

1.本发明属于大地土壤电性结构建模及变压器偏磁电流计算领域，特别涉及一种三维非对称结构土壤模型下基于场路耦合的偏磁电流计算方法，主要用于直流偏磁的计算。


背景技术：

2.高压直流输电具有输送容量大，经济性高等特点，在我国能源与负荷中心不平衡分布的背景下发展迅速。随着输送容量，电压等级的提高，直流接地极对周边电力系统的直流偏磁影响也愈发突出。所以，在直流接地极选址与设计阶段，必须对周边电力系统的偏磁影响进行评估。而评估的精准度直接受所采用的土壤模型影响。
3.现有的土壤电性结构建模主要有两种方法，一是将土壤在深度方向分层，认为在单个层深范围内的土壤电阻率是均匀的；二是先将土壤在深度方向分层，然后在单个层深范围内认为土壤电阻率关于垂直于层深方向的轴对称，电阻率自注流点呈射线状变化。上述两种建模方法均均忽略了土壤电阻率各向异性特性，致使其计算的地表电位分布存在较大误差，最终导致对直流偏磁影响的评估失准，不利于电力系统的安全运行。
4.现有的土壤电性建模中，主要有水平分层土壤模型和复合分层土壤模型。水平分层结构将某一深度区域电阻率变化较小的土壤近似处理为一个电阻率分布均匀的分层结构，一般将整个土壤模型分为三至七个分层；主要缺点在于只考虑了土壤电阻率在深度方向的变化，且分层数较少。复合分层结构在深度方向的分层结构与水平分层结构类似，在垂直于深度方向认为电阻率以接地极为原点沿射线方向变化，其最终的分层结构关于深度方向呈轴对称分布；然而实际的土壤因为岩层和土体的组合方式的不同，在横、纵方向上呈现出不同的电阻率，且可能具有很大差异，如在横向上表现为高阻特性，而在纵向上表现为低阻特性。所以，上述的两种土壤模型均不能真实的反映土壤的电性结构，最终将使得偏磁电流评估失准，危害电网的安全运行和增大经济投入。
5.现阶段的偏磁电流主要是根据水平分层或复合分层土壤模型正演得到的地表电位进行计算。现有的关于偏磁电流计算方法的专利，如zl201510093440.0“多直流接地极不同运行方式下直流偏磁电流影响站点的预测方法”，提供了一种多直流接地极不同运行方式下直流偏磁电流影响站点的预测方法，通过监测装置获取变压器中的偏磁电流数据，采用这些数据对假想的水平分层模型进行修正，由此得到较为准确的地表电位分布情况，进而对偏磁电流的分布进行预测。这一方法虽然对土壤模型进行了修正，但其采用的水平分层土壤模型将真实的电性结构做了很大程度的简化，只考虑了土壤电阻率在深度方向上的变化，与实际情况差别较大，所以使得其偏磁电流计算结果仍存在较大误差。


技术实现要素：

6.为此，本发明提供一种三维非对称结构土壤模型下基于场路耦合的直流偏磁计算方法，通过建立高度模拟真实电性结构的土壤模型，充分考虑土壤电阻率的各向异性特性，提高偏磁电流的计算精度。
7.本发明采取的技术方案为：
8.一种三维非对称结构土壤模型下基于场路耦合的直流偏磁计算方法，通过需要计及的直流偏磁影响范围确定土壤模型的尺寸，将模型离散为边长2m的小立方体；通过实际的测点数据，确定模型的区块划分；将反演的三维电阻率映射至各区块，非三维电阻率数据通过转换后映射；通过注流试验数据确定注流点坐标，施加激励，确定边界条件，划分网格，进行地表电位计算；确定观测路径，与注流试验结果进行对比，修正土壤模型；获取接地极周边电力系统接线图，坐标等参数，搭建直流电路模型；输入各节点电势，进行直流偏磁计算。
9.一种三维非对称结构土壤模型下基于场路耦合的直流偏磁计算方法，包括以下步骤：
10.步骤1：根据需要评估的直流偏磁影响范围，确定待建立的土壤模型尺寸，并对土壤模型进行离散化操作；
11.步骤2：根据浅层大地电阻率和深层大地电阻率测点数据，对土壤模型进行区块划分；
12.步骤3：读取电阻率数据，判断是否为三维电阻率数据，若不是则将其转换为三维电阻率；
13.步骤4：将三维电阻率数据映射到相应的土壤模型区块中；
14.步骤5：读取注流试验相关设置参数，确定注流点坐标及注流试验测点路径；
15.步骤6：确定土壤模型边界条件，施加激励，将模型离散为四面体有限元网格；
16.步骤7：计算地表电位分布；
17.步骤8：取与注流试验一致的测点路径，并将其地电位与注流试验结果对比，判断|v
数
‑
v
注
|<δv；
18.步骤9：当判断结果为否时，对土壤电阻率进行修正，并再次进行计算；
19.步骤10：当判断结果为是时，读取周边电力系统接线数据，确定各变压器坐标及直流电路模型，然后获取相应坐标地电位分布；
20.步骤11：将获取的地电位分布映射至直流电路模型中，并进行直流偏磁计算；
21.步骤12：输出结果。
22.步骤1中，根据需要评估的直流偏磁影响范围，确定待建立的土壤模型尺寸，并对土壤模型进行离散化操作，土壤模型离散化程度可以根据计算精度需求自由控制。
23.步骤2中，根据浅层大地电阻率和深层大地电阻率测点数据，将浅层土壤中电阻率变化较小的区域合并为一个区块，减小深层土壤的离散程度，将其划为一个较大的区块。
24.大地的电性结构与土壤中导电离子和含水量相关，在接地极的广域范围内可能存在湖泊、河流、堰塘、金属矿床等电阻率均匀的区域，在步骤1的离散化过程中，是将整个求解域内的土壤离散为若干个小土壤块，当然也包括上述电阻率较为均匀的区域。而在进行地表电位计算时，会重复计算这些电阻率均匀的小立方体的边界方程使得求解的矩阵数量增多，计算时间增长。同时根据文献(邱立,等.直流接地极注流试验参数对地表电位分布的影响[j].三峡大学学报(自然科学版),2017,39(1):79
‑
83.)，深层土壤及远离接地极区域的土壤对地表电位的影响较小。综上两点，为了进一步优化所述方法
的计算时间，故将浅层土壤中电阻率变化较小的区域合并为一个区块，减小深层及距离接地极较远区域土壤的离散程度，将其划为一个较大的区块。
[0025]
步骤3中，读取电阻率数据，判断是否为三维电阻率数据，若不是则将其通过线性插值或二次函数插值的方式转换为三维电阻率。
[0026]
步骤4中，三维电阻率数据为：由大地电磁法或大地音频电磁法经三维反演得到的电阻率数据。
[0027]
在对接地极深层电阻率勘测方面，现工程上普遍采用大地电磁法，该方法也运用于物理勘探领域，主要涉及地球物理学。大地电磁法反演根据地表实测的视电阻率、相位等数据来求取大地深部电导率结构，该电导率结构的正演响应能极好地拟合视电阻率、相位等实测数据。反演算法的目的是寻找一个合理且与某种地球物理观测结果相符合的地球物理模型，其原理是通过相应的数学运算，利用实际测得的地表电磁场响应，如视电阻率、阻抗相位、倾子等，获得一个符合实际情况的地电模型。根据反演维度可以分为一维反演、二维反演、三维反演。三维反演的主要优点在于，其考虑了大地电性在三个维度上的变化，较一维反演具有分辨率高，对异常体判断准确的优点，更加适合用于立方体模型计算。现阶段三维正演方面的研究已趋于成熟，随着三维正演的发展，mt的三维反演研究也日趋升温，反演方法众多，主要有共轭梯度法极大似然反演、非线性共轭梯度反演、快速松弛反演和人工神经网络反演等。在接地极深层电阻率勘测中主要采用的是一维反演方法，现将其他领域的最新成果引入到接地极地表电位计算领域，以使得参与地表电位计算的初始数据更加准确。
[0028]
步骤8中，当判断结果为否时，对土壤电阻率进行修正，当v
数
<v
注
时，增大所取路径上的土壤电阻率；当v
数
>v
注
时，减小所取路径上的土壤电阻率，然后再次进行计算。
[0029]
与现有最好技术相比，本发明一种三维非对称结构土壤模型下基于场路耦合的直流偏磁计算方法，的优点在于：
[0030]
1.本发明所建立的土壤模型为三维模型，且土壤块电阻率赋值为张量，更能真实反映土壤电阻率的各向变化，减小了模型误差，提高了直流偏磁的评估精度；
[0031]
2.本发明所建立的土壤模型可以模拟土壤中存在电阻率异常岩土体的情况，并可以分析由此带来的地表电位分布的影响；
[0032]
3.本发明所提出的三维非对称结构土壤模型下基于场路耦合的直流偏磁计算方法，可以为直流接地极的选址与设计工作提供有效的参考。
[0033]
4.本发明提供了一种三维非对称结构土壤模型下基于场路耦合的直流偏磁计算方法，从建立真实土壤模型出发，考虑了土壤电阻率的各向异性特性，减小了地表电位分布的计算误差，提高了直流偏磁影响的评估精度。
附图说明
[0034]
图1为步骤1)中的离散化过程示意图。
[0035]
图2为本发明方法的计算流程图。
[0036]
图3是本发明所提出的土壤模型几何结构示意图。
[0037]
图4是本发明所提出的经离散化处理后的土壤模型几何结构示意图。
[0038]
图5是某接地极浅层大地电阻率测点位置分布示意图。
[0039]
图6是本发明所提出的直流电路模型示意图。
[0040]
图7是本发明所提出的划分网格后土壤模型几何结构示意图及单个四面体网格示意图。其中：图7
‑
a为单个四面体网格的几何示意图。
[0041]
图8是采用本发明方法计算得到的某接地极广域范围内直流偏磁结果图。
具体实施方式
[0042]
一种三维非对称结构土壤模型下基于场路耦合的直流偏磁计算方法，通过需要计及的直流偏磁影响范围确定土壤模型的尺寸，将模型离散为边长2m的小立方体；通过实际的测点数据，确定模型的区块划分；将反演的三维电阻率映射至各区块，非三维电阻率数据通过转换后映射；通过注流试验数据确定注流点坐标，施加激励，确定边界条件，划分网格，进行地表电位计算；确定观测路径，与注流试验结果进行对比，修正土壤模型；获取接地极周边电力系统接线图，坐标等参数，搭建直流电路模型；输入各节点电势，进行直流偏磁计算。本发明从初始模型出发，提高了直流偏磁的计算精度，在直流工程投运后减少了保护误动或拒动的可能，提高了电力系统运行的稳定性。
[0043]
一种三维非对称结构土壤模型下基于场路耦合的直流偏磁计算方法，包括以下步骤：
[0044]
步骤1：根据需要评估的直流偏磁影响范围，确定待建立的土壤模型尺寸，并对土壤模型进行离散化操作。
[0045]
进行离散化操作的目的是将整个求解区域的土壤块离散为小的土壤块，以减少求解难度。
[0046]
具体操作如图1所示，根据土壤模型尺寸及求解精度，首先在x轴一端建立一个小立方体(元)，然后通过复制操作在x轴上形成一列小立方体(线)，然后再通过复制操作将这列小立方体在xy平面上形成一面域内小立方体(面)，最后通过复制操作将这一面域内的小立方体在xyz直角坐标系中形成以小立方体为基的大立方体(体)。上述过程只是离散化的一种方法，但不局限于该方法。
[0047]
步骤2：根据浅层大地电阻率和深层大地电阻率测点数据，对土壤模型进行区块划分；
[0048]
步骤3：读取电阻率数据，判断是否为三维电阻率数据，若不是则将其转换为三维电阻率；
[0049]
步骤4：将三维电阻率数据映射到相应的土壤模型区块中；
[0050]
步骤5：读取注流试验相关设置参数，确定注流点坐标及注流试验测点路径；
[0051]
步骤6：确定土壤模型边界条件，施加激励，将模型离散为四面体有限元网格；
[0052]
步骤7：计算地表电位分布；
[0053]
步骤8：取与注流试验一致的测点路径，并将其地电位与注流试验结果对比，判断|v
数
‑
v
注
|<δv；
[0054]
步骤9：当判断结果为否时，对土壤电阻率进行修正，并再次进行计算；
[0055]
步骤10：当判断结果为是时，读取周边电力系统接线数据，确定各变压器坐标及直流电路模型，然后获取相应坐标地电位分布；
[0056]
步骤11：将获取的地电位分布映射至直流电路模型中，并进行直流偏磁计算；
[0057]
步骤12：输出结果。
[0058]
步骤1中，根据需要评估的直流偏磁影响范围，确定待建立的土壤模型尺寸，并对土壤模型进行离散化操作，土壤模型离散化程度可以根据计算精度需求自由控制。
[0059]
步骤2中，根据浅层大地电阻率和深层大地电阻率测点数据，将浅层土壤中电阻率变化较小的区域合并为一个区块，减小深层土壤的离散程度，将其划为一个较大的区块。
[0060]
步骤3中，读取电阻率数据，判断是否为三维电阻率数据，若不是则将其通过线性插值或二次函数插值的方式转换为三维电阻率。
[0061]
步骤4中，三维电阻率数据为：由大地电磁法或大地音频电磁法经三维反演得到的电阻率数据。
[0062]
步骤8中，当判断结果为否时，对土壤电阻率进行修正，当v
数
<v
注
时，增大所取路径上的土壤电阻率；当v
数
>v
注
时，减小所取路径上的土壤电阻率，然后再次进行计算。
[0063]
本发明所述的三维非对称结构下基于场路耦合的直流偏磁计算方法，其计算流程图如图2所示，主要分为数据预处理部分，几何建模部分，地表电位计算部分，直流偏磁计算部分。
[0064]
本发明所述计算方法由以下步骤组成：
[0065]
1)根据需要评估的直流偏磁影响范围：根据dl/t5224
‑
2014，应对极址周围地电位升大于3v的变电站，以及与其有直接电气联系的变电站进行仿真计算。现随着直流输电技术的发展，直流接地极的入地电流也逐步增大，其产生偏磁影响的范围也随之增大，所以建议将150km范围内的变电站、电厂等纳入计算网络。
[0066]
确定待建立的土壤模型尺寸：为了保证直流偏磁的计算精度，所建立的土壤模型要大于需要计算偏磁电流的计算网络。所以，土壤模型的尺寸确定为400km，如图3所示。将接地极作为点元处理，其坐标设为(0,0,0)，根据模型尺寸建立电流场计算几何模型，即为原始的土壤模型，如图3所示。其中原点(0,0,0)为直流电流注入点，除地表外其余五面为零电势面。并将模型离散为边长为2m的小立方体，离散后的土壤模型如图4所示。
[0067]
2)根据浅层大地电阻率测点数据，如表1所示，即为某接地极的浅层大地电阻率数据，所有测点位置分布如图5所示。现以此为例进行浅层大地电阻率的区块划分，但不限于此种划分方式。
[0068]
a，首先确定在深度方向上的区块层深，依次为2m、3m、5m、5m、
……
、300m，共计13层；
[0069]
b，确定在地表平面上的区块个数为49个，故设置每个区块的长、宽分别为150m、150m。
[0070]
综上，浅层土壤区块划分完毕。其中部分区块无对应的电阻率数据，采用对应层深的视电阻率填充。获取分层坐标(x
i
,y
i
,z
i
)，对接地极(0,0,0)邻近区域进行区块(x
i
‑1,x
i
,y
i
‑1,y
i
,z
i
‑1,z
i
)划分。
[0071]
表1某接地极浅层大地电阻率数据
[0072][0073]
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。，深层区块划分。
[0074]
表2某接地极深层大地电阻率数据
[0075][0076]
4)读取电阻率数据，判断是否为三维电阻率数据。如表3所示，为某局部区域的三维电阻率数据，表3中层厚即为该土壤层在深度方向的厚度，x、y、z为指定的直角坐标系对应的三个坐标轴。实际的大地因为土壤电阻率各向异性特性所体现出的电性在各方向上并
不一致，如表3所示，各轴上电阻率的值及变化各不相同，且同一轴的正负半轴其值及变化也不相同，这体现了所建模型的非对称特性。若不是则将其通过线性插值或二次函数插值转换为三维电阻率。
[0077]
表3某三维电阻率数据
[0078][0079]
上述表格为本发明所建议采用的三维电阻率数据格式。实际工程中所采用的勘测方法为一维大地电磁法及一维反演方法，使得所得到的电阻率数据不能通过插值等方式实现三维电阻率转换。但针对三维大地电磁法及三维反演所得到的电阻率数据则可以通过插值等方式实现三维转换。附录一中列举了一种插值方法，但不限于该方法。
[0080]
5)将三维电阻率(ρ
xi
,ρ
yi
,ρ
zi
)数据映射到相应的土壤模型区块(x
i
‑1,x
i
,y
i
‑1,y
i
,z
i
‑1,z
i
)中。该部分的实质为对划分区块后的土壤块进行材料属性设置，该步骤可直接在相关软件中实现。
[0081]
6)读取注流试验相关设置参数。注流试验是将一试验电极埋设在极址中心区域，在离开试验电极一定距离处埋设辅助电极，在二电极间用电流线将直流电源、电流表串联，并与大地构成试验回路。通过改变测量电极的位置测得不同位置处的地表电位，注流试验法模拟了接地极运行情况，获得的参数真实可靠。注流试验参数主要有：注入电流大小、注流深度、测量路径。确定注流点坐标(x
m
,y
m
,z
m
)及注流试验测点路径(x
m
‑1,y
m
‑1,z
m
‑1；x
m
,y
m
,z
m
)。
[0082]
7)确定土壤模型边界条件：所谓边界条件即指在接地极真实运行的情况下，由直流电流散流所形成的恒定电流场，在土壤中产生的电位分布情况。因此，易知在无穷远边界处的电位为零。即如图3所示的，除地表外的其余五面为零电势面。
[0083]
施加激励：所谓激励即为接地极运行时注入大地的直流电流，即在如图3所示原点(0,0,0)处，设置注入电流为6250a。
[0084]
将模型离散为四面体有限元网格：与步骤1)中离散化类似，当求解整个模型时使十分困难的，为了简化求解的难度，通过将整个模型剖分为若干个小的单元，分别求解这样单元是比较简单的。因此，完成网格剖分后的模型如图7所示，其中图7
‑
a为单个四面体网格的几何示意图。
[0085]
8)对以下电流场方程进行数值计算，求得地表电位分布；
[0086][0087]
第一类边界条件：
[0088][0089]
第二类边界条件：
[0090][0091]
在上式中，(x,y,z)为单位元的坐标，v为该单位元的电势，ρ为该单位元的电阻率，f(x,y,z)为该单位元的电流函数。
[0092]
9)取与注流试验一致的测点路径(x
m
‑1,y
m
‑1,z
m
‑1；x
m
,y
m
,z
m
)，并将其地电位v
数
与注流试验结果v
注
对比，判断|v
数
‑
v
注
|<δv。
[0093]
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)进行计算。
[0094]
11)当判断结果为是时，读取周边电力系统接线数据：当接地极运行时，强大的直流电流改变了地表电位的分布，使得不同位置的地表电位各不相同。因此使得周边通过中性点接地的电力设施间存在电位差，从而使得在交流电网产生一定的直流分量，而该直流分量受到交流系统参数的影响。所以主要包括周边厂、站的地理接线图即相对于接地极的坐标，厂、站变压器接地电阻，线路分裂数、回数、距离。确定各变压器坐标(x
n
,y
n
,z
n
)及直流电路模型：如图6所示，为某接地极周边两变电站构成的直流电路模型。同上所述，因为a、b两站间存在电位差，使得两站通过输电线路及大地构成了直流通路，由此产生了直流分量，即偏磁电流。偏磁电流的大小与回路中的电阻参数直接相关。
[0095]
然后获取相应坐标的地电位v
n
。
[0096]
12)将获取的地电位分布v
n
映射至直流电路模型中，按如下电路公式进行直流偏磁计算；
[0097][0098]
其中a为接地点，b为非接地点。
[0099]
13)、输出计算结果：如图8所示，为采用本发明方法计算得到的某接地极广域范围内的直流偏磁影响结果。通过采用本发明所提出的计算方法，可以直接求得相关厂、站的偏磁电流大小，可以直接为评估偏磁影响提供基础数据。如图8中第一行所示，2.75a即为从站点1流入大地的偏磁电流，而在相关设计中，该站允许通过的最大中性点电流为2.6a。通过
计算，可知在接地极正常运行的情况下，该站的偏磁电流已经超标6％，这将严重影响该站变压器的正常运行，并且威胁电力系统的正常运行。所以，通过本发明所提出方法，可以在接地极选址阶段对周边的直流偏磁影响进行有效的评估，避免由接地极投运后偏磁电流超标产生的经济损失，并对接地极的优化设计提供有效的参考。
[0100]
附录一：
[0101]
插值拟合:
[0102]
将该问题用数学方法描述，即已知某一区域(r3是三维欧式空间)内n个点(x
i
,y
i
,z
i
)的测量值ρ
i
(i＝1,2,
…
,n)，对需要求出该点的值ρ。可以采用四维数据插值方法求解这一数学问题，进行土壤电阻率数据的拟合。
[0103]
(1)二次样条函数边界条件的确定
[0104]
记i＝[a,b]，
△
:a＝x0＜x1＜
…
＜x
n
＝b为i的一个分划。若给定型值序列及m0，m
n
，则存在唯一的半节点二次插值样条s(x)，它适合
[0105]
s(x
i
)＝y
i
(i＝0,1,
…
,n)
[0106]
s
′
(x0)＝m
0 s
′
(x
n
)＝m
n
[0107]
若记m
i
＝s
′
(x
i
),m
i
＝s
″
(x
i
)(i＝0,1,
…
,n)，则在上s(x)可以表示成
[0108][0109]
由于
[0110][0111]
所以二次插值样条s(x)的二阶导数在各内半节点处有跳跃。如果利用最小二乘法，使得内半节点处二阶导数的变化最小，即
[0112][0113]
这相当于保证了曲率的变化最小。
[0114]
s(x)的m
‑
关系式为
[0115]
_
i
m
i
‑1+3m
i
+λ
i
m
i+1
＝d
i
(1≤i≤n
‑
1)
ꢀꢀꢀ
(3)
[0116]
其中，
[0117]
λ
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
[0118]
由_1m0+3m1+λ1m2＝d1[0119]
得m1+a1m2＝b1+c1m0ꢀꢀꢀ
(4)
[0120]
其中，a1＝λ1/3,b1＝d1/3,c1＝
‑
_1/3。
[0121]
由_2m1+3m2+λ2m3＝d2[0122]
再将式(4)代入，得m2+a2m3＝b2+c2m0ꢀꢀꢀ
(5)
[0123]
其中，a2＝λ2/(3
‑
_2a1),b2＝(d2‑
_2b1)/(3
‑
_2a1),c2＝
‑
_2c1/(3
‑
_2a1)。
[0124]
反复运用上述方法推下去，可得一般式
[0125]
m
i
+a
i
m
i+1
＝b
i
+c
i
m0(1≤i≤n
‑
1)
ꢀꢀꢀ
(6)
[0126]
其中，a
i
＝λ
i
/(3
‑
_
i
a
i
‑1),b
i
＝(d
i
‑
_3b
i
‑1)/(3
‑
_
i
a
i
‑1),c
i
＝
‑
_
i
c
i
‑1/(3
‑
_
i
a
i
‑1)。
[0127]
若令
[0128]
e
i
＝(3
‑
_
i
a
i
‑1)
‑1,a0＝b0＝0,c0＝1
[0129]
则a
i
＝λ
i
e
i
,b
i
＝(d
i
‑
_
i
b
i
‑1)e
i
,c
i
＝
‑
_
i
c
i
‑1e
i
，特别地，由
[0130]
m
n
‑1+a
n
‑1m
n
＝b
n
‑1+c
n
‑1m0[0131]
得m
n
‑1＝f
n
‑1m0+g
n
‑1m
n
+h
n
‑1ꢀꢀꢀ
(7)
[0132]
其中，f
n
‑1＝c
n
‑1,g
n
‑1＝
‑
a
n
‑1,h
n
‑1＝b
n
‑1。
[0133]
由
[0134]
m
n
‑2+a
n
‑2m
n
‑1＝b
n
‑2+c
n
‑2m0[0135]
及式(7)，得
[0136]
m
n
‑2＝f
n
‑2m0+g
n
‑2m
n
+h
n
‑2[0137]
其中，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，反复运用这种方法，可得一般式
[0138]
m
i
＝f
i
m0+g
i
m
n
+h
i
(0≤i≤n)
ꢀꢀꢀ
(8)
[0139]
其中，f
i
＝c
i
‑
a
i
f
i+1
,g
i
＝
‑
a
i
g
i+1
,h
i
＝b
i
‑
a
i
h
i+1
，并规定：f
n
＝h
n
＝0,g
n
＝1，由式(8)，式(2)可写成
[0140][0141]
若令可得关于m0，m
n
的线性代数方程组
[0142][0143]
其中，其中，
[0144]
由柯西不等式，即得
[0145]
ad
‑
b2≥0
[0146]
此不等式等号仅在{f
i
‑
f
i
‑1}与{g
i
‑
g
i
‑1}(1≤i≤n)线性相关时，等号成立。
[0147]
当ad
‑
b2≠0时，从方程组(9)可以解出
[0148][0149]
再由式(8)就可以解得m1,
…
,m
n
‑1。
[0150]
由二次样条的i型边界条件，可得
[0151][0152]
即
[0153][0154]
式(12)即为边界条件。
[0155]
(2)三元二次样条函数
[0156]
定义：设在三维直角坐标系o
‑
xyz中的立方体区域上给定一个立方体网格分割其中
[0157]
△
x
:a＝x0＜x1＜
…
＜x
l
＝b,
[0158]
△
y
:c＝y0＜y1＜
…
＜y
m
＝d,
[0159]
△
z
:e＝z0＜z1＜
…
＜z
n
＝f
[0160]
在r上满足下述两个条件的函数t(x,y,z)称为三元二次样条函数，简称三二次样条。
[0161]
在每个子立方体
[0162]
上关于x,y,z都是二次多项式函数，即
[0163][0164]
其中半节点并约定
[0165][0166]
在整个r3上，连续，简记作t(x,y,z)∈c
i,i,i
(r3)。
[0167]
此外，给定一组数据{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)为三二次插值样条函数。
[0168]
在x轴上关于
△
x
的二次样条函数的全体组成l+3维的线性空间s2(x；
△
x
)，它的基样条{h
r
(x)}(r＝0,1,
…
,l+2)满足条件
[0169][0170]
在y轴上关于分割
△
y
的二次样条函数的全体组成m+3维线性空间s2(y；
△
y
)；在z轴上关于分割
△
z
的二次样条的全体组成n+3维的线性空间s2(z；
△
z
)，它们的基样条{j
s
(y)}与{e
t
(z)}也像{h
r
(x)}一样，在节点处适合类似的条件。
[0171]
对r3上固定分割由定义的三二次样条的全体构成一个(l+3)(m+3)(n+3)维的线性空间s2(x,y,z；
△
)，即
[0172][0173]
并且三个一元二次基样条的直积
[0174]
{h
r
(x)j
s
(y)e
t
(z)}(0≤r≤l+2；0≤s≤m+2；0≤t≤n+2)
[0175]
恰好构成它的一组基底，称之为三二次基样条。于是s2(x,y,z；
△
x
)中任一个三二次样条函数可以表示为
[0176][0177]
上式共有(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型边界给出的，所以这里很自然地应用下面形式的边界条件
[0178]
立方体r3的六个边界平面上的所有节点处的一阶法向偏导数
[0179][0180]
其中，t＝0,l；u＝0,m；v＝0,n；0≤i≤l；0≤j≤m；0≤k≤n。
[0181]
立方体r3的12条边界棱边上的所有节点处的二阶混合偏导数
[0182][0183]
t,u,v,i,j,k的取值范围与式(14)中相同。
[0184]
立方体r3的8个顶点处的三阶混合偏导数
[0185]
s
tuv
＝t
′″
xyz
(x
t
,y
u
,z
v
)
ꢀꢀꢀ
(16)
[0186]
t＝0,l；u＝0,m；v＝0,n。
[0187]
在r上关于给定分割
△
的任意一个三二次样条t(x,y,z)，总可以表示为式(13)，直接代入插值条件(3)，并使之满足边界条件(14)、(15)和(16)，则从基样条的性质就可以得到
[0188][0189]
上面最后四式中i,j,k分别规定如下
[0190][0191]
上述八式的右端无重复地出现了式(13)中全部系数
[0192]
a
rst
(0≤r≤l+2；0≤s≤m+2；0≤t≤n+2)，并且函数t(x,y,z)的偏导数的连续性可由式(13)中每个一元基样条是c
i
连续而保证。因此设给定了直角坐标系o
‑
xyz中的立方体区域r3及其立方体网格分割
△
，并且给出(l+3)(m+3)(n+3)个常数t
ijk
,p
tjk
,q
iuk
,r
ijv
,pp
iuv
,qq
tjv
,rr
tuk
,s
tuv
，
[0193]
(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)为边界条件。
[0194]
(3)三二次样条函数边界条件的确定
[0195]
按定义进行三二次插值需具体给定边界条件(14
‑
16)。利用插值条件(3)来确定边界条件式(14
‑
16)，方法如下：
[0196]
设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)中获得，它在两端点x0,x
l
处的一阶导数值
[0197]
t
′
x
＝(x0,y
j
,z
k
)＝p
0jk
[0198]
t
′
x
＝(x
l
,y
j
,z
k
)＝p
ljk
[0199]
可按下式求出
[0200]
[0201]
其中m
ijk
＝t
″
xx
(x
i
,y
j
,z
k
)，
[0202][0203][0204]
而f
i
,g
i
,h
i
满足下列递推公式
[0205][0206]
且f
l
＝h
l
＝0,g
l
＝1,a0＝b0＝0,c0＝1，m
1jk
＝f1m
0jk
+g1m
ljk
+h1,m
l
‑
1,jk
＝f
l
‑1m
0jk
+g
l
‑1m
ljk
+h
l
‑1。这样由式(18)确定了边界条件式(14)中的第一个条件，同样方法可以确定另外两个条件。
[0207]
固定x＝x
i
,z＝z
v
(0≤i≤l；v＝0,n)，t
′2(x
i
,y,z
v
)是关于y的二次样条函数，它在节点y
j
处的函数值t
′2(x
i
,y
j
,z
v
)可由刚确定的边界条件式(14)的第三个条件获得。它在两端点y0,y
m
处的一阶导数值
[0208]
t
″
y2
(x
i
,y0,z
v
)＝pp
i0v
[0209]
t
″
y2
(x
i
,y
m
,z
v
)＝pp
imv
[0210]
可由前面求边界条件式(14)的方法类似地求出，这样便可求得式(15)中的第一个条件，同理可求出另外两个条件。
[0211]
固定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)中第一个条件获得。
[0212]
它在两端点x0,x
l
处的一阶导数值
[0213]
t
′″
xyz
(x0,y
u
,z
v
)＝s
0uv
[0214]
t
′″
xyz
(x
l
,y
u
,z
v
)＝s
luv
[0215]
可用求边界条件式(14)的方法类似地求出，这样就确定了边界条件式(16)，至此确定了全部边界条件。
[0216]
(4)三元二次样条计算方法
[0217]
对三元二次样条函数可以表示成
[0218]
[0219]
其中系数a
rst
可由插值条件(3)及边界条件(2)
‑
(4)完全确定。给定一点(x
*
,y
*
,z
*
)∈r3，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)，则
[0220][0221]
其中，分别表示基样条ψ
s
(y),σ
t
(z)依次在x
l
,y
m
,z
n
处的一阶导数值，可由二次样条的m
‑
连续性方程组求出。而f0,f1,g0,g1为混合函数，代入式(20)可以求出t(x
*
,y
*
,z
*
)。
[0222]
(5)四维超曲面图形表示
[0223]
三元函数t＝t(x,y,z)∈c(r3)对应高程值w(常数)的等值面图形绘制方法的基本思想如下：
[0224]
1)将r3作网格分割与三二次样条定义域分划类似；
[0225]
2)在每个平面z＝z
k
(k＝0,1,
…
,n)上，等值面方程t(x,y,z
k
)＝w实质上是关于变量x,y的等值线方程。求出平面z＝z
k
上矩形网格上的全部等值点，再作正等轴测投影变换，得到它们在xoy平面上的投影点，然后按一定次序逐点连接这些投影点，绘出等值线的正等轴测投影图形。
[0226]
3)按k＝1,2,
…
,n的次序，绘出所有平面z＝z
k
上矩形区域内的等值线图形，就形成了一张三元函数t＝t(x,y,z)对应于高程值w的等值面图形。
[0227]
其中等值线的绘制方法如下：
[0228]
4)矩形域的网格划分
[0229]
设其中对[a,b],[c,d]分别作均匀网格分划：
[0230]
△
x
:a＝x0＜x1＜
…
＜x
l
＝b
[0231]
△
y
:c＝y0＜y1＜
…
＜y
m
＝d
[0232]
得到的分划若记子矩形的横边与纵边分别为inc,ind，则
[0233]
x
i
＝a+i
·
inc,(0≤i≤l)
[0234]
y
j
＝c+j
·
ind,(0≤j≤m)
[0235]
记f
i,j
＝f(x
i
,y
j
)，网格点(x
i
,y
j
)记为p
i,j
(0≤i≤l；0≤j≤m)。
[0236]
5)等值点确定
[0237]
等值线与网格横边或纵边的交点称为等值点。为绘出等值线，必须先求出所有等值点。假设ind,inc充分小，等值点可用线性插值求出。检查每一网格横边或纵边是否有等值点的具体方法如下：
[0238]
[1]当(w
‑
f
i
‑
1,j
)(w
‑
f
i,j
)＜0时，在网格横边上有一等值点，其投影坐标为
[0239][0240]
[2]当(w
‑
f
i,j
‑1)(w
‑
f
i,j
)＜0时，在网格横边上有一等值点，其投影坐标为
[0241][0242]
记
[0243][0244]
规定当(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)分别称为相对横、纵坐标。它们分别是横、纵边上有等值点的标志。
[0245]
[3]等值点的追踪
[0246]
建立规则，有次序地将等值点逐点连成等值线。
[0247]
[4]等值线起点和终点的搜索
[0248]
二元函数对应于高程值的等值线可能有若干条分支，其中有的是开曲线，有的是闭曲线。对每条支曲线，如果能找到其起点，按照等值点的追踪方法，求出其上所有等值点，并逐段连接相邻两等值点，直到其终点，则完成该支曲线的绘制。对于开曲线，其线头在r
xy
的边界上，线尾也必在边界上。而对于闭曲线，其上任何一个等值点均可作为起点，该点也同时为终点。
[0249]
按上述方法用matlab编制程序可实现从大地电磁法三维反演模型到三维非对称结构土壤电阻率模型的数据拟合。