CN106199064A - 用于管道内轴向流场成像的迭代超声层析成像方法 - Google Patents

Abstract

一种用于管道内轴向流场成像的迭代超声层析成像方法，通过在管道上下游两个截面安装具有大发射角的超声波传感器构建密集声道网络；激发并接收超声波信号，通过计算超声波沿各声道顺逆流传播的时间获得管道轴向流场沿各声道的平均流速，进而获得若干原始投影数据；对投影数据按平行声道分组，内插并细分，扩充投影数据数量；离散化待建的管道轴向流速场图像，在先验约束条件下，基于投影数据进行迭代层析成像，实现管道轴向流场重建；本发明可以在不干扰原始流态的情况下实现对管道内部轴向流场分布的探测，进而实现对管道轴向流场的高精度重建。

Description

用于管道内轴向流场成像的迭代超声层析成像方法

技术领域

本发明属于信息科学技术领域，应用于管道输送系统中对内部流体流速进行测量，特别涉及一种用于管道内轴向流场成像的迭代超声层析成像方法。

背景技术

在日常生活和工业生产中，为实现物质的输送和能量的转换，常构建一些明渠和管道等对自然界的流动现象主动加以应用。对管道内轴向流场进行定量表征具有多方面的重要意义，如实现高精度的流量计量，研究流体力学理论，解决与流体力学相关的工程问题等。

早期的流场测量方法通过对空间中的不同位置点进行流速测量发展而来，也称为基点法，典型的有毕托管测速仪和热丝风速仪。两者同属于插入式流场测量方法，在测量时均会对流场产生干扰。利用多普勒效应的超声多普勒测速仪和激光多普勒测速仪是两种典型的非插入式基点法流场测量仪器。多普勒测速仪通过利用流体中自然存在的杂质或人为添加粉尘、氧化铁等散射微粒，利用超声或激光等测量微粒的多普勒频移来确定微粒速度，间接得到流场中某些点的流速【刘友,杨晓梅,马修真.基于激光多普勒测速的流场测试技术.激光与光外,2012,42(1):18-21.】。与基点法不同，全场式流场测量方法能够在一次测量中直接获得二维或三维流速分布剖面，包括传统的阴影法、纹影法、干涉法、莫尔偏折法等和目前应用较广泛的以粒子图像测速仪为代表的基于示踪粒子的移动进行流场测量的方法【王浩,曾理江.二维及三维流场的光学测量方法.光学技术,2001,27(2):139-142.】【Dabiri D.Digital particle image thermometry/velocimetry:a review.Exp Fluids,2008,46:191-241.】。全场式流场测量方法还包括电容层析成像【Fuchs A,BrandstatterB,Watzenig D,et al.Flow profile estimator for closed pipes based onelectrical capacitance tomography techniques.Proc.21st IEEEInstrum.Meas.Technol.Conf.,IEEE Press,2004,2326–2631.】、X射线层析成像【Choi YJ,Mccarthy K L,Mccarthy M J.Tomographic techniques for measuring fluid flowproperties.J Food Sci 2002,67(7):2718–2724.】、核磁共振层析成像【Heindel TJ,GrayJN,Jensen TC.An X-ray system for visualizing fluid flows.Flow MeasInstrum,2008,19(2):267–278.】和超声波层析成像【廖光洪,朱小华,林巨等.海洋声层析应用与观测研究综述.地球物理学进展2008,23:1782-1790.】等。超声波层析成像方法利用流体流动对超声波的调制作用来测量沿声道方向的流体平均流速并基于层析成像原理由声道平均流速重建流场分布。基于超声波层析成像的流场测量方法可以弥补光学测量技术在非透光条件下的使用限制，无需像粒子图像测速仪等方法一样添加示踪粒子，并且具有成本低、通用性强和测速精度较高等优点，具有重要的研究意义和应用价值。

