超声层析成像理论与实现

超声层析成像的理论与实现•英国从事超声成像的专家P.N.TWells在2000年的文章《超声成像技术的现状与未来》一文中指出：“在最近的十几年里，有关超声成像技术的研究在医学成像领域至少占百分之二十五以上的份额，并且这种趋势还在继续增长。”•Wells还指出：“目前成功地应用于医学领域的超声成像设备大都是基于反射波，且其成像也只是定性的，根据超声散射波的信息，定量地生成人体内部的结构图，是超声应用技术的研究者追求的新目标。”“未来的超声成像技术应该是制造出不需成像专家或医学专家才能识别的反映客观现实真实图像的超声成像设备，即使是这种设备是不完美的。”主要内容•一.超声层析成像技术的发展历史•二.超声层析成像技术的基本模型及方法•三.问题的不适定性及其正则化•四.模型噪声的判断方法——Picard准则•五.静态正则化技术在超声层析技术中的应用•六.迭代正则化技术在超声层析技术中的应用•七.总结与展望一.超声层析成像的发展历史•1.折射系数层析成像方法•2.衰减系数层析成像方法•3.射线跟踪方法•4.透射式衍射层析成像及反射式衍射层析成像方法•5.基于精确场描述的层析成像方法1.折射系数层析成像方法[1(,)]BwAnxydsctRefractive-indextomography2.超声衰减系数层析成像Attenuationtomography•衰减系数(,,)xyf•综合衰减系数0(,)xyf0(,)BAxyds21ln[/]EE3.射线跟踪方法RayTracingMethod4.透射式衍射层析成像及反射式衍射层析成像方法物体傅里叶变换频域空域入射波前向散射场从不同方向照射物体时,前向散射场数据的傅里叶变换5.基于精确场描述的层析成像方法二.超声层析成像技术的基本模型及方法非齐次亥姆霍兹方程(HelmholtzEquation)220()()()()kprorpr2220()()1()/()1ornrcrcr0sprprpr•1.波动方程及其解全场方程(TotalFieldEquation)(第二类Fredholm积分方程)0SprprGrrorprdr散射场方程(ScatteringFieldEquation)探测器方程(DetectorEquation)sSprGrrorprdr2.积分方程的离散化─矩量法''0()()jjiiijjprprorprc''()()jjsmmjjprorprd2010201000()22()()2ijmjimjjkaHkaji=jcdjkaimjJkaHkR若若向量形式:()()()tintPPCOP()()stPDOP3.波动方程的近似•①Born近似0sprpr()()tinPP()()sinPDOPBorn逆解O应满足的条件:/4an②Rytov近似()()rpre0()()()srrr2()()or2()2snr2（）应满足的条件:''''001()()()()sGrrprordrpr4.基本方法•Born迭代算法(BI)•Levenberg-Marquardt和Newton-Kantorovich方法•变形Born迭代方法(DBI)Born迭代算法(BI)求Born逆解O由全场方程确定全场()()()tintkPPCOP()tkP由散射场方程求散射场,并计算()()stkkkPDOP()()()ssskkPPP()skP由方程求改变量()()stkkkPDOPkO求1kkkOOO求Born逆解O由全场方程确定全场()()()tintkPPCOP()tkP由散射场方程求散射场,并计算()()stkkkPDOP()()()ssskkPPP变形Born迭代算法(DBI)根据最新求得的Ok改变散射方程的系数矩阵D()skP求1kkkOOO由方程求改变量()()stkkkPDOPkOLevenberg-Marquardt和Newton-Kantorovich方法()()()tintPPCOP()()stPDOP代入()()()[]sintPDOPCOP三.问题的不适定性及其正则化•适定性问题是指：对于连续算子方程Kx=y，如果解x满足：(1).存在；(2).唯一；(3).连续地依赖于数据y。否则，即上述三个条件有一个不满足，则称其为不适定的（Ill-posed）。离散不适定问题(DiscreteIll-PosedProblem)2minbAx若:(1).矩阵A的条件数非常大，或者说矩阵A的最大奇异值和最小奇异值之比非常大；(2).矩阵A的奇异值逐渐下降趋于零。对于线性方程组Ax=b或最小二乘问题:Tikhonov正则化222022min()AxbLxxL=In，x0=0时，称为Tikhonov正则化的标准形式,其解可表示为:2221,()niiiiiiubxv四.模型噪声的判断方法:Picard准则•离散Picard准则：若方程组Ax=b的傅里叶系数趋于零的速度在平均意义下快于矩阵A的奇异值趋于零的速度的话，则称该方程组满足离散Picard准则(条件)。1,nilsiiiubxv,iub2221,()niitikiiiiubxv最小二乘解:Tikhonov正则化解:受噪声污染和无噪声污染的Picard图污染严重污染较轻A对比度为30％时对比度为20％时对比度为10％时五.静态正则化技术•1.截断奇异值分解正则化方法TruncatedSingularValueDecomposition(TSVD)•2.截断完全最小二乘正则化方法TruncatedTotalLeastSquares(TTLS)1.截断奇异值分解正则化方法(TSVD)对于线性方程组Ax=b或最小二乘问题2minbAx1,nilsiiiubxv2221,()niitikiiiiubxv最小二乘解:Tikhonov正则化解:1,kitsvdiiiubxvTSVD正则化解:正则化参数的选取方法•离差原理(DiscrepancyPrinciple)方法•广义交叉验证(GCV)方法•L曲线(L-Curve)方法10logbAx10logx减小时增加时由L曲线方法确定k’采用一维搜索的方法确定更精确的kTSVD方法的数值仿真结果BACDE对比度为10％时对比度为20％时对比度为30％时原始图像迭代过程的相对误差和相对残差曲线迭代过程的相对误差和相对残差曲线2.截断完全最小二乘正则化方法满足：2minAxbAxb最小二乘问题：min(,)(,)FAbAb完全最小二乘问题：满足：Axb截断完全最小二乘的步骤1.首先，计算增广矩阵(A，b)的奇异值分解：11(,)nTiiiiAbUVuv2．确定截断参数k≤min(n，rank(A，b))使得：1kk221,11,1(,)0nknnVvv3.记q=n-k+1,将矩阵分块11122122VVVVV4.则完全最小二乘问题的解为：212221222222TkxVVVVVTTLS方法的数值仿真结果原始图像对比度为10％时对比度为20％时对比度为30％时迭代过程的相对误差和相对残差曲线六.迭代正则化技术1.求解最小二乘问题的共轭梯度方法(cgls)2.LSQR方法1.求解最小二乘问题的共轭梯度方法(cgls)将共轭梯度法应用于法方程TTAAxAb10000(;),,,kTTTTTTTkAAArspanArAAArAAAr相当于在Krylov子空间：产生的序列xk，使得：22*2()TkkkAAfxxxAxbcgls方法的解可表示为:2211()Tnikikiiiibuxqv1()kq2的k-1次多项式，其系数的确定是其中：依赖于：(1).方程的右侧项b的特征；(2).矩阵A的奇异值的分布；(3).迭代的次数迭代次数增加，残差变化不大，但解的范数受影响较大正则化参数对迭代的影响cgls方法的数值仿真结果BACDE对比度为10％时原始图像对比度为20％时对比度为30％时迭代过程的相对误差和相对残差曲线图5.15采用clgs方法，五种不同图像在对比度为30％时的相对残差(RRE)曲线，cgls迭代次数为10LSQR迭代方法Lanczos三对角过程LanczosTAA应用于将矩阵A双对角化Golub和Kahan(1965)Paige和Saunders(1982)2minbAx线性方程组Ax=b和应用于LSQR方法的优点:1.速度快2.对不适定性问题数值稳定3.从迭代过程很容易求得数值分析的数值BACDE原始图像LSQR方法的数值仿真结果对比度为10％时对比度为20％时对比度为30％时七.总结与展望•首先利用Picard理论，分析了超声层析成像问题的中的模型噪声问题，给出了入射波的确定方法、以及正则化方法的适用范围的判断方法。•采用了两类四种正则化方法对超声层析成像问题中的不适定性问题进行了研究，通过对正则化参数选择的修正，完成了较大对比度物体的成像问题。•结论：静态正则化方法数值稳定，但速度慢；迭代正则化方法速度快，但数值稳定性不如静态方法。今后需要进一步研究的工作1.前向散射问题的研究（波动方程的精确程度）2.离散化方法有限元法、边界元法矩量法中基函数的确定3.正则化问题基于非对成方程的Krylov子空间方法谢谢！