分治法在管道涡流检测阻抗解析中的应用

王志刚，罗清旺，师奕兵

（电子科技大学自动化工程学院 成都 611731）

·自动化技术·

分治法在管道涡流检测阻抗解析中的应用

王志刚，罗清旺，师奕兵

（电子科技大学自动化工程学院 成都 611731）

介绍了一种铁磁性管道涡流阻抗模型的数值解析方法。该方法利用阻抗模型中贝塞尔函数在大宗量（分治法产生的一个子区间）时的渐近性，在整个广义积分区间采用分治法简化积分函数，降低了对阻抗模型广义积分的计算难度。基于对计算量与计算准确度的折中，讨论了该解析方法中分治点的选择原则。将该方法的阻抗模型计算结果与通过物理检测设备测试实际管道的值进行比较，验证了分治法解阻抗模型的可行性。该方法对于解析柱坐标系下管道涡流检测阻抗模型的应用，具有简单、快速和高精度的优点。

贝塞尔函数； 大宗量； 分治法； 阻抗模型； 数值解析

在管道涡流无损检测研究中，常利用Maxwell方程及某种边界条件，列出电磁波在铁磁性管道中传播的模式方程，然后根据该模式方程，研究电磁波在铁磁性管道中的传播阻抗特性，并从传播阻抗特性中提取出管道物理信息［1-3］，最后基于这些物理信息判别铁磁性管道的损伤情况。上述针对铁磁性管道的检测方法，涉及到对电磁场中传播阻抗的数值解析。由文献［4-5］可知，管道涡流阻抗模型是关于修正贝塞尔函数的一种复杂广义积分形式，该积分形式决定了它难于直接利用贝塞尔函数的积分性质进行计算。虽然也有一些关于贝塞尔广义积分函数的研究性论文发表，但是它们多是针对特定的贝塞尔函数形式［6-9］，如较常见且采用的有复杂递推法［7］以及针对特定形式的贝塞尔函数积分［8］方法。前者只能用于特定的递推形式，后者应用于有限的积分区间，文献［9］提出了修正贝塞尔函数的数值积分方法的数学演绎。对于本文所需解析的涡流阻抗模型中具有的函数复杂性、函数高速振荡衰减及在广义区间积分等特征，它们均不适用。

另一方面，函数的渐近性被广泛应用于解析复杂积分函数，并具有很好的解析效果［10-13］。可以利用贝塞尔函数渐进性，采用分治法处理复杂积分函数。分治法的思想是将一个难以直接解决的大问题，分割成一些规模较小的相同问题，以便分而治之。依此将阻抗模型的广义积分区间分治为小宗量积分区间和大宗量积分区间，并在大宗量积分区间应用贝塞尔函数的渐近性，选定逼近函数简化阻抗函数。一般的数列排序分治法［14］，通常将最后的数值点作为分治点，它的选取不影响计算结果。考虑到需对涡流阻抗模型中积分区间［0，∞］分治，所以不能取最后数点作为分治点，这样对分治点的选取时需要考虑逼近函数对原函数的逼近精度，以及阻抗函数的数值计算准确度。

本文以分治法求解阻抗模型的过程包含：分治积分区间及阻抗模型积分函数的衰减特性和逼近特性分析；积分函数计算和结果准确度分析两个阶段。依靠分治法和贝塞尔函数逼近的特性的简化，可在保证计算结果精度的条件下大大减少计算量。

1 铁磁性管道涡流阻抗模型

图1所示，激励及接收线圈均处于管道内部（管道无限长）。a、b、c分别为激励、接收线圈、管道的内半径，线圈轴向距离为s，s

图1 管道电磁传播阻抗模型

在柱坐标系中，各区域中任一点的磁位量设为A（r，z；r′，z′）其为包含复变量的贝塞尔方程：

式中，I1（xnr）,K1（xnr）分别为1阶第一类和第二类修正贝塞尔函数；Cn（x）、Dn（x）分别代表在区域n的函数系数； xn为介质传播参数，可表示为，

式中，下标n=1，2，3，4分别对应区域I，II，III，IV，区域I和II为空气；x=x1=x2,μ0=μ1=μ2。对式（1）求解，并设a=b＜c，经复杂推导可得铁磁性管道电磁涡流检测阻抗模型为［15］：

其中，式（3）为空气中线圈组传输阻抗，式（4）表示由于铁磁性管道存在所导致阻抗函数变化的增量。因此，阻抗函数可写为：Z=Zd+Λ。由式（3）、式（4）的积分特性可知，被积函数是衰减函数，且计算难度较大。另外，由管道趋肤特性可知，|x3c|＞＞1，根据修正贝塞尔函数的性质，式（6）成立：

将式（6）代入式（5）简化，有：

由于管道具有较大的电导率，可对式（2）化简，得到：