目前，将计算机层析成像与时差式超声波法相结合来进行流场测量的方法在海洋洋流观测、大气风场监测和管道内流场测量等很多场合得到了研究和应用。将超声层析成像方法用于管道轴向流场测量时，超声波的传播路径与在其他应用场合下相比更短，这对飞行时间测量精度提出了更高的要求。正因为如此，基于时差式超声波法的管道轴向流场层析成像技术起步较晚。然而，伴随着近年来电子电路技术和计算机技术的发展，该方面的研究受到国内外学者越来越多的关注。1991年，Teerawatanachai等人设计了一种可旋转的超声波探测管段【Teerawatanachai S,Komiya K,Sasamoto H.Estimating the velocityprofile of air flow by means of the ultrasonic translation time computedtomography method.Proc.1991Int.Conf.Ind.Electron.Control Instrum.,IEEE Press,1991:2379–2384.】。管段上安装有1个发射器和7个接收器，共构成7条扇束形声道。由于探测管段可旋转到的角度位置有8个，所以总共可获得56个声道上的投影数据。这项研究首次实现了对管道轴向流场的超声层析重建。Sejong等人将可旋转探测管段工艺与超声波流量计中常用的平行式声道布置方式相结合，实现对流场的平行束式扫描【Chun S,Yoon B-R,Lee K-B.Diagnostic flow metering using ultrasound tomography.J MechSciTechnol2011,25(6):1475-1482.】。与扇形束扫描相比，这种布置方案下超声波收发换能器一一对应，能够保证超声波信号具有较高的信噪比，但相应的扫描效率有所降低。Rychagov等人设计了一种与之类似的平行束式流场扫描方式，区别在于是将可旋转管段与超声波流量计中的折射式声道布置方式相结合【Rychagov M N,Ruchkin S V,Tereshchenko S A,etal.Imaging of fluid flow by tomographic reconstruction using enhancedmultipath ultrasonic measurements.Proc IEEE Ultrason.Symp,IEEE Press,2003:803–806.】。在这种方式下，位于同一侧的发射接收探头构成V字形声道，能够有效地补偿横流对测量的影响。Kurniadi等人设计了一种无需管道进行旋转的超声波传感器布置方案【Kurniadi D,Trisnobudi A.A multi-path ultrasonic transit time flow meterusing a tomography method for gas flow velocity profile measurement.PartPartSystCharact,2006,23(3):330–338.】。该方案分别将探测管段上下游的两个横截面作为发射平面和接收平面，其上各布置有8个发射或接收换能器，共构成64条声道。由于超声波在各声道方向上单向传播，因而在计算流速时需要以零流状态下的超声波传播速度作为标定数据，所以测量结果受流体温度的影响较大。另外，该布置方案的测量精度受管道横流的影响较大并且管道轴心处的成像分辨率较低。

发明内容

为了克服上述现有技术的缺点，本发明的目的在于提供一种用于管道内轴向流场成像的迭代超声层析成像方法，在不干扰原始流态的前提下，无盲区地监控管道内轴向流场分布，实现高精度的二维流场重建。

为了实现上述目的，本发明采用的技术方案是：

一种用于管道内轴向流场成像的迭代超声层析成像方法，通过在管道上下游两个截面安装具有大发射角的超声波传感器构建密集声道网络；激发并接收超声波信号，通过计算超声波沿各声道顺逆流传播的时间获得管道轴向流场沿各声道的平均流速，进而获得若干原始投影数据；对投影数据按平行声道分组，内插并细分，扩充投影数据数量；离散化待建的管道轴向流速场图像，在先验约束条件下，基于投影数据进行迭代层析成像，实现管道轴向流场重建；本发明可以在不干扰原始流态的情况下实现对管道内部轴向流场分布的探测，进而实现对管道轴向流场的高精度重建。

具体包括以下步骤：

步骤一：在管道的上下游两个横截面上分别安装N个具有大发射角和接收角的超声波传感器，每一个传感器均能够与另一截面上所有传感器之间互相收发，共形成N(N-1)条可利用声传播路径，即声道(注：在管道横截面上投影相同的两传感器之间构成的声道对轴向流场成像而言无意义)。

