0 引言
不同于交流输电系统，直流输电系统中单极性直流电场的方向不随时间变化，使得电场中的自由电荷沿电场线方向定向移动并积聚在绝缘子表面，造成局部电场畸变，严重时甚至会引发绝缘子沿面闪络，大大降低绝缘子的耐受性能[1-4]。为了抑制直流电场下绝缘子表面电荷积聚，提高直流电气设备的绝缘性能，有必要弄清楚直流电场下电荷的聚散行为及形成机制。因此，首先要对绝缘子表面电荷分布进行准确的测量。
目前表面电荷分布的主要测量手段有“粉尘图”法、基于Pockels 效应的光学测量法及无源/有源静电探头法[5-7]。其中，有源静电探头法因其测量精度高、稳定性高和对测量对象要求较低的优势，是目前应用最为广泛的一种方法[8-10]。但是，静电探头法只能直接测量绝缘子表面的电位分布，要得到电荷密度分布，还需要对所测的表面电位分布与电荷密度分布之间的关系进行推演，即进行表面电荷的反演计算。
D.K.Davies 最初于1967 年提出静电探头法的同时也提出了线性标度法对表面电荷进行换算[11]，认为探头所测电压与待测表面实际电荷密度分布呈线性关系[11-13]。而事实上，对于一个测量点数为N的测量系统，每一个微元内的电位均是由绝缘子表面电荷在该微元内的线性叠加所得。可见，电位与电荷之间并不是简单的线性关系，而需要通过电位-电荷转换矩阵H 来进行表达[14]。基于此，λ 函数法[15]、φ 函数法[16]以及视在电荷法[17]被提出。但是，这些方法均需要进行大规模的矩阵运算或对矩阵进行直接求逆，不仅消耗大量的计算资源且计算精度较差。此外，研究人员通常会通过增加测量点数来提高表面电荷的测量精度，这将导致电位-电荷转换矩阵具有非常高的维数（N×N）。而大规模矩阵的求逆往往伴随着不适定问题，即使转换矩阵或电位矩阵中存在着很小的扰动，在反演计算中也会引入较大的误差[18]。
针对这一问题，研究者们采用二维傅里叶变换将空间域中向量的卷积运算转换为频域中的乘积运算，并构建维纳滤波器抑制噪声干扰，计算表面电荷密度分布[19-20]。潘子君等在结合二维傅里叶变换的基础上引入了约束最小二乘方滤波器来抑制噪声干扰，并通过迭代手段自适应获取滤波系数的最优解[21]。这些方法有效地改善了表面电荷反演计算的精度和稳定性，但仅适用于平板或圆柱形绝缘子这类“平移不变系统”。而盆式或锥形绝缘子属于“平移变化系统”，激励施加的位置会影响该系统的响应[22]，因此对于这类“平移变化系统”，上述方法不再适用。针对这一问题，研究人员采用基于标准Tikhonov 正则化的维纳滤波器对噪声进行抑制，减弱矩阵的病态性。但维纳滤波器的滤波系数需要通过预先实验测定噪声和信号的功率密度谱来获得，而在实际情况中它们很难通过实验得到[23]。因此，有研究者提出通过点扩散函数（Point Spread Function,PSF）理论，在一个微元内设置单位电荷后，仿真计算不同正则化参数下该微元内表面电荷的反演结果，借助Gaussian 分布拟合及傅里叶变换得到电荷在空间频域下的分布特性后，再根据幅频特性的截止频率选择正则化参数，由此得到在协方差矩阵HTH 最大特征值θmax 的0.03%～1%之间的常数均可作为正则化参数的值[18,24-25]。但为了获取更好的表面电荷反演结果，该方法在实际使用时往往需要人为地多次调节参数，选择最优的空间分辨率，其参数选择标准及反演精度都难以控制。此外，标准吉洪诺夫正则化方法对于解决较大规模的不适定问题而言效率较低，而常用的解决大规模不适定问题的迭代正则化方法在实际应用中由于迭代终点的不精确选择，往往存在半收敛问题[26]。因此，有必要对正则化方法进行改进，研究正则化参数的自适应选择，提高算法效率及精度。
综上所述，本文考虑迭代方法和Tikhonov 方法的优点，引入迭代正则化技术，提出了基于Lanczos 双对角化（Lanczos Bidiagonalization,LBD）和Tikhonov 正则化算法混合方法（以下简称混合Lanczos-Tikhonvo 算法）的绝缘子表面电荷反演算法，并采用自适应加权广义交叉验证法（Adaptive Weighted Generalized Cross Validation,A-WGCV）自动选择正则化参数，完成绝缘子表面电荷密度分布反演。本文首先阐释了该混合方法的计算原理、正则化参数的选择方法及迭代停止条件；然后与视在电荷法、维纳滤波反演算法的计算精度进行比较；同时结合实验验证了算法的准确性和可信度；最后通过该方法获取了针电极下施加不同电压时的绝缘子表面电荷分布，进一步验证了本算法的可靠性。
1 表面电荷测量
1.1 表面电位测量平台
本文以图 1 所示的具有旋转对称结构的绝缘子为研究对象。该绝缘子以Al2O3 为填料掺杂环氧树脂，根据工业生产的浇注工艺制作。实验时，可将绝缘子中心金属柱连接高压电极，外圈环形电极接地。
图1 绝缘子尺寸及剖面图
Fig.1 Dimensions of insulator used in experiment
将表面清理干净的待测绝缘子安装于图2 所示的实验系统中进行实验，实验过程中外部罐体封闭，将罐内抽成真空后再充入定量气体可以得到所需的稳定气体环境。实验所用电源为大连泰思曼科技有限公司TCM6002P60-150-N 型号的直流高压电源，输出电压范围为-60～60 kV。直流高压电源的输出电压则通过罐体内与套管相连的可升降金属导杆施加于绝缘子。测量绝缘子表面电位时，关闭直流电源并升高金属导杆，夹具上的测量探头在预置程序的控制下，由位于罐体附腔处的两轴三维运动装置带动，移动至绝缘子表面上方3 mm 处，沿绝缘子中心向外的径向方向进行移动测量。探头每向外移动一次，绝缘子在与中心柱相连的外部电机带动下 旋转一圈，协同完成绝缘子表面电位的测量。实验所用测量探头为Trek-341B 型静电位计的配套高压探头，探头距离待测绝缘子表面2～4 mm 时测量精度可达到满量程的±0.1%以内。
图2 绝缘子表面电位测量系统
Fig.2 Insulator surface potential measurement system
1.2 表面电位与表面电荷密度分布的关系
目前，研究者们在进行表面电荷反演计算时普遍认为绝缘子表面电位为待测表面所有电荷共同作用的结果[14-17]。如果将待测表面分为N 个相同的小单元，那么整个待测表面的电位φ 和电荷密度σ 的关系可以表示为[17,27-28]
式中，hij 为表面第j 个测量单元中的电荷密度对第i 个测量单元电位值的贡献，与电场和电荷条件无关，由采用模型的几何结构及材料决定，因此当模型确定后即可通过计算得到。
从式（2）可以看出，当σN×1 中的σi=1 μC/m2，其他元素都为0 时，有
式中，H i为H 矩阵的第i 列，因此，此时的表面电位值就是H 矩阵的第i 列。依照上述思路，可依次将计算模型的第1～N 个测量单元的电荷密度设为1 μC/m2，其他单元的电荷密度设为0，分别计算表面电位值以得到H 矩阵的1～N 列。设绝缘子表面沿径向分为r 圈，每圈分为相等的d 个小单元。由于本文选择的研究对象具有旋转对称结构，对于每一圈测量单元，当对第一个测量单元设置单位电荷密度并计算其电位分布后，其余单元格对应的电位分布通过对第一次结果进行旋转位移即可得到。
因此，若将绝缘子计算模型沿表面自内而外等分为50 圈，每圈每1°设置一个单元格，共360 个小单元，即绝缘子表面被分为18 000（50×360）个测量单元。选择沿径向的一列单元格（50 个）依次设置单位电荷密度，每设置一次后计算所得的电位分布即分别为H 矩阵的第1+360(r-1)列，其中r =1,2,…,50。根据旋转对称性，对各列数据进行旋转位移后得到完整的H 矩阵。通过实验测得绝缘子表面电位分布φ 后，即可根据式（1）计算表面电荷分布。
以图3a 所示划分测量单元，并依前述步骤将50次电位计算结果导出，如图3b 所示，再分别旋转位移179 次后得到H 矩阵，此时H 矩阵是一个规模为18 000×18 000 的大型方阵，其条件数为
图3 计算模型的测量单元划分
Fig.3 The measurement unit division of the calculation model
因此，H 矩阵为病态矩阵，即在求解式（1）时，即使H 矩阵和测量的电位数据中仅存在极小的扰动，也会使计算出的电荷密度分布产生较大误差[18]。而在实际实验中，这样的小扰动通常是不可避免的，因此无法直接对H 求逆来完成表面电荷反演计算。
2 基于混合Lanczos-Tikhonvo 算法的表面电荷反演算法
对于式（1）所示问题，引入标准Tikhonov 正则化方法，即等价为求解最小二乘问题[29]
将式（5）展开，即求解方程
得到正则化解为
式中，λ 为正则化参数，该参数的选择对反演计算的精度有着重要影响； σλ 为正则化解，是边界约束和残差范数的加权组合的最小值；I 为单位矩阵。
如式（7）所示，采用标准Tikhonov 正则化方法求解时，虽然矩阵的不适定问题得到缓解，但仍需进行18 000×18 000 维度的矩阵求逆运算。因此，标准Tikhonov 正则化的性能受到了一定的限制。
2.1 迭代LBD 与Tikhonov 正则化算法混合方法
基于以上分析，本文引入迭代LBD 与Tikhonov正则化算法的混合方法进行表面电荷密度分布的反演。利用LBD 方法将大规模问题投影到低维子空间中，然后对低维问题使用Tikhonov 正则化方法求解，再将低维子空间的问题解进行一定变换即可得到原大规模问题的解[26]，这一过程可大大降低矩阵的计算量，提高反演计算的效率。
不同于标准Tikhonov 正则化，迭代正则化方法需要解决的最小二乘问题为
对H 矩阵进行迭代LBD 变换，对于k（k ≪N）维的两个正交子空间，定义其正交基分别为P 和Y，则有
式中，B k为一个低维双对角矩阵，表示为
在迭代计算B 矩阵的过程中，各参数需满足[30]
将迭代法则以矩阵形式表示，有
式中，e1 为单位向量；ek+1 为k+1 维单位矩阵的最后一列；θk+1 为矩阵Bk+1 的第k+1 个对角线元素；Pk、Yk 分别为
每次迭代计算后，引入Tikhonov 正则化对所得更低维度的矩阵Bk 进行反问题求解。
经k 步迭代后，原问题的残差表示为
式中，f 为子空间问题的解；Pk+1 为标准正交矩阵。因此，原本的最小二乘问题近似为
式中，σ 为矩阵Yk 的列向量。
由式（15）可以发现，LBD 的每一次迭代都需要解决一个涉及双对角矩阵Bk 的最小二乘问题。而Bk 的维度相比H 矩阵而言小得多，因此可以在每次迭代中使用Tikhonov 正则化方法。根据式（15）计算得到每次迭代后的正则化解为
2.2 A-WGCV 正则化参数选择
正则化方法的关键在于正则化参数的选择，但目前对该参数优劣的评价并无统一的标准，这也成为了限制正则化方法性能的另一个重要因素。
广义交叉验证（Generalized Cross Validation,GCV）方法为目前Tikhonov 正则化方法常用的正则化参数选择方法之一，该方法认为一个好的正则化参数λ 应该可以预测一组数据的缺失值[31]。因此，将初始数据分为两个部分，一部分用以计算近似解，另一部分验证计算结果。令φ(i)为φ 去掉第i 个元素后得到的向量；H(i)为矩阵H 去掉第i 行之后所得的矩阵；H(i,:)为矩阵H 第i 行；σλ(i)是将φ(i)和H(i)作为已知量，正则化参数取λ 时的计算解，那么被删去的元素φi 应可以通过H(i,:)σλ(i)估计得到。而所寻找的λ 应满足所有元素的期望误差最小的要求，即
依次忽略测量电位向量φ 的每一个元素，并根据式（17）寻找使期望误差最小的λ 值，这一过程可用GCV 函数表示为
式中，trace(A)为矩阵 A 的迹。因此正则化解为σλ=
。将H 矩阵的奇异值分解（Singular Value Decomposition,SVD）H =UΣV T代入式（18）即可得到
式中，U、V 矩阵中的列向量分别为H 矩阵的左、右奇异向量；ui 为矩阵U 的第i 列；bi 为H 矩阵的奇异值；Σ 为一个对角线元素均为H 的奇异值的对角线矩阵。
对于进行LBD 迭代的H 矩阵，矩阵Bk 的SVD为
式中，Rk+1、Wk 分别为包含Bk 左、右奇异向量的正交矩阵；Δk=diag(δ1, δ2,…, δk)，其对角线元素δi 为矩阵Bk 的奇异值，并按照δ1≥δ2≥…≥δk≥0 的顺序排列。那么此时GCV 函数应写作
理想情况下，随着迭代的进行，解的误差会趋于稳定，但实际情况并不总是如此，在混合Lanczos-Tikhonvo 算法中运用GCV 时，其计算得到的正则化参数总是存在过度平滑解的问题[19]，因此本文采用A-WGCV 方法进行正则化参数的选择。在原本GCV 方法的基础上，引入一个权重参数ω。对于投影后的问题，加权广义交叉验证法（Weighted Generalized Cross Validation,WGCV）函数写作
从式（22）可以看出，ω 的选择是混合方法成功的关键。令 ηk,opt为第k 次迭代时的最优正则化参数，假设在第k 次迭代时，已知 ηk,opt，那么就可以通过最小化WGCV 函数求得ω，即
但实际上ηk,opt为未知数，因此可以计算ηk,opt=δmin (Bk)时对应的
（δmin (Bk)为矩阵 B k最小的奇异值）。但当k 较大时，由于矩阵的病态性，δmin (Bk)更趋近于0，会导致这一方法计算出不合适的参数。因此，选择 
，用问题中早期的良好条件成分来帮助稳定较小奇异值带来的有害影响。如此，即可自适应地选择最合适的ω。
2.3 迭代停止准则
迭代停止的判断准则也是混合方法求解中至关重要的一步。A-WGCV 要求迭代停止时需要满足式（21）最小，在理想条件下，此时求得的正则化参数ηk 应收敛到了一个定值，那么式（21）也应收敛到一个定值tol。因此得到LBD 迭代停止的标准为
式中，tol 为一个定值，此处取为1×10-12。
但在处理实际问题时，LBD 有可能无法完全稳定迭代而出现半收敛现象，表现为
满足式（24）要求后又增大。因此，提出了第二个停止条件。
当迭代次数k0 满足式（25）时，即停止迭代。
2.4 算例分析
2.4.1 算法复原性验证
综合2.1、2.2 和2.3 节所述，给出了基于混合Lanczos-Tikhonov 算法的反演计算流程，如图4所示。
基于图4，为了分析混合Lanczos-Tikhonov 算法的计算精度，借助仿真算例模拟实际测量中的完整计算过程。
图4 基于混合Lanczos-Tikhonov 算法的反演计算流程
Fig.4 The flow chart of inversion calculation based on Lanczos-Tikhonov hybrid method
以图1 所示的绝缘子试样为例，采用COMSOL Multiphysics 在试样表面设置一圈圆环状表面电荷带，如图5a 所示。所设电荷密度从中间向周围逐渐减小，最大电荷密度为1 C/m2，最小电荷密度为0.6 C/m2，计算此时试样表面的电位分布，结果如图5b 所示。为了模拟实际测量电位时的噪声干扰，图5b 的表面电位分布图谱上叠加了均值为0、标准差分别为表面电位φ 最大值的0.1%、0.5%和2%的高斯噪声。
图5 模拟设置的表面电荷及电位分布结果
Fig.5 Surface charge distribution setting and surface potential distribution result of the simulation
采用混合Lanczos-Tikhonov 算法对此算例进行反演计算，计算基于CPU 为i7-7700K、运行内存为32G 的硬件系统完成，迭代开始至计算结束平均用时246 s。得到不同外加噪声水平下电荷反演结果如图6 所示。为了直观地与图5a 中的预置电荷密度分布情况进行比较，将后续各电荷反演结果的色条范围均设置为0～1 C/m2。电荷密度大小超出色条范围的测量单元颜色用色条极值的颜色表示。
图6 基于混合Lanczos-Tikhonov 算法的电荷反演结果
Fig.6 Charge inversion results based on Lanczos-Tikhonov hybrid method
由图6 可知，算法反演结果与模拟设置的电荷密度分布吻合良好。在本算例中，分析各图像导出数据发现，当没有外加噪声时，反演所得电荷密度在-0.03～1.01 C/m2 范围内变化，与预设电荷分布基本一致。随着叠加噪声水平的升高，反演所得图像开始出现噪声残留，当外加噪声等级为0.5%时，反演所得电荷密度的变化范围为-0.19～1.18 C/m2，仍与预设表面电荷分布较为一致。当外加噪声等级为2%时，反演所得电荷密度在-0.20～1.08 C/m2 范围内变化，有较好的反演效果。本文计算了1 000 次迭代条件下，A-WGCV 方法在无噪声时的收敛特性（由相对误差反映），如图7 所示。可见A-WGCV 方法具有较好的收敛稳定性，达到收敛后相对误差不再发生改变。
图7 A-WGCV 方法的收敛特性
Fig.7 Convergence properties of A-WGCV method
2.4.2 精度分析
为了更客观地评价本算法反演计算的精度，引入信噪比（Signal-Noise Ratio,SNR）和峰值均方误差（Peak Mean Square Error,PMSE）的二次方根进行讨论，SNR 和
的计算式分别为
式中，σi 为表面电荷实际值σtrue 中的元素；
为表面电荷反演值
中的元素；A 为表面电荷密度实际值的最大值。
利用SNR 和
对前述算例中不同噪声水平下的电荷反演结果进行计算并对比，结果见表1。当没有噪声时，混合方法反演结果的信噪比为 43.63 dB，其峰值均方误差的二次方根为0.66%；随着噪声等级的提升，信噪比逐渐减小，峰值均方误差的二次方根增大，反演计算的精度下降；但在2%噪声干扰下，其信噪比仍有23.84 dB。可见，本算法具有优异的反演精度。
表1 不同噪声水平下反演电荷的SNR 和
Tab.1 SNR and 
of charge inversion under
different noises
2.5 与其他电荷反演方法的对比
目前，对于平移时变系统，视在电荷法和基于维纳滤波的表面电荷反演算法使用较为广泛。
图8 给出了前述算例在不同噪声下基于视在电荷法的表面电荷反演结果。分析图8 导出数据，在噪声水平为0%时，视在电荷法能够得到较好的反演结果；但当信号中叠加0.1%噪声时，电荷密度在 -0.42～1.16 C/m2 的范围内变化，反演结果开始出现失真现象；当信号中叠加0.5%噪声时，电荷密度在-2.03～2.49 C/m2 的范围内变化，图像失真现象更为明显；当信号中叠加2%噪声时，电荷密度变化范围为-7.47～8.44 C/m2，反演结果已经严重失真。
图8 基于视在电荷法的电荷反演结果
Fig.8 Inversion results of charge based on the apparent charge method
进一步地，计算了利用视在电荷法所求前述算例在不同噪声下反演结果的信噪比和峰值均方误差的二次方根，结果见表2。可知，视在电荷法在无噪声下的反演效果较好，信噪比高，峰值均方误差小；但稍有噪声（0.1%噪声等级）时，其信噪比就发生了明显下降，峰值均方误差大幅度升高，与反演图谱表现一致。
表2 不同噪声下反演电荷的SNR 和
(基于视在电荷法反演结果)
Tab.2 SNR and 
of charge inversion under different noises (based on the apparent charge method inversion result)
可见，微小的干扰对于不适定问题的反问题有着重要影响，严重破坏了反演计算的稳定性。在实际测量中，测量的电压信号均包含了一定的噪声信息。因此，视在电荷法仅适用于对测量及计算精度要求不高时的表面电荷反演计算，对于本文的表面电荷测量系统已不再适用。
基于维纳滤波的表面电荷反演算法本质上仍是基于标准Tikhonov 正则化构造的，以解的二次函数为正则化项，限制其解为平滑解。利用文献[24]提出的PSF 方法选择正则化参数，在绝缘子表面径向中间位置的微元内设置单位电荷，计算得到不同正则化参数下该微元内的电荷分布及其分布的空间幅频特性如图9 所示。
图9 通过PSF 方法反演所得微元内电荷分布及幅频特性
Fig.9 The charge distribution and amplitude-frequency characteristics obtained by PSF
随着正则化参数的减小，反演计算所得到的表 面电荷密度逐渐增大。定义幅频特性曲线衰减到一半时所对应的频率为截止频率，截止频率的倒数为系统的空间分辨率。由图9b 可见，当正则化参数的二次方为HTH 矩阵最大特征值θmax 的0.000 8%～0.03%时，截止频率对应的空间分辨率均可以满足测量要求，但是最优的正则化参数还需要进行多次人工测试后获取。在计算了不同正则化参数的反演结果后，对于本文的测量系统，选定0.03%θmax 的二次方根作为维纳滤波法的正则化参数。
由此，利用基于维纳滤波的反演算法计算了前述算例在不同噪声下的反演结果，该计算基于相同的硬件系统完成，平均用时724 s，相比于混合Lanczos-Tikhonov 算法，计算效率明显降低。基于维纳滤波的电荷反演结果如图10 所示，将该计算结果与本文算法所得的反演结果对比可知，基于维纳滤波的表面电荷反演算法也能良好地复原电荷分布图谱。
图10 基于维纳滤波的电荷反演结果
Fig.10 Charge inversion results based on Wiener filtering
分析基于维纳滤波的表面电荷反演算法计算数据可知，无噪声时，反演所得电荷密度在0～0.98 C/m2的范围内变化，反演效果较好；当信号中叠加0.1%噪声时，电荷密度在-0.07～1.04 C/m2 的范围内变化；随着叠加噪声的等级升高至0.5%，反演所得电荷密度的变化范围为-0.31～1.28 C/m2，反演图像开始出现明显失真；而当信号中叠加2%噪声时，电荷密度大小在-1.22～1.99 C/m2 的范围内变化，所得表面电荷分布图谱相比于混合Lanczos-Tikhonov 算法产生了更多的寄生纹波，并且失真现象更为严重。
为了直观地展示基于混合Lanczos-Tikhonov 算法和维纳滤波法在不同噪声水平下反演效果的优劣，将两种方法的电荷反演结果（即图6、图10 所示结果）绘制折线图如图11 所示。选择图3a 中每一圈测量单元的第一个单元格内的电荷密度代表该圈电荷密度，将有预置电荷的圈，即第21～30 圈的数据进行放大观察。由图11 可见，在无噪声（0%）和高噪声（2%）情况下，混合Lanczos-Tikhonov 算法的反演结果均与预置电荷分布更加吻合。
图11 两种方法反演结果的数据对比
Fig.11 Comparison of inversion results of two the methods
同样地，计算了利用维纳滤波法所求前述算例在不同噪声下反演结果的信噪比和峰值均方误差，结果见表3。维纳滤波法在低噪声下的反演效果较为稳定，信噪比较高，峰值均方误差较小；而在无噪声和高噪声情况下基于混合Lanczos-Tikhonov 算法具有更高的信噪比和更小的峰值均方误差。
表3 不同噪声下反演电荷的SNR 和
(基于维纳滤波反演算法反演结果)
Tab.3 SNR and 
of charge inversion under different noises (based on the inversion results of Wiener filtering inversion algorithm)
可见，基于维纳滤波的反演算法可以良好地复原电荷分布图谱，在测量信号中含有较低水平噪声信号时比较有优势。但其最优正则化参数需要通过反复测试才能选定，且选择标准人为拟定，参数选择效率较低，选择标准不统一，在实际使用中难以控制计算精度。本文基于迭代正则化，利用 AWGCV 方法自动地选择迭代过程中投影正则化问题的正则化参数和权重系数，相较于维纳滤波法具有更高的运算精度和计算效率。
3 实验验证
为了验证本文算法的有效性，采用针电极对绝缘子表面注入电荷，针电极布置方式如图12 所示。
图12 电晕充电针电极布置方式
Fig.12 Arrangement of needle electrode for corona discharge
3.1 结合粉尘图法的算法验证
本文选用粉尘图法[32]表征的绝缘子表面真实电荷分布作为反演算法的对比验证。对针电极施加-10 kV 直流电压10 min 后，测量锥形绝缘子表面电位，然后将带电碳粉与该绝缘子置于同一密闭空间，利用鼓风机吹拂带电有色碳粉于空气中后，碳粉微粒在重力作用下自由下落并被绝缘子表面的异极性电荷吸附，由此表征出绝缘子表面电荷分布的基本轮廓。所选黑色碳粉来自激光打印机中的硒鼓，平均粒径为5 μm，具有单一极性（正极性）。所得绝缘子粉尘图及对应表面电位分布如图13 所示。
图13 粉尘图法及对应表面电位分布
Fig.13 Lichtenberg figure and the surface potential distribution
由图13a 可见，负极性电晕放电使得带正极性电荷的碳粉在绝缘子表面呈“圆斑状”聚集，而在绝缘子中心导体附近无碳粉聚集。在图13b 中可见，绝缘子表面负电位呈椭圆斑状分布，椭圆斑处电位由中心向外逐渐减小。被测电位极性主要为负极性，与针电极极性相同，中心导体附近出现了正电位分布，与粉尘图法所得结果基本一致。
基于该表面电位的测量结果，图14 给出了采用混合Lanczos-Tikhonov 算法、维纳滤波法、视在电荷法三种不同反演算法所得的表面电荷分布。其中，混合Lanczos-Tikhonov 算法在迭代正则化参数时根据A-WGCV 方法自动选择，维纳滤波法的正则化参数仍然选择为HTH 矩阵最大特征值θmax 的0.03%的二次方根。
图14 不同算法表面电荷反演结果
Fig.14 Surface charge inversion results of different algorithms
由图14 可见，尽管维纳滤波法也还原了表面电荷分布，但是其图像中存在着较多寄生纹波，图像边缘模糊，对于噪声的抑制也有限，无法清晰地分辨电荷分布的轮廓。视在电荷法的反演结果中则存在着大量噪声，湮没了电荷信号，对于本文中的大规模矩阵不适定问题而言并不适用。而混合Lanczos-Tikhonov 算法能够较好地还原表面电荷分布，其图像轮廓与粉尘图像基本一致，相比前两个方法可以更准确地反演出本文所用绝缘子的表面电荷密度分布。
3.2 负极性电晕放电下的绝缘子表面电荷分布
在0.1 MPa 空气中进行不同电压下、通过电晕放电使绝缘子表面电荷积聚的实验，并结合粉尘图法和混合Lanczos-Tikhonov 算法，从仿真计算与实验的角度同时观察不同条件下绝缘子表面电荷分布情况，进一步验证本算法在本文实验系统下的可靠性。
对针电极施加-10 kV、-15 kV 和-20 kV 直流电压15 min，测量结果如图15 所示。当电压较低时，如-10 kV，针电极正对的区域主要积聚圆斑状的负极性电荷，其外圈被少量环状正电荷环绕；随着电压升高至-15 kV，“圆斑状”负电荷积聚区域变大，电荷密度升高，同时在负电荷与中心导体之间出现了正极性电荷分布；当电压达到-20 kV 时，针电极正对区域积聚正极性电荷，且被“月牙状”的负极性电荷包围。
图15 不同施加电压下绝缘子表面电荷测量结果
Fig.15 Measurement results for surface charge of the insulator withstand different applied voltage levels
4 结论
针对盆式或锥形绝缘子这类平移变化系统，本文提出了一种基于 A-WGCV 的混合 Lanczos-Tikhonov 算法的绝缘子表面电荷反演算法，得到主要结论如下：
1）本文引入混合Lanczos-Tikhonov 算法解决了电荷-电位转换矩阵的不适定问题，采用 AWGCV 方法选择正则化参数，自适应地求取权重系数ω，有效地减小了矩阵计算量，提高了表面电荷反演计算的效率。
2）本文算法对比视在电荷法在有较高噪声水平时明显具有更高的精度。而采用PSF 方法选择正则化参数时，基于维纳滤波的反演算法在较高噪声水平下的结果出现了更严重的失真。视在电荷法和维纳滤波法在反演效果、SNR 和
方面的表现均在一定程度上弱于本文所用算法。
3）结合粉尘图法的实验结果对比发现，视在电荷法的反演结果中存在大量噪声，维纳滤波法可以还原电荷分布的大致形状，但图像中寄生纹波更重，边缘模糊。两种方法的反演效果与粉尘图像的吻合度均低于本算法。
4）利用针电极放电，结合粉尘图法和算法反演结果发现，随着施加电压的升高，针电极正对区域积聚的电荷发生极性反转，原本的负极性“圆斑”变为“月牙状”包围正极性电荷斑。粉尘图像与反演结果吻合良好，进一步证明了本算法的有效性。
[1]唐炬,潘成,王邸博,等.高压直流绝缘材料表面电荷积聚研究进展[J].电工技术学报,2017,32(8): 10-21.Tang Ju,Pan Cheng,Wang Dibo,et al.Development of studies about surface charge accumulation on insulating material under HVDC[J].Transactions of China Electrotechnical Society,2017,32(8): 10-21.
[2]张贵新,李大雨,王天宇.交流电压下气固界面电荷积聚与放电特性研究进展[J].电工技术学报,2022,37(15): 3876-3887.Zhang Guixin,Li Dayu,Wang Tianyu.Progress in researching charge accumulation and discharge characteristics at gas-solid interface under AC voltage[J].Transactions of China Electrotechnical Society,2022,37(15): 3876-3887.
[3]田汇冬,李乃一,吴泽华,等.直流GIL 柱式绝缘子表面电荷积聚特性[J].电工技术学报,2018,33(22): 5178-5188.Tian Huidong,Li Naiyi,Wu Zehua,et al.Accumulation characteristics of surface charge for post insulators of DC enclosed transmission line[J].Transactions of China Electrotechnical Society,2018,33(22): 5178-5188.
[4]齐波,张贵新,李成榕,等.气体绝缘金属封闭输电线路的研究现状及应用前景[J].高电压技术,2015,41(5): 1466-1473.Qi Bo,Zhang Guixin,Li Chengrong,et al.Research status and prospect of gas-insulated metal enclosed transmission line[J].High Voltage Engineering,2015,41(5): 1466-1473.
[5]王涵,薛建议,陈俊鸿,等.SF6/N2 混合气体中金属微粒对GIL 盆式绝缘子表面电荷积聚特性的影响[J].电工技术学报,2018,33(20): 4663-4671.Wang Han,Xue Jianyi,Chen Junhong,et al.Influence of metal particles on surface charge accumulation characteristics of GIL spacer in SF6/N2 mixture[J].Transactions of China Electrotechnical Society,2018,33(20): 4663-4671.
[6]Liang Zuodong,Lin Chuanjie,Liang Fangwei,et al.Designing HVDC GIS/GIL spacer to suppress charge accumulation[J].High Voltage,2022,7(4): 645-651.
[7]Zhang Lei,Lin Chuanjie,Li Chuanyang,et al.Gassolid interface charge characterisation techniques for HVDC GIS/GIL insulators[J].High Voltage,2020,5(2): 95-109.
[8]穆海宝,张冠军,郑楠,等.Pockels 效应表面电荷测量中电荷反演算法的研究[J].中国电机工程学报,2011,31(13): 135-141.Mu Haibao,Zhang Guanjun,Zheng Nan,et al.Charge inversion algorithm in surface charge measurement based on pockels effect[J].Proceedings of the CSEE,2011,31(13): 135-141.
[9]薛建议,王涵,李科峰,等.直流GIL 中盆式绝缘子表面电荷分布特性研究[J].中国电机工程学报,2018,38(20): 6164-6172.Xue Jianyi,Wang Han,Li Kefeng,et al.Research on charge distribution characteristics on spacer surface in DC GIL[J].Proceedings of the CSEE,2018,38(20): 6164-6172.
[10]Pedersen A.On the electrostatics of probe measurements of surface charge densities[C]// Gaseous Dielectrics V: Proceedings of the Fifth International Symposium on Gaseous Dielectrics,Knoxville,Tennessee,USA,1987: 235-241.
[11]Davies D K.The examination of the electrical properties of insulators by surface charge measurement[J].Journal of Scientific Instruments,1967,44(7): 521-524.
[12]汪沨,胡红利,张乔根,等.用于绝缘子表面电荷测量的电容探头法的标度问题[J].高压电器,2002,38(2): 40-43.Wang Feng,Hu Hongli,Zhang Qiaogen,et al.Calibration of capacitive probe measurement for insulator surface charge[J].High Voltage Apparatus,2002,38(2): 40-43.
[13]刘文静,汪沨,印峰,等.固体绝缘材料空间及表面电荷测量方法[J].绝缘材料,2006,39(6): 59-61,64.Liu Wenjing,Wang Feng,Yin Feng,et al.Measurement of space charge and surface charge accumulated in solid dielectrics[J].Insulating Materials,2006,39(6): 59-61,64.
[14]Zhang Boya,Zhang Guixin.Interpretation of the surface charge decay kinetics on insulators with different neutralization mechanisms[J].Journal of Applied Physics,2017,121(10): 105105.
[15]Rerup T O,Crichton G C,McAllister I W.Using the λ function to evaluate probe measurements of charged dielectric surfaces[J].IEEE Transactions on Dielectrics and Electrical Insulation,1996,3(6): 770-777.
[16]Faircloth D C,Allen N L.High resolution measurements of surface charge densities on insulator surfaces[J].IEEE Transactions on Dielectrics and Electrical Insulation,2003,10(2): 285-290.
[17]Ootera H,Nakanishi K.Analytical method for evaluating surface charge distribution on a dielectric from capacitive probe measurement-application to a cone-type spacer in ±500 kV DC-GIS[J].IEEE Transactions on Power Delivery,1988,3(1): 165-172.
[18]Zhang Boya,Gao Wenqiang,Qi Zhe,et al.Inversion algorithm to calculate charge density on solid dielectric surface based on surface potential measurement[J].IEEE Transactions on Instrumentation and Measurement,2017,66(12): 3316-3326.
[19]Kumada A,Okabe S,Hidaka K.Resolution and signal processing technique of surface charge density measurement with electrostatic probe[J].IEEE Transactions on Dielectrics and Electrical Insulation,2004,11(1): 122-129.
[20]邓军波,王涵,薛建议,等.基于静电探头法的表面电荷密度分布及电场分布的改进反演算法研究[J].中国电机工程学报,2018,38(4): 1239-1247,1301.Deng Junbo,Wang Han,Xue Jianyi,et al.Improved reverse algorithm of surface charge density distribution and electric field distribution based on electrostatic probe[J].Proceedings of the CSEE,2018,38(4): 1239-1247,1301.
[21]潘子君,潘成,唐炬,等.基于图像复原技术与约束最小二乘方滤波器的绝缘子表面电荷反演算法[J].电工技术学报,2021,36(17): 3627-3638.Pan Zijun,Pan Cheng,Tang Ju,et al.Inversion algorithm for surface charge on insulator based on image restoration technology and constrained least square filter[J].Transactions of China Electrotechnical Society,2021,36(17): 3627-3638.
[22]Bracewell R N.Fourier analysis and imaging[M].New York: Springer-Science+Business Media,2003.
[23]张博雅,张贵新.直流GIL 中固-气界面电荷特性研究综述Ⅰ: 测量技术及积聚机理[J].电工技术学报,2018,33(20): 4649-4662.Zhang Boya,Zhang Guixin.Review of charge accumulation characteristics at gas-solid interface in DC GIL,part Ⅰ: measurement and mechanisms[J].Transactions of China Electrotechnical Society,2018,33(20): 4649-4662.
[24]张博雅,张贵新.直流GIL 中固-气界面电荷特性研究综述Ⅱ: 电荷调控及抑制策略[J].电工技术学报,2018,33(22): 5145-5158.Zhang Boya,Zhang Guixin.Review of charge accumulation characteristics at gas-solid interface in DC GIL,part Ⅱ: charge control and suppression strategy[J].Transactions of China Electrotechnical Society,2018,33(22): 5145-5158.
[25]薛建议,王涵,李科峰,等.圆台形绝缘子的表面电荷密度反演算法[J].高电压技术,2019,45(11): 3554-3561.Xue Jianyi,Wang Han,Li Kefeng,et al.Inverse algorithm for surface charge of cone-type spacer[J].High Voltage Engineering,2019,45(11): 3554-3561.
[26]闫纪源,梁贵书,段祺君,等.等离子体表面阶跃型梯度硅沉积对环氧树脂闪络性能的影响[J].电工技术学报,2022,37(9): 2377-2387.Yan Jiyuan,Liang Guishu,Duan Qijun,et al.Effect of plasma surface step gradient silicon deposition on flashover properties of epoxy resin[J].Transactions of China Electrotechnical Society,2022,37(9): 2377-2387.
[27]胡琦,李庆民,刘智鹏,等.基于表层梯度电导调控的直流三支柱绝缘子界面电场优化方法[J].电工技术学报,2022,37(7): 1856-1865.Hu Qi,Li Qingmin,Liu Zhipeng,et al.Interfacial electric field optimization of DC tri-post insulator based on gradient surface conductance regulation[J].Transactions of China Electrotechnical Society,2022,37(7): 1856-1865.
[28]李传扬,林川杰,陈庚,等.高压直流盆式绝缘子气-固界面电荷行为研究综述[J].中国电机工程学报,2020,40(6): 2016-2026.Li Chuanyang,Lin Chuanjie,Chen Geng,et al.Review of gas-solid interface charging phenomena of HVDC spacers[J].Proceedings of the CSEE,2020,40(6): 2016-2026.
[29]Hansen P C.Regularization tools: a Matlab package for analysis and solution of discrete ill-posed problems[J].Numerical Algorithms,1994,6(1): 1-35.
[30]易丽娅.图像复原的Bregman 迭代正则化方法研究[D].武汉: 华中科技大学,2011.
[31]杨艳飞.一般形式正则化的大规模离散线性不适定问题算法的研究[D].北京: 清华大学,2018.
[32]薛乃凡,李庆民,刘智鹏,等.微纳粉尘运动行为与微弱放电探测技术研究进展[J].电工技术学报,2022,37(13): 3380-3392.Xue Naifan,Li Qingmin,Liu Zhipeng,et al.Research advances of the detection technology for kinetic behavior and weak discharge of the micro-nano dust[J].Transactions of China Electrotechnical Society,2022,37(13): 3380-3392.
of the data results obtained by the first two algorithms are analyzed.Analysis shows that the accuracy and computational efficiency of the proposed algorithm in this paper are outstanding.Finally,in order to further verify the correctness of the inversion algorithm,the Lichtenberg figures is used to characterize the real charge distribution on the insulator surface.After negative polarity voltage was applied to the insulator surface by needle electrode,the surface potential of the insulator was measured and inverted by hybrid Lanczos-Tikhonov algorithm,Wiener filtering method and apparent charge method to obtain the charge distribution,which was compared with the Lichtenberg figures.It is proved that the hybrid Lanczos-Tikhonov algorithm can well restore the surface charge distribution of the insulator,and the image contour obtained by the hybrid Lanczos-Tikhonov algorithm is basically consistent with the Lichtenberg image.On this basis,the needle electrode is further used to apply different negative polarity voltage to the insulator surface,and the inversion results of hybrid Lanczos-Tikhonov algorithm are compared with the Lichtenberg image to verify the effectiveness of the algorithm.The experimental results also show that,with the increase of the applied voltage,the polarity of the accumulated charge in the insulator area aligned with the needle electrode will be reversed,and original negative polarity "round spot" charge will become "crescent shape" surrounding the positive polarity charge spot.