观察式（4）～式（8），由于贝塞尔函数特性及广义的积分区间，化简之后仍然很难直接积分计算或利用贝塞尔函数性质进行计算。

2 分治法解阻抗模型

分治法解阻抗模型，是将广义积分区间分治为小宗量和大宗量区间；在大宗量区间，以简化近似函数代替原函数进行计算；而在小宗量区间，则进行定积分或求和计算。

2.1 贝塞尔函数的Melin-Barnes逼近

由分治法思想可知，在积分区间采用分治法，若要对积分函数化简，其在大宗量区间时必须有高精度逼近的近似函数。由文献［13］可知，汉克函数在大宗量时具有良好的逼近函数，并且，修正贝塞尔函数可由汉克函数线性表示。所以，如果已知汉克函数在大宗量时的高精度逼近函数，就可以通过两类贝塞尔函数的线性关系，求得其在大宗量时修正贝塞尔函数的高精度逼近式，该线性关系为：

Borel求和法和Melin-Barnes数值法对贝塞尔函数均可实现良好逼近［13］，而两者比较，Borel求和法计算量大，本文采用Melin-Barnes数值法对大宗量区间的修正贝塞尔函数进行逼近。

由于柯西留数定理在一些特殊实积分，如反常积分、广义积分是很好的分析手段，而Melin-Barnes数值法的基本思想是：在形如式（11）的一般幂级数中应用柯西留数定理，有：

在一定条件下［13］，在式（11）中应用柯西留数定理，可将 S（N，z）表示成Melin-Barnes积分，有：

式中，N-1＜c＜N 。由于汉克函数的完全等式是形如式（11）所示的无穷级数的幂级数，不利于直接数值计算，所以应用柯西留数定理将汉克函数的幂级数形式转换成形如式（12）的Melin-Barnes积分，得到汉克函数的高精度逼近式，以便于达到快速数值计算汉克函数的目地，具体过程参见文献［13］。

在大宗量（λ＞＞1）时，可由Melin-Barnes数值法［13］逼近的汉克函数简化为：

式中，α表示阶数。式（13）、式（14）提供了大宗量区间精确的汉克函数逼近公式，再结合式（9）、式（10）可以得到修正贝塞尔函数在大宗量区间内高精度的逼近式。然而，这种逼近不可避免存在一定误差，它与分治点的选取有直接关系。

2.2 选取分治点

在分治法解析阻抗模型的过程中，分治点λ的选取决定了渐近函数引入起点。理论上，分治点λ越大，大宗量积分区间越小，计算结果越准确，但是计算量就会越大；反之，分治点λ越小，计算量越小，大宗量区间的逼近误差就会越大。合理选取分治点，使得计算量和结果准确性得以折中。

为了便于分治点的选取，首先需要对阻抗模型被积函数的变化趋势进行分析，它是一个稳定系统的输出（该稳定系统的输出是一个电压信号）。由稳定系统特性可知，其输出信号必然收敛，阻抗模型在［0，∞）内广义积分被积函数必然衰减，这为大宗量区间内被积函数的逼近提供了可能，函数的衰减特性导致大宗量区间内的积分量占积分总量较少。图2所示为阻抗Zd、Λ的被积函数在有限区间内的分布情况。

由图2可知， Zd的被积函数呈震荡衰减的趋势，Λ的被积函数的实部、虚部均呈快速衰减的趋势，其值在λ＞4时基本趋于0。阻抗模型的被积函数变化趋势为分治算法提供了可行性，亦为分治算法中分治点的起点提供一个大致范围，如λZd≥10，λΛ≥4。为获得误差小的渐近贝塞尔函数分治起点，结合文献［16］，修正贝塞尔函数渐近曲线的误差为：

式中，α为修正贝塞尔函数的阶数；λ为分治点。由式（15）可知λ值越大，修正贝塞尔函数的渐近误差越小，Melin-Barnes数值法的逼近度越高，但是会增加数值计算量。

图2 涡流阻抗模型的被积函数分布

分治点的选取通过观察函数计算结果的稳定度并进一步确定阻抗模型数值结果的准确度为原则。结合图2和式（15）就可为分治算法提供误差小，数值大小合适的分治点。

2.3 分治法解涡流阻抗模型

在阻抗模型中，物理参数满足如下关系式，0

实际测量中通常取 a=b，而a＞b和a＜b情况下的阻抗函数衰减速度均大于 a=b情况下函数的衰减速度。因此，若该方法对 a=b的情况适用，那么对a＞b或a＜b也适用。在式（3）中，选取分治点 λ0分治积分区间，在大宗量区间进行渐近性分析，将式（9）、式（10）、式（13）、式（14）代入式（3）的被积函数，当xa=xb≥λ0，a=b＜c时，有：

所以，式（3）可化为：