步骤二：分别测量超声波沿各声道顺流和逆流传播的时间，由下式计算管道流体沿各声道的轴向平均流速：

${\stackrel{\‾}{V}}_{path}=\frac{L}{2\cos \mathrm{\φ}}\left(\frac{t_{B\→A}-t_{A\→B}}{t_{A\→B}{t}_{B\→A}}\right)$

其中，t_B→A和t_A→B分别为超声波顺流和逆流传播的时间，φ为声道与管道轴线方向的夹角，L为声道长度；

步骤三：将上一步获得的投影相同声道的轴向平均流速数据取平均，然后分别与各声道在管道横截面上的投影长度相乘，得到原始流场沿各声道的个投影数据，将投影数据按照是否位于同一平行声道组进行分类，共N组；

投影相同声道指的是空间中的若干条不同的但在管道横截面上具有重合的投影的声道。

步骤四：对投影数据进行预处理操作，设定管壁处的投影值为0，为更好地反映投影数据的连续性，选用三次样条曲线进行原始数据拟合，内插函数为：

$\begin{array}{c}p=B\left(t\right)=g_i\left(t\right)=\frac{M_{i-1}}{6u_{i-1}}{(t_i-t)}^3+\frac{M_i}{6u_{i-1}}{(t-t_{i-1})}^3+(\frac{p_{i-1}}{u_{i-1}}-\frac{M_{i-1}}{6}{u}_{i-1})(t_i-t)\\ +(\frac{p_i}{u_{i-1}}-\frac{M_i}{6}{u}_{i-1})(t-t_i),t_{i-1}<t<t_i,(i=1,2,\mathrm{...},n)\end{array}$

式中，(t_i,p_i)为区间[a,b]上的原始数据点，a＝t₀<t₁<…<t_n＝b，u_i-1＝t_i-t_i-1，M_i-1＝B″(t_i-1)，M_i＝B″(t_i)。t₀是内插采用的第一个原始数据点横坐标，t_n是内插采用的最后一个即第n个原始数据点横坐标，内插共采用了n个原始数据点，B″(t_i)是内插函数在第i个投影数据点处的二阶导数。

步骤五：对投影数据的内插曲线进行等间隔离散获得N组细分的投影数据；

步骤六：将包含了重建目标即管道轴向流场的某个正方形图像离散化为J＝n×n个像素，图像向量X＝(x₁,x₂,...,x_J)^T和投影向量P＝(p₁,p₂,...,p_I)^T之间存在下面的关系：

AX＝P

式中，A＝(a_i,j)_I,J是长度矩阵或称为系统矩阵，a_i,j为直线L_i被第j个像素区域截取的长度，代表第j个像素对第i条射线投影的贡献，j＝1,2,3,…,J,i＝1,2,...,I,J为所有像素的个数，I为所有投影数据的个数，图像向量X中x_j为待建流场在第j个像素点区域内流速的平均值，投影向量P中p_i为直线L_i上的投影数据，图像向量X在管壁外侧处的值始终限制为0；

步骤七：执行迭代，将线性方程组中的每个方程看作一个超平面，对于给定的初始值X⁽⁰⁾，将其投影到第一个方程所代表的超平面得到X⁽¹⁾，然后将X⁽¹⁾投影到第二个方程所代表的超平面得到X⁽²⁾，依次类推得到第I个超平面上的投影X^(I)，这样完成了一次迭代。然后，以X^(I)为初始值进行第二次迭代，最终得到投影X₂ ^(I)。迭代过程可表示为：

$\begin{array}{c}{X_k}^{\left(i\right)}={X_k}^{(i-1)}+{\mathrm{\&lambda;}}^{\left(i\right)}\frac{p_i-{X_k}^{(i-1)}\&CenterDot;A_i}{A_i\&CenterDot;A_i}{A}_i\\ ={X_k}^{(i-1)}+{\mathrm{\&lambda;}}^{\left(i\right)}Q(P-{{\mathrm{AX}}_i}^{(i-1)}),i=1,2,\mathrm{...},I\end{array}$

