本发明专利涉及无损检测超声相控阵、数据稀疏优化及图像处理领域，特别是一种基于阵元稀疏优化及新边缘指导插值算法的超声相控阵稀疏全聚焦成像方法。
背景技术
超声相控阵检测是按照一定的延迟时间通过控制传感阵列中的各个阵元，有规律地激励和接收超声波，实现超声波束的偏转和聚焦，对工件内部的缺陷进行无损检测。在航天、核电、复合材料等领域关键设备的无损检测，超声相控阵以其检测速度快、灵敏度高、对象适应能力好被广泛使用。目前采用全矩阵采集(fullmatrixcapture,fmc)及全矩阵数据的全聚焦(totalfocusmethod,tfm)成像方法，实现了被测区域内任意位置的聚焦成像，大大提高了图像分辨率。如申请公开号cn106802323a，申请公开日2017年6月6日的专利文献公布了一种基于全矩阵数据的超声全聚焦成像系统，都采用超声相控阵探头，使用全矩阵采集所有阵元的超声数据，可实现高精度成像，提高图像分辨率。但由于全聚焦成像采用的全矩阵存在数据量大及数据冗余，使得计算很费时，成像效率低，因此一定程度上限制了其在工业无损检测上的实际使用。
为解决此问题，可以通过改进硬件结构来提高处理速度，如采用并行计算、多fpga等，如申请公开号为cn106651740a，申请公开日2017年5月10日的专利文献公布了一种基于fpga的超声全数据聚焦快速成像方法及系统，采用fpga提高了全聚焦成像的速度，但这种方法需要的硬件成本较高，难以在众多工业检测领域推广。
技术实现要素：
本发明针对现有相控阵全聚焦成像法中存在全矩阵数据量大、计算费时、成像效率低的问题，提出一种基于阵元稀疏优化及新边缘指导插值的超声相控阵稀疏全聚焦成像方法，在保证高质量成像的同时，极大地缩减了计算时间，提高超声相控阵全聚焦的成像速度。
为了实现上述技术目的，本发明的技术方案是，一种基于阵元稀疏优化及新边缘指导插值算法的超声相控阵稀疏全聚焦成像方法，包括以下步骤：
步骤一、换能器阵列阵元稀疏优化，减少相控阵换能器的有效发射/接收阵元数量，结合几乎差集与遗传算法，优化稀疏阵列单元位置。优化目标为最小化稀疏阵列的旁瓣峰值(peakside-lobe,psl)及主瓣宽度(main-lobewide,mlw)，其值由稀疏阵列指向性函数得到的其阵列方向性图中测得，具体关键步骤如下：
步骤1、种群初始化，根据所设计目标稀疏阵列总阵元数n，有效阵元个数k，从ads库中选择对应参数的几乎差集(n,k,λ,t)-ads及其相应的二值序列，对二值序列进行循环移位可产生它的n-1个转换类型，并根据这n个几乎差集对应的二值序列的psl值对它们排序，其中λ,t分别表示几乎差集的n阶阿贝尔群中非零元素的个数以及非零元素在复合集中出现的次数，选择排序靠前的p/2(p为迭代过程中每代种群的个体总数)个个体作为初始种群，即：ρp＝{bp(n)＝w(p)(n)；n＝0,…,n-1},1≤p≤p/2，式中，p为每代种群中某个个体的编号，ρp表示第p个个体对应的二值序列，n表示二值序列的位数，w(p)(n)为第p个几乎差集对应的二值序列中第n位的数值，若n∈d(σ)，即循环位移数为σ，且n为几乎差集转换类型d(σ)中的某个元素时，w(p)(n)＝1，否则w(p)(n)＝0；rp(n)为随机数0或1，bp(n)表示第p个种群个体中第n位的基因值，在bp(n)＝1处放置阵元，否则，不放置阵元，剩余p/2个初始个体随机产生，即：ρp＝{bp(n)＝rp(n)；n＝0,...,n-1},p/2≤p≤p式中rp(n)为随机数0或1；
步骤2、适应度函数评估。根据适应度函数fit＝k1*mlw+k2*psl对上述步骤1所得旁瓣峰值(peakside-lobe,psl)及主瓣宽度(main-lobewide,mlw)进行评估k1、k2为根据不同优化目标选定的系数值；
步骤3、选择、交叉、变异，对初始群体p(k)按适应度函数值进行排序，去掉阵列性能较差的一半群体，预先设置一定的交叉几率，将剩下的m/2个个体两两配对，并按照这一几率对各个组合的某些基因位进行调换，得到新的m/2个个体，完成该过程后群体中恢复为m个个体。再依据各变异几率将种群中个体的某些基因替换为相应的等位基因，即基因值由0变异到1或者由1变异到0，设稀疏率为v，b(m)、b(n)分别表示种群个体中第m、n位上的基因值，基因的位变异概率pbm(n)由式(1)计算：
步骤4.重复步骤2、3，直到循环完所设置的迭代次数，记录下每次循环中的最优二值序列，并将其还原为对应的阵列结构，得到最优稀疏阵列；
步骤二、稀疏阵列有效孔径修正，根据有效孔径一致的稀疏阵列与全阵列的在远场区域具有相同的点扩散函数为条件，选择合适的发射阵元权重函数ωt和接收阵元权重函数ωr，对稀疏后的各发射/接收阵元进行加权，使得稀疏阵列与全阵列的有效孔径相同；
步骤三、工件成像区域进行成像点稀疏并预成像，设工件成像区域原有成像点设置为2l×2s，2l和2s分别表示像素点的行数和列数，首先将成像区域的横向和纵向成像点数均减少一半，即参与稀疏全聚焦计算的成像点数为l×s，然后在所述步骤一和所述步骤二的基础上进行预成像，如对稀疏全聚焦预成像图ilr(x,z)取像素点l×s个，l和s分别表示像素点的行数和列数，依次激励第i(i＝1,2,...,n)个有效阵元，其它有效阵元接收信号,每次接收的发射－接收对回波信号数据保存为e1，1(t)，e1，2(t)，...，en，n(t)，如图2所示，将一维线性相控阵换能器放置于被测介质的表面，建立直角坐标系oxz，x轴沿阵列方向平行于介质表面，z轴垂直于阵列并指向试件内部，原点位于换能器中心位置。针对某一目标聚焦点(x,z)，将换能器中所有有效发射-接收阵元组合的超声回波信号在目标聚焦点(x,z)处叠加，获得表征该点信息的幅值i(x,z)。依次计算被测区域内每一个聚焦点的幅值，根据式(2)得到整个被测区域内的预成像：
式中(xtx,ztx)表示被激励阵元的坐标，(xrx,zrx)为接收阵元的坐标，nt、nr分别表示稀疏阵列中发射、接收阵元数目，ωt、ωr分别表示发射、接收阵列的权重函数，ei,j(t)是第i个阵元发送第j接收的信号，c表示超声波在介质中的传播速度；
步骤四、对步骤三得到的低分辨率的预成像图像进行新边缘指导插值成像，对稀疏全聚焦预成像图(大小为l×s)进行2×2插值计算，相应的新图像大小为2l×2s，原始像素点对应的像素值保持不变，插值按以下两个阶段进行：
第一阶段计算4个原始像素点插值得到的像素值。如图3所示，图中黑点表示已知像素点，白点表示待插值点，待插点的像素值为：
式中，i(2l+1,2s+1)表示待插值像素点像素强度，i(2(l+p),2(s+q))表示待插值像素点周围4个最近邻的像素点，p,q为待插值点距中心点(l,s)的横向、纵向偏移点数，α2p+q表示插值权系数。α2p+q根据α＝(cct)-cy，
第二阶段计算原始像素点和第一阶段所得像素点插值得到的像素点值，如图4所示，白点表示已知像素点，黑点表示待插值点，其像素点值为：
i(2l+1,2s)＝α0i(2l,2s)+α1i(2l+1,2s-1)+α2i(2l+2,2s)+α3i(2l+1,2s+1)，
按照上述类似的方法可得到插值权系数，最终求得待插点的像素值。
附图说明
图1为本发明的流程图；
图2为相控阵全聚焦成像示意图；
图3为新边缘指导插值的第一阶段示意图；
图4为新边缘指导插值的第二阶段示意图；
图5为本发明超声相控阵检测实验示意图
图6为试块全聚焦成像图；
图7为未插值的试块稀疏全聚焦成像图；
图8为本发明得到的试块稀疏全聚焦成像图；
具体实施方式
下面结合附图对本发明作进一步说明。
以标准相控阵b型试块上圆弧型分布的侧边孔为检测对象，如图5所示，检测区域在图中通过白色矩形框标记出来，该区域共有18个直径为2mm侧边孔，由上到下依次编号为1～18，实验使用广州多浦乐电子科技有限公司生产的5l32-0.6×10型超声相控阵换能器，换能器内阵元总数为32，各阵元间距为0.6mm，阵元宽度为0.5mm，换能器的中心频率为5mhz，超声波在试块中的声速为5900m/s。以100mhz的采样频率，按照全矩阵数据采集方式采集并存储全矩阵数据。成像时所用pc机的处理器为corei5，运行内存4gmb，成像区域为60mm×60mm。设置横向、纵向成像点间隔为0.4mm，将成像区域划分为151×151像素网格，对此成像区域进行基于阵元稀疏优化及新边缘指导插值的超声相控阵稀疏全聚焦成像，具体实施步骤如下：
步骤一、换能器阵列阵元稀疏优化，使用32个阵元中的8个作为有效发射阵元，接收使用全部的32个阵元，对8个有效阵发射阵元优化的具体关键步骤如下：步骤1、种群初始化，根据所设计目标稀疏阵列总阵元数32，有效阵元个数8，从ads库中选择对应参数的几乎差集(32,8,λ,t)-ads及其相应的二值序列，对二值序列进行循环移位可产生它的31个转换类型，并根据这32个几乎差集对应的二值序列的psl值对它们排序，其中λ,t分别表示几乎差集的32阶阿贝尔群中非零元素的个数以及非零元素在复合集中出现的次数，选择排序靠前的200个个体作为初始种群，即：ρp＝{bp(n)＝w(p)(n)；n＝0,...,31},1≤p≤200，式中，p为每代种群中某个个体的编号，ρp表示第p个个体对应的二值序列，n表示二值序列的位数，w(p)(n)为第p个几乎差集对应的二值序列中第n位的数值，若n∈d(σ)，即循环位移数为σ，且n为几乎差集转换类型d(σ)中的某个元素时，w(p)(n)＝1，否则w(p)(n)＝0；rp(n)为随机数0或1，bp(n)表示第p个种群个体中第n位的基因值，在bp(n)＝1处放置阵元，否则，不放置阵元，剩余200个初始个体随机产生，即：ρp＝{bp(n)＝rp(n)；n＝0,...,31},200≤p≤400，式中rp(n)为随机数0或1；
步骤2、适应度函数评估。根据适应度函数fit＝k1*mlw+k2*psl对上述步骤1所得旁瓣峰值psl及主瓣宽度mlw进行评估，k1、k2为根据不同优化目标选定的系数值，此时需同时优化阵列的主瓣特性和旁瓣特性，令k1、k2均为1；
步骤3、选择、交叉、变异，对初始群体p(k)按适应度函数值进行排序，去掉阵列性能较差的一半群体，预先设置交叉几率为0.3，将剩下的200个个体两两配对，并按照这一几率对各个组合的某些基因位进行调换，得到新的200个个体，完成该过程后群体中恢复为400个个体。再依据各变异几率将种群中个体的某些基因替换为相应的等位基因，即基因值由0变异到1或者由1变异到0，设稀疏率为v，b(m)、b(n)分别表示种群个体中第m、n位上的基因值，基因的位变异概率pbm(n)由式(1)计算；
步骤4.重复步骤2、3，直到循环完所设置的迭代次数，记录下每次循环中的最优二值序列，并将其还原为对应的阵列结构，得到最优稀疏阵列为“10000010000010011100010000000001”，其中数字1代表是有效阵元；
步骤二、稀疏阵列有效孔径修正，根据有效孔径一致的稀疏阵列与全阵列的在远场区域具有相同的点扩散函数为条件，选择合适的发射阵元权重函数ωt和接收阵元权重函数ωr，对稀疏后的各发射/接收阵元进行加权，使得稀疏阵列与全阵列的有效孔径相同；
步骤三、工件成像区域进行成像点稀疏并预成像，设工件成像区域原有成像点设置为302×302，首先将成像区域的横向和纵向成像点数均减少一半，即参与稀疏全聚焦计算的成像点数为151×151，然后在所述步骤一和所述步骤二的基础上进行预成像，如对稀疏全聚焦预成像图ilr(x,z)取像素点151×151，依次激励第i(i＝1,2,...,32)个有效阵元，其它有效阵元接收信号,每次接收的发射－接收对回波信号数据保存为e1，1(t)，e1，2(t)，...，e32，32(t)，如图1所示，将一维线性相控阵换能器放置于被测介质的表面，建立直角坐标系oxz，x轴沿阵列方向平行于介质表面，z轴垂直于阵列并指向试件内部，原点位于换能器中心位置。针对某一目标聚焦点(x,z)，将换能器中所有有效发射-接收阵元组合的超声回波信号在目标聚焦点(x,z)处叠加，获得表征该点信息的幅值i(x,z)。依次计算被测区域内每一个聚焦点的幅值，根据式(2)得到整个被测区域内的预成像：
式中(xtx,ztx)表示被激励阵元的坐标，(xrx,zrx)为接收阵元的坐标，ωt、ωr分别表示发射、接收阵列的权重函数，ei,j(t)是第i个阵元发送第j接收的信号，c表示超声波在介质中的传播速度；
步骤四、对步骤三得到的低分辨率的预成像图像进行新边缘指导插值成像，对稀疏全聚焦预成像图进行2×2插值计算，相应的新图像大小为302×302，原始像素点对应的像素值保持不变，插值按以下两个阶段进行：
第一阶段计算4个原始像素点插值得到的像素值。如图2所示，图中黑点表示已知像素点，白点表示待插值点，待插点的像素值为：
式中，i(2l+1,2s+1)表示待插值像素点像素强度，i(2(l+p),2(s+q))表示待插值像素点周围4个最近邻的像素点，p,q为待插值点距中心点(l,s)的横向、纵向偏移点数，α2p+q表示插值权系数。α2p+q根据α＝(cct)-cy，
第二阶段计算原始像素点和第一阶段所得像素点插值得到的像素点值，如图3所示，白点表示已知像素点，黑点表示待插值点，其像素点值为：
i(2l+1,2s)＝α0i(2l,2s)+α1i(2l+1,2s-1)+α2i(2l+2,2s)+α3i(2l+1,2s+1)，
按照上述类似的方法可得到插值权系数，最终求得待插点的像素值。
图6为采用全聚焦算法的试块成像结果，图7为步骤三得到的预成像结果，即未插值的试块稀疏全聚焦成像图，图8为本发明得到的试块稀疏全聚焦成像图。与全聚焦成像结果相比，本发明方法的成像时间由17.2分钟减小到1.08分钟，成像效率提高了约16倍，而与全聚焦方法的精度误差仅相差很小，18个孔的平均定量误差仅1.09％。