离散正则化（共7篇）
离散正则化 篇1
1 引言及主要结论
随着科技日新月异的发展,不适定问题[1]广泛出现在图像恢复、数学物理等科学领域[2,3],
实际上,上述问题转化为线性离散不适定问题,即求解
其中b = btr ue+ e是实际观测到的带有噪音的右端项,btr ue是理想状态下的输出结果,e是观测过程中无法消除的噪音数据。而我们真正需要求解的线性系统为:
因此,求解的难度就在于如何由已知的观测系统(1.1)来求解未知的xtr ue。
2011 年,Hochstenbach[5]和Reichel提出了Fractional Tikhonov正则化方法。本文在此基础上提出基于大规模系数矩阵A的投影Fractional Tikhonov正则化方法,进一步考虑广义的Arnoldi Fractional Tikhonov正则化方法,并提供了大量数值试验结果,数据表明新方法的有效性。
2 Fractional Tikhonov正则化方法简单描述及Arnoldi Fractional Tikhonov正则化方法
首先,我们考虑这样一个极小化问题:
其中定义,这里矩阵W是对称半正定矩阵。一般地,矩阵W取如下形式:
其中 α > 0 。若0 < α < 1 ,则矩阵W可以借助AAT的广义逆来表示。参数 α 的介入,会使得(2.1)式的解xλ,α更接近线性系统(1.1)的精确解。我们把上述方法称作Fractional Tikhonov正则化方法或Weighed Tikhonov正则化方法。对系数矩阵进行奇异值分解有
那么Fractional Tikhonov正则化解表示如下:
其中x = Vy 。Fractional Tikhonov正则化方法的滤波因子为
参数 λ 利用广义交叉准则确定。
Fractional Tikhonov正则化方法建立在对系数矩阵A进行奇异值分解的基础上,然而对于大规模矩阵,它的奇异值分解需要极大的计算量,因此我们选择先将大规模问题投影到维数较低的Krylov子空间上,对小规模问题进行求解。Lewis和Reichel在2008年提出Arnoldi Tikhonov方法[6],详细地介绍了该类方法。若考虑对投影后的小规模问题采用Fractional Tikhonov正则化方法,即本文提出的Arnoldi Fractional Tikhonov正则化方法,简记成算法AFT。算法过程如下:
算法3.1 Arnoldi Fractional Tikhonov方法
1)选初始向量x0,计算残量r0=b-Ax0,;
3) 利用Fractional Tikhonov正则化方法求解最小二乘问题,记解为yλ。
4) 计算近似正则解xλ= x0+ Vmyλ。
3 广义Arnoldi Fractional Tikhonov正则化方法
考虑如下广义形式Tikhonov正则化
这里L称为一个正则化算子(矩阵)。这里的范数是向量2范数。特别地,当L取单位矩阵,即L = I时,回归到标准形式Tikhonov正则化。
这里假设L的选择使得
其中N(A) 表示矩阵A的零空间。于是其正规方程组为
从而得到其唯一解
其中正则化参数 λ 决定了解对误差的敏感度和对精确解逼近的准确度。
通常使用中,正则化矩阵L的选择分别表示一阶导数算子和二阶导数算子,有如下两种形式
如果矩阵A和L是中小型规模矩阵,利用矩阵对{A,L}的广义奇异值分解(GSVD:Generalized Singular Value Decomposition)直接求解。这里我们考虑的是A和L为大规模矩阵,以至于不能进行广义SVD分解的情形。
2012 年, Reichel等针对大规模线性离散不适定问题,将Tikhonov正则化方法作进一步的研究,提出了矩阵{A,L} 的广义Arnoldi约化,建立了一种在广义Krylov子空间的Tikhonov正则化方法[7]。通过广义Arnoldi约化,大规模不适定问题转化成下述小规模问题
对矩阵HL(1:βk,1:k)作QR分解,有
令,代入(4.2)就可以得到Tikhonov正则化方法的标准形式:
对极小化问题(3.2)采用经典的Tikhonov正则化方法求解,得到正则解wλ(k)。再由得到原问题的解,这是已有的Generalized Arnoldi Tikhonov正则化方法[8]。
同样,对极小化问题(3.2)采用Fractional Tikhonov正则化方法求解,得到正则解w(k)λ。再由得到原问题的解,即得新方法—Generalized Arnoldi Fractional Tikhonov正则化方法。试验中选取的矩阵L∈Rn×n,
4 数值试验
本节对提出的所有改进算法都进行了数值试验,并与已有的算法进行了数值比较,以显示其优越性。本节数值试验从matlab(2010b)工具箱“regularization tool”[9]中得到系数矩阵A ,右端向量b和精确解x . 算例中的 α 的选取与[10]文中的取法类似。
例1 求解方程Ax = b ,其中系数矩阵A分别取baart、foxgood、gravity、i_laplace(2)、phillips、shaw等6类矩阵,矩阵阶数N取1000,投影到的子空间维数m = 10 ,取参数 α = 0.2,0.4 ,分别取右端扰动项 ν = 1%,5%,10% ,其中。比较Arnoldi-Fractional Tikhonov正则化方法(简记为AFT)与Arnoldi Tikhonov正则化方法(简记为AT)的计算解与精确解相对误差范数结果如表1、2所示:
由表1、2 可以看出,对于大规模线性离散不适定问题,采取Arnoldi-Fractional Tikhonov正则化方法求得的解的精度要比Arnoldi Tikhonov正则化方法高,数值试验说明,该方法是有意义的。
例2求解方程Ax=b,其中系数矩阵A分别取baart、fox⁃good、gravity、i_laplace(2)、phillips、等5类矩阵,矩阵阶数N取1000,投影到的子空间维数m=5,取参数α=0.2,0.5,分别取右端扰动项ν=1%,5%,10%,其中。比较General⁃ized Arnoldi Fractional Tikhonov正则化方法(简记为GAFT)与Generalized Arnoldi Tikhonov正则化方法(简记为GAT)的计算解与精确解相对误差范数结果如表3所示:
由表3 的结果知,GAFT的计算结果精度比GAT的计算结果高,但是参数 α 的选取会直接影响到求解精度,并且不同矩阵所取的最合适的 α 不一样。
根据以上三组数值试验的结果,可以看出:在Fractional Tikhonov正则化方法基础上给出的Arnoldi-Fractional Tikhonov和Generalized-Arnoldi- Fractional Tikhonov正则化方法在求解大规模线性离散不适定问题时是具有一定优势的,计算结果的精度更高一些。
5 结束语
对于大规模不适定问题的求解难度在于其系数矩阵的奇异值分解计算量过大。首先将其投影到小规模子空间上,再对投影问题的求解采用Fractional Tikhonov正则化方法,并推广到广义情形。数值试验表明新方法是具有一定优势的。但其中参数 α 的值影响解的精度,如何确定合适的 α 是下一步需要解决的问题。
参考文献
[1]黄光远,刘小军.数学物理反问题[M].济南:山东科学技术出版社,1993.
[2]Kirsch A.An Introduction to the Mathematical Theory of In-verse Problems[M].New York:Springer-Verlag,1996.
[3]王彦飞.反演问题的计算方法及其应用[M].北京:高等教育出版社,2007.
[4]Tikhonov A N.Solution of incorrectly formulated problems andthe regularization method[J].Soviet.Math.Dokl.,1963(4):1035-1038.
[5]Hochstenbach M E,Reichel.L.Fractional Tikhonov regular-ization for linear discrete ill-posed problems[J].BIT,2011(51):197-215.
[6]Lewis B,Reichel L.Arnoldi-Tikhonov regularization methods[J].J.Comput.Appl.Math.,2009,226:92-102.
[7]Wilkinson J H.The Algebraic Eigenvalue problem[J].OxfordUniversity Press.,1965:104-108.
[8]Reichel L,Sgallari F,Ye Q.Tikhonov regularization based ongeneralized Krylov subspace methods[J].Applied NumericalMathematics.,2012(62):1215-1228.
[9]肖庭延,于渗根,王彦飞.反问题数值解法[M].北京:科学出版社,2003.
[10]Morikuni K,Reichel L,Hayami K.FGMRES for linear dis-crete ill-posed problems[J].Appl.Numer.Math.,2014(75):175-187.
离散正则化 篇2
图像的全变差模型是Rudin等人于1992年提出的, 最初是用于图像去噪处理的［3］。1998年, Chan等人首次将全变差模型运用到了图像的去模糊复原处理中, 提出了一种基于正则化的全变差图像去模糊方法［4］。2008年Money等人将shock滤波器与全变差模型结合在一起, 提出了一种半盲的模糊图像复原方法［5］。该方法在对模糊图像进行Shock滤波器的预处理之后, 将预处理后的图像作为初始值代入到基于正则化的全变差模型中进行复原, 并运用了不同的范数对图像进行复原。2009年, Beck等人提出了一种基于图像梯度约束的全变差图像去模糊复原算法［6］。2010年, Chantas等人提出了一种基于权重全变差的贝叶斯图像复原方法［7］。该方法可等价于各向异性扩散的某种形式, 具有较好的抗噪性和边缘保护能力。2010年Chen等人提出了一种交互式的最小化策略来实现非局部的图像复原［8］。同年, Zhang等人也提出了一种非局部的图像复原方法［9］。2010年, Afonso等人利用变量分裂策略对全变差模型进行了有效的求解, 提出了一种快速的图像去模糊复原方法［10］。
虽然全变差模型能够有效地保护图像的边缘, 但是在现有的方法中, 全变差模型都只是等价地利用了图像水平方向和垂直方向, 两个方向的梯度信息［4—8，10］, 因此, 并不能充分利用图像中所有方向上的梯度信息, 不能将图像边缘的稀疏特性完全、准确地表现出来, 从而降低复原图像的质量。
针对现有方法存在的缺陷, 本文提出了一种权重的稀疏性正则化约束图像复原方法。本论文首先运用差分算子, 得到图像中四个方向上的梯度信息; 然后, 利用提取的梯度信息, 计算出图像边缘四个方向上各自的权重, 提出了一种多方向的权重的全变差模型对复原图像进行权重的稀疏性正则化约束; 最后, 运用了一种布雷格曼迭代 ( Bregman iteration, BI) 策略对提出的方法进行最优化求解。实验结果表明, 较近几年的一些具有代表性的图像复原方法相比, 不仅主观的视觉效果得到了较为明显的改进, 而且客观的信噪比增量也增加了0. 3 ~2. 5 dB。
1权重的稀疏性正则化图像复原方法
图像的模糊退化模型可由下式表示:
式 ( 1) 中, y为观察到的模糊噪声图像, 大小为p × q, h为已知的模糊退化函数 ( 也称为点扩散函数PSF) , x为p × q大小的原始清晰图像, n为零均值、 方差为 σ2的加性高斯噪声, * 表示卷积运算。本论文的目的就是在已知y和h的情况下, 求解出原始清晰图像x。
首先, 本论文采用一阶差分算子来获得图像四个方向上的一阶梯度:
式 ( 2) 中, i = 1, 2, …, p, j = 1, 2, …, q。xv、xh、x135和x45分别表示图像的水平方向、垂直方向、135°方向和45°方向四个方向。由公式 ( 2) 可知, 本论文几乎利用了图像所有方向上的一阶梯度, 与只利用水平方向和垂直方向两个方向的梯度相比, 显然能够获得更加丰富的梯度信息。
定义四个大小为pq × pq的矩阵: Dh、Dv、D135和D45, 因此图像四个方向上的一阶梯度可用矩阵的形式表示为:
众所周知, 图像的边缘是稀疏且连续的, 而且, 清晰图像的边缘比模糊图像的边缘更加的稀疏。利用图像边缘的这一稀疏特性, 再结合公式 ( 2) , 本论文提出了一种权重的稀疏性正则化约束图像复原模型:
式 ( 4) 中, ‖·‖22表示L2范数的平方, λ 为正则化参数。因此最小化L1范数得到的解是最稀疏的［11］, 因此, φ ( v) 定义为:
式 ( 5) 中, ‖·‖L1表示L1范数, v被定义为:
W表示图像边缘的稀疏性权重:
式 ( 7) 中, ε 为一个足够小的正数, 避免分母为零。
由公式 ( 4) 和公式 ( 5) 可知, 提出的方法利用了图像四个方向上的梯度信息, 较传统的全变差模型相比, 能够更加充分地利用图像中所有方向上的梯度信息。另一方面, 权重的引入则能够对图像边缘的稀疏性进行更加准确的表示。
2布雷格曼迭代策略
由公式 ( 5) 可知, 提出的复原模型引入了L1范数的最优化问题, 为了能够有效地求解这类模型, 本文运用了布雷格曼迭代的最优化策略［12］。
基于分裂的准则［12—15］, 引入一个辅助变量b1, 从而产生一个约束项: b1= Wv, 那么提出的模型便转化成了一种约束的最优化问题:
约束条件: b1= Wv。
再利用二次性的惩罚函数将约束项合并到公式 ( 4) 中, 于是约束的最优化问题 ( 8) 又转化成了一个非约束的最优化问题:
式 ( 9) 中, b1为辅助变量, t1为惩罚变量。
运用一种交互式的方法来迭代地求解辅助变量b1, 惩罚变量t1和复原图像x 。
2. 1固定b1和t1, 求解x
由公式 ( 5) 、式 ( 6) 和式 ( 7) 可知, Wv是对图像x进行一种矩阵的运算, 因此, 我们用一种矩阵运算D来代替Wv, 即Dx = Wv, 代入公式 ( 10) 中得到:
利用快速的傅里叶变换得到:
2. 2固定x, 求解b1和t1
利用二维的shrinkage函数来对变量b1和t1进行更新:
在提出的方法中, 初始化变量: x0= y, 同时辅助变量和惩罚变量的初始值均设为零, 即: b11= t11= 0。
3实验及分析
为了验证提出方法的性能, 在实验及分析部分, 本论文分别在256 × 256大小的“Lena”, “Camera- man”和“Barbara”三幅经典灰度级测试库图像上进行了实验, 并将提出的方法与近几年提出的一些具有代表性的图像复原方法进行了比较［6，10，14］. 本论文采用了增加的信噪比[ISNR = 10 lg ( ‖y - x‖2/ ‖x' - x‖2) ]来衡量实验中涉及到的复原方法的性能, 其中y为观察到的模糊噪声图像, x为原始的清晰图像, x'为复原图像。提出的方法是在windows7的操作系统上运行的, 实验所用的计算机具有4 G内存和2. 3 GHz的双核Intel处理器。在所有的实验中, 迭代的最大次数均被设为100次, 同时, 所有图像的像素值都被归一到0 ~1之间。
3. 1实验参数的配置
由公式 ( 9) 可知, 提出的方法一共引入了2个参数。通过实验我们发现: 参数 λ 决定了复原图像中边缘细节的锐化程度, 如果 λ 过小, 复原图像会变得过于锐化, 同时还会造成严重的噪声增强干扰; 如果 λ 过大, 则会出现严重的过度平滑现象, 平滑掉复原图像的边缘细节。经过反复的实验, 在本论文的实验中, λ 被设置为: λ ∈[0. 000 5, 0. 005]。
由公式 ( 12) 和式 ( 14) 可知, 参数 γ 越大, 每次迭代对变量b1的跟新就越少, 从而导致复原图像x的跟新速度就越慢, 算法迭代次数增加。因此, 为了在不影响算法性能的条件下尽可能地提高算法的效率, 通过反复的实验得到, 在本论文的所有实验中, 我们对参数 γ 进行如下的设置: γ = 10λ。
3. 2提出方法的性能比较实验
接下来, 本论文提出的方法将分别与文献[6]、 文献[10]和文献[14]的图像复原方法进行了比较。 为了保证实验的公平性, 各种方法中所涉及到的参数都是严格按照相应文献中作者的要求进行设定的。图1所示为实验所用的“Lena”, “Cameraman” 和“Barbara”三幅灰度级测试库图像; 图2所示为实验所用的三种不同的模糊退化函数。
图3、图4和图5从主观的视觉效果上对提出的方法与文献[6], 文献[10]和文献[14]的方法进行了比较。图3 ( b) 运用了均匀模糊 ( UB) 对Cam- eraman标准测试图像进行模糊; 图4 ( b) 运用了高斯模糊 ( GB) 对Lena标准测试图像进行模糊; 5 ( b) 运用了任意的运动模糊 ( RMB) 对Barbara标准测试图像进行模糊。所有的模糊图像都同时加入了BSNR = 45 dB的高斯噪声。
如图3 ~ 图5所示, 本论文提出的方法与文献[6], 文献[10]和文献[14]的方法相比, 均能得到更加清晰的复原图像, 能够复原出更加丰富的图像边缘和细节, 同时还能够有效地避免文献[10]方法中出现的假边缘瑕疵[如图4 ( c) 和图4 ( d) 所示]。
接下来, 增加的信噪比指标 ( ISNR) 将客观地量化以上四种图像复原方法的性能。表1列出了在Lena, Cameraman和Barbara灰度级标准测试图像上, 所有实验中所得到的ISNR值。可以很明显地看到, 本论文提出的方法明显优于其余的三种图像复原方法, 在所有的试验中均获得了最高的ISNR值。较近几年提出的一些较好的图像复原方法相比, 提出方法的信噪比增量增加了0. 3 ~2. 5 dB。
4结论
本论文针对现有方法存在的缺陷, 提出了一种权重的稀疏性正则化约束图像复原方法。提出的方法充分利用了图像四个方向上的梯度信息, 计算出了图像边缘四个方向上各自的权重, 提出了一种多方向的权重的全变差模型对复原图像进行权重的稀疏性正则化约。运用了一种布雷格曼迭代的最优化策略来得到清晰的复原图像。大量的实验结果表明, 提出的方法不仅能够复原出更加丰富的图像边缘和细节, 同时还能够有效地消除假边缘瑕疵, 得到更高质量的复原图像。较近几年提出的一些具有代表性的图像复原方法相比, 不仅主观的视觉效果得到了较为明显的改进, 而且客观的信噪比增量也增加了0. 3 ~2. 5 dB。将提出的方法运用到模糊视频序列的复原和图像的超分辨率重建中是接下来研究工作的重点。
离散正则化 篇3
同步电机是电力系统中的重要部件,其运行行为影响到电力系统的各个方面,而掌握精确的电机参数,对准确分析和计算其动态行为有重要的意义。在实际工作过程中,电机的实际参数值并不是一成不变的,而是随着环境和工况的不断变化在一定范围内变化,如温度变化引起的集肤效应,会影响电机定、转子的电阻值,磁场饱和程度不同也会影响电感参数等。 因此同步电机参数的辨识一直是电力系统研究的重要内容[1]。
传统的参数辨识算法中,最小二乘法是比较常用的算法,具有算法简单、易于理解、易于实现等优点,因此被广泛应用[2,3,4]。 但最小二乘法存在一定的局限性,没有考虑到系统的病态性问题。 所谓病态性问题就是系统数据微小的变化引起解的巨大变化[5],当病态性严重时,算法会存在收敛性和多值性的问题,结果将偏离真实值。
由于同步电机也是高维非线性系统,其病态性是参数辨识过程中无法回避的问题。 文献[6-8]都提到了同步电机系统的病态性,并采用子集选择法来克服系统的病态,但子集选择法有它的局限性,即需要一些先验知识来帮助确定哪些参数是固定的。 本文将参数辨识看成是一种非线性反问题,反问题具有不适定性,也就是病态性问题,其求解过程就是解决病态性问题的一个过程[9]。 在反演问题理论中,正则化是解决病态问题的基本思路,本文将经典的Tikhonov正则化方法引入到同步电机的参数辨识中,并通过在仿真中设置多个场景,证明了该方法能克服系统的病态性并有效地进行参数辨识。
1 病态性分析及其度量
一般而言,对于一个系统的模型,如果原始数据的微小变化引起解的巨大变化,则称该模型为病态的,反之则称为良态的。 病态与良态,是模型本身固有的属性,它表征了模型抗干扰性的强弱,即稳定性的好坏[5]。
通常用条件数K来度量病态性的严重程度。 统计应用经验表明:若0<K<100,则认为没有病态;若100 < K< 1 000,则认为存在中等程度的病态;若K >1 000,则认为存在严重病态[5]。
引理1 设, 若对Cn ×n上的某一矩阵范数‖·‖有,则非齐次线性方程组Ax=y与(A+δA)(x+δx)= y +δ y的解满足:
其中,‖·‖v为Cn上与矩阵范数‖·‖相容的向量范数,证明过程见文献[10];δ 为表征输出误差的参数。
由式(1)可以知道,数据的误差对逆矩阵和求解线性方程组解的影响与‖A‖‖A-1‖的大小有关,当‖A‖‖A-1‖较大时,近似逆矩阵或线性方程组的相对误差可能较大,因此‖A‖‖A-1‖可作为影响求解线性方程组解的大小的一种度量。
定义1 A,则称
为矩阵A的条件数。 一般地,如果系数矩阵A的条件数大就称A对于求逆或求解线性方程组是病态的,否则称为良态。
2 反问题与Tikhonov正则化方法
2.1 反问题与参数辨识
反问题是相对于正问题而言的,一个先前被研究的相对充分或完备的问题称为正问题,而与此相对应的另一个问题称为反问题。 从实际应用中来看,可以概括地说,有2 种动机驱动着反问题的研究:想了解物理过程过去的状态或辨识参数;想了解如何通过干预当前的状态或调整某些参数去影响或控制该系统,以使其在未来到达人们所期望的状态[9]。
图1描述了反问题的基本原理,而参数辨识是指在输入和输出数据的基础上,从给定系统的数学模型中确立系统模型参数,因此参数辨识实际上就是一种典型的反问题。反问题求解面临2个根本困难:
a. 用于反问题求解的原始数据可能不属于该问题的精确解所对应的数据集合,因而,在经典意义下的近似解可能不存在;
b. 近似解的不稳定性,即原始数据小的观测误差(这个在工程中是不可避免的)会导致近似解与真解的严重偏离。
这是反问题求解中要面对的2 个难点和关键所在,即所谓的反问题的不适定性。 其中,a为反问题解的存在性问题,对于参数辨识而言,即是参数的可辨识问题[11];b是关于解的唯一稳定问题,实际上也是病态性问题。 反问题求解主要是解决病态性问题。而在反演理论中,正则化方法是解决病态问题的基本思路。
2.2 Tikhonov正则化方法
假设反问题可以用一个抽象的算子方程(3)来描述,其中x代表系统的未知量,y代表系统的输出,A为系统算子。 反问题为:已知A和y来求未知量x。当A为线性算子,称其为线性反问题,否则为非线性的反问题:
求解反问题不适定性的普遍方法是:用一组与原不适定问题相“邻近”的适定问题的解去逼近原问题的解,这种方法称为正则化方法。 如何建立有效的正则化方法是反问题领域中不适定问题研究的重要内容。 解决不适定性的典型的方法是变分正则化方法,又称为Tikhonov正则化方法[9,12,13]。
通常测量值都是存在误差的,当y=yδ(yδ为输出测量值)时,问题式(3)的准确解为x=A+yδ;按照广义逆的定义来求A+,在数值上是不稳定的,换言之,求解极小化问题式(4)在Hadamard意义下是不适定的:
按照正则化思想,可以用一系列与问题式(3)相邻近的适定问题来近似,例如用下述带有参数 α(α≥0)的极小化问题来近似:
称Mα[x,y,A]为Tikhonov泛函,α≥ 0 为正则参数,易见式(5)的欧拉方程为:
因此式(3)的极小解xα为正则解:
对于任何 α≥0 而言,其解存在、唯一,且连续依赖于A、y和 α。 则余下的工作是如何选取合适的正则参数 α 的问题。 总体而言,正则参数 α 的选取要兼顾近似解的数值稳定性和与原问题的好的逼近程度这2 个要求。
这里正则参数 α 主要采用MOROZOV偏差原理来获得,假设,具体步骤如下[9],其中 αn为第n次迭代中的正则参数值,xαn,δ为第n次迭代中的x值。
a.给定初始正则参数α0≥0,令n=0。
b.解方程
c.对步骤b中的方程求导得到方程:
求解该方程得。
d . 分别由表达式,计算出F (αn)和F′(αn)。
e. 令 αn+1= αn- F(αn) / F′ (αn),若小于某指定精度,则计算终止,否则进入步骤f。
f. 令n = n + 1,转步骤b。
3 同步电机模型及病态性分析
同步发电机参数的计算依赖于数学模型的建立,模型不同参数也有所不同,本文选用同步发电机在dq旋转坐标系下的稳态方程(8)作为数学模型,同时忽略饱和、磁滞和涡流的影响,并且忽略阻尼绕组[14,15]。
其中,id、iq和ud、uq分别为定子绕组d、q轴的电流和电压;if、uf分别为励磁绕组的电流和电压;Rs、Rf分别为定子绕组和励磁绕组电阻;Ld、Lq分别为定子绕组d、q轴上的自感;Lf为励磁绕组自感;Lmd、Lmq分别为定子绕组d、q轴与励磁绕组间的互感。
由式(8)可知,本文主要识别的参数为Rs、Rf、Lq、Ld、Lf、Lmd、Lmq,因此式(8)可以写成如下形式:
其中,A为利用观察值建立的矩阵;输出y = [uduquf]T;参数矩阵x =[LdLqLfLmdLmqRs]T,所以本文中的参数识别的反问题即为:已知A和y,求x。
同步发电机的数学模型的病态性主要表现为矩阵A的病态性,可以通过计算法矩阵ATA的条件数来度量系统病态的严重程度,条件数的计算可以采用式(2)。
本文将通过下一节的同步发电机实例的测试数据来计算法矩阵的条件数。 由于矩阵A需要用到电流的微分量,所以必须进行离散化处理:i′ = (i(k) -i(k - 1)) / Ts,其中Ts为采样周期。 对同步电机的运行电流进行采样,就能确定矩阵A,然后通过式(2)就能计算得到条件数K。 当Ts= 1 × 10-4s时,计算得K = 4.687 × 105。 由此可知,系统的病态性严重。
4 参数辨识的仿真
4.1 仿真模型
为了获得同步发电机的测量数据,采用MATALB的Simulink平台[16]搭建了同步发电机模型,其中电机采用Simulink自带的同步电机模型,反问题的Tikhonov正则化方法采用S-function来编写,并作为一个模块嵌入到同步发电机系统仿真环境中,主要实现对发电机参数的识别。 电机的参数为某航空独立交流电源中主发电机的实际参数,如表1所示。
4.2 仿真分析
为了真实地模拟观测数据,并验证Tikhonov正则化方法求解病态问题的能力以及对参数辨识的有效性,分别对测量数据加入10%、20%、30% 的Gauss白噪声,即:
由式(9)可知,反问题的输入为电机模型输出的iq、id、if、uq、ud、uf,如图2 所示(图2 所示为t = 0.1 s时机械功率Pm阶跃增大时的发电机电流、电压波形),然后通过反问题的正则化方法,计算出同步电机的参数Rs、Rf、Lq、Ld、Lf、Lmd、Lmq,其中Ld= L1s+Lmd,Lq= L1s+Lmq,Lf=L1fd+Lmd。
本文采用MOROZOV偏差原理来求解正则化参数 α。 选取正则参数 α 必须非常小心,如果 α 太大,则新得到的问题对原问题的逼近程度太差;相反如果α 太小,则问题的不适定性并没有克服,数值计算仍然很不稳定。
a. 加入10 % 白噪声:表2 是在测量数据中加入10 % 白噪声后的辨识结果,经过迭代最终得到正则参数 α=1.36×10-3,迭代次数为8。 由表中数据可以看出辨识结果较好。
b. 加入20 % 白噪声:表3 是在测量数据中加入20 % 白噪声的辨识结果,经过迭代最终得到正则参数 α=8.65 × 10-3,迭代次数为16。 由表中数据可以看出辨识结果虽然不太精确,但可以接受。
c. 加入30 % 白噪声:表4 是在测量数据中加入30 % 白噪声的辨识结果,经过迭代最终得到正则参数 α =3.64 × 10-2,迭代次数为20。 30 % 白噪声代表比较严重的工况,由表中数据可以看出辨识结果开始偏离真实值。
4.3 与传统最小二乘法的比较
如何建立有效的正则化方法是反问题领域中病态问题研究的重要内容。 通常的正则化方法有Tikhonov正则化方法、信赖域法、正则化的内积法等,Levenberg-Marquardt[17]也是一种特殊的正则化方法,可以看作是对非线性问题作先线性化后正则化的过程,文献[9]对它的正则化进行了证明。 本文选择的是经典的Tikhonov正则化方法。 传统的参数辨识方法中,使用最多的是最小二乘辨识法以及一些改进的最小二乘法,如增广二乘法、广义最小二乘法、加权最小二乘法,与它们相比,本文方法的本质区别在对系统病态性能力的克服上,这里通过Tikhonov正则化方法与最小二乘辨识法的比较来进行证明。 对于系统Y=AX,X的最小二乘解是x=(ATA)-1ATy,而正则解为xα=(ATA + α I)-1ATye。 系统的病态表现为矩阵A的病态,即法矩阵ATA的病态性。 与最小二乘法相比,正则化方法增加了 α I一项,这一项的引入使法方程的病态性得到改善,因而能得到好的估计值。
表5 显示了对测量数据加入10 %、20 %、30 %的白噪声,并采用最小二乘法进行辨识的同步电机参数值。 可以看出,在10 % 的白噪声污染下,最小二乘法的结果已经开始偏离,但勉强可以接受,而在强噪声(30% 的白噪声污染)的工况下,最小二乘法的辨识结果已经完全偏离真实值。 由此可知,Tikhonov正则化方法克服病态性的能力优于最小二乘法。
5 结论
离散正则化 篇4
用迭代方法处理各种反问题已有悠久的历史。但是研究表明, 使用迭代方法求解反问题, 有时会出现所谓的“半收敛”现象, 即在迭代的早期阶段, 近似解可稳定地得到改进, 展现出“自正则化”效应, 但当迭代次数超过某个阈值后便会趋向于发散。因而, 使用迭代法求解的关键是要寻找一个恰当的终止原则, 在迭代次数和原始数据误差水平之间找到平衡值。研究表明, 迭代指数, 即迭代步数正好起到正则化参数的作用, 而这个终止准则对应着正则化参数的某种选择方法。并且使用迭代方法求解还有很多优点, 因此, 在正则化问题求解中通常选用迭代的方法, 常用的迭代方法有:Landweber迭代法、Van Cittert迭代方法、最速下降方法和迭代Tikhonov正则化的求解方法, 以及正则化方法的快速数值实现。
2 基于解空间分解的GMRES算法及图像复原应用
2.1 正则化模型与图像复原
设F和U分别表示度量空间, 度量为ργ和ρμ, 算子A:F到U映F到U, 则该问题变为线性反问题 (当A为线性算子时) , 或非线性反问题 (当A为非线性算子时) 。“不适定性” (病态性) 是所有反问题所具有的一个共同的特性。一般情况下, 不适定性是反问题本身的固有特征:如果问题的先验信息是未知的, 那么就无法得到理想的结果。因此, 我们应该尽可能多的收集先验信息, 最大限度的复原原问题。通常, 人们将求解反问题 (不适定问题) 的理论和方法称为正则化方法。对于图像处理问题, 由于涉及到大规模的方程组求解, 法方程的维数太大, 此时再应用代数方法求解就会遇到一些难以实现的技术问题, 而选用正则化方法不但可以克服上述缺点, 还具有某些优点, 当问题从无穷维度变到有限维度时, 迭代求解不会影响系数结构, 而且能够起到节约运算空间的效果。这些优势在大规模计算中非常有利。
对于图像恢复的病态性问题, 利用正则化思想进行图像复原时, 需要利用先验信息, 构造某种约束条件, 使用数理统计方法, 将图像复原这一不适定问题转变成适定问题, 进而使得近似解满足适定性的三项约束, 这也是正则化方法的优势所在。
2.2 解空间分解的广义极小残量算法
在对线性方程组Ax=b, A紧算子, 进行求解时, 为了尽可能减少存储空间和计算开销, Krylov子空间迭代法是求行之有效的方法。当系数矩阵A对称正定, 共扼梯度法 (CG) 或预共辘梯度法 (PCG) 可快速准确求解该方程组的近似解;当A对称但不正定时, 极小残量法或预极小残量法则能有效求解方程组。对于一般的非对称矩阵, 常采用广义极小残量法、共扼梯度法来求解。GMRES算法利用Arnold过程产生Krylov子空间Kj的正交基, Arnold过程中每次迭代运算, 都要调用所有前面的迭代所产生的正交基来生成下一个正交基。
2.3 光学图像复原结果
对于方程Ax=b, 利用基于解空间分解的加速GMRES算法迭代求解。计算步骤如下:
Step1.置初始值x0=0, 并令δ=10-8;
Step2.用解空间分解的加速GMRES算法迭代求解式Ax=b, 在第j步的值为xj;
Step3.若终止迭代;否则置j=j+1, 继续进行Step2。在迭代运算中, 正则化参数αj=15, 随着迭代的进行自动更新。
图像复原实验中处理的是256*256尺寸的0-255灰度级的liftingbody图像。用改进信噪比来衡量算法的复原性能。从复原之后的对比效果看, 共扼梯度法 (CG) 并不能有效的抑制模糊退化, 复原结果仍然比较模糊, 图像边缘有振铃波纹出现。解空间分解的加速GMRES算法复原结果的边界纹路比较清晰, 很好的显示出原图像边缘细节部分, 与此同时, 振铃波纹因为加窗处理得到有效抑制, 整体视觉效果很好。
3 线性代数方法与图像复原应用
在涉及到复杂矩阵和向量的离散图像复原模型中, 可以从线性代数方法中得到一种效率较高的求解方法, 常用的方法是约束最小二乘法。对约束最小二乘法进行改进, 根据先验信息, 把正则化思想和约束最小二乘法等有机结合在一起, 并将其运用到离散图像复原中, 得到的约束最小二乘的空域迭代法可以出色的抑制噪声, 而且在噪声很强是也可以得到很好地复原结果。
将正则化思想与约束最小二乘法相结合, 继而复原退化图像。通过对噪声能量的限制来使用正则化理论, 运用空域迭代时很好的抑制了噪声放大现象, 同时克服了病态性, 而且计算速度得到了提升。实验数据表明, 本方法更适合复原污染程度较大的图像, 但不适合复原模糊程度较大的图像。
4 总结与展望
图像复原近年来受到了越来越广泛的关注, 正则化方法理论的发展也越来越得到完善, 许多学者从模型上、理论上、应用上分别展开了对于正则化的图像复原的深入研究, 本文的研究虽然力求有较强的实用性, 但是由于受到多方面的限制, 在理论和工程应用等方面仍存在很多的待丰富和改进之处, 需要在以后的工作中继续深入研究。
首先, 要加深对观测图像的先验信息的挖掘, 因为如果能够有效利用先验信息, 就能极大的改善估计精度以及问题的病态性。要充分的利用各种先验知识, 构造更加精确地目标泛函, 设计出更加优良的算法, 同时要充分分析现有正则化参数的选择方法, 结合各个方法的优缺点, 构造出更加高效的正则化算子;要注意噪声扩大与图像复原的平衡, 充分利用成像时的分段平滑性质, 去除图像边缘模糊和振铃现象。其次, 要注重正则化图像复原方法与其他图像复原方法的有机结合。现阶段, 神经网络、小波分析和遗传算法等新式的算法在图像复原方面取得了极大的进展, 如果能够将这些理论结合在一起, 形成优势互补, 一定能得到性能更好的图像处理算法。
参考文献
[1]肖庭延, 于慎根, 王延飞.反问题的数值解法[M].北京:科学出版社, 2003.
[2]邹谋炎.反卷积和信号复原[M].北京:国防工业出版社, 2001.
[3]苏超伟.偏微分方程逆问题的数值解法及其应用[M].西安:西北工业大学出版社, 1995.
离散正则化 篇5
关键词：稀疏表示,缺失数据填补,分类属性
1 引言
稀疏表示理论[1]是机器学习领域近几年出现的新方法,其应用最小化l1范数[2]的优化方法获得基于过完备字典特征表示的稀疏向量,是获取、表示数据的有效工具。在现有的分类研究应用中,稀疏表示获得了比传统方法更好的分类性能,已成功应用于人脸识别、语音识别等信号和图像识别领域[3]。另一方面,传统的数据挖掘理论研究需要数据完整而确定,但在实际应用中,由于数据测量误差、获取限制、存储介质故障等原因,人们所能收集到的数据往往存在缺失现象。应对数据缺失现象的常规做法是寻求合适的算法进行缺失数据填补。相较于能够直接进行数值计算的数值型属性,在处理分类属性时,由于其不具备直接进行数值计算的原理,需要进行相应处理后方可进行填补[4]。
基于上述情况,Shao,et.al提出了基于K最近邻局部约束稀疏表示的分类属性缺失数据填补方法CSRimpute(Categorical Sparse Representation imputation)[5]。该方法针对分类属性缺失数据的特点,利用完整数据设计字典,在保留局部结构特征的同时改善分类属性缺失数据的填补效果。
2 CSRimpute算法介绍
CSRimpute算法是在局部约束稀疏表示的基础上,结合K最近邻算法设计字典,力图解决缺失数据的填补问题。该算法可以适用于包含一个缺失值或被推广到包含多个缺失值的数据对象上。为了方便说明,需要定义一些概念如下:
X=[x1,x2,…,xi,…,xn]∈Cm×n表示一个包含n个数据对象的分类属性数据集。
列向量xi∈Cm×1表示第个数据对象:
第i个数据对象在第j个属性上的缺失值成为缺失属性值,记做:X(j,i)=xi(j)=xij=“*”。分类数据集共有m个属性行,每个属性行分别有cj种取值,且c1+c2+…+cm=M。
该算法的具体过程如下:
输入:含有缺失数据对象的数据矩阵X=[x1,x2,…,xi,…,xn]∈Cm×n;正则化约束参数λ1,λ2>0;字典包含原子数据对象数量k;
输出,填补后的完整数据集X;
过程:
(1)将原始数据集X转化为二进制矩阵A。在xi的第j个分类属性行所对应的cj行中,仅在代表其取值的属性行取值为1,其他取值为0;若属性缺失,则cj行取值均为0;
(2)将A划分为A=[ACAM]两部分,其中完整数据集AC=[a1,a2,…,ac,…,anc]∈Cm×nc,缺失数据集AM=[anc+1,anc+2,…,am,…,anc+nm]∈Cm×nm,假设A中前nc个数据对象都是完整的;
(3)应用K最近邻作为字典构造方法,针对AM中的每个缺失数据对象am分别构造字典AN(m)=[aN(1),…,aN(k),…,aN(K)]∈CM×K,重复步骤4至步骤7;
(4)将am和AN(m)在所有am非缺失属性上进行投影得到am*和A*N(m),即去除am中的所有的缺失属性并在AN(m)中移除相应的属性;
(5)计算am*和A*N(m)中每个数据对象的欧几里得距离,根据公式,
得到距离向量d;
(6)针对缺失数据集AM中的每个缺失数据对象am求解公式,
对应的局部约束稀疏表示优化问题,得到稀疏表示系数向量αm*∈RK。其中稀疏表示系数向量为am*在A*N(m)的重构表示;d为新投影后的缺失数据对象xl*与投影后的字典B*中每个原子数据对象的距离构成的向量。
(7)针对xi包含的每一个缺失属性值xij,对于它在重构数据中对应的cj个重构值,选取其中最大重构值在原始分类数据集X中所对应的属性值进行填补。其中;
(8)算法结束,输出填补后的数据矩阵X。
3 实验分析
本实验从UCI机器学习数据库中选择了4个经典的分类属性数据集(Soybean,ZOO,SPECT Heart,Chess)。为了将原始数据值和填补估计值进行对比,针对每个数据集,首先删除其中包含缺失属性值的数据对象,得到完整数据集。然后,随机选取数据集中的部分数据对象构成缺失对象数据集,对于每个被选取的数据对象,从中随机选择一定数量的属性值,人为地将这些属性值设定为缺失属性值。
设计本试验的目的是测试正则化参数l1和l2的敏感性。为了便于操作,将缺失属性比率和缺失数据集的比率均设定为20%。对于设定的正则化参数,分别选取它们的2n倍计算缺失填补正确率,结果如表1、表5所示。
从前3张表中可以看出,在4个数据集中,无论如何变化,算法的填补正确率相对都比较平稳,没有出现大幅度增加或减少的情况。两个正则化参数λ1和λ2的变化使得CSRimpute的填补效果略有浮动,但浮动的范围较小,并未过度偏离最优结果。最后的统计表显示,Soybean和Chess数据集的极差和标准差较小,而ZOO和SPECT Heart数据集的极差和标准差相对较大。这在一定程度上说明正则化参数λ1和λ2对填补效果的影响在较大数据集的情况下反而较小,这是因为在较大数据集中能够更容易找到与目标数据对象更相似的数据,从而能够得到更加理想的稀疏表示。
4 结论
本文针对CSRimpute算法中的正则化参数该如何选择的角度出发,通过4个数据集验证分析了λ1和λ2对算法填补效果的影响。实验结果表明,缺失数据的填补效果随着正则化参数在最优值附近较为平稳的变化,即CSRimpute算法受正则化参数λ1和λ2变化的影响并不明显,在实际应用中能够比较容易确定。
参考文献
[1]Wright J,Yang A Y,Ganesh A,et al.Robust Face Recognition vi Sparse Representation[J].Pattern Analysis and Machine Intelligence,IEEE Transactions on,2009,31(2):210-227.
[2]Candès E J,Romberg J,Tao T.Robust Uncertainty Principles:Exact Signal Reconstruction from Highly Incomplete Frequency Information[J].Information Theory,IEEE Transactions on,2006,52(2):489-509.
[3]Duan C H,Chiang C K,Lai S H.Face Verification with Local Sparse Representation[J].Signal Processing Letters,IEEE,2013,20(2):177-180.
[4]Shekhar S,Patel V M,Nasrabadi N M,et al.Joint Sparse Representation for Robust Multimodal Biometrics Recognition[J].Pattern Analysis and Machine Intelligence,IEEE Transactions on,2014,36(1):113-126.
离散正则化 篇6
翁爱华[13]曾采用著名的OCCAM[14]反演法直接对感应电动势进行反演,OCCAM反演法采用线性搜索来自适应优求取每一反演迭代步中的最佳正则化因子,使得在每一迭代步中数据目标函数极小(充分拟合),且在此条件下的模型粗糙度函数极小(模型极光滑)。这一措施使得OCCAM具有极高的稳定性和模型分辨率,但由于每次迭代要额外多次求解反演方程和正演计算,导致运算量过大、耗时较长。本文采用陈小斌[15]在大地电磁测深法中使用的自适应正则化反演法(adaptive regularized inversion algorithm,ARIA)直接对感应电动势进行反演拟合。相较于OCCAM反演,ARIA法主要优点是采用简单有效的正则化因子求取算法,其收敛时间更快,拟合效果基本相同。
1反演方法基本原理
按照正则化反演思路,即要寻求一个模型,它既可以尽量地拟合观测数据,同时又满足模型目标函数所限定的模型性质。因此正则化反演目标函数可以表示为
式(1)中Φd(m)为观测数据目标函数,Φm(m)是模型先验约束条件目标函数,m为模型向量,λ为正则化因子。Φd(m)由式(2)给出。
式(2)中Δd为观测数据与理论响应差,Wd为数据加权矩阵,Wd=diag{1/σ1,1/σ2,…,1/σj,…,1/σm},σj为数据标准差。
构建模型约束目标函数常用的三种是模型参数的平方和最小、模型参数导数的平方和最小以及模型参数二阶导数的平方和最小;又分别对应简称为最小模型约束、最平缓模型约束和最光滑模型约束。由于本文采用的是层状介质模型,上述三种约束离散化表示为
式(3)中R表示粗糙度,表示为
正则化因子λ为先验模型目标函数与数据目标函数之间的加权系数,其作用是平衡两个目标函数,决定了反演的拟合对象。λ→∞表示完全以先验信息作为反演目标,数据拟合程度差;λ=0则表示完全以拟合数据作为反演目标,模型约束差。因此如何选取正则化因子对反演结果影响很大,一方面要让数据充分拟合,同时也要放得到的模型满足先验信息。参照陈小斌提出的自适应正则化算法,正则化因子λ有两种自适应调整方案:
式中k表示第k次反演迭代。由于电磁法的一维反演常以均匀半空间作为初始模型,第一次迭代m计算值为0,这会造成MD方案的正则化因子计算无意义,因此采用CMD方案。
将式(1)对模型求导,根据极小化原则,可得到反演方程
令A=JkTWdTWdJK+λRTR,b=JkTWdTWdΔdk-λRTRmk,式(4)可简化为AΔm=b。式中,Δdk为第k次迭代中观测数据与当前模型的理论响应向量差,JK为当前模型的灵敏度矩阵,这里一维反演的灵敏度矩阵采用差分法求解。在A非奇异的情况下,模型修正向量可表示为Δm=A-1b,从而可得到下一次迭代的初始模型m=mk+Δm。
反演的终止条件为相对拟合差(rms)小于一个给定的值。
2数值结果
讨论地电断面被剖分为30层的等差网格,反演深度由线圈大小、采样时间和断面平均电阻率确定,初始模型为100Ω·m的均匀半空间,选取最平缓模型约束,迭代要求至少5%拟合差或者最多15次的迭代。
2.1理论数据
为了便于比较分析,计算了一个类似积水采空区的H型地电模型[图1(b)],从上往下三层电阻率分别为400、20和400Ω·m,积水层埋深120 m、厚度50 m,采用自适应正则化反演经过10次迭代得到了反演结果。图1(a)中的大图是由H模型计算的感应电动势(实线)和反演结果模型响应(虚线)对比曲线,小图是在反演迭代过程中的观测数据目标函数即拟合差随正则化因子的变化情况。
可以看出反演稳定收敛,正则化因子在前5次迭代中变化较大,之后趋于稳定,拟合误差快速减小;感应电动势曲线拟合效果好(第10次迭代的收敛误差为8.63×10-4),只是在晚期几毫秒以后有微小的偏离。从图1(b)中的反演模型对比中可以看到低阻层的反应很好,反演结果是相当可靠的。
2.2实测数据
实测数据来源于一次对某煤矿区的瞬变电磁测深,已知地质资料显示该区煤层、采空积水区与围岩之间存在较为明显电性差异,表现为相对低阻异常。采用中心回线装置,选取工区中某测线作分析,该测线有点距10 m的测点共400个。图2(a)先是根据晚期视电阻率公式[7],从感应电动势计算得到的晚期视电阻,再对晚期视电阻率采用目前解释中常用的“烟圈反演”得到的断面图,图2(b)是对各测点采用自适应正则化反演法直接对感应电动势进行反演得到的断面图,图2中的红线是根据地震勘探资料推断划分出的煤层(红实线)或采空区积水层(红虚线)。可以看到,同地震勘探资料划分出的异常区的深度和厚度相比较,采用自适应正则化反演要比“烟圈反演”结果对应得更加准确。
选取测线中位置在1 200和4 000的两个测点分析,图3是两个测点的反演拟合情况。可以看出两个测点的感应电动势曲线拟合都很好(1 200位置测点第10次拟合差为4.65×10-2,4 000位置测点第10次拟合差为4.96×10-2),正则化因子刚开始变化较大,两个目标函数也随之变化大,但在5次迭代后反演就趋于稳定。
3结论
离散正则化 篇7
80年代初,Johnson将微波技术中矩量法引入超声成像领域[1],矩量法不再局限于Born近似的限制,可以方便地引入被测物体的某些先验知识(这些先验知识是克服任何逆散射问题的不适定性和对噪声敏感性的关键),可以对大物体、高对比度、强散射的物体进行成像,因此适应范围更广[2]。
不适定性是图像重建逆问题求解中非常普遍的情况,本文采用Tikhonov和TSVD正则化方法解决其不适定性,两种正则化算法的成功应用依赖于正则化参数λ的合适选取。其中容差原理法确定正则化参数需要知道数据项中噪声大小,而广义交叉验证法操作较为繁琐且结果不稳定。本文基于L-曲线法,提出以解的范数和残差变化量的加权形式作为确定正则化参数的依据,在迭代过程根据问题不适定性程度,自适应地调整搜索范围。通过两种正则化技术和优化策略来实现在空间域内的图像重构。
1 UCT图像重建原理
当用超声波从某一方向作用于物体的被成像截面时,超声波与物体内部介质相互作用产生的全场满足非齐次亥姆霍兹方程[3],其解可用Lippmann-Schwinger积分方程表示:
undefined
式中:p(undefined)和pinundefined分别为全场和入射场;undefined为自由空间的格林函数;undefined为物体内部介质声学特性的未知函数。
在实际应用中,我们只能测到物体的散射场和入射场,将式(1)转化为散射场形式:
undefined
通过矩量法将上述两个卷积方程转化为代数方程组的形式:
P(t)=P(m)+COP(t). (3)
P(s)=DOP(t). (4)
方程(3)、(4)仅描述了物体在某一方向入射波作用下接受到的散射波,为获得足够多的信息,需要从K个方向发射超声波,由此得到K×M个方程,为使上述方程组可解,通常KM>N,即要求方程组为超定的。
这个超定方程组的求解具有较强的非线性。这里采用Born迭代方法[4]。具体描述如下:
①先假设全场等于入射场P(t)=P(in),由方程(4)求得未知函数的初始解O0,常称其为Born逆解;
②然后将Ok代入方程(3)求更接近实际的物体内部的全场Pundefined,k表示迭代次数;
③由第②步求得的全场Pundefined,代入方程(4),计算散射场Pundefined=DOPundefined,并求计算散射场与测量散射场的差ΔPundefined=Pundefined-P(s),若‖ΔPundefined‖小于事先给定的δ,则停止迭代,否则,转④;
④根据散射场的改变量ΔPundefined,由方程ΔPundefined=DΔOkPundefined求未知函数的改变量ΔOk,然后赋Ok+1=Ok+ΔOk作为未知函数的新值;
⑤由新的未知函数Ok+1(rn),返回步骤②,进一步求更近似的全场。
在迭代过程中,将遇到不适定性逆散射方程的求解问题,而对不适定性问题的正则化技术的选择是否得当,将直接影响迭代算法的稳定性。所以,正则化技术是图像重建问题的关键。
2 正则化方法及参数的选取
影响Born重建效果的一个关键因素是散射方程ΔPundefined=DΔOkPundefined的求解。先将其简化为Ax=b,最小二乘解min‖Ax-b‖2可表示为:undefined。其系数矩阵A具有两个特点:(a) 条件数较大;(b) 奇异值迅速下降。方程组是不适定的,利用Tikhonov和TSVD正则化方法并选择L曲线准则来确定最佳的正则化参数[5]。
算法1:用Tikhonov方法计算Ax=b的算法如下[6]。
①矩阵A的奇异值分解:=svd(A),其中S=diag(σ1,σ2,…σn)
②对σn≤λ≤σ1,
计算解的范数:undefined;
残差的范数:undefined。
整数λ称为正则化参数。
③作图(lg‖xλ‖2),lg‖Axλ-b‖2)),确定曲线拐角点,拐点对应的xλ是正则解xtik。
④采用Tikhonov法求解正则解undefined。
算法2:用TSVD方法计算Ax=b的算法如下[7]。
①矩阵A的奇异值分解:=svd(A)。
②对1≤k≤T,计算解的范数:
undefined,整数k称为截断参数;
残差的范数:undefined,其中T为奇异值σi的维数。
③作图(lg(‖xλ‖2),lg(‖Axλ-b‖2)),确定曲线拐角点,拐点对应的xk是正则解xtsvd。
④采用TSVD法求解正则解undefined。
研究表明,对离散不适定问题,lg(‖xk‖2),lg(‖Axk-b‖2)的关系曲线总是呈现L-曲线形状,该关系曲线有一个明显的拐角,最优正则化参数位于L-曲线拐点附近[8]。因此,可通过寻找L-曲线拐点作为最优正则化参数。
3 重建结果比较
实验装置我们采用圆环形结构,发射器(同时也是接收器)等间隔地放在围绕物体的圆环上,半径取20×λ,背景介质为水。采用F=200 kHz的超声波作为入射波,超声波在水中的波速约为C0=1 500 m/s,波长约为7.5 mm,对应的波数为k0=0.837 758 mm-1。按照矩量法,我们对包含物体的正方形区域进行均匀采样,采样间隔均为d=λ/10=0.75 mm,采样点的个数为35×35,划分后单元的总数为N=1 225,且在每个小单元上做内接圆,其半径为a=d/2=0.375 mm。换能器共有50个,可工作于发射、接收两种模式,每次仅有一个换能器发射,其他换能器接收来自物体的散射波。
如图1所示,实验选用的样品如下:对比度为20%的人工绘制的图形,轮廓信息较为丰富。对于已知像函数O(undefined),先由全场方程(3)求出ROI中的全场分布,然后通过散射场方程(4)求出物体外的散射场P(s)分布,再使用BI进行图像重建。两种方法的重建结果如图2和图3所示:
由相对误差可以判别仿真图像与目标原图像的接近程度。Tikhonov正则化方法的结果与目标原函数相比,其相对误差为:undefined,当采用TSVD正则化方法的结果与目标原函数相比,其相对误差为:undefined。相对误差越小,重建图像越接近目标原图像。
由图4可以看出采用Tikhonov方法时,相对误差在第4次相对误差最小,当超过第4次后,相对误差变大。图5可以看出采用TSVD方法迭代到第18次时最小误差恒定不变,因此在第18次中止迭代。对比两种方法,TSVD方法明显好于的Tikhonov方法重建的结果。
4 结论
本文运用Born迭代算法解决了方程的非线性问题,对于离散不适定问题,采用Tikhonov和TSVD两种正则化方法分别对方程进行了求解。经实验验证:TSVD正则化方法明显优于Tikhonov正则化方法。超声逆散射成像技术是一个比较复杂的系统工程,算法的稳定性及速度只是问题的一个方面。目前尚没有成熟的理论和技术可用于产品的自动检测,对该方面的一些关键技术进行研究具有一定的学术价值和广阔的应用前景。
参考文献
[1]Tracy M L,Johnson S A.Inverse Scattering Solutions bya Sinc Basis,Multiple Source,Moment Method-Part II:Numerical Evaluations[J].Ultrasonic Imaging,1983,5(4):376-392.
[2]陶进绪,张东文,叶寒生.频域法超声逆散射成像[J].信号处理,2009,25(2):169-173.
[3]刘超,汪元美.超声层析成像的理论与实现[D].浙江大学,生物医学工程,博士论文,2003.
[4]Chew W C,Wang Y M.Reconstruction of Two-Dimen-sional Permittivity Distribution Using the Distorted BornIterative Method[J].IEEE Transactions on Medical Ima-ging,1990,9(2):218-225.
[5]韩旭,刘杰,李伟杰.时域内多源动态载荷的一种计算反求技术[J].力学学报,2009,41(4):595-601.
[6]Hansen P C,O'Leary D P.The Use of the L-Curve in theRegularization of Discrete Ill-Posed Problems[J].SIAMJournal on Scientific Computing,1993,14(6):1487-1503.
[7]Hansen P C,Jensen TK.An Adaptive Pruning Algorithmfor the Discrete L-curve Criterion[J].Journal of Computa-tional and Applied Mathematics,2007,198(2):483-492.