式中，X_k ⁽ⁱ⁾代表经过第k轮迭代的第i次投影后得到的图像向量，A_i为系统矩阵A的第i个行向量，λ⁽ⁱ⁾∈(0,1)为松弛因子，Q是定义的算子；

步骤八：按上式进行若干轮迭代，直到迭代终止的条件满足，即某轮迭代后的图像向量与前一轮迭代后的图像向量X_s-1 ⁽ⁱ⁾之间满足以下关系：

X_s ⁽ⁱ⁾-X_s-1 ⁽ⁱ⁾<ε

式中，ε为给定的某一个极小量，由此，可获得管道内轴向流场成像图。

所述步骤一中，超声波传感器布置在管道的上下游两个横截面上，每个横截面上安装N个，且沿轴向一一对应。

所述步骤一中，超声波传感器的发射角和接收角至少满足：

上游任一超声波传感器发射信号能够被下游任一超声波传感器接收；

下游任一超声波传感器发射信号能够被上游任一超声波传感器接收。

所述步骤四中，存在n-1个约束条件，即：

g′_i(t_i)＝g′_i+1(t_i)(i＝1,2,...,n-1)

且人为限定边界点处的二阶导数为零，即：

M₀＝B″(t₀)＝0，M_n＝B″(t_n)＝0。

所述步骤七中，为了在保证图像重建质量的同时加快迭代速度，最初迭代时步长λ⁽ⁱ⁾应满足下式：

${\mathrm{\&lambda;}}^{\left(i\right)}<\frac{2<{{\mathrm{QAE}}_k}^{\left(i\right)},{E_k}^{\left(i\right)}>}{<{{\mathrm{QAE}}_k}^{\left(i\right)},{{\mathrm{QAE}}_k}^{\left(i\right)}>}=\frac{2<{{\mathrm{AE}}_k}^{\left(i\right)},P-{{\mathrm{AX}}_k}^{\left(i\right)}>}{<{{\mathrm{AE}}_k}^{\left(i\right)},{{\mathrm{AE}}_k}^{\left(i\right)}>}.$

在此基础上，最初迭代时步长λ⁽ⁱ⁾再尽量满足下式：

${\widehat{\mathrm{\&lambda;}}}^{\left(i\right)}=\frac{<Q(P-{{\mathrm{AX}}_k}^{\left(i\right)}),P-{{\mathrm{AX}}_k}^{\left(i\right)},>}{<Q(P-{{\mathrm{AX}}_k}^{\left(i\right)}),Q(P-{{\mathrm{AX}}_k}^{\left(i\right)})>}=\frac{<E_k^{\left(i\right)},P-{{\mathrm{AX}}_k}^{\left(i\right)},>}{<E_k^{\left(i\right)},E_k^{\left(i\right)}>}$

且后期迭代λ⁽ⁱ⁾需要选择较小的常量以保证较好的重建质量。

与现有技术相比，本发明旨在利用超声波法非接触无损测量的优点，通过合适的传感器设计和声道布置方式，形成对管道内轴向流场的网格化密集探测。通过选用大收发角的超声波传感器，使得其区别于传统方法中一对传感器构成一条声道的特点，声道数目大大扩展；通过投影数据的平行重排和内插扩充，克服了实际获取数据在数量上的局限性；采用添加了先验知识的迭代层析成像算法实现了管道轴向流场的成像。本发明为实现在复杂流态下对管道内部轴向流场分布的监控提供了条件。

附图说明

图1是超声波传感器布置方案总体示意图，其中图1(a)是传感器布置方式和收发形式，图1(b)是声道在管道横截面上的投影图。

图2是平行声道组示意图，其中图2(a)是包含N/2个声道的弦声道组，图2(b)是包含个声道的直径声道组。

图3是待建图像离散模型示意图。

图4是消除漩涡影响的设计原理，其中图4(a)是轴向流速矢量分解和横向漩涡矢量分解图，图4(b)是横向漩涡示意图。

具体实施方式