同理，在式（4）中选取分治点 λ1，在大宗量区间进行渐近分析。当1xaλ≥时，有：

将式（18）代入式（4），得：

由式（1）、式（2）可知，dx表示磁位量 A（r， z；r′， z′）在介质传播点x的微积分， x3可由式（8）化简。式（17）所示的积分式在数值计算软件中可直接采用函数命令调用的方式得到数值结果。而式（19）在积分区间［0，λ1/a］的表达式仍很难直接函数调用得到积分结果。由图2可知，Λ的被积函数衰减很快，故在［0，λ1/a］分治区间内，可用积分定义将式（19）变化为：

式（20）中，在分治区间［0，λ1/a］中累加和的步进xΔ在一定程度上影响结果的精度，本文选取Δx=0.005，以保证计算结果的准确性。利用MATLAB软件函数库及高效的数值计算和符号计算功能，对于可直接计算的贝塞尔函数数值积分，采用函数命令调用的形式实现，如式（17）。对于难于积分的复杂贝塞尔函数形式，采用积分定义的方法积分，如式（19）在［0，λ1/a］的积分。计算中，约定在计算量和计算结果精度满足要求前提下，当3次结果相对误差低于0.5%时，即采用第2次的计算结果，其对应的截断点λ判为合适。算法过程如下：

1） 参数初始化。包括：c ＞ a = b ＞ 0，ω ＞ 0，μ0＞ 0，σ3＞0,ε3＞0,λ=1,λstep=2,x=0，xstep=0.005等。

2） 若x≤λ / a ，λ=λ+λstep

3 结果验证

为了验证分治法解析阻抗模型的正确性，选取合适的管道参数值，进行数值计算。依实际情况取参数：a=b=0.018 9 m，s=0.063 5m，c=0.067 2 m ，f =20kHz ，μ0=4π ×10-7H/m，μ3=80×4π ×10-7H/m，σ3=5.95×106S/m。根据图2和式（15）取：λ0≥10， λ1≥4 。由分治点λ与参量的关系式可知，在大宗量区域：λ≤xa。应用分治法计算Zd和Λ的数值结果如表1、表2所示。

表1 λ0不同时　Zd的值

表2 λ1不同时Λ的值

由表1、表2可知，在分治点较小时（如 λ0=1，λ1=1时）Zd与Λ的分治结果与后面稳定的分治结果有一定的差距，正如前文分析，高精度的贝塞尔逼近函数须在大宗量区间引入，即分治点较大时引入。当分治点较小时引入，会造成分治结果不正确。而在 λ0≥10和 λ1≥4的范围，Zd与Λ的分治结果基本趋于稳定，其值的稳定性体现了应用分治法解阻抗模型的正确性。由表1可知 λ0=15作为Zd的分治点；由表2知 λ1=5作为Λ的分治点最合适，既减小了计算量，又保证了结果值的准确度。为了体现分治点对积分结果的影响（因为在 λ0≥10的范围， Zd的值趋于稳定，呈现水平的特性，所以分治点对分治结果的影响主要体现在0≤λ0≤10范围内），选取 Zd～λ0（λ0≤12）绘图，如图3所示。

由图3可知，分治点值越大，Zd积分结果越稳定。在分治点 λ0≤8的区域内，其对Zd积分结果影响很大， λ0越小，影响效果越大。表1、表2和图3的积分结果的准确性是通过积分结果渐近误差的连续稳定性来保证的，即

图3 Zd积分值与分治点取值关系

为了更好地验证分治法解阻抗模型的正确性，本文通过实际管道阻抗检测值来验证其一致性。

对于式（7）所示的算子，在电磁涡流近场的阻抗模型中反映了管道带来的耦合影响，对其可做如下变形，以获得更加简洁的表述：

式中，

首先，运用分治法计算获得式（3）、式（4）的L值，再通过管道测试试验，获得实际涡流阻抗的L值，将上述结果与去耦度ζ做变化曲线进行比较，曲线如图4所示。图4从电学试验角度展示了分治算法解阻抗模型的可行性和正确性。由图4可知运用分治算法解电磁涡流近场阻抗模型的结果与实际物理检测值具有很好的一致性（个别实测异常点除外，如2c=0.153 7 m，ζ=13.1579时）。另外，由于分治算法的快速、资源占用少的特点，其在解阻抗模型中更具有可操作性。

图4 分治法L值与实测L值比较图

4 结 论

通过对铁磁性管道近场涡流传输阻抗模型的分析，利用分治法将该阻抗模型的广义积分区间分治为小宗量积分区间和大宗量积分区间，并在大宗量积分区间应用贝塞尔函数的渐近性，选定汉克函数逼近简化阻抗函数，完成复杂阻抗函数的解析，同时，通过计算结果的逼近误差稳定性来确定结果的准确性。电学对比实验的结果表明，该方法不仅可以快速有效的解析输出，且与实际测试结果具有很好的一致性，充分证明了该方法的可行性。

