BPT算法对SIRT算法的一种加权研究.docVIP

BPT算法对SIRT算法的一种加权研究 摘要：对声波层析成像中的SIRT算法提出了一种加权算法。先通过BPT算法得到慢度分布，以之作为SIRT算法的初始值，以空间各点的慢度值与平均慢度值之差作为权系数对SIRT算法中射线矩阵进行加权。算例表明，该改进算法保持了SIRT算法收敛特点的同时，能够加快SIRT算法的收敛速度，提高计算效率，尤其能提高计算精度，是一种可行的加权方法。 关键词：层析成像；SIRT；反演；加权 中图分类号：TU528 引 言 在声波层析成像检测技术（CT）中为得到比较清晰的图像，一般采用迭代重建的方法。迭代重建方法通常又分为矩阵求解方法和代数求解方法两种形式。其中，代数求解方法具有原理简单，容易编程实现的优点，从而应用较广。 代数类求解方法主要是代数重建法(ART)及联合迭代重建法(SIRT)。ART算法具有计算速度快、计算效率高的优点，但缺点是对测量误差比较敏感，对初始值的依赖性比较高，且它有一最佳迭代次数，超过一定次数则不收敛；SIRT算法具有收敛性好、对初值的要求比较低、对测量误差不敏感的特点，从而在层析成像中应用较广，但缺点是计算效率比ART算法低得多[1-3]。 另外还有一种非迭代的反投影法（BPT）T在一侧某点发射，接收换能器R在另一侧接收，使每个网格都有1条及以上的测线通过，利用测量得到的走时通过反演技术获得测区各部分的波速分布图，从而确定异常区域的位置、尺寸以及异常区域本身的波速等。如图1所示[5-6] 图1 测区网格划分示意图 设图像的整个区域划分成个成像单元，测区有条测线通过，为第条测线的总长度，则第条测线从激发点T到接收点R的走时为： （1） 式中：—第个成像单元的波速； —第个成像单元的慢度（波速的倒数）。 假定成像单元足够小，可将每个成像单元的视为常数，设第条射线在第个成像单元内的线段长度为，则式(1)可离散成如下级数形式： （2） 2 SIRT算法的改进 从数学角度看，式(2)实际上是一个线性方程组，并可以写成矩阵方程的形式： （3） 其中矩阵为走时矩阵；A为射线矩阵，是第条射线在第个成像单元内的线段长度的集合；为所要反演单元的慢度矩阵。 由于被检测对象的结构、环境等客观条件的限制，采集得到的投影数据往往存在严重缺失，所以方程(3)具有稀疏性、欠定性、不相容性等特点，在数学上直接对其求解时非常困难，一般采用迭代的方法，根据测量得到的射线走时，寻求使取得最小值的向量。目前使用的较为广泛的代数类求解算法是联合迭代重建法(SIRT)。 在使用代数类求解算法之前，一般先采用BPT算法生成迭代所需的慢度初始值，这样可以加快计算速度。BPT反投影技术的计算公式为[8]： （5） 在得到慢度初值的基础上，SIRT算法利用某一成像单元所有通过它的射线的修正值来确定此单元的平均修正值，第个成像单元的慢度迭代修正增量表示为[9]： （6） 式中为矩阵A中第列元素中非零元素的个数，即表示穿过第个单元格的射线的条数，为迭代计算次数。为第条射线次迭代后的走时误差。 虽然SIRT算法总是收敛的，但其收敛速度偏慢，计算精度有待提高。如果能对射线矩阵进行加权，在速度正常区域添加较小的权系数，而在速度异常区域添加较大的权系数，就有可能加快收敛速度，提高计算精度。 本文选择权系数对每个成像单元进行加权： （7） 式中表示第个成像单元的慢度值，表示所有成像单元慢度的平均值。两者通过BPT算法得到。 对于待测量的对象来说，速度异常单元通常只占一小部分，所有成像单元慢度的平均值与正常速度区域的慢度值相接近，而与异常速度区域的慢度值相差较大。（7）式表明在速度异常单元添加了较大的权重，在每次迭代过程中的修改量变大，而在正常速度单元添加了较小的权重，迭代过程中的修改量较小，从而有可能使算法的收敛速度加快。 加权后SIRT算法慢度的迭代修正增量公式为： （8） 慢度的修正公式为 （9） 方程(8)和(9)即为改进算法的基本迭代公式。 由于涉及到权系数的计算，因此与常规SIRT算法在初始值的要求上有所不同。改进后SIRT算法的具体计算过程可概括为以下三步： 第一步，根据（5）式，由测量得到的值利用BPT算法产生迭代所需的初始慢度。 （10） 由计算出所有成像单元慢