下面结合附图和实施例详细说明本发明的实施方式。

一种用于管道轴向流场成像的迭代超声层析成像方法，包括以下步骤：

步骤一：参照图1，以管径500mm的管道为例，在上下游间距为500mm的两个横截面上加工传感器安装孔。考虑到尺寸限制，出于折中考虑，每个横截面上分别安装12个收发一体式超声波传感器。该布置方法下探头的个数可以增删，具体数目结合实际应用场合下的管道直径、传感器尺寸和成像分辨率的要求进行改动。

步骤二：加工制作具有大声束角的超声波传感器，用于发射和接收超声波信号。传感器的空间指向和声束角均满足一定的条件，使得其可以向另一横截面上的所有传感器收发超声波信号，如图1(a)所示。在本实施例情况下，传感器的空间指向与管道轴线的夹角可设计为22.45°，半指向角设计为大于24.47°。

步骤三：将传感器置于管道上的安装孔中，管道空间中可构成144条声道，其中132条可用声道。由图1(b)所示，此132条声道在管道横截面上的投影线共66条，即每条投影线实际代表空间的两条声道。这样的设计可以消除流场内部漩涡对测量的影响。

如图4所示，探头A、B之间的声道和探头C、D之间的声道是空间中的两条不同声道，但在管道横截面上具有重合的投影。假定位于管道横截面上的流场横向漩涡如图4(b)所示，则由上游探头A向下游探头B发射的超声波速度会因此横流而减慢，导致流速被低估，由上游探头D向下游探头C发射的超声波速度则因此横流而加快，导致流速被高估。漩涡的影响可通过将两声道的流速测量值取平均来抵消。

步骤四：通过一定的时序控制方式完成所有超声波传感器之间信号的收发，利用信号处理方法处理采集的超声波波形，获取超声波沿各声道顺流和逆流传播的时间，由下式计算管道流体沿各声道的轴向平均流速：

${\stackrel{\&OverBar;}{V}}_{path}=\frac{L}{2\cos \mathrm{\&phi;}}\left(\frac{t_{B\&RightArrow;A}-t_{A\&RightArrow;B}}{t_{A\&RightArrow;B}{t}_{B\&RightArrow;A}}\right)---\left(1\right)$

其中，t_B→A和t_A→B分别为超声波顺流和逆流传播的时间，φ为声道与管道轴线方向的夹角，L为声道长度；

步骤五：将上一步获得的投影相同声道的轴向平均流速数据取平均，以消除横向漩涡对速度的影响。将此平均值与声道在管道横截面上的投影长度相乘，可得到原始流场在图1(b)所示投影线上的66个投影数据。将投影数据按照是否位于同一平行声道组进行分类，共12组，有图2所示的两种形式，分别是图2(a)所示的包含N/2个声道的弦声道组以及图2(b)所包含个声道的直径声道组。

步骤六：对投影数据进行预处理操作，设定管壁处的投影值为0。为更好地反映投影数据的连续性，选用三次样条曲线对每组投影数据拟合，内插函数为：

$\begin{array}{c}p=B\left(t\right)=g_i\left(t\right)=\frac{M_{i-1}}{6u_{i-1}}{(t_i-t)}^2+\frac{M_i}{6u_{i-1}}{(t-t_{i-1})}^3+(\frac{p_{i-1}}{u_{i-1}}-\frac{M_{i-1}}{6}{u}_{i-1})(t_i-t)\\ +(\frac{p_i}{u_{i-1}}-\frac{M_i}{6}{u}_{i-1})(t-t_i),t_{i-1}<t<t_i,(i=1,2,\mathrm{...},n)\end{array}---\left(2\right)$

式中，(t_i,p_i)为区间[a,b]上的原始数据点，a＝t₀<t₁<…<t_n＝b。另外，u_i-1＝t_i-t_i-1，M_i-1＝B″(t_i-1)，M_i＝B″(t_i)。

由于原始数据点处的左右导数相同，所以存在n-1个约束条件，即：

g′_i(t_i)＝g′_i+1(t_i)(i＝1,2,...,n-1) (3)