［1］ DODD C V， DEEDS W E. Analytical solutions to eddy -current probe-coil problems［J］. J Appl Phys， 1968， 39（6）：2829-2838.

［2］ DODD C V， CHENG C C， DEEDS W E. Induction coils coaxial with an arbitrary number of cylindrical conductors［J］. J Appl Phys， 1974， 45（2）： 638-647.

［3］ DESJARDAINS D P R， KRAUSE T W， GAUTHIER N. Analytical modeling of the transient response of a coil encircling a ferromagnetic conducting rod in pulsed eddy current testing［J］. NDT and E International， 2013， 60：127-131.

［4］ VASIC D， BILAS V. Model-based stochastic inversion of coil impedance for determination of tube inner radius and electromagnetic properties［J］. AUTOMATIKA： Journal for Control， Measurement， Electronics， Computing and Communications， 2009， 50（3）： 227-233.

［5］ VASIC D， BILAS V. Validation of a coil impedance model for simultaneous measurement of electromagnetic properties and inner diameter of a conductive tube［J］. IEEE Trans Instrum Meas， 2006， 55（5）： 337-341.

［6］ CHEN Ru-yun， AN Cong-pei. On the evaluation of infinite integrals involving Bessel functions［J］. Applied Mathematics and Computation， 2014， 235： 212-220.

［7］ LIN Q G. Infinite integrals involving Bessel functions by contour integration［J］. Integral Transforms and Special Functions， 2013， 24（10）： 783-795.

［8］ CHEN Ru-yun， AN Cong-pei. On evaluation of Bessel transforms with oscillatory and algebraic singular integrands［J］. Journal of Computational and Applied Mathematics， 2014， 264： 71-81.

［9］ BARICZ A， POGANY T. On a sum of modified Bessel functions［J］. Mediterranean Journal of Mathematics， 2014，11（2）： 349-360.

［10］ NICULESCU C， POPOVICI F. The asymptotic behavior of integrable functions［J］. Real Analysis Exchange， 2012，38（1）： 157-168.

［11］ XIANG Shu-huang. Laplace transforms for approximation of highly oscillatory Volterra integral equations of the first kind［J］. Applied Mathematics and Computation， 2014， 232：944-954.

［12］ GIULIANO R， MACCI C. Large deviation principles for sequences of maxima and minima［J］. Communications in Statistics-Theory and Methods， 2014， 43（6）： 1077-1098.

［13］ KOWALENKO V. Exactification of the asymptotics for Bessel and Hankel functions［J］. Applied Mathematics and Computation， 2002， 133（2）： 487-518.

［14］ CORMEN T H， LEISERSON C E， RIVEST R L， et al. Introduction to algorithms［M］. 3rd ed. Beijing： China Machine Press， 2013： 16-21.

［15］ 王志刚， 罗清旺， 师奕兵. 铁磁性管道内涡流线圈耦合分析与管道参数检测［J］. 仪器仪表学报， 2014， l35（12）：2843-2851. WANG Zhi-gang， LUO Qing-wang， SHI Yi-bing. Analysis of eddy current coil coupling in ferromagnetic pipe and parameter detection［J］. Chinese Journal of Scientific Instrument， 35（12）： 2843-2851.

［16］ 仇海亮. 复宗量贝塞尔函数的数值解法及在测井解释中的应用［J］. 中国石油和化工标准与质量， 2013（17）： 128. QIU Hai-liang. Numerical solution of complex argument Bessel function and its application in the log interpretation［J］. China Petroleum and Chemical Standard and Quality， 2013（17）： 128.

编 辑 漆 蓉

Divide-and-Conquer PAippp e rEodadchy ACpuprrlieendt tToe Ismtpedance Analysis of

WANG Zhi-gang， LUO Qing-wang， and SHI Yi-bing
（School of Automation engineering， University of Electronic Science and Technology of China Chengdu 611731）

This paper proposes a unique analytical method for eddy current impedance model of ferromagnetic pipe. Based on the asymptotic behavior of Bessel function in model under bulk quantity occasion，the method exploits the divide-and-conquer approach in the generalized integral interval to simplify model calculation. The calculation results of the model are compared with the results of practical pipe eddy current test，verifying the correctness and feasibility of calculation on the divide-and-conquer approach. The analytical method has the advantages of conciseness and high precision for impedance model applications under cylindrical coordinate systems.

Bessel function； bulk quantity； divide-and-conquer approach； impedance model；numerical analysis

TE1

A

10.3969/j.issn.1001-0548.2016.02.011

2014 - 07 - 21；

2015 - 11 - 25

“十二五”国家科技重大专项（2011ZX05020-006-005）

王志刚（1967 - ），男，副教授，主要从事检测技术方面的研究.