基于贝叶斯弹性网正则化方法的电流元三维反演方法转让专利
申请号 : CN201910363304.7
文献号 : CN110276094A
文献日 : 2019-09-24
发明人 : 胡军 , 赵根 , 何金良 , 王善祥 , 欧阳勇 , 王中旭 , 曾嵘 , 庄池杰 , 张波 , 余占清
申请人 : 清华大学
摘要 :
一种基于贝叶斯弹性网正则化方法的电流元三维反演方法，通过待求解参数的先验分布信息及已知的测量数据。其有益效果是：贝叶斯正则化方法克服了离散电流元计算过程中的病态性，具有较高的计算鲁棒性和准确度。另一方面，其参数选择的方式更为温和，可以设定电流元的先验分布，输出参数有一定的置信区间。
权利要求 :
1.一种基于贝叶斯弹性网正则化方法的电流元三维反演方法，通过待求解参数的先验分布信息及已知的磁场测量数据，进行空间体电流或线电流分布的反演重建，其特征在于，步骤1、将磁场传感器测量结果整理为磁场强度H向量；
步骤2、以贝叶斯定理表示待求电流元分布I的后验概率 ；
步骤3、引入超参数σ，由先验分布及超先验分布表示步骤2中后验概率分布的各个参数；
步骤4、按照贝叶斯弹性网正则化框架设计分层模型中的电流元先验概率分布，得到基于弹性网正则化的分层贝叶斯模型；
步骤5、根据贝叶斯公式进行后验概率分布的求解，采用吉布斯采样方法，根据所述分层贝叶斯模型计算得到后验概率分布值；
步骤6、根据已有磁场数据和待求离散电流元网格的位置坐标，进行上式中超参数及正则化参数的求解；
根据毕奥-萨伐尔定律，磁场可表示为电流元的积分形式，对于待求空间离散网格中的电流元，磁场与电流元间的关系可以矩阵形式表示。
2.根据权利要求1中所述的基于贝叶斯弹性网正则化方法的电流元三维反演方法，其特征在于，认为磁场传感器测量结果中引入了误差ε： （1）
其中H为磁场强度，I为待求电流元分布，G为引导场矩阵。
3.根据权利要求1中所述的基于贝叶斯弹性网正则化方法的电流元三维反演方法，其特征在于，在引导场矩阵G确定的情况下，根据贝叶斯定理，待求电流元分布I表示为：公式2中， 为是似然函数， 为已知电流元的先验概率分布，为边缘似然，在此为常数。
4.根据权利要求1中所述的基于贝叶斯弹性网正则化方法的电流元三维反演方法，其特征在于，所述步骤4中，基于弹性网正则化框架，对应电流元的先验概率分布为：公式（4）中，λ1，λ2为正则化系数。
5.根据权利要求1中所述的基于贝叶斯弹性网正则化方法的电流元三维反演方法,其特征在于，步骤6中，根据蒙特卡洛最大期望算法进行边缘似然的估计，计算需要求解的正则化参数，在每次的循环迭代过程中更新正则化参数值，在E步（Expectation）中，以当前参数计算对数似然的期望值，通过M步（Maximization）找到能使E步产生的似然期望最大化的参数值，再将参数代回E步中进行迭代计算，直到计算收敛，实际计算时采用交叉验证的方法确定式（4）中的两个正则化参数。
6.根据权利要求3中所述的基于贝叶斯弹性网正则化方法的电流元三维反演方法,其特征在于，所述公式（2）满足正态高斯分布，即磁场传感器测量误差满足正态高斯分布：公式（3）中，L为单位矩阵。
7.根据权利要求4中所述的基于贝叶斯弹性网正则化方法的电流元三维反演方法,其特征在于，基于弹性网正则化的分层贝叶斯模型为：。
8.根据权利要求7中所述的基于贝叶斯弹性网正则化方法的电流元三维反演方法,其特征在于，步骤5中，采用吉布斯采样的方法，对公式（4）进行改写：其中 是归一化常数，
在原分层模型的基础上增加一层，第一层为高斯分布，表示为 ，第二层为截断伽玛分布，表示为 ，相对应地，增加超参数τ，得到改进后的分层贝叶斯模型为 ：
基于以上分层模型，由下至上依次求解，得到似然函数，再根据公式（2）求解，最终得到待求电流元的分布情况，完成电流元的三维反演。
说明书 :
基于贝叶斯弹性网正则化方法的电流元三维反演方法
技术领域
[0001] 本发明涉及电气电磁反演应用领域，特别是一种基于贝叶斯弹性网正则化方法的电流元三维反演方法。
背景技术
[0002] 近年来，伴随着先进的传感和测量技术及其对应的控制方法，智能电网及能源互联网快速发展。智能电网的重要基石是对电力系统的全景实时状态信息进行深度感知。能源互联网作为智能电网的拓展，是“能量流、信息流、业务流”的高度融合，呈现出信息能源基础设施深度融合的趋势。智能电网中的电气信号具有频带宽、范围广、数据量大等特征，对实时信号的采集、传输、存储及分析成为智能电网及能源互联网的关键技术之一。
[0003] 针对电力系统中的磁场及电流，在不同的电压等级下及不同的电气设备中有不同的量测要求，从而对测量用传感器的参数要求各不相同，目前在电力系统中稳态电流的量测范围为5mA 6kA，小电流诸如泄漏电流的量测范围为100μA 10mA，大电流诸如暂态电流的~ ~量测范围为10kA 100kA，不同电流的幅值跨越8 9个数量级，频率由直流到50MHz左右，对应~ ~
测量磁场强度范围与之类似。
测量磁场强度范围与之类似。
[0004] 在测量及监控系统中，传感器作为监测各类电气状态信息量的主要工具，成为智能电网研究的重要方向。针对电网电流及磁场信号的特点及具体的量测需求，对传感器的要求包括测量准确度高、量程范围大、易于安装、成本低等，从而为电网状态信息的辨识及深度感知提供坚实的数据支撑。
[0005] 目前基于磁场信号量测的电流传感器包括霍尔传感器、各类巨磁阻传感器以及光电流互感器等，其中巨磁阻及隧穿磁阻传感器因其高灵敏度、宽频带、小温漂、高可集成度、小体积、廉价等特点，具备广阔的发展前景，是智能电网中磁场及电流监测的优选方案。如前所述，在电力系统应用场景中，磁场及电流的幅值和频率动态变化范围巨大，对基于巨磁电阻及隧穿磁阻效应的传感元件进行材料选择及结构设计，以实现各种特性参数的调整，适用智能电网广泛的测量需求。
[0006] 在电网实际量测应用场景中，存在各种直接测量电流、电压量较为困难的情况，需要对其产生的磁场、电场进行测量，并通过对所测场量的反演，实现复杂电磁环境下源的定位成像及空间场的分布重建。针对电力系统中的各类电磁现象，例如针对局部放电，目前已经开展了相应的反演研究，对于系统的状态感知及参数辨识具有重要意义。
[0007] 局部放电是指电气设备的绝缘部分区域在足够强的电场作用下局部范围内发生放电的现象，这种放电未形成固定的放电通道。局部放电将使设备的绝缘劣化，并引起系统的电能损失和无线电干扰，因此对局放进行反演的主要目的是实现故障识别及放电定位，发现故障点，从而保证电气设备运行可靠性。
[0008] 局部放电的测量主要包括侵入式和非侵入式测量两种。侵入式测量是在电路中串入耦合电容测量出现的局放信号，通过神经网络、模糊理论、支持向量机等高级算法，对信号进行去噪处理，并进行特征提取与识别，从而获得准确全面的故障识别反演结果。非侵入式测量主要是针对局放产生的声、光、电信号，包括特高频天线、光学传感器、超声波传感器等测量手段，主要实现局放的反演定位。特高频段电磁波信号在均匀介质中匀速传播，同一放电产生的电磁波由于传播距离的不同，使得不同位置处的传感器放电信号的到达时间存在差异。在采用特高频天线法测量时，通过设置传感器阵列，对阵列中不同传感器对同一放电脉冲测量的时间差进行计算，利用空间几何关系得到放电点坐标。在这种定位方法中，时间差直接影响了测量精度，最初的研究通过信号初始峰值获得时间差，但这种方法受到多次折反射、电磁干扰噪声的影响，因此准确度较低。采用能量累积法得到特高频传感器阵列中各测量点的时间差，实现了变压器内局部放电的定位，定位误差为cm级。采用脉冲峰值及相关系数法得到时间差，通过矩形及Y形两种传感器的阵列摆放方式，可以实现在三维空间内对变压器内局部放电及变电站现场15m范围内的局部放电定位通过小波去噪对信号进行处理得到时间差，利用二维双曲线定位模型获得局放信号位置。而采用光学传感器及超声波传感器进行测量时，基于光学传感器测量结果，通过信号模式识别及支持向量机的处理方法实现了对放电信号的检测及定位。通过超声传感器阵列搭建放电定位系统，根据多平台测向及全局搜索的方法有效提高了定位精度。
[0009] 电气工程其它领域的电磁反演，也已经有了一定研究基础。在工频电场反演问题的研究中，基于全局正则化建立测量电场与架空输电线等效电荷间的反演问题，采用阻尼高斯-牛顿法实现反演问题求解，通过有限的电场数据重建得到线路周围大范围区域的电场分布，有效减少了电场测量的工作量。在接地网故障的诊断中，通过测量接地网电阻抗、磁场分布、节点电压值等参数对接地网的故障位置进行反演，评估接地网工作状况及故障情况，类似方法应用于大地异常体的检测中。在变压器的状态评估中，测量变压器状态信息量，包括绝缘试验数据、电气试验数据等，通过一段时间内数据的变化情况，对变压器的状态进行反演，根据评估结果实现故障的检测。
[0010] 由磁场到磁场源的反演问题，通常具有病态、不适定性的特性，具体而言，其解具有Hadamard意义下的不存在性、不唯一性及不连续性。在实际问题中，引导场矩阵通常条件数较大，呈现病态性，同时通常选择重建的分布式电流元数目N远远大于测量磁场强度数目M，使得方程高度欠定（Underdetermined）。具体而言，电流元反演问题的不适定性主要表现于，由于磁场测量结果列向量H通常带有测量误差（噪声），较小的测量误差将对求解结果产生很大的影响，因此问题的解不具有连续性，且解通常不唯一。
发明内容
[0011] 本发明目的：设计一种基于贝叶斯弹性网正则化方法的电流元三维反演方法，通过空间磁场分布得到待求电流元分布。将待求解的线电流及体电流抽象为离散电流元的网格模型，针对空间中远多于传感器个数的待求解区域离散网格，通过正则化方法处理，在求解过程中引入先验信息，将其耦合到反演问题中，得到满足方程的近似解逼近真实电流元分布，使问题转化为切合实际要求的适定性问题。
[0012] 发明思路：在统计学中的一个主要学派是贝叶斯学派，在该学派的思想中，假设模型参数是随机变量，且依赖于先验概率，通过固定的观测数据以及先验概率信息得到参数的后验概率分布。
[0013] 本发明的目的是为了解决上述问题，设计了一种基于贝叶斯弹性 网正则化方法的电流元三维反演方法。具体设计方案为：一种基于贝叶斯弹性网正则化方法的电流元三维反演方法
电流元模型：
一段载有均匀同相的时变电流导线称为电流元，其直径远小于长 度，长度远小于波长及测量距离，设坐标系原点处有一段电流为I，长 度为l的电流元，方向向量为(m,n,p)，周围空间是均匀线性且各 向同性的媒质。当磁场频率较低(通常小于100MHz)，且在空气中进 行测量时，通常满足测量距离R远远小于磁场波长λ的条件，此时磁场 作为磁准静态(Magnetoquasistatic)场处理，空间中各处位移电流 密度远小于传导电流密度，全电流方程的表达形式简化为：
公式(8)中，JC为空间中的传导电流， 为空间中由电场变化 引起的位移电流。
电流元模型：
一段载有均匀同相的时变电流导线称为电流元，其直径远小于长 度，长度远小于波长及测量距离，设坐标系原点处有一段电流为I，长 度为l的电流元，方向向量为(m,n,p)，周围空间是均匀线性且各 向同性的媒质。当磁场频率较低(通常小于100MHz)，且在空气中进 行测量时，通常满足测量距离R远远小于磁场波长λ的条件，此时磁场 作为磁准静态(Magnetoquasistatic)场处理，空间中各处位移电流 密度远小于传导电流密度，全电流方程的表达形式简化为：
公式(8)中，JC为空间中的传导电流， 为空间中由电场变化 引起的位移电流。
[0014] 空气介质中电流元产生的空间磁场强度为：公式（9）中，r为测量场点矢量，r'为源点矢量，Idl'为线电流C区域内的每一段电流元，J(r')为与C等效的区域V内各源点处的体电流密度。
[0015] 于线电流I为标量，体电流J(r')为矢量，两者在公式（7）中等价。在源点r'处的体电流密度J(r')可以表示为附近区域内N个体积元素jk(t)的矢量和。在计算过程中，考虑体电流密度的时变特性，需要将J(r',t)进行变量分离，变为时间及空间乘积的形式，如公式：公式（10）中，jk(t)为不同点处的时变体电流密度，Φk(r')为与空间量相关的部分。
[0016] 将进行变量分离后的体电流密度J(r',t)代入公式（9），得到公式：从而磁场强度被分解为时间及空间量乘积的形式，将其中空间相关量的分项进行整理：
公式（11）、公式（12）整理成矩阵及标量形式，得到N处源点的体电流在M个测量点处产生的磁场强度为：
电流元反演的一般表示：
根据公式（13），分布式电流元的反演问题转化为根据磁场强度列向量H(t)对引导场矩阵G或电流值列向量J(t)进行求解。将电流求解区域进行离散处理，并整理引导场矩阵G，在磁场列向量H(t)及引导场矩阵G已知的情况下，设定合适的惩罚函数求解电流值列向量J(t)，从而得到整个区域内的电流分布。以上模型的推导过程中，采用了体电流的表达形式，由于体电流及线电流在毕奥-萨伐尔定律中可以相互转化，因此上述分析同样适用于线电流形式。
公式（11）、公式（12）整理成矩阵及标量形式，得到N处源点的体电流在M个测量点处产生的磁场强度为：
电流元反演的一般表示：
根据公式（13），分布式电流元的反演问题转化为根据磁场强度列向量H(t)对引导场矩阵G或电流值列向量J(t)进行求解。将电流求解区域进行离散处理，并整理引导场矩阵G，在磁场列向量H(t)及引导场矩阵G已知的情况下，设定合适的惩罚函数求解电流值列向量J(t)，从而得到整个区域内的电流分布。以上模型的推导过程中，采用了体电流的表达形式，由于体电流及线电流在毕奥-萨伐尔定律中可以相互转化，因此上述分析同样适用于线电流形式。
[0017] 根据公式（12）引导场矩阵G中的元素需要通过积分方式求得，过程过于复杂，因此对其进行线性化处理。设电流元起始位置坐标为（xs, ys, zs），方向向量为（m, n, p），长度为l，第i个测量点的坐标为（x0i, y0i, z0i），则根据公式（9）得到第i个测量点处的磁场强度矢量：公式（14）中，x = xs+dl·m，y = ys+dl·n，z = zs+dl·p
即：
在对求解区域及测量区域进行合理规划的基础上，得到与磁场测量结果对应的引导场矩阵G。由此，反演问题转化为一个线性方程的求解，即：针对公式（13）的矩阵形式，在已知磁场测量结果H(t)及与空间位置量相关的引导场矩阵G的情况下，求解对应的电流值列向量J(t)或者I(t)。
即：
在对求解区域及测量区域进行合理规划的基础上，得到与磁场测量结果对应的引导场矩阵G。由此，反演问题转化为一个线性方程的求解，即：针对公式（13）的矩阵形式，在已知磁场测量结果H(t)及与空间位置量相关的引导场矩阵G的情况下，求解对应的电流值列向量J(t)或者I(t)。
[0018] 一种基于贝叶斯弹性网正则化方法的电流元三维反演方法，通过待求解参数的先验分布信息及已知的磁场测量数据，进行空间体电流或线电流分布的反演重建，其特征在于，步骤1、将磁场传感器测量结果整理为磁场强度H向量；
步骤2、以贝叶斯定理表示待求电流元分布I的后验概率 ；
步骤3、引入超参数σ，由先验分布及超先验分布表示步骤2中后验概率分布的各个参数；
步骤4、按照贝叶斯弹性网正则化框架设计分层模型中的电流元先验概率分布，得到基于弹性网正则化的分层贝叶斯模型；
步骤5、根据贝叶斯公式进行后验概率分布的求解，采用吉布斯采样方法，根据所述分层贝叶斯模型计算得到后验概率分布值；
步骤6、根据已有磁场数据和待求离散电流元网格的位置坐标，进行上式中超参数及正则化参数的求解。
步骤2、以贝叶斯定理表示待求电流元分布I的后验概率 ；
步骤3、引入超参数σ，由先验分布及超先验分布表示步骤2中后验概率分布的各个参数；
步骤4、按照贝叶斯弹性网正则化框架设计分层模型中的电流元先验概率分布，得到基于弹性网正则化的分层贝叶斯模型；
步骤5、根据贝叶斯公式进行后验概率分布的求解，采用吉布斯采样方法，根据所述分层贝叶斯模型计算得到后验概率分布值；
步骤6、根据已有磁场数据和待求离散电流元网格的位置坐标，进行上式中超参数及正则化参数的求解。
[0019] 根据毕奥-萨伐尔定律，磁场可表示为电流元的积分形式，对于待求空间离散网格中的电流元，磁场与电流元间的关系可以矩阵形式表示。
[0020] 认为磁场传感器测量结果中引入了误差ε： （1）。
[0021] 其中H为磁场强度，I为待求电流元分布，G为引导场矩阵。
[0022] 在引导场矩阵G确定的情况下，根据贝叶斯定理，待求电流元分布I表示为：（2）
公式2中， 为是似然函数， 为已知电流元的先验概率分布，
为边缘似然，在此为常数。
公式2中， 为是似然函数， 为已知电流元的先验概率分布，
为边缘似然，在此为常数。
[0023] 所述步骤4中，基于弹性网正则化框架，对应电流元的先验概率分布为：公式（4）中，λ1，λ2为正则化系数。
[0024] 步骤6中，根据蒙特卡洛最大期望算法进行边缘似然的估计，计算需要求解的正则化参数，在每次的循环迭代过程中更新正则化参数值，在E步（Expectation）中，以当前参数计算对数似然的期望值，通过M步（Maximization）找到能使E步产生的似然期望最大化的参数值，再将参数代回E步中进行迭代计算，直到计算收敛。实际计算时采用交叉验证的方法确定式（4）中的两个正则化参数。
[0025] 所述公式(2)满足正态高斯分布，即磁场传感器测量误差满足正态 高斯分布：p(H|G,I)～N(GI,σ2L) (3)
公式(3)中，L为单位矩阵
基于弹性网正则化的分层贝叶斯模型为：
步骤5中，采用吉布斯采样的方法，对公式(4)进行改写：
其中C(λ1,λ2,σ2)是归一化常数，
在原分层模型的基础上增加一层，第一层为高斯分布，表示为 N(0,σ2(t-1)/(λ2t))，第二层为截断伽玛分布，表示为TG(1/2,8λ2σ2/λ12,(1,∞))。 相对应地，增加超参数τ，得到改进后的分层贝叶斯模型为：
基于以上分层模型，由下至上依次求解，得到似然函数，再根据公 式(2)求解，最终得到待求电流元的分布情况，完成电流元的三维反演。 通过本发明的上述技术方案得到的基于贝叶斯弹性网正则化方法的电 流元三维反演方法，其有益效果是：
贝叶斯正则化方法克服了离散电流元计算过程中的病态性，具有较 高的计算鲁棒性和准确度。另一方面，其参数选择的方式更为温和，可 以设定电流元的先验分布，输出参数有一定的置信区间。
公式(3)中，L为单位矩阵
基于弹性网正则化的分层贝叶斯模型为：
步骤5中，采用吉布斯采样的方法，对公式(4)进行改写：
其中C(λ1,λ2,σ2)是归一化常数，
在原分层模型的基础上增加一层，第一层为高斯分布，表示为 N(0,σ2(t-1)/(λ2t))，第二层为截断伽玛分布，表示为TG(1/2,8λ2σ2/λ12,(1,∞))。 相对应地，增加超参数τ，得到改进后的分层贝叶斯模型为：
基于以上分层模型，由下至上依次求解，得到似然函数，再根据公 式(2)求解，最终得到待求电流元的分布情况，完成电流元的三维反演。 通过本发明的上述技术方案得到的基于贝叶斯弹性网正则化方法的电 流元三维反演方法，其有益效果是：
贝叶斯正则化方法克服了离散电流元计算过程中的病态性，具有较 高的计算鲁棒性和准确度。另一方面，其参数选择的方式更为温和，可 以设定电流元的先验分布，输出参数有一定的置信区间。
附图说明
[0026] 图1是本发明所述电流元模型的示意图；图2是本发明所述基于贝叶斯弹性网正则化方法、Lasso正则化方法和Ridge正则化方法对应的先验概率分布；
图3是本发明所述分布式电流元反演计算模型结构示意；
图4是本发明所述分布式电流元反演计算模型中呈现电流元的形态示意图，其中：
图4（a）是本发明所述分布式电流元反演计算模型中呈现单个部分线电流元的形态示意图
图4（b）是本发明所述分布式电流元反演计算模型中呈现两个部分线电流元的形态示意图；
图4（c）是本发明所述分布式电流元反演计算模型中呈现单个体电流元的形态示意图；
图4（d）是本发明所述分布式电流元反演计算模型中呈现两个体电流元的形态示意图。
图3是本发明所述分布式电流元反演计算模型结构示意；
图4是本发明所述分布式电流元反演计算模型中呈现电流元的形态示意图，其中：
图4（a）是本发明所述分布式电流元反演计算模型中呈现单个部分线电流元的形态示意图
图4（b）是本发明所述分布式电流元反演计算模型中呈现两个部分线电流元的形态示意图；
图4（c）是本发明所述分布式电流元反演计算模型中呈现单个体电流元的形态示意图；
图4（d）是本发明所述分布式电流元反演计算模型中呈现两个体电流元的形态示意图。
[0027] 图5是本发明所述磁场强度误差±0.5%时单个部分线电流元的反演结果示意图，其中：图5（a）是本发明所述采用Ridge正则化方法的计算结果的单个部分线电流元的反演结果示意图；
图5（b）是本发明所述采用贝叶斯弹性网正则化方法的计算结果的单个部分线电流元的反演结果示意图；
图5（c）是本发明所述采用Lasso正则化方法的计算结果的单个部分线电流元的反演结果示意图；
图6是本发明所述磁场强度误差±0.5%时两个部分线电流元的反演结果示意图，其中：
图6（a）是本发明所述采用Ridge正则化方法的计算结果的两个部分线电流元的反演结果示意图；
图6（b）是本发明所述采用贝叶斯弹性网正则化方法的计算结果的两个部分线电流元的反演结果示意图；
图6（c）是本发明所述采用Lasso正则化方法的计算结果的两个部分线电流元的反演结果示意图；
图7是本发明所述磁场强度误差±0.5%时单个体电流元的反演结果示意图，其中：
图7（a）是本发明所述采用Ridge正则化方法的计算结果的单个体电流元的反演结果示意图；
图7（b）是本发明所述采用贝叶斯弹性网正则化方法的计算结果的单个体电流元的反演结果示意图；
图7（c）是本发明所述采用Lasso正则化方法的计算结果的单个体电流元的反演结果示意图；
图8是本发明所述磁场强度误差±0.5%时两个体电流元的反演结果示意图，其中：
图8（a）是本发明所述采用Ridge正则化方法的计算结果的两个体电流元的反演结果示意图；
图8（b）是本发明所述采用贝叶斯弹性网正则化方法的计算结果的两个体电流元的反演结果示意图；
图8（c）是本发明所述采用Lasso正则化方法的计算结果的两个体电流元的反演结果示意图。
图5（b）是本发明所述采用贝叶斯弹性网正则化方法的计算结果的单个部分线电流元的反演结果示意图；
图5（c）是本发明所述采用Lasso正则化方法的计算结果的单个部分线电流元的反演结果示意图；
图6是本发明所述磁场强度误差±0.5%时两个部分线电流元的反演结果示意图，其中：
图6（a）是本发明所述采用Ridge正则化方法的计算结果的两个部分线电流元的反演结果示意图；
图6（b）是本发明所述采用贝叶斯弹性网正则化方法的计算结果的两个部分线电流元的反演结果示意图；
图6（c）是本发明所述采用Lasso正则化方法的计算结果的两个部分线电流元的反演结果示意图；
图7是本发明所述磁场强度误差±0.5%时单个体电流元的反演结果示意图，其中：
图7（a）是本发明所述采用Ridge正则化方法的计算结果的单个体电流元的反演结果示意图；
图7（b）是本发明所述采用贝叶斯弹性网正则化方法的计算结果的单个体电流元的反演结果示意图；
图7（c）是本发明所述采用Lasso正则化方法的计算结果的单个体电流元的反演结果示意图；
图8是本发明所述磁场强度误差±0.5%时两个体电流元的反演结果示意图，其中：
图8（a）是本发明所述采用Ridge正则化方法的计算结果的两个体电流元的反演结果示意图；
图8（b）是本发明所述采用贝叶斯弹性网正则化方法的计算结果的两个体电流元的反演结果示意图；
图8（c）是本发明所述采用Lasso正则化方法的计算结果的两个体电流元的反演结果示意图。
具体实施方式
[0028] 下面结合附图对本发明进行具体描述。
[0029] 一种基于贝叶斯弹性网正则化方法的电流元三维反演方法电流元模型：
图1是本发明所述电流元模型的示意图，如图1所示，一段载有均匀同相的时变电流导线称为电流元，其直径远小于长度，长度远小于波长及测量距离，设坐标系原点处有一段电流为I，长度为l的电流元，方向向量为（m, n, p），周围空间是均匀线性且各向同性的媒质。
当磁场频率较低（通常小于100MHz），且在空气中进行测量时，通常满足测量距离R远远小于磁场波长λ的条件，此时磁场作为磁准静态（Magnetoquasistatic）场处理，空间中各处位移电流密度远小于传导电流密度，全电流方程的表达形式简化为：
（8）。
图1是本发明所述电流元模型的示意图，如图1所示，一段载有均匀同相的时变电流导线称为电流元，其直径远小于长度，长度远小于波长及测量距离，设坐标系原点处有一段电流为I，长度为l的电流元，方向向量为（m, n, p），周围空间是均匀线性且各向同性的媒质。
当磁场频率较低（通常小于100MHz），且在空气中进行测量时，通常满足测量距离R远远小于磁场波长λ的条件，此时磁场作为磁准静态（Magnetoquasistatic）场处理，空间中各处位移电流密度远小于传导电流密度，全电流方程的表达形式简化为：
（8）。
[0030] 公式（8）中，JC为空间中的传导电流， 为空间中由电场变化引起的位移电流。
[0031] 空气介质中电流元产生的空间磁场强度为： （9）
公式（9）中，r为测量场点矢量，r'为源点矢量，Idl'为线电流C区域内的每一段电流元，J(r')为与C等效的区域V内各源点处的体电流密度。
公式（9）中，r为测量场点矢量，r'为源点矢量，Idl'为线电流C区域内的每一段电流元，J(r')为与C等效的区域V内各源点处的体电流密度。
[0032] 于线电流I为标量，体电流J(r')为矢量，两者在公式（7）中等价。在源点r'处的体电流密度J(r')可以表示为附近区域内N个体积元素jk(t)的矢量和。在计算过程中，考虑体电流密度的时变特性，需要将J(r',t)进行变量分离，变为时间及空间乘积的形式，如公式： （10）
公式（10）中，jk(t)为不同点处的时变体电流密度，Φk(r')为与空间量相关的部分。
公式（10）中，jk(t)为不同点处的时变体电流密度，Φk(r')为与空间量相关的部分。
[0033] 将进行变量分离后的体电流密度J(r',t)代入公式（9），得到公式：（11）
从而磁场强度被分解为时间及空间量乘积的形式，将其中空间相关量的分项进行整理：
（12）
公式（11）、公式（12）整理成矩阵及标量形式，得到N处源点的体电流在M个测量点处产生的磁场强度为：
（13）
电流元反演的一般表示：
根据公式（13），分布式电流元的反演问题转化为根据磁场强度列向量H(t)对引导场矩阵G或电流值列向量J(t)进行求解。将电流求解区域进行离散处理，并整理引导场矩阵G，在磁场列向量H(t)及引导场矩阵G已知的情况下，设定合适的惩罚函数求解电流值列向量J(t)，从而得到整个区域内的电流分布。以上模型的推导过程中，采用了体电流的表达形式，由于体电流及线电流在毕奥-萨伐尔定律中可以相互转化，因此上述分析同样适用于线电流形式。
从而磁场强度被分解为时间及空间量乘积的形式，将其中空间相关量的分项进行整理：
（12）
公式（11）、公式（12）整理成矩阵及标量形式，得到N处源点的体电流在M个测量点处产生的磁场强度为：
（13）
电流元反演的一般表示：
根据公式（13），分布式电流元的反演问题转化为根据磁场强度列向量H(t)对引导场矩阵G或电流值列向量J(t)进行求解。将电流求解区域进行离散处理，并整理引导场矩阵G，在磁场列向量H(t)及引导场矩阵G已知的情况下，设定合适的惩罚函数求解电流值列向量J(t)，从而得到整个区域内的电流分布。以上模型的推导过程中，采用了体电流的表达形式，由于体电流及线电流在毕奥-萨伐尔定律中可以相互转化，因此上述分析同样适用于线电流形式。
[0034] 根据公式（12）引导场矩阵G中的元素需要通过积分方式求得，过程过于复杂，因此对其进行线性化处理。设电流元起始位置坐标为（xs, ys, zs），方向向量为（m, n, p），长度为l，第i个测量点的坐标为（x0i, y0i, z0i），则根据公式（9）得到第i个测量点处的磁场强度矢量： （14）
公式（14）中，x = xs+dl·m，y = ys+dl·n，z = zs+dl·p
即：
（15）
在对求解区域及测量区域进行合理规划的基础上，得到与磁场测量结果对应的引导场矩阵G。由此，反演问题转化为一个线性方程的求解，即：针对公式（13）的矩阵形式，在已知磁场测量结果H(t)及与空间位置量相关的引导场矩阵G的情况下，求解对应的电流值列向量J(t)或者I(t)。
公式（14）中，x = xs+dl·m，y = ys+dl·n，z = zs+dl·p
即：
（15）
在对求解区域及测量区域进行合理规划的基础上，得到与磁场测量结果对应的引导场矩阵G。由此，反演问题转化为一个线性方程的求解，即：针对公式（13）的矩阵形式，在已知磁场测量结果H(t)及与空间位置量相关的引导场矩阵G的情况下，求解对应的电流值列向量J(t)或者I(t)。
[0035] 通过待求解参数的先验分布信息及已知的磁场测量数据，进行空间体电流或线电流分布的反演重建，其特征在于，步骤1、将磁场传感器测量结果整理为磁场强度H向量；
步骤2、以贝叶斯定理表示待求电流元分布I的后验概率 ；
步骤3、引入超参数σ，由先验分布及超先验分布表示步骤2中后验概率分布的各个参数；
步骤4、按照贝叶斯弹性网正则化框架设计分层模型中的电流元先验概率分布，得到基于弹性网正则化的分层贝叶斯模型；
步骤5、根据贝叶斯公式进行后验概率分布的求解，采用吉布斯采样方法，根据所述分层贝叶斯模型计算得到后验概率分布值；
步骤6、根据已有磁场数据和待求离散电流元网格的位置坐标，进行上式中超参数及正则化参数的求解。
步骤2、以贝叶斯定理表示待求电流元分布I的后验概率 ；
步骤3、引入超参数σ，由先验分布及超先验分布表示步骤2中后验概率分布的各个参数；
步骤4、按照贝叶斯弹性网正则化框架设计分层模型中的电流元先验概率分布，得到基于弹性网正则化的分层贝叶斯模型；
步骤5、根据贝叶斯公式进行后验概率分布的求解，采用吉布斯采样方法，根据所述分层贝叶斯模型计算得到后验概率分布值；
步骤6、根据已有磁场数据和待求离散电流元网格的位置坐标，进行上式中超参数及正则化参数的求解。
[0036] 根据毕奥-萨伐尔定律，磁场可表示为电流元的积分形式，对于待求空间离散网格中的电流元，磁场与电流元间的关系可以矩阵形式表示。
[0037] 认为磁场传感器测量结果中引入了误差ε： （1）。
[0038] 其中H为磁场强度，I为待求电流元分布，G为引导场矩阵。
[0039] 在引导场矩阵G确定的情况下，根据贝叶斯定理，待求电流元分布I表示为：（2）
公式2中， 为是似然函数， 为已知电流元的先验概率分布，
为边缘似然，在此为常数。
公式2中， 为是似然函数， 为已知电流元的先验概率分布，
为边缘似然，在此为常数。
[0040] 所述步骤4中，基于弹性网正则化框架，对应电流元的先验概率分布为： （4）
公式（4）中，λ1，λ2为正则化系数。
公式（4）中，λ1，λ2为正则化系数。
[0041] 步骤6中，根据蒙特卡洛最大期望算法进行边缘似然的估计，计算需要求解的正则化参数，在每次的循环迭代过程中更新正则化参数值，在E步（Expectation）中，以当前参数计算对数似然的期望值，通过M步（Maximization）找到能使E步产生的似然期望最大化的参数值，再将参数代回E步中进行迭代计算，直到计算收敛。实际计算时采用交叉验证的方法确定式（4）中的两个正则化参数。
[0042] 所述公式（2）满足正态高斯分布，即磁场传感器测量误差满足正态高斯分布： （3）
公式（3）中，L为单位矩阵
基于弹性网正则化的分层贝叶斯模型为：
（5）。
公式（3）中，L为单位矩阵
基于弹性网正则化的分层贝叶斯模型为：
（5）。
[0043] 步骤5中，采用吉布斯采样的方法，对公式(4)进行改写：其中C(λ1,λ2,σ2)是归一化常数，
在原分层模型的基础上增加一层，第一层为高斯分布，表示为 N(0,σ2(t-1)/(λ2t))，第二层为截断伽玛分布，表示为TG(1/2,8λ2σ2/λ12,(1,∞))。 相对应地，增加超参数τ，得到改进后的分层贝叶斯模型为：
基于以上分层模型，由下至上依次求解，得到似然函数，再根据公 式(2)求解，最终得到待求电流元的分布情况，完成电流元的三维反演。
实施例1
图3是本发明所述分布式电流元反演计算模型结构示意，如图3 所示，根据上述算法推导，对基于贝叶斯弹性网正则化方法的电流元三 维成像结果进行分析，构建如图3所示结构的计算模型。其中点代表三 维磁场强度的测量点，三种测量点所在平面分别平行于xy，yz和xz平 面，每个平面的大小为20m×20m，同一平面的测量点间距为1m，测量点 总数为
264个。图中中心立方体区域为计算的电流元区域，该区域边长 为4.5m，并设定该区域内每个独立的电流元为边长为0.5m的立方体，从 而在该区域内共有729个待求电流元值。该模型中，通过264个三维磁 场强度得到729个分布式电流元的值及其方向，同时对空间磁场进行重 建，是典型的不适定及欠定问题。设置电流元立方体与测量距离的比值 小于0.1，可以采用上述近似线性模型求解。
在原分层模型的基础上增加一层，第一层为高斯分布，表示为 N(0,σ2(t-1)/(λ2t))，第二层为截断伽玛分布，表示为TG(1/2,8λ2σ2/λ12,(1,∞))。 相对应地，增加超参数τ，得到改进后的分层贝叶斯模型为：
基于以上分层模型，由下至上依次求解，得到似然函数，再根据公 式(2)求解，最终得到待求电流元的分布情况，完成电流元的三维反演。
实施例1
图3是本发明所述分布式电流元反演计算模型结构示意，如图3 所示，根据上述算法推导，对基于贝叶斯弹性网正则化方法的电流元三 维成像结果进行分析，构建如图3所示结构的计算模型。其中点代表三 维磁场强度的测量点，三种测量点所在平面分别平行于xy，yz和xz平 面，每个平面的大小为20m×20m，同一平面的测量点间距为1m，测量点 总数为
264个。图中中心立方体区域为计算的电流元区域，该区域边长 为4.5m，并设定该区域内每个独立的电流元为边长为0.5m的立方体，从 而在该区域内共有729个待求电流元值。该模型中，通过264个三维磁 场强度得到729个分布式电流元的值及其方向，同时对空间磁场进行重 建，是典型的不适定及欠定问题。设置电流元立方体与测量距离的比值 小于0.1，可以采用上述近似线性模型求解。
[0044] 图4是本发明所述分布式电流元反演计算模型中呈现电流元的形态示意图，如图4所示，立方体中的电流元形态示意如图3所示，设置四种电流元分布形式，在图中以矢量线段表示。其中，图4(a)，图4(b)中电流元呈现部分线电流的形态，用于模拟非长直线电流形式，分别以情况1和情况2表示，图4(c)，图4(d)中电流元呈现簇状形态，用于模拟不规则的体电流形式，分别以情况3和情况4表示。在每种分布形式的计算中分别采用Ridge，弹性网和Lasso正则化方法进行电流元分布的三维反演成像计算，分析不同误差等级下反演计算得到分布与真实分布的差异，同时计算磁场分布的重建误差。
[0045] 实施例2在实施例1的基础上，对于情况1中的电流元分布，当磁场强度理论值分别叠加±0.5%之间的随机相对误差时，在每种分布形式的计算中分别采用Ridge，弹性网和Lasso正则化方法进行电流元分布的三维反演成像计算，分析不同误差等级下反演计算得到分布与真实分布的差异，同时计算磁场分布的重建误差，采用三种正则化方法的计算结果如图5所示。
[0046] 实施例3在实施例1的基础上，对情况2在磁场强度理论值叠加±0.5%之间的随机相对误差时，在每种分布形式的计算中分别采用Ridge，弹性网和Lasso正则化方法进行电流元分布的三维反演成像计算，分析不同误差等级下反演计算得到分布与真实分布的差异，同时计算磁场分布的重建误差，采用三种正则化方法的计算结果如图6所示。
[0047] 实施例4在实施例1的基础上，对情况3在磁场强度理论值叠加±0.5%之间的随机相对误差时，在每种分布形式的计算中分别采用Ridge，弹性网和Lasso正则化方法进行电流元分布的三维反演成像计算，分析不同误差等级下反演计算得到分布与真实分布的差异，同时计算磁场分布的重建误差，采用三种正则化方法的计算结果如图7所示。
[0048] 实施例4在实施例1的基础上，对情况4在磁场强度理论值叠加±0.5%之间的随机相对误差时，在每种分布形式的计算中分别采用Ridge，弹性网和Lasso正则化方法进行电流元分布的三维反演成像计算，分析不同误差等级下反演计算得到分布与真实分布的差异，同时计算磁场分布的重建误差，采用三种正则化方法的计算结果如图8所示。
[0049] 对比例1对实施例1-4的计算结果进行横向对比：
每种情况下，Ridge正则化方法得到的分布结果最为密集，电流元的方向性较好，Lasso正则化方法得到分布结果最为稀疏，在电流元定位方面的效果较好，弹性网正则化方法则具备两种方法的优点。相对于代数正则化方法，贝叶斯Ridge正则化方法的结果中有效电流元个数明显减小，而Lasso正则化方法与代数正则化方法的定位结果基本一致，但整体分布较为平均，从而使得结果更为温和。
每种情况下，Ridge正则化方法得到的分布结果最为密集，电流元的方向性较好，Lasso正则化方法得到分布结果最为稀疏，在电流元定位方面的效果较好，弹性网正则化方法则具备两种方法的优点。相对于代数正则化方法，贝叶斯Ridge正则化方法的结果中有效电流元个数明显减小，而Lasso正则化方法与代数正则化方法的定位结果基本一致，但整体分布较为平均，从而使得结果更为温和。
[0050] 对比例2对实施例1-4的计算结果进行综合对比：
弹性网正则化方法仍然在电流元分布的计算方面取得了较好的效果，对于情况2，贝叶斯正则化方法的结果在定位方面要优于代数正则化在此种情况下的求解结果。采用贝叶斯正则化方法较代数正则化方法可以将求解权重更均匀地分配至各分布电流元，更符合实际求解需求。在实际求解过程中，可以根据各方法的特点，综合使用各类正则化方法，根据弹性网正则化方法计算源的位置，结合Ridge正则化方法求解源的方向，从而得到更为准确的电流元分布。
弹性网正则化方法仍然在电流元分布的计算方面取得了较好的效果，对于情况2，贝叶斯正则化方法的结果在定位方面要优于代数正则化在此种情况下的求解结果。采用贝叶斯正则化方法较代数正则化方法可以将求解权重更均匀地分配至各分布电流元，更符合实际求解需求。在实际求解过程中，可以根据各方法的特点，综合使用各类正则化方法，根据弹性网正则化方法计算源的位置，结合Ridge正则化方法求解源的方向，从而得到更为准确的电流元分布。
[0051] 实施例5在实施例1-4的基础上，即在对电流元分布的反演结果进行分析的基础上，进一步对该电流元周围磁场分布进行重建计算，并与真实的磁场强度分布结果进行对比，分析反演算法在磁场重建计算方面的可行性。
[0052] 设定磁场重建计算区域为中心位于坐标原点的28m×28m×28m的立方体表面，计算点均匀分布，共384个。采用电流元反演结果计算各点处三维磁场合成值，得到其与真实值间的误差如表1、表2所示。
[0053] 其中表1表示情况1在三种误差水平下，采用三种正则化方法进行磁场重建的平均计算误差，表2表示在采用弹性网正则化方法对四种电流元分布情况进行重建时，在三种误差水平下磁场重建的平均计算误差。
[0054] 由表1表2可知，贝叶斯正则化方法在磁场重建计算方面仍有较高的准确度和稳定性，当叠加误差水平发生变化时，同一种正则化方法下的重建计算误差增大不明显，可以认为这是由于在贝叶斯正则化方法中求取最大后验概率分布值，使得在求解过程中对测量误差的敏感性变小而造成。而在表3.9的四种电流元分布情况中，重建计算误差均较小，情况2和情况4分别相对于情况1和情况3而言，计算误差增大，表明电流元模型的复杂程度将增大计算误差。由此，在磁场的重建计算方面，代数正则化和贝叶斯正则化方法效果类似，均能得到准确和稳定的计算结果。考虑到实际测量中误差的复杂性和不可预测性，在计算中优先采用贝叶斯正则化方法。
[0055] 上述技术方案仅体现了本发明技术方案的优选技术方案，本技术领域的技术人员对其中某些部分所可能做出的一些变动均体现了本发明的原理，属于本发明的保护范围之内。