另外，人为限定边界点处的二阶导数为零，即：

M₀＝B″(t₀)＝0 (4)

M_n＝B″(t_n)＝0 (5)

根据原始投影数据的分布和数值特点，可将投影数据改写为投影角度θ_k(即声道倾角)和声道旋转角β_t的函数p(θ_k,β_t)。图2示出了投影角度θ_k、声道旋转角β_t和声道位置r_t。由于β_t为等间隔的，因此相比于以r_t为横坐标，以β_t为横坐标进行样条拟合能够改善两侧声道处的插值效果。

步骤七：对投影数据进行横坐标转换(依据两者之间的关系r_t＝Rcosβ_t)，得到离散化后的等间隔声道(r_t的变化为等距的)投影数据。以每组投影数据扩充为101个为例，总共可获得1212个投影数据。

步骤八：将包含了重建目标即管道轴向流场的某个正方形图像离散化为J＝n×n(例：200×200)个像素，如图3所示。图像向量X＝(x₁,x₂,...,x_J)^T和投影向量P＝(p₁,p₂,...,p_I)^T之间存在下面的关系：

AX＝P (6)

式中，A＝(a_i,j)_I,J是长度矩阵或称为系统矩阵，a_i,j为直线L_i被第j个像素区域截取的长度，代表第j个像素对第i条射线投影的贡献，j＝1,2,3,…,J,i＝1,2,...,I,J为所有像素的个数，I为所有投影数据的个数。另外，图像向量X中x_j为待建流场在第j个像素点区域内流速的平均值，投影向量P中p_i为直线L_i上的投影数据。

步骤九：在迭代之前和迭代过程中，添加先验约束条件：始终设定图像向量X在管壁外侧处的值为0；将真实声道投影数据设置为严格约束条件；添加平滑算法满足流场连续分布的特点等。

步骤十：执行迭代。迭代过程可以表示为：

$\begin{array}{c}{X_k}^{\left(i\right)}={X_k}^{(i-1)}+{\mathrm{\&lambda;}}^{\left(i\right)}\frac{p_i-{X_k}^{(i-1)}\&CenterDot;A_i}{A_i\&CenterDot;A_i}{A}_i\\ ={X_k}^{(i-1)}+{\mathrm{\&lambda;}}^{\left(i\right)}Q(P-{{\mathrm{AX}}_k}^{(i-1)}),i=1,2,\mathrm{...},I\end{array}---\left(7\right)$

式中，X_k ⁽ⁱ⁾代表经过第k轮迭代的第i次投影后得到的图像向量，A_i为系统矩阵的第i个行向量，λ⁽ⁱ⁾∈(0,1)为松弛因子,Q是定义的算子。

为了在保证图像重建质量的同时加快迭代速度，最初迭代时步长λ⁽ⁱ⁾应在满足式(8)条件的基础上尽量满足式(9)。后期迭代λ⁽ⁱ⁾需要选择较小的常量以保证较好的重建质量。

${\mathrm{\&lambda;}}^{\left(i\right)}<\frac{2<{{\mathrm{QAE}}_k}^{\left(i\right)},{E_k}^{\left(i\right)}>}{<{{\mathrm{QAE}}_k}^{\left(i\right)},{{\mathrm{QAE}}_k}^{\left(i\right)}>}=\frac{2<{{\mathrm{AE}}_k}^{\left(i\right)},P-{{\mathrm{AX}}_k}^{\left(i\right)}>}{<{{\mathrm{AE}}_k}^{\left(i\right)},{{\mathrm{AE}}_k}^{\left(i\right)}>}---\left(8\right)$

${\widehat{\mathrm{\&lambda;}}}^{\left(i\right)}=\frac{<Q(P-{{\mathrm{AX}}_k}^{\left(i\right)}),P-{{\mathrm{AX}}_k}^{\left(i\right)},>}{<Q(P-{{\mathrm{AX}}_k}^{\left(i\right)}),Q(P-{{\mathrm{AX}}_k}^{\left(i\right)})>}=\frac{<E_k^{\left(i\right)},P-{{\mathrm{AX}}_k}^{\left(i\right)},>}{<E_k^{\left(i\right)},E_k^{\left(i\right)}>}---\left(9\right)$

