一种柱面波成像方法

1.本发明涉及探测与成像技术领域，尤其是一种柱面波成像方法。


背景技术：

2.平面波成像(plane wave imaging,pwi)，通过激励单元一次发射平面波，覆盖整个成像区域，同时接收单元采集整个成像区域的数据，实现高帖率成像，被广泛应用于医学成像和无损检测等领域。2009年，montaldo等人又提出了复合平面波成像，用于提高图像的对比度和分辨率。
3.但是，对于工业管道检测、测井、人体食管和血管的医学成像等应用领域，往往采用环形阵列激励和检测，其激励的波阵面近似为柱面波。目前，尚无针对该类检测方式的类似于平面波成像的高帖率成像算法。


技术实现要素：

4.有鉴于此，本发明提出一种通用的柱面波成像算法，通过一次激励产生柱面波，覆盖整个成像区域，实现高帖率成像，主要应用于工业管道检测、测井、人体食管和血管的医学成像等应用领域。
5.本发明的技术方案为：一种柱面波成像方法，包括如下步骤：
6.步骤1，通过激励单元产生柱面波，所述柱面波是通过环形阵列激励单元、圆弧形阵列激励单元、柱面阵列激励单元之一同步激励产生，或者通过相控阵激励单元利用波束成形技术产生；柱面波成像中的激励波阵面是柱面，或者为位于柱面内的不规则形状，柱面波传播范围覆盖整个成像区域；
7.柱面波成像中的波阵面的传播方向可以向内，也可以向外。波阵面向内传播适合成像区域在初始波阵面内部的情况，例如工业管道检测；波阵面向外传播时适合成像区域在初始波阵面外侧的情况，例如测井、人体食管和血管检测。
8.步骤2，通过多个接收单元采集波形原始信号，所述接收单元与激励单元复用，或者与激励单元分离、根据成像需求布置。
9.步骤3：将成像区域剖分为若干成像目标点。根据成像区域大小、维度、成像分辨率等要求将成像区域剖分为若干网格，网格尺寸小于或等于成像分辨率要求，每个网格为一个成像目标点。
10.步骤4：选择数学模型计算每个成像目标点的成像值，完成整个成像区域的柱面波成像。
11.如果成像区域为柱面波传播范围内的任意体或者面，成像区域内单个成像目标点的成像值的数学模型由公式(1)和(2)表示：
[0012][0013]
其中，成像区域剖分为m个成像目标点，(xj,yj,zj)的表示第j个成像目标点的坐标，j＝1,
…
,m；共采用n个接收单元，(xi,yi,zi)表示第i个接收单元的坐标，i＝1,
…
,n；r表示初始波阵面的半径，c为波速，t
ij
表示柱面波经过第j个成像目标点散射到达第i个接收单元需要的传播时间；s(i,t
ij
)表示第i个接收单元在第j个成像目标点的成像值，通过对第i个接收单元采集的信号在t
ij
时刻插值计算得到；将n个接收单元在第j个成像目标点的成像值叠加，得到第j个成像目标点的成像值vj。
[0014]
计算第i个接收单元在第j个成像目标点的成像值s(i,t
ij
)时，具体插值方法与激励波形特征、采样率等相关，一般采用邻近线性插值。根据实际需要，可以提高插值阶数，选择二阶拉格朗日插值、或者三次样条插值等多种插值方法。
[0015]
当成像区域为一个二维水平断面时，成像目标点的成像值的数学模型由公式(1)、(2)和(3)表示：
[0016]
zj＝z0ꢀꢀꢀ
(3)
[0017]
公式(3)对所有j＝1,
…
,m成立。数学模型(1)、(2)和(3)是数学模型(1)和(2)的一个特例，对成像目标点增加约束(3)，表示所有成像目标点位于高度为z0的二维水平断面内。需要注意，接收单元并不要求位于成像断面内，可以根据实际情况布置。特别地，如果接收单元也位于z0水平断面内，则公式(2)可以进一步简化为：
[0018][0019]
公式(1)、(2)和(3)表示的数学模型适合接收单元布置在成像水平断面外的二维水平断层成像；公式(1)、(3)和(4)表示的数学模型适合接收单元位于成像水平断面内的二维水平断层成像。
[0020]
当成像区域为一个过柱面中心轴的二维竖直断面时，作为公式(1)和(2)表示的数学模型的一个特例，成像目标点的成像值的数学模型由公式(1)、(2)和(5)表示：
[0021][0022]
其中tg-1
表示反正切函数，公式(5)对所有j＝1,
…
,m成立，表示所有成像目标点位于扇角为的二维竖直断面内。同样地，接收单元并不要求位于成像断面内，可以根据实际情况布置。
[0023]
当成像区域为一个圆柱面时，作为公式(1)和(2)表示的数学模型的一个特例，成像目标点的成像值的数学模型由公式(1)、(2)和(6)表示：
[0024]
[0025]
公式(6)对所有j＝1,
…
,m成立，表示所有成像目标点位于半径为r0的柱面内。同样地，接收单元可以根据实际情况布置。
[0026]
成像目标点的成像值由各接收单元在该成像目标点的成像值直接叠加得到，或者对各接收单元在该成像目标点的成像值施加数学变换后再叠加得到，即前述数学模型的公式(1)可以由公式(7)替换：
[0027][0028]
其中，f{s(i,t
ij
)}表示对第i个接收单元在第j个成像目标点的成像值施加数学变换f。
[0029]
所述数学变换之一为归一化变换，用公式(8)-(11)表示：
[0030][0031][0032][0033][0034]
其中，公式(8)和(9)为有极性归一化，即归一化后的成像值仍保留正负极性；公式(10)和(11)为无极性归一化，即归一化后的成像值失去正负极性，只保留幅值特征。根据实际情况选择具体使用哪一种规一划变换。
[0035]
所述数学变换的另一种方式为系数加权法，用公式(12)表示：
[0036][0037]
其中，公式(12)中的w
ij
表示第i个接收单元在第j个成像目标点的成像值的加权系数，根据成像目标点与接收单元的之间相对位置关系、接收单元的接收口径与灵敏度函数、以及柱面波场的物理特征等确定。通常，将接收单元对成像目标点所张的立体角作为权系数；或者根据接收单元的指向角确定，成像目标点位于接收单元指向角范围内，权系数为1；成像目标点位于接收单元指向角范围外，权系数为0。
[0038]
上述两类数学变换可以同时使用，或者使用其中之一，或者均不使用。
[0039]
有益效果
[0040]
本发明提出的柱面波成像适用于声波、超声波、以及微波多种探测方法，采用一次发射柱面波，覆盖整个成像区域，成像帖率高，能够满足工业、石油、矿产、医学等多个领域采用环形阵列激励和检测的高帖率成像的需求。
附图说明
[0041]
图1：柱面波成像原理示意图；
[0042]
图2：二维水平断面成像实施例模型示意图；
[0043]
图3：二维水平断面成像实施例激励单元、接收单元与成像区域分布示意图；
[0044]
图4：二维水平断面成像实施例计算单个接收单元在单个成像目标点的成像值的方法示意图；
[0045]
图5：二维水平断面成像实施例成像结果，每个接收单元在成像目标点的成像值不施加归一化变换；
[0046]
图6：二维水平断面成像实施例成像结果，每个接收单元在成像目标点的成像值施加归一化变换(8)；
[0047]
图7：二维水平断面成像实施例成像结果，每个接收单元在成像目标点的成像值施加归一化变换(9)；
[0048]
图8：二维水平断面成像实施例成像结果，每个接收单元在成像目标点的成像值施加归一化变换(10)；
[0049]
图9：二维水平断面成像实施例成像结果，每个接收单元在成像目标点的成像值施加归一化变换(11)。
具体实施方式
[0050]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅为本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域的普通技术人员在不付出创造性劳动的前提下所获得的所有其他实施例，都属于本发明的保护范围。
[0051]
根据本发明的一个实施例，提出一种柱面波成像方法，包括如下步骤：
[0052]
步骤1，通过激励单元产生柱面波，所述柱面波是通过环形阵列激励单元、圆弧形阵列激励单元、柱面阵列激励单元同步激励产生、或者通过相控阵激励单元利用波束成形技术产生；柱面波成像中的激励波阵面为柱面，或者为位于柱面内的不规则形状，柱面波传播范围覆盖整个成像区域；
[0053]
柱面波成像中的波阵面的传播方向可以向内，也可以向外。波阵面向内传播适合成像区域在初始波阵面内部的情况，例如工业管道检测；波阵面向外传播时适合成像区域在初始波阵面外侧的情况，例如测井、人体食管和血管检测。
[0054]
步骤2，通过多个接收单元采集波形原始信号，所述接收单元与激励单元复用，或与激励单元分离、并根据成像需求布置。
[0055]
步骤3：将成像区域剖分为若干成像目标点。根据成像区域大小、维度、成像分辨率等要求将成像区域剖分为若干网格，网格尺寸小于或等于成像分辨率要求，每个网格为一个成像目标点。
[0056]
步骤4：选择数学模型计算每个成像目标点的成像值，完成整个成像区域的柱面波成像。
[0057]
如果成像区域为柱面波传播范围内的任意体或者面，成像区域内单个成像目标点的成像值的数学模型由公式(1)和(2)表示：
[0058][0059]
其中，成像区域剖分为m个成像目标点，(xj,yj,zj)的表示第j个成像目标点的坐标，j＝1,
…
,m；共采用n个接收单元，(xi,yi,zi)表示第i个接收单元的坐标，i＝1,
…
,n；r表示初始波阵面的半径，c为波速，t
ij
表示柱面波经过第j个成像目标点散射到达第i个接收单元需要的传播时间；s(i,t
ij
)表示第i个接收单元在第j个成像目标点的成像值，通过对第i个接收单元采集的信号在t
ij
时刻插值计算得到；将n个接收单元在第j个成像目标点的成像值叠加，得到第j个成像目标点的成像值vj。
[0060]
计算第i个接收单元在第j个成像目标点的成像值s(i,t
ij
)时，具体插值方法与激励波形特征、采样率等相关，一般采用邻近线性插值。根据实际需要，可以提高插值阶数，选择二阶拉格朗日插值、或者三次样条插值等多种插值方法。
[0061]
当成像区域为一个二维水平断面时，成像目标点的成像值的数学模型由公式(1)、(2)和(3)表示：
[0062]
zj＝z0ꢀꢀꢀ
(3)
[0063]
公式(3)对所有j＝1,
…
,m成立。数学模型(1)、(2)和(3)是数学模型(1)和(2)的一个特例，对成像目标点增加约束(3)，表示所有成像目标点位于高度为z0的二维水平断面内。需要注意，接收单元并不要求位于成像断面内，可以根据实际情况布置。特别地，如果接收单元也位于z0水平断面内，则公式(2)可以进一步简化为：
[0064][0065]
公式(1)、(2)和(3)表示的数学模型适合接收单元布置在成像水平断面外的二维水平断层成像；公式(1)、(3)和(4)表示的数学模型适合接收单元位于成像水平断面内的二维水平断层成像。
[0066]
当成像区域为一个过柱面中心轴的二维竖直断面时，作为公式(1)和(2)表示的数学模型的一个特例，成像目标点的成像值的数学模型由公式(1)、(2)和(5)表示：
[0067][0068]
其中tg-1
表示反正切函数，公式(5)对所有j＝1,
…
,m成立，表示所有成像目标点位于扇角为的二维竖直断面内。同样地，接收单元并不要求位于成像断面内，可以根据实际情况布置。
[0069]
当成像区域为一个圆柱面时，作为公式(1)和(2)表示的数学模型的一个特例，成像目标点的成像值的数学模型由公式(1)、(2)和(6)表示：
[0070]
[0071]
公式(6)对所有j＝1,
…
,m成立，表示所有成像目标点位于半径为r0的柱面内。同样地，接收单元可以根据实际情况布置。
[0072]
成像目标点的成像值由各接收单元在该成像目标点的成像值直接叠加得到，或者对各接收单元在该成像目标点的成像值施加数学变换后再叠加得到，即前述数学模型的公式(1)可以由公式(7)替换：
[0073][0074]
其中，f{s(i,y
ij
)}表示对第i个接收单元在第j个成像目标点的成像值施加数学变换f。
[0075]
所述数学变换之一为归一化变换，用公式(8)-(11)表示：
[0076][0077][0078][0079][0080]
其中，公式(8)和(9)为有极性归一化，即归一化后的成像值仍保留正负极性；公式(10)和(11)为无极性归一化，即归一化后的成像值失去正负极性，只保留幅值特征。根据实际情况选择具体使用哪一种规一划变换。
[0081]
所述数学变换的另一种方式为系数加权法，用公式(12)表示：
[0082][0083]
其中，公式(12)中的w
ij
表示第i个接收单元在第j个成像目标点的成像值的加权系数，根据成像目标点与接收单元的之间相对位置关系、接收单元的接收口径与灵敏度函数、以及柱面波场的物理特征等确定。通常，将接收单元对成像目标点所张的立体角作为权系数；或者根据接收单元的指向角确定，成像目标点位于接收单元指向角范围内，权系数为1；成像目标点位于接收单元指向角范围外，权系数为0。
[0084]
上述两类数学变换可以同时使用，或者使用其中之一，或者均不使用。
[0085]
根据本发明的又一实施例，如图1所示为柱面波成像原理示意图，柱面波往中心方向传播，当柱面波碰到散射体后，散射体成为次波源，向外散射波，被接收单元采集。
[0086]
本实施例中，采用超声波。参考图2，采用环形超声换能器阵列同步激励以产生柱面波，环形阵列半径为6.5cm，超声换能器等间距排列。成像目标为圆筒型实物模型，成像目
标中心与超声换能器环阵中心重合，外半径为3cm，内半径为2.5cm。接收单元与激励单元复用超声换能器。采用同步激励产生柱面波，柱面波向内传播，传播速度为1500m/s。接收单元同时采集超声数据。
[0087]
成像区域为二维水平断面，并且接收单元与成像区域位于同一断面内。因此，成像目标点的成像值的数学模型由公式(1)和(4)表示，其中，成像区域剖分为m个成像目标点，在本实施例中，m＝160000；(xj,yj)的表示第j个成像目标点的坐标，j＝1,
…
,m；共采用n个接收单元，在本实施例中，n＝128；(xi,yi)表示第i个接收单元的坐标，i＝1,
…
,n；r表示初始波阵面的半径，在本实施例中，r为6.5cm；c为波速，在本实施例，c为1500m/s。
[0088]
具体地，如图3所示，成像区域为中心位于环阵圆心、边长为6cm的正方形区域。成像区域剖分为400x400网格，每个网格边长为0.15mm，每个网格即为一个成像目标点，共m＝160000个成像目标点。
[0089]
坐标系如图4所示，坐标原点位于环阵中心。成像目标点从左往右、从下往上编号，超声换能器以x坐标轴为起点逆时针开始编号，图中标出了第1个成像目标点与第1个超声换能器。根据公式(4)，分别取i＝1和j＝1，第1个成像目标点坐标为(-0.03,-0.03)，第1个接收单元坐标为(0.065,0)，计算得到，柱面波经过第1个成像目标点散射到达第1个接收单元的传播时间，
[0090][0091]
采用同样的方法，可以依次计算柱面波经过第j个成像目标点散射到达第i个接收单元的传播时间t
ij
，然后采用邻近线性插值方法计算第i个接收单元采集的原始波形在t
ij
时刻的值，即为第i个接收单元在第j个成像目标点的成像值s(i,t
ij
)。
[0092]
将第i个接收单元在第j个成像目标点的成像值s(i,t
ij
)施加系数加权的数学变换公式(12)，加权系数根据接收单元的指向角确定，成像目标点位于接收单元指向角范围内，权系数为1，成像目标点位于接收单元指向角范围外，权系数为0。然后叠加，成像结果如图5所示，真实的模型位置在图中用白色虚线表示。从图5可以看出，成像结果很好地反映实物模型的位置及内外边界。
[0093]
将第i个接收单元在第j个成像目标点的成像值同时施加归一化变换(8)和系数加权变换(12)，然后叠加，成像结果如图6所示。
[0094]
将第i个接收单元在第j个成像目标点的成像值同时施加归一化变换(9)和系数加权变换(12)，然后叠加，成像结果如图7所示。
[0095]
将第i个接收单元在第j个成像目标点的成像值同时施加归一化变换(10)和系数加权变换(12)，然后叠加，成像结果如图8所示。
[0096]
将第i个接收单元在第j个成像目标点的成像值同时施加归一化变换(11)和系数加权变换(12)，然后叠加，成像结果如图9所示。
[0097]
从图5-图9可以看出，不施加规一划变换或者施加四种规一划变换之一，成像结果
均与实物模型相符。本实施例中，同时施加归一化变换(10)和系数加权变换(12)的成像结果对比度最佳。但是，实际应用时，需要比较不同归一化变换的结果确定哪一种最佳。
[0098]
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，且应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。