按式(7)进行若干轮迭代，直到迭代终止的条件满足，即某轮迭代后的图像向量与前一轮迭代后的图像向量X_s-1 ⁽ⁱ⁾之间满足以下关系：

X_s ⁽ⁱ⁾-X_s-1 ⁽ⁱ⁾<ε (10)

式中，ε为给定的某一个极小量。由此，可获得管道内轴向流场成像图。

Claims (6)

1.一种用于管道内轴向流场成像的迭代超声层析成像方法，其特征在于，包括以下步骤：

步骤一：在管道的上下游两个横截面上分别安装N个具有大发射角和接收角的超声波传感器，每一个传感器均能够与另一截面上所有传感器之间互相收发，共形成N(N-1)条可利用声道；

步骤二：分别测量超声波沿各声道顺流和逆流传播的时间，由下式计算管道流体沿各声道的轴向平均流速：

${\stackrel{\&OverBar;}{V}}_{path}=\frac{L}{2\cos \mathrm{\&phi;}}\left(\frac{t_{B\&RightArrow;A}-t_{A\&RightArrow;B}}{t_{A\&RightArrow;B}{t}_{B\&RightArrow;A}}\right)$

其中，t_B→A和t_A→B分别为超声波顺流和逆流传播的时间，φ为声道与管道轴线方向的夹角，L为声道长度；

步骤三：将上一步获得的投影相同声道的轴向平均流速数据取平均，然后分别与各声道在管道横截面上的投影长度相乘，得到原始流场沿各声道的个投影数据，将投影数据按照是否位于同一平行声道组进行分类，共N组；

步骤四：对投影数据进行预处理操作，设定管壁处的投影值为0，为更好地反映投影数据的连续性，选用三次样条曲线进行原始数据拟合，内插函数为：

$\begin{array}{c}p=B\left(t\right)=g_i\left(t\right)=\frac{M_{i-1}}{6u_{i-1}}{(t_i-t)}^3+\frac{M_i}{6u_{i-1}}{(t-t_{i-1})}^3+(\frac{p_{i-1}}{u_{i-1}}-\frac{M_{i-1}}{6}{u}_{i-1})(t_i-t)\\ +\left(\frac{p_i}{u_{i-1}}-\frac{M_i}{6}{u}_{i-1}\right)\left(t-t_i\right),t_{i-1}<t<t_i,\left(i=1,2,...,n\right)\end{array}$

式中，(t_i,p_i)为区间[a,b]上的原始数据点，a＝t₀<t₁<…<t_n＝b，u_i-1＝t_i-t_i-1，M_i-1＝B”(t_i-1)，M_i＝B”(t_i)，t₀是内插采用的第一个原始数据点横坐标，t_n是内插采用的最后一个即第n个原始数据点横坐标，内插共采用了n个原始数据点，B″(t_i)是内插函数在第i个投影数据点处的二阶导数；

步骤五：对投影数据的内插曲线进行等间隔离散获得N组细分的投影数据；

步骤六：将包含了重建目标即管道轴向流场的某个正方形图像离散化为J＝n×n个像素，图像向量X＝(x₁,x₂,...,x_J)^T和投影向量P＝(p₁,p₂,...,p_I)^T之间存在下面的关系：

AX＝P

式中，A＝(a_i,j)_I,J是长度矩阵或称为系统矩阵，a_i,j为直线L_i被第j个像素区域截取的长度，代表第j个像素对第i条射线投影的贡献，j＝1,2,3,…,J,i＝1,2,...,I,J为所有像素的个数，I为所有投影数据的个数，图像向量X中x_j为待建流场在第j个像素点区域内流速的平均值，投影向量P中p_i为直线L_i上的投影数据，图像向量X在管壁外侧处的值始终限制为0；

步骤七：执行迭代，迭代过程表示为：

$\begin{array}{c}{X_k}^{\left(i\right)}={X_k}^{\left(i-1\right)}+{\mathrm{\&lambda;}}^{\left(i\right)}\frac{p_i-{X_k}^{\left(i-1\right)}\&CenterDot;A_i}{A_i\&CenterDot;A_i}{A}_i\\ ={X_k}^{\left(i-1\right)}+{\mathrm{\&lambda;}}^{\left(i\right)}Q\left(P-{{\mathrm{AX}}_k}^{\left(i-1\right)}\right),i=1,2,...,I\end{array}$

式中，X_k ⁽ⁱ⁾代表经过第k轮迭代的第i次投影后得到的图像向量，A_i为系统矩阵A的第i个行向量，λ⁽ⁱ⁾∈(0,1)为松弛因子，Q是定义的算子；

步骤八：按上式进行若干轮迭代，直到迭代终止的条件满足，即某轮迭代后的图像向量与前一轮迭代后的图像向量X_s-1 ⁽ⁱ⁾之间满足以下关系：

X_s ⁽ⁱ⁾-X_s-1 ⁽ⁱ⁾<ε

式中，ε为给定的某一个极小量，由此，可获得管道内轴向流场成像图。

2.根据权利要求1所述用于管道内轴向流场成像的迭代超声层析成像方法，其特征在于，所述步骤一中，超声波传感器布置在管道的上下游两个横截面上，每个横截面上安装N个，且沿轴向一一对应。

3.根据权利要求1所述用于管道内轴向流场成像的迭代超声层析成像方法，其特征在于，所述步骤一中，超声波传感器的发射角和接收角至少满足：

上游任一超声波传感器发射信号能够被下游任一超声波传感器接收；

下游任一超声波传感器发射信号能够被上游任一超声波传感器接收。

4.根据权利要求1所述用于管道内轴向流场成像的迭代超声层析成像方法，其特征在于，所述步骤四中，存在n-1个约束条件，即：

g′_i(t_i)＝g′_i+1(t_i)(i＝1,2,...,n-1)

且人为限定边界点处的二阶导数为零，即：

M₀＝B”(t₀)＝0，M_n＝B”(t_n)＝0。

5.根据权利要求1所述用于管道内轴向流场成像的迭代超声层析成像方法，其特征在于，所述步骤七中，为了在保证图像重建质量的同时加快迭代速度，最初迭代时步长λ⁽ⁱ⁾应满足下式：

${\mathrm{\&lambda;}}^{\left(i\right)}<\frac{2<Q{{\mathrm{AE}}_k}^{\left(i\right)},{E_k}^{\left(i\right)}>}{<Q{{\mathrm{AE}}_k}^{\left(i\right)},Q{{\mathrm{AE}}_k}^{\left(i\right)}>}=\frac{2<{{\mathrm{AE}}_k}^{\left(i\right)},P-{{\mathrm{AX}}_k}^{\left(i\right)}>}{<{{\mathrm{AE}}_k}^{\left(i\right)},{{\mathrm{AE}}_k}^{\left(i\right)}>}.$

6.根据权利要求5所述用于管道内轴向流场成像的迭代超声层析成像方法，其特征在于，最初迭代时步长λ⁽ⁱ⁾应尽量满足下式：

${\widehat{\mathrm{\&lambda;}}}^{\left(i\right)}=\frac{<Q(P-{{\mathrm{AX}}_k}^{\left(i\right)}),P-{{\mathrm{AX}}_k}^{\left(i\right)},>}{<Q(P-{{\mathrm{AX}}_k}^{\left(i\right)}),Q(P-{{\mathrm{AX}}_k}^{\left(i\right)})>}=\frac{<E_k^{\left(i\right)},P-{{\mathrm{AX}}_k}^{\left(i\right)},>}{<E_k^{\left(i\right)},E_k^{\left(i\right)}>}$

且后期迭代λ⁽ⁱ⁾需要选择较小的常量以保证较好的重建质量。