范文一：基于迭代方法的软阈值估计小波去噪

第 31卷 第 1期 系统工程与电子技术

Vol. 31 No. 12009年 1月

Systems Engineering and Electronics Jan. 2009

文章编号 :1001-506X(2009) 01-0036-04

收稿日期 :2007-09-10; 修回日期 :2007-11-15。 基金项目 :国家自然科学基金资助课题 (60572135)

() 基于迭代方法的软阈值估计小波去噪

孙 蕾 , 谷德峰 , 罗建书

(国防科技大学理学院 , 湖南 长沙 410073)

摘 要 :基于小波变换的阈值去噪方法依赖于信号的长度和噪声的方差 。 在大多数情况下 , 噪声的方差是未 知的 , 需要对其进行估计 。 提出一种基于软阈值函数估计阈值大小的算法 。 该算法是一种迭代算法 , 证明了该算 法的收敛性 , 并分析了该算法运算量 。 数值实验表明该算法估计阈值的准确性和 M AD 方法相当 , 同时计算量远 远小于 MAD 方法 。

关键词 :小波 ; 去噪 ; 软阈值 ; 阈值估计 中图分类号 :T N 911. 72 文献标志码 :A

Wavelet denoising to estimate the soft -threshold based on iterative methods

SUN Lei, GU De -feng, LUO Jian -shu

(Coll. of S ciences, National Univ. of Def ense Technology , Changsha410073, China)

Abstract:Thresholding based on w avelet depends on the length of the s ignal and the level of the noise. In

m ost cases, the variance of the noise is unknow n, so it needs to be estimated. T his paper proposes an iterat ive m ethod to est im ate the threshold based on soft -thresholding. The convergence of t he iterat ive algorithm is proved, and an analysis of t he computat ional time is given. Num erical experiment show s that the algorithm is compet it ive to the M AD met hod, and the comput ational t ime is m uch less than M AD met hod.

Keywords:w avelet; denoise; soft -threshold; threshold est imation

0 引 言

在信号的获取和传输过程中 , 不可避免的会引入噪声 , 并对后续的信号处理将产生不利的影响。由于小波良好的 时频局部化特性 , 在信号去噪中 , 取得了广泛的应用。因为 一般噪声分解后的小波系数幅值都比较小 , 所以可以利用 小波系数阈值收缩法来区分噪声和信号 , 去除较小的小波 系数 , 得到去噪后重构图像的小波系数 , 再经过小波逆变换 得到重构图像。

1994年 , Donoho 和 Johnstone 首次提出了小波阈值去 噪方法 (VisuShrink 方法 ) 。此方法在信号属于 Besov 空间 时 , 在大量风险函数下获得近似理想的去噪风险 , 取得了优 于其他的线性估计的效果 [1]

。该方法在理论上非常理想 , 但在实际应用中 , 随着小波系数个数增加 , 相应的阈值随之 增大 , 有趋向将所有小波系数置零 , 产生 /过扼杀 0系数的现 象。 Donoho 和 Johnstone 随后提出 SU REShrink 方法 [2]

, 该准则是基于均方差准则的无偏估计 , 在实际应用中能获 得较为满意的效果 , 但该方法也存在过保留系数的特点。 此后针对 Donoho 和 Johnstone 的小波阈值估计方法 , Na -son [3]提出的交叉验证法 (cross validation, CV) , 以及 Weg -rich 和 Warchola

[4]

提 出 的广 义交 叉验 证 法 (generalized

cross validation, GCV) 都对 Donoho 的算法有所改进。但

是这些方法是通过对信号进行交叉处理和小波变换来确定 阈值 , 其精度虽有所提高 , 计算量却大大增加。李旭超和朱 善安 [5]指出尽管 Donoho 的 Visushrink 方法是利用信号的 小波细节的统计值作为阈值 , 其精度稍有不足 , 但由于计算 量小 , 因而很实用。关新平等

[6]

也提出一种新的模糊小波

阈值去噪方法 , 根据不同的模糊隶属度在 VisuShrink 阈值 与 SUREShrink 阈值之间选取去噪阈值 , 并用软阈值函数 消减 , 取得了良好的效果。所以本文针对 VisuShrink 方法 提出一种新的阈值估计方法 , 在运算速度和阈值估计精度 上都有极大的改进。

Donoho 和 Johnstone 提出的小波阈值去噪方法 (Vis -ushrink 方法 ) 阈值为 T =R 2ln N , 这个阈值是依赖于信号 长度 N 和噪声标准差 R 的。在实际应用中 , 由于 R 是未知 的 , 因此需要估计 R 的大小。 Donoho 给出了噪声信号方差 的一个估计方法 MAD (median absolut e deviation) , 该方法 能较准确的估计出被高斯白噪声污染的信号的去噪阈值。

第 1期 孙蕾等 :基于迭代方法的软阈值估计小波去噪 #37#

小波阈值去噪准则最常见的有硬阈值准则和软阈值准

则两 种。 Alexandre Azzalini, M arie Farge, Kai Schne-i

der [7]提出了一种基于硬阈值去噪准则的迭代算法估计去

噪阈值。在小波系数足够稀疏的表示下 , 该方法估计阈值

接近理论阈值。且计算量比 MAD 小很多。但是硬阈值函

数去噪容易产生振铃和伪吉布斯效应 , 软阈值去噪整体的

连续性较好 , 能使重构信号保持原信号的光滑性 [8], 特别是

在图像信号处理中 , 软阈值去噪的重构图像的视觉效果更

好 , 所以人们更倾向于使用软阈值函数去噪 [9]。因此本文

提出一种了基于软阈值函数的小波去噪方法。这里基于软

阈值函数并不是指套用别人的迭代方法求出一个阈值 , 然

后用该阈值和软阈值函数进行去噪 , 而是指整个迭代的过

程设计都用到了软阈值函数 , 然后根据我们的迭代算法求

出的阈值进行软阈值函数去噪。由于阈值函数选择的不

同 , 该算法的迭代函数 , 迭代停止的条件 , 迭代收敛性的证

明都和文献 [7]有所不同 , 需要重新设计算法和证明其收敛

性。本文首先给出采用迭代估计软阈值的算法 , 然后证明

了该迭代方法的收敛性 , 分析了该算法运算量 , 并通过数值

实验验证了该算法的有效性。

1软阈值去噪

1. 1阈值函数

进行小波域信号去噪时 , 首先将被噪声污染的信号 Y

进行小波分解得到小波系数 W , 然后对 W 采用适当的去噪

准则进行处理 , 得到处理后的系数 W ~, 再对 W ~进行小波逆变

换 , 得到重构信号 Y ~。

由于软阈值函数去噪能较好的保持信号原有的连续性

和光滑性 , 因而得到了广泛的应用。故本文提出的阈值估

计方法将基于软阈值函数。

设 T 阈值 , 软阈值去噪函数的定义为

G T (x) =

x -T, x >T

0, |x |[T

T +x, x <>

设源信号为 s i , 被具有噪声标准差为 R 的高斯白噪声 污染。考虑污染后的离散信号

y i =s i +e i , i =1, , , N

其中 , 信号长度为 N =2J , e i 独立同分布于 N (0, R 2) 。经过 小波变换 , 得到小波系数

W ={w-1, 0}G {w j , k B j =0, , , J -1; k =0, , , 2j -1}其中 , w -1, 0表示最粗尺度下的低频系数 , w j , k (j=0, , , J -1, k =0, , , 2j -1) 表示第 j 尺度下第 k 个小波系数。 由于噪声是高斯形式的 , 它的正交变换也将服从高斯 分布。本文提出的算法的主要思想是 :先假设小波变换的 最精细的尺度 J -1下所有的信号都是噪声 , 估计初始方差

R 20=

2J -1

E 2J-1-1

k=0

|w J-1, k |2, 求解出初始阈值 T 0=R 02ln N , 然 后通过迭代逐步将阈值逼近理论上的阈值。

下面给出基于软阈值函数去噪的阈值 T 的估计算法。 算法 1

步骤 1初始化。给定 E 足够的小。在最精细的尺度

J -1下 , 粗略估计 R 20, R 20=

M

E M -1

k=0

|w J -1, k |2, 其中 M =2J -1为该尺度下的小波系数的个数 , T 0=R 02ln N 。假设所有 的系数都是噪声 , 令噪声信号的个数 N e =M 。

步骤 2令 T =T n , 计算绝对值小于 T n 的小波系数的 个数 N e , 认为他们是噪声信号。

步骤 3计算新的方差 R 2n +1=

M

E M-1

k=0

|G T (w J -1, k ) -w J -1, k |2, T n +1=R n +12ln N 。

步骤 4n =n +1如果 |T -T n |>E , 转到步骤 2。否则 算法停止。

1. 2算法收敛性证明

考虑迭代过程 :

I B +y +使得 T n +1=I(T n )

其中

I(T) =

M

E

k=0

|G T (wJ-1, k ) -w J-1, k |2=

M

E

|w J-1, k |\T

|G T (wJ-1, k ) -w J -1, k |2+E

|w J-1, k |<>

|G T (wJ-1, k ) -w J-1, k |2=

M

E

|w J -1, k |\T

|T |2+E

|w

J-1, k

|<>

|w J -1, k |2(1)

定理 1由算法 1计算得到的序列 {Tn }n I +是收敛的。

证明 证明分 3步。

(1) 首先证明函数 I B +y +是递增的。

对任意的 T a , T b I +, 不妨设 T a [T b 。

|I(T b ) |2-|I(T a ) |2=

M

(E

|w , k |\T b

|T b |2+E

|w J -1, k |

|w J-1, k |2-

E

|w J -1, k |\T a

|T a |2-E

|w J-1, k |

|w J-1, k |2) =

M

E

|w J-1, k |\T b

(|T b |2-|T a |2) +E

|w J-1, k |[T a

(|w J-1, k |2-

|w J-1, k |2) +E

T

a

|w

J-1, k

|<>

b

(|T b |2-|w J -1, k |2) \0故 I(T b ) \I(T a )

#38 #

系统工程与电子技术

第 31卷

(

2) {Tn }n I +

是关于 n I

+

递减的。

T 20

-|I(T 0) |2

=

M

E

M-1

k=0

|w J , k |2-M

E

|w J-1, k |\T 0

|T 0|2-

E

|w J-1,k |

|w J-1, k |2

=M E |w J-1,k |\T 0

(|w J-

1, k |2-|T 0|2) \0

故 T 0\I(T 0) 。由于函数 I B

+

y

+

是递增的 , 故

I(T 0) \I . I(T 0)

以此类推 , I . I(T 0) \I . I . I (T0), , , I n (T 0) \I n +1(T 0), 即 T n \T n +1。

(

3) 序列 {T n }n I

+

是收敛的。

由于 T n \0, 所 以 序 列 {T n }n I

+

有 下 界。又 由 于

{Tn }n I

+

是递减的 , 故 {T n }n I

+

是收敛的。

1. 3 收敛速度的估计

设算法 1收敛时的迭代次数为 l, 得到的阈值为 T l 。 下面说明该算法估计阈值收敛速度快 , 而且该阈值是非常 接近理论上的阈值的。

由于小波变换的正交性和归一性 , 高斯白噪声的变换 系数仍是高斯白噪声。最精细层的高频小波系数刻画了信 号的非平稳特征 , 如边缘 , 尖峰 , 断点等信息。对于含噪声 的信号而言 , 这一层主要被噪声所控制。故

R 2

0=

M

E M-1

k=0

|

w J-1, k |

2

是噪声方差 R 2的一个很好的估计 , 即 R 0U R 。在迭代算法 的第一步迭代 , 我们有 T 0=(2ln N ) 1/2

R 0U (2ln N)

1/2

R =T ,

这说明第一次迭代得到的阈值已经非常接近理论上的阈

值了。

Berm an 证明 [10], 噪声的最大幅度有一个刚好低于 T =R , 即

lim N y +]Pr T -ln N

[max i=1, , , N

|e i |[T =1

对于最精细层的高频小波系数而言 , 即 lim M y +]

Pr T -ln M

[

max k=0, , , 2J

-1

|w J-1, k |[T =1(2)

这个结论表明当 M 充分大时 , T 是一个噪声最大模的一个 很好的估计。由式 (1) 知

I (T0) =

M

E

|w J -1, k |>T 0

|T 0|2+

E

|

w J -

1, k

|[T 0

|w J-1, k |2

T 0非常接近理论阈值 T , 由式 (2) 知 , T 0非常接近 max k =0, , , 2J

-1

|w J -1, k |, 满足 |w J -1, k |>T 0的 w J -1, k

的个数也非常的少 ,

因而有

I(T 0) =

M

E

|w J-1, k |>T 0

|T 0|2+

E

|w J-1, k |[T 0

|w J-1, k |

2

U

M

E

k=0

|w J-1, k |2=T 0

这表明阈值 T 0已经很接近 I(#) 的不动点 , 故迭代将很快

的收敛。

2 数值计算

选取 4种 M ATLAB 中典型的源信号对该算法作测 试。原始信号上加上高斯白噪声 N (0, 1) 。样本的数量 N =213=8192, 此时根据文献 [1]知 , 理论上的阈值应该为 4. 245。通过 M AD 方法得到的阈值为 T m 。

T m =

0. 6745m ed 0[k<>

J-1

(|w J-1, k |)

重建信号的质量可以由信噪比 SNR 客观的评价 , SNR 的定义为

SNR =10lg

E

N

i =1|s i |

2

E

i =1

|s i -t ~

i |2

其中 , s i 为原始信号 , t ~

i 为去噪后的信号 , N 为信号长度。

表 1中对 M AD 估计的阈值也采用软阈值函数去噪。 从表 1来看 , 本文提出的阈值估计方法估计的阈值非常的 接近理论阈值。同时除了 Block 信号外 , 本文算法重建信 号的客观评价和其他方法质量相当。对于 Block 信号这种 情况的出现 , 根据文献 [11], 如果觉得软阈值函数对系数幅 值减少了 T 的作用太强 , 可以将软阈值减少一半。在这个 例子 中 , 将阈 值减 少一 半后 , 本文 算法 的 SNR 增 加到 32. 116dB 。图 1是信号测试图。

表 1 四种测试的结果比较

迭代次数 估计阈值 SNR 文献 [7]

软阈值 M AD 文献 [7]软阈值 M AD 文献 [7]软阈值 Block 344. 2534. 1784. 27329. 69932. 98329. 677Bumps 224. 3584. 2444. 25331. 49831. 33031. 498H eaviSin e 244. 0764. 1854. 20334. 60535. 39434. 549Doppler

2

2

4. 147

4. 224

4. 233

24. 575

24. 434

24. 575

第 1期 孙蕾等 :基于迭代方法的软阈值估计小波去噪 #39#图 1测试结果的比较 :(a) ~(d) 分别为四种原始信号 , 污染后信号以及去噪后的信号

虽然本文算法的 SNR 在 Block 信号上和基于硬阈值 估计的比有一定的差距 , 但是从信号重构的图像上看 , 我们 可以很明显的看到在信号的不连续点处 , 硬阈值函数的重 构信号有很大的振荡 , 噪声的影响依然很大 , 而采用本文的 算法 , 不连续点处的振荡相对小了很多。对于其他三个信 号 , 本文的算法和硬阈值算法的去噪效果是相当的。 在计算量上 , 由于 M AD 方法要用到排序 , 最快的排序 法的运算量级为 O(N log 2N ) 。本文设计的迭代算法的运 算量为 O(Nl) , 且 1. 3节的理论分析和本节的数值试验都 表明迭代次数 l 是非常小的。为了验证本文算法在速度上 的优势 , 同时采用 /冒泡 0排序和快速排序两种排序法作为 时间参照。 /冒泡 0排序思想简单 , 运算量为 O(N 2) 。快速 排序的运算量级为 O(N log 2N) , 就排序的平均时间而言 , 快速排序被认为是目前最好的一种排序方法。由表 2可以 看出 , 在速度上 , 本算法也占有很大的优势。

表 2算法时间比较

冒泡排序 /s 快速排序 /s 迭代算法 /s Block 0. 4380. 3750. 016 Bump s 0. 4220. 4060. 016 HeaviS ine 0. 4220. 3750. 016 Doppler 0. 4220. 3750. 016

3结 论

本文提出的基于软阈值函数估计阈值大小的迭代算法 经证明是收敛的 , 分析表明该算法估计的非常接近理论阈 值。 数值实验表明该算法估计阈值的准确性和 MAD 方法 (下转第 224页 )

#224#系统工程与电子技术 第 31卷

当 az 2随时间变化时 , 状态值 x 的 h 2分量逐渐偏离正 常值 , 可见基于权值选优粒子滤波器能很好地跟踪系统状 态的变化情况。为了说明基于权值选优粒子滤波器的故障 预测算法在故障概率预测方面的优越性 , 在图 5中实线和 虚线分别表示基于权值选优粒子滤波器的故障预测算法和 SIR 故障预测算法的故障预测概率。从图 5可见 , SIR 故障 预测算法预测的概率具有较大的振荡 , 因而不能准确地估 计系统将要发生的故障。结合表 1和图 5, 根据式 (13) 可 知 :在 t =70$t 的时候 , 系统已经进入故障状态。基于权值 选优粒子滤波器的故障预测算法在 t =66$t 时的预测概率 超过 0. 5, 且逐渐增加 , 直到 t =72$t 时故障预测概率为 1, 可见基于权值选优粒子滤波器的故障预测算法能及时准确 地预测系统即将发生的故障。

3结 论

SIR 粒子滤波算法在实际应用中易受样本贫化问题的 影响 , 对于估计那些较长时间维持不变的量时样本贫化问 题的影响尤为突出 , 而系统故障预测正是属于此种类型。 针对上述问题 , 本文提出了一种基于权值选优粒子滤波器 的故障预测算法 , 按照粒子权值的大小 , 从大量粒子中选择 出比较好的粒子用于滤波 , 增加样本的多样性 , 从而缓解样 本贫化问题 , 提高故障预测算法的跟踪能力。仿真结果表 明该算法是可行有效的 , 通过与 SIR 算法的比较说明该算 法稳定性较好 , 而且具有较强的跟踪能力。但由于大量备 选粒子参与更新 , 使得该算法运算量比较大 , 如何减少算法 的计算量将是下一步的研究方向。 参考文献 :

[1]陈敏泽 , 周东华 . 动态系统的故障预 报技术 [J ].控制 理论与应 用 . 2003, 20(6) :819-820.

[2]W illiam Ng, Jack Li, S imon Godsill, et al. A review of recen t results in multiple targ et track ing[C ]M Pr oc. of the 4th Inter-national S ymp osiu m on Imag e and S ignal P rocessin g and A nal -y sis , 2005:40-45.

[3]Ch eng Chan g, Ras hid Ans ari. Kernel Particle Filter for visu al tracking [J]. IE EE S ig nal Pr oc essing letter s, 2005, 12(3) :242 -245.

[4]T ao Wei, Huang Yu fei, Philip Chen. Particle filterin g for adap -tive sensor fault detection and identification [C]M Pr oc. of the IE EE I nter national Conf er ence on Robotics an d Au tomation Or-land o, Flor id a, 2006:3807-3812.

[5]Sanjeev Arulampalam , Simon M askell, Neil Gordon, et al. A tutorial on particle filters for online non -linear/n on -gauss ian Bayesian trackin g [J ]. I EE E T rans. on S ignal Pr oce ssing , 2002, 50(2) :174-188.

[6]Djuric c M Petar, Kotecha H Jayesh, Zhang Jianqui, et al. Part-i cle filtering [J]. IE EE S ig nal Pr oc essing M ag az ine , 2003, 9: 19-38.

[7]Cody K, Dieter F, M arina M. Rea-l time particle filters [J]. P roc. of th e IE EE , 2004, 92(3) :471-472.

[8]周东华 , 叶银忠 . 现代故障诊断与容错控制 [M ]. 北京 :清华大 学出版社 , 2000:265-267.

[9]Ch en M Z, Zh ou D H. Particle filtering based fault prediction of nonlinear systems [C ]M IF A C S ymp osium Pr oc. of S af e Pr oce ss. Washing ton :E lsev ier Science , 2001:2971-2977.

(上接第 39页 )

相当 , 同时计算量远远小于 MAD 方法。和迭代估计硬阈 值去噪方法比 , 整体去噪水平差不多。虽然该算法在某些 算例中 , SNR 会稍低一些 , 但是该算法能够较好的消除硬 阈值去噪在不连续点的振荡现象。

参考文献 :

[1]Don oho D L, J oh nstone I M. Ideal spatial adaptation via w avelet s hrinkage[M ]. B iometrika , 1994, 425-455.

[2]Don oho D L, John stone I M. Adapting to u nknow n smoothness via w avelet sh rinkage[J ]. J our nal of the A mer ic an S tatistical A ssocia tion , 1995, 12(90) :1200-1224.

[3]Nason G. Wavelet shrink age u sing cross -validation[J]. Jour nal of the Royal S tatistic al S ocietey , 1996, 58(6) :463-479. [4]W eyrich N, Warh ola G. W avelet s hrinkage and generalized cross validation for image denoising[J]. I EE E T rans. on I ma ge P ro -c essing , 1998, 7(1) :82-90.

[5]李旭超 , 朱善安 . 小波域图 像降噪 概述 [J]. 中国 图像图 形学报 , 2006, 11(9) :1201-1209.

[6]关新平 , 刘冬 , 唐英干 , 等 . 基于局部方差的模糊小波阈值图像去

噪 [J ]. 系 统工 程与电 子技 术 , 2006, 28(5) :650-653. (Gu an Xinping, Liu Dong, T an g Yingan, et al. Fuzzin g w avelet thr esh olding b ased on local variance for im age den ois ing [J ]. S y stems En gineering and Electronic s, 2006, 28(5) :650-653. ) [7]Azz alin i A, Farg e M , S chneider K. Nonlin ear w avelet threshol -ding:A recu rsive m ethod to determine the optimal denoising thr esh old[J]. A pp l. Compu t. H ar mon. Anal , 2005, 18:177 -185.

[8]Donoho D L. De -noising by s oft -thres holding[J]. I EE E Tr ans. on I nf ormation The or y , 1995, 41:613-627.

[9]Chang S G, Yu B, Vetterli M. Adaptive w avelet thresholding for image d enoisin g and compr ess ion [J]. IE EE T rans. on I m -a ge P rocessing , 2000, 9(9) :1532-1546.

[10]Berman S M. S ojou rners an d extremes of stochastic process es [M ]. W adsw orth , R ead ing , M A , 1989.

[11]M allat S. A wavelet tour of sign al processin g[M]. S econd ed-i tion, Ac ademic P ress , 1999.

范文二：压缩感知的冗余字典及其迭代软阈值实现算法

第18卷第1期 电路与系统学报 V01(18 No(1 2013年2月 JOURNAL OF CIRCUITS AND SYSTEMS February，2013

01(0059(06 文章编号:1007—0249(20131

压缩感知的冗余字典及其迭代软阈值实现算法+

赵慧民?， 倪霄2

(1(广东技术师范学院电子与信息学院，广东广州510665;2(冠微科技(深圳)有限公司，广东深圳518125) 摘要一冗余字典的信号稀疏分解是一种新的信号表示理论，采用超完备的冗余函数系统代替传统的正交基函数，

为信号自适应地稀疏扩展提供了极大的灵活性。本文研究了压缩感知理论下的冗余字典、测量矩阵及其限制等容特性

(RIP，Restricted lsometry Property)，并给出了RIP、字典大小、稀疏度和测量次数的关系，提出了一种新的迭代软

阈值(IST)算法，与正交匹配追踪(OMP)算法和迭代硬阈值(IHT)算法相比较，实验结果表明了IST算法具有 高的信号恢复率。 更 关键词t压缩感知;冗余字典;迭代阈值算法;限制等容特性;测量矩

TN911(72 中圈分类号t 文献标识码:A 阵

1 引言 压缩感知CS(Compressed Sensing)的稀疏信号表示理论通常基于正交线性变换，但许多信号是

各种自然现象的混合体，这些混合信号在单一的正交基变换中不能非常有效地表现出来。例如，一个

含有脉冲和正弦波形的混合信号，既不能用单一的脉冲基函数，也不能用单一的正弦基函数有效地表

示。为此，CS理论采用一种新的信号稀疏表示技术，称为冗余字典，来表示这类混合信号。其基函数

用称之为字典的超完备的冗余函数系统取代，字典的选择尽可能好地符合被逼近信号的结构，其构成

可以没有任何限制，而字典中的元素被称为原子。从字典中通过测量矩阵找到具有最佳线性组合的m

项原子来表示一个信号，称作信号的稀疏逼近或高度非线性逼近。所以，测量矩阵的选择是关键，而 压缩感知测量矩阵满足的充分必要条件就是RIP。因此，RIP为CS信号的稀

测量和精确重构提供了理论保证。 疏 CS理论的三个组成要素是信号

的 稀疏变换(目前的稀疏变换有离散

余弦 变换DCT，小波Wavelet，

Curvelet，过 (a)原始图像 (b)噪声污染80,图像 (c)CS恢复图像 atom完备原子分解(Overcomplete 图l CS重构恢复的原始信号( Decqmposition)等;稀疏信号的非相干 测量及稀疏信号的重建算法 1。CS主要思想是对一类具有稀疏先验的信号，先经小部分非线性测量矩

阵进行采样，其包含足够信息良好逼近信号，再通过一定类型的线性或非线性解码机制就可高概率精

确重建原始信号(见图1)。其中，非相干测量矩阵的RIP界定和快速有效的迭代阈值实现算法是CS 理 论实用化的关键因素之一，也是CS的主要研究内容之一。 2 相关知识 CS的基本问题是确定线性非自适应测量矩阵，使稀疏信号测量样本数据玎最小化，最终使s个

非 零元素的信号x?R。能够稳定恢复。其中以的每一个矢量值都可作为稀疏信号z的内积，所有矢 作为矩阵甲的行，吵=(吵(，y:，(((，少。)，妒?R”。，即有S=败。基于上面的概念，所谓的CS理量值

论为:在

2012-06-11 +收稿日期t 修订日期:2012-08—26

基金项目:国家自然科学基金资助项目(61272381);广东省自然科学基金项目($2012010008639);广东省科技计划项目 (20128010100035)

万方数据

电路与系统学报 第18卷

适当选择的稀疏基y上具有s稀疏描述的以个采样信号石，可通过它在另一组非相干基 m(S兰m<>，i?{1 2((，m)获得精确重建‘扪。即 痧=(卵，杉，(((，《 )上ffl y:m:卯rs。其中，咖为测量基，其构成的m×玎维矩阵为测量矩阵。咖妒 7为传感系统矩阵，J，为

测量向量。从S中恢复x的方法，主要是凸优化算法和贪婪算法。凸优化算法是基于三，最小范数模型

实现的，主要包括迭代阈值和Bregman(布雷格曼)算法。贪婪算法是基于迭代计算工的支撑，包括正 Pursuit)[5，61。这些Pursuit)和基追踪BP算法(Basis 交匹配追踪算法OMP(OrthogonalMatching c{l'(((，方 法对测量矩阵西都有一定的严格要求，即咖必须满足RIP条件。精确地说，设支撑A d}，由以

索引的列组成矩阵罗的子矩阵y以，其局部RIP常数以=以(罗)为最小值，且满足不等式:

(1) (1一以)||工||;?IIy。xll;-<(1+ 一)iixll;="">

而对于全局RIP常数则为: (2)S?N 喀=菇(咖:=sup8A(e)， l,lJ=S 基于上述理论及其条件，本文主要解决的问题是，当信号不能用正交基稀疏表示时，采用冗余字典

D?R“‘表示，并在满足RIP要求下如何选择合适的测量矩阵咖?R“4使传感系统矩阵咖P7具有更少

非零值【”，同时运用迭代阈值算法概率恢复原始的信号。也就是说，如何从字典中D中找到具有的 最佳 线性组合的m项原子来表示一个信号，实现信号的稀疏逼近或高度非线性逼近。

3 非相干字典、测量矩阵与RIP界定问题

从非线性逼近的角度来讲，高度非线性逼近包含两个层面:一是根据目标函数从一个给定的基库

中挑选好的或最好的基;二是从这个好的基中拣选最好的m项组合。利用贪婪算法和凸优化自适应追 踪，从一个冗余函数系统中进行m项逼近或阈值逼近也属此例。下面从非相干字典与RIP界定方法研 此问题。 究 3(1非相干字典 要想在冗余字典中获得高度非线性逼近的建设性结果，必须首先将注意力集中在某 些特别的字典

上【8】。许多研究人员瞄准了非相干字典， 也就是说相干系数口小于某个常数的字典。相干系数的定义

为:

(3),z=supI(识，?， )| 暂

当相干系数较大时，原子间的相互关联也较强。如果,2=1，则意味着字典中至少包含了两个一模

一样的原子。反之，当相干系数较小时，就称字典是非相干的，正交基的相干系数为零。相干系数为

字典的冗余性提供了另一种可能的测度手段。 (3)式说明当?充分小时，字典D(虽然可能是超完珞 的)接近一个正交基。 3(2测量矩阵 为了实现最佳的逼近并重构稀疏信号，Candes和Tao给出并证明了传感矩阵多妒7必须满足限制

容性条件RIP[21。Baraniuk在文献[9】中给RIP的等价条件是测量矩阵咖和稀疏表示的基妒不相 等

关，即 要求咖的行痧，不能由y的列{6，i稀疏表示，且y的列缈;不能由痧的行妒，稀疏表示。满足RIP

条件的测量 矩阵有高斯随机矩阵、贝努利矩阵、一致球测量矩阵、二值随机矩阵、局部哈达玛矩阵 (Toeplitz)矩阵等。对于大小为聆×d的随机矩阵咖，其满足的条件是【?】: 以及托普利兹 n(2 (4)P(]ll西D』2一忪吲?占怕n?2e一 2，占?(0,1, 3)

其中，常数c>O，任意0?R4。 3-3 RIP界定问题

万方数据

61 第1期 赵慧民等:压缩感知的冗余字典及其迭代软阈值实现算法

设D?R出K为大小为K的字典，具有限制性等容常数蠡(D)，S?N。对于满足(4)式的随机测量

矩阵西?R“4，当测量次数为:

(5)n?C8—2(Slog(K,S)+log(2e(1+12,8))+f)， 常数C?9,c，万?(0，1)，f>0 则传感系统矩阵如具有限制等容性常数: (6)蟊(吐移)?菇(D)+万(1+瓦(D))

对于固定的万和计数器t，测量次数可以更精确地通过字典大小K和稀疏度s表示为:

(7) 门?删l。 g了K

(1)字典是一个正交基情况 这时字典D是正交基，那么6(D)=0。这时对于随机测量矩阵西， RIP的条件就是公式(1)的不

等式。如果己知字典D的相干系数?和Babel函数关系??(K): 2川m吨a刨x蚤|<谚，办>I，那么RIP问题表示 为: (8)蟊??1(S一1)?(S一 1)M

这里，?可通过字典大小和信号维数界定为,2>?三芸,，等价于?, c,?万，也就是说字典越大相干

系数越小。这样，对字典的稀疏度界定为s<?d,c。 (2)字典由规范正交基(Dirac基)和DCT基联合组成情况 设Dirac基和DCT基联合组成字

典D的维数为d=22川(P为控制字典搜索收敛速度因子，P越

小字典收敛越快)，则D的相干系数,2=?2肛=2,，D的原子数为K=22川。这时，在BP 算法下，当 信号稀疏度S<2p一6时，根据(6)式传感矩阵彻具有限制等容性常数为万><1,3。把s、k代入(5) 式，测量次数为:="" (9)n?cl(4s(2pl092一logs)+c2+f)="">

Basis)基下稀疏，则测量次数为:如果信号仅在Dirac(Canonical

(10)11?c,(4S(2pl092一l092S)+c :+t)

4 信号恢复的迭代软阈值算法

假设D为dxK的字典，A是大小11×d的随机高斯测量矩阵，“为字典D的Babel函数相干

Y Ih=1。则信号在解码恢复时，设接收的信号Y=D。x的支撑为,1，且IAI=S，归一化有 |1系数。 (8)式RIP的界定，推导出信号X成功恢复的阈值条件为: 根据(6), l’一 (11) HI(x)=导掣>“(s)+,2l(S-1) |l工11? 1)可知，字典下的阈值门限取决于字典的相干系数和支撑的大小，称这里，I x,I=mini薯I。由(1

Soft 迭代软阈值算法(IST算法，Iterative Thresholding Algorithm)。 (1)测量次数11和字典大为 小K以及稀疏度S的关系 1)，可以从玎维测量矢量(即n次测量样本)z=匆当信号残差为占时，通过迭代阈值条件(1 中高概率恢复原始信号z: =ADx (12)"?C(e)(1092K+t1

z7(6885。其中，c(占)=4Cl,。2+2C26"一，CI和C2为信号收敛调节因子，并有C1={ ;?2( 5044，C2=?8e ?6厅 残差占可由字典迭代更新过程中的原子间最佳误差估计得到【?]:

万方数据

第18卷 62电路与系统学报

r如班ImMax 州E如么)Is=卿I(赢，谚)I- cbmin =I‰砂丽1 ，

l—II工Il。?1(s一1)一||xII 。?1(sp ?百—鲁 -(IXmin

(13)代入由于(12)式中，c(占)=4Cl占一2十2C2占一1?c3占-2 C3?4Cl+2C2?25(40。把公式(11)、

(12)式的c(e1有:

， c?鲻小船【涮叫c 沪州?，j

因此， (12)式中n的大小可写为:

(15) 娩25(4s(1+绯弘1))【制叫 (沪麒 p1)](1092K+t)

这时Babel函数相干 特殊地，当字典D为一个正交基字典时，信号总可以满足恢复的RIP限定条件，

系数??岱)2?-(s一1)=o，信号恢复的测量次数 n降为j，

(16)

脱5小s，(黔卜 2nD

(2)信号X恢复的概率计算 根据Hilbert卒闻

卜范辑的定理[12】(有关系: (17)ra刚inI(砂，彳破)|>ma— xI(砂，彳么)J kaA

P(?辫I(匆，爿谚)|>m登I(砂，4焱 )I)<>

七E^ +可=I(_:以,七?,)pI 警-, I(弘织)l+主]

P(嘶?删?卿?州I一主)?P必 (圳州似删一詈 )} ?，

’ lE 、一l。o 2 ，?驴A(埔)一(砂以)l?争<2l小冲卜看‰ )="">

P(=?删?詈?枷詈H小圳州 (训鲥+詈)}?，

I(刁 l e即(一疗高]?三P(I(4y，彳以 )一 (y，纯)陲詈)?2 ’ 、 一 ’ 七E^

e冲卜高)?，噌枷纠?芎枷椰 cxD-，疗高卜2I

叫 砷e冲卜者勃卜ed 一”者舞]

万方数据

63第1期 赵慧民等:压缩感知的冗余字典及其迭代软闽值实现算法

e(m辨l

@2’

小2细卜高 ] 可见，当聆一定时，信号恢复的概率和字典的迭代更新误差以及字典的大小有关。因此，迭代阈

值算法能有效减少每次迭代中的误差，在一个可达到的最佳误差估计的一个常数内被保证。事实上，

这个算法需要一个固定的迭代数目，其仅仅依赖于测量矩阵的稀疏度S的一种计算形式。 (3)IST算 法的实现过程如下:

1)初始信号设置为Xo=0，初始残差为占=Y，计算器(即迭代次数初值)t=-0。

2)计算疗=C(e)(1092K+，)，xt+l=xt+D7(y一三啊)。

3)选择字典D的最大K个元:X川=H。(x川)。其中，H。是阈值算子，即保持能够恢复信号 数量 大小的最大K个元;设置其它元时为零。

满足条件t时(可选择一定的迭代次数)迭代停止，否则扣t+1，返回第2)步。 5 实验结果

一 实验用Dirac基和DCT基构成 大小d=256的字典D，则其相干

zo(0884。本文试数?:?l,128 系

验6 个大小为n×d的高斯随机 测量矩 (a)Dirac基和DCT基构成字典D (b) 一个正交基(ONB) 阵爿，其中珂的变化范围为64,224， 图2 支撑,1和测量次数门作为OMP算法函数时信号的恢复率 测量元素仍，』,N(O，玎一)。稀疏度S

变化范围为4,64，在字典大小和稀疏

度防(爿 5之

撑4。对本身稀疏 ? 间随机选择 (葺)j。A可能的支 , 1,1 的0—1随机信号进 行

(a)Dirac基和DCT基构成字典D (b)一个正交基(ONB) 测量，本文主要 图3 支撑和测量次数n作为IHT算法函数时信号的恢复比较了本文IST算

法与迭代硬阈值算 率 法(IHT)[13】和

OMP算法[51恢 信号时测量次复 数和

5010 20 30 40 支撑之间的关系，

,1及其在此关系下信 (a)Dirac基和DCT基构成字典D (b) 一个正交基(ONB) 号的恢复率，如图 图4支撑,1和测量次数，l作为IST算法函数时信号的恢复率 2,4所示。

由图2,4可见，在充分测量次数下，各种算法均可高概率恢复信号。但在测量次数一定时，本

文 提出的IST算法具有更高的信号恢复率。这说明，本文的IST算法实现速度比IHT和OMP算法 实验中，IST算法也有一点缺陷，即测量次数取决于最大和最小基系数的比率。 更快。

6 总结和展望

万方数据

第18卷电路与系统学报

通过分析CS理论中RIP条件，深入研究了一个正交基ONB与Dirae基和DCT基构成的

冗余字 典，确定了随机测量矩阵对信号测量次数的要求，以及测量次数与字典大小、信号稀疏度的

关系。在 此基础上，结合凸优化算法和贪婪算法算法的优却点，提出了一种快速有效新的迭代软阈 现算法。该算法在测量次数玎一定时，信号恢复的概率和字典的迭代更新误差以及字典的大 值(IST)实 小有关。 本文认为，虽然此算法不是最优的，但可以为信息安全以及数字水印数据的稀疏生成提供研

究工具。 本文进一步地研究方向是，寻找更优的测量矩阵，发展更有效和低复杂度的图像,视频压缩 其信号的重构算法。 感知系统及

参考文献; Trans Inform 306[1】D Donoho(Compressed sensing【J1(IEEE Theory，2006，51(4):1289—1 E Tao(Robust reconstruction from [21 Candes，J Romberg，T uncertainty principles:Exact signal highly incomplete frequency information【J 】J IEEE Trans Inform，2006，52(2):489-509( M reconstruction convex relaxation:Fourier and Gaussian Proe(CISS [3】 Rudelson，R Vershynin(Sparse by measurements【A】(In 2006(40th on Annual Conference Information Sciences and Systems)[C】(2006( 4 S J Kim，K Koh，M method for with Lustig，S Boyd，D Gorinevsky(A interiorpoint large—scale L1-regularized least-squares problems and signal processing statistics【J】(JournalofMachineLearningResearch，2007，7(8):1519—1555( applicationsin 5 J from random measurements via Tropp，A Gilbert(Signal recovery orthogonal matching pursuit【J】(IEEE Tram(Inform(Theory，2007， 53(12):4655-4666(

6 杨海蓉，张成，丁大为，韦穗(压缩传感理论与重构算法[J】:电子学报，201 l，39(1):142—148 (7 Rauhut H，Sehnass P and redundant Transactions on Information K，Vandergheynst Compressed sensing dictionaries【J】IEEE Theory，

2008，54(5):2210—2219( 8 S Chen，D Donoho，M Saunders(Atomic basis Sci decomposition by pursuit【J】(SIAMJ Comput，1999，20:33—61( R9 Baraniuk G(Compressivesensing[J]IEEESignalProcessing Magazine，2007，24(4):118-121( 1 andviaD Needell，RVershynin(Uniformuncertaintyrecoveryprinciple signal regularized orthogonal matching pursuit【J】Found(Comput 1 7-334( Math，2009，9(3):3 X Bai(Minimum error L l Mei，H Wu，E Blasch，L bounded efficient tracker with occlusion IEEE Conference on [11 】Ling，Y detection[A】(in Vision and Pattern 1(Computer Recognition【CVPR)[C](201 112】 杨必成(算子范数与Hilbert型不等式[M](北京:科学出版社，2009， T E Davies(Iterative for [13 】Blumensath，M of Fourier and thresholding sparse approximations【J】Journal Analysis Applications，2008， 14(5—6、:629?654( 作者简介:赵?民(1966()，男，广东省图像图形学会理事，广东技术师范学院教授，硕士研宄生导师，2001年获 中山大学信号与信息处理专业博士学位，主要研究方向为数字信号处理与信息安全技术。 Redundant dictionaries of and an compressed sensing application algorithm

of iterative SOft threshold

ZHAOXia02Hui(min91，NI Normal (1(Guangdong Polytechnic University，Guangzhou 510665，China; 5 18 1 Technology(Shenzhen)Limitted，Shenzhen 25，China) 2(Grandvio

of redundant dictionaries is a new for Abstract:Signal decomposition can sparse theory signal representation(The theory

a flexible method for extension via redundant function of adaptively provide sparsity instead signal using overcomplete

orthonormal-basis function(Based on for redundant dictionaries and measurement matrix as well investigation conventionalas restricted for a novel of iterative soft isometry property compressed sensing，the paper presents algorithm IST with the iterative hard Comparing proposed oahogonal matching threshold(IST)( pursuits(OMP)and

results show that the IST has ratio(threshold(IHT)respectivelyexperiment， algorithm higher signal recovery threshold Key words:compressed sensing;redundant dictionaries;iterative algorithm;Restricted Isometry Property; measurement matrix

万方数据

范文三：迭代阈值法

数字图像处理的目的之一是图像识别, 而图像分割是图像识别工作的基础。图像分割是指把图像分解成具有特性的区域并提取出感兴趣目标的技术和过程，是计算机视觉领域的一个重要而且基本的问题，分割结果的好坏将直接影响到视觉系统的性能。因此从原理，应用和应用效果的评估上深入研究图像分割技术具有十分重要的意义。

本课题主要介绍了图像分割的基本知识。图像分割的算法有阈值分割法，边缘检测法，区域分割等，本设计重点介绍了基于最小点阈值方法，基于最优阈值分割方法，基于迭代图像分割方法，最大类间方差法（OTSU）的图像分割法的原理和他们的MATLAB的实现代码与运行结果。

关键词：图像分割； MATLAB；阈值分割；

1 课程设计目的.......................................................................................................... 3

2 课程设计要求.......................................................................................................... 3

3 相关知识.................................................................................................................. 3

3.1 图像分割的概述........................................................................................... 3

3.2 阈值分割的基本原理................................................................................... 4

3.3 阈值分割方法的分类................................................................................... 5

3.3.1 基于点的全局阈值方法.................................................................... 6

3.3.2 基于区域的全局阈值方法................................................................ 6

3.3.3 局部阈值法和多阈值法.................................................................... 6

4 程设计分析.............................................................................................................. 6

4.1 基于迭代的方法实现图像切割................................................................... 6

4.2 最大类间方差的方法实现图像切割........................................................... 7

5 程序设计.................................................................................................................. 7

5.1 程序简单介绍............................................................................................... 7

5.2 程序代码....................................................................................................... 8

6 结果与分析............................................................................................................ 11

结束语.......................................................................................................................... 13

参考文献...................................................................................................................... 14

迭代阈值法

1 课程设计目的

本设计的课题任务是掌握图像阈值分割算法研究，实现对图像的分割。了解图像分割的应用及基本方法，理解阈值化图像分割原理，理解三类典型的阈值化分割算法，并利用之进行图像分割，给出实验结果并做出分析。

2 课程设计要求

⑴查阅相关资料；

⑵理解基于各像素值的阈值分割算法，基于区域性质的阈值分割算法, 基于坐标位置的阈值分割算法；软件编程实现利用基于各像素值的阈值分割算法进行图像分割，要求完成如下内容：包括极小值点阈值、最优阈值、迭代阈值，基于最大方差的阈值，基于最大熵的阈值等方法，利用之实现图像分割，这里的图像可以针对核磁共振图像（MRI）

⑶用MATLAB实现，并观察各算法之间的区别。 [2][1]

3 相关知识

3.1 图像分割的概述 在对图像的研究和应用中，人们往往仅对图像中的某些部分感兴趣，这些部分称为目标或前景(其他部分称为背景)，他们一般对应图像中特定的、具有独特性质的区域。为了辨识和分析目标，需要将他们分离提取出来，在此基础上才有可能对目标进一步利用。图像分割就是指把图像分成格局特性的区域并提取出感兴趣目标的技术和过程。这里特性可以是象素的灰度、颜色、纹理等，预先定义的目标可以对应单个区域，也可以对应多个区域。现有的图像分割算法有：阈值分割、边缘检测和区域提取法。本文着重研究基于阈值法的图像分割技术。 所谓图像分割是指根据灰度、彩色、空间纹理、几何形状等特征把图像划分成若干个互不相交的区域，使得这些特征在同一区域内，表现出一致性或相似性，而[3]

在不同区域间表现出明显的不同。简单的讲，就是在一幅图像中，把目标从背景中分离出来，以便于进一步处理。图像分割是图像处理与计算机视觉领域低层次视觉中最为基础和重要的领域之一，它是对图像进行视觉分析和模式识别的基本前提。同时它也是一个经典难题，到目前为止既不存在一种通用的图像分割方法，也不存在一种判断是否分割成功的客观标准。

阈值法是一种传统的图像分割方法，因其实现简单、计算量小、性能较稳定而成为图像分割中最基本和应用最广泛的分割技术。已被应用于很多的领域，例如，在红外技术应用中，红外无损检测中红外热图像的分割，红外成像跟踪系统中目标的分割；在遥感应用中，合成孔径雷达图像中目标的分割等；在医学应用中，血液细胞图像的分割，磁共振图像的分割；在农业工程应用中，水果品质无损检测过程中水果图像与背景的分割。在工业生产中，机器视觉运用于产品质量检测等等。在这些应用中，分割是对图像进一步分析、识别的前提，分割的准确性将直接影响后续任务的有效性，其中阈值的选取是图像阈值分割方法中的关键技术。 [4][3]

3.2 阈值分割的基本原理

图像阈值化分割是一种最常用，同时也是最简单的图像分割方法，它特别适用于目标和背景占据不同灰度级范围的图像。它不仅可以极大的压缩数据量，而且也大大简化了分析和处理步骤，因此在很多情况下，是进行图像分析、特征提取与模式识别之前的必要的图像预处理过程。图像阈值化的目的是要按照灰度级，对像素集合进行一个划分，得到的每个子集形成一个与现实景物相对应的区域，各个区域内部具有一致的属性，而相邻区域布局有这种一致属性。这样的划分可以通过从灰度级出发选取一个或多个阈值来实现。

阈值分割法是一种基于区域的图像分割技术，其基本原理是：通过设定不同的特征阈值，把图像像素点分为若干类．常用的特征包括：直接来自原始图像的灰度或彩色特征；由原始灰度或彩色值变换得到的特征．设原始图像为f(x，y)，按照一定的准则在f(x，y)中找到特征值T，将图像分割为两个部分，分割后的图像为 : [5]

g x,y = b0,f x,y

若取b0=0（黑）,b1=1（白），即为我们通常所说的图像二值化。

一般意义下，阈值运算可以看作是对图像中某点的灰度、该点的某种局部特

性以及该点在图像中的位置的一种函数，这种阈值函数可记作

T=T x,y,n x,y ,f x,y (3.2.2)

式中f(x，y)是点(x，y)的灰度值；n(x，y)是点(x，y)的局部邻域特性．根据对T的不同约束，可以得到3种不同类型的阈值，即

(1)点相关的全局阈值

T=T f x,y (3.2.3)

只与点的灰度值有关

(2)区域相关的全局阈值

T=T n x,y ,f x,y (3.2.4)

与点的灰度值和该点的局部邻域特征有关

(3)局部阈值或动态阈值

T=T x,y,n x,y ,f x,y (3.2.5)

它与点的位置、该点的灰度值和该点邻域特征有关

因此本文分三大类对阈值选取技术进行综述：

1) 基于点的全局阈值方法；

2) 基于区域的全局阈值方法

3) 局部阈值方法和多阈值方法

3.3 阈值分割方法的分类

全局阈值法指利用全局信息对整幅图像求出最优分割阈值，可以是单阈值，也可以是多阈值；局部阈值法是把原始的幅图像分为几个小的子图像，再对每个子图像应用全局阈值法分别求出最优分割值。其中全局阈值法又可分为基于点的阈值法和基于区域的阈值法。阈值分割法的结果很大程度上依赖于阈值的选择，因此该方法的关键是如何选择合适的阈值。由于局部阈值法中仍要用到全局阈值

法，因此本文主要对全局阈值法中基于点的阈值法和基于区域的阈值法分别进行了研究。根据阈值法的原理可以将阈值选取技术分为三大类

3.3.1基于点的全局阈值方法 基于点的全局阈值算法与其他几大类方法相比，算法时间复杂度较低，易于实现，适合应用于在线实时图像处理系统。

3.3.2基于区域的全局阈值方法

对一幅图像而言，不同的区域，比如说目标区域或背景区域，同一区域内的象素，在位置和灰度级上同时具有较强的一致性和相关性。

3.3.3局部阈值法和多阈值法 局部阈值(动态阈值) 当图像中有如下一些情况：有阴影，照度不均匀，各处的对比度不同，突发噪声，背景灰度变化等，如果只用一个固定的全局阈值对整幅图像进行分割，则由于不能兼顾图像各处的情况而使分割效果受到影响。有一种解决办法就是用与象索位置相关的一组阈值(即阈值使坐标的函数)来对图像各部分分别进行分割。这种与坐标相关的阈值也叫动态阈值，此方法也叫变化阈值法，或自适应阈值法。这类算法的时间复杂性和空间复杂性比较大，但是抗噪能力强，对一些用全局阈值不易分割的图像有较好的效果。

多阈值法很显然，如果图像中含有占据不同灰度级区域的几个目标，则需要使用多个阈值才能将他们分开。其实多域值分割，可以看作单阈值分割的推广。 [7][6]4 程序设计分析

4.1 基于迭代的方法实现图像切割

迭代法是基于逼近的思想，其步骤如下：

(1)求出图象的最大灰度值和最小灰度值，分别记为Zmax和Zmin，令初始阈值 T0= Zmax+Zmin ÷2[2] (4.1.1)

(2)根据阈值Tk将图象分割为前景和背景，分别求出两者的平均灰度值Z0和Zb；

(3)求出新阈值:

Tk+1= Z0+Zb ÷2[2] (4.1.2)

(4)若Tk=Tk+1，则所得即为阈值；否则转到第二步，迭代计算。

迭代所得的阈值分割的图象效果良好。基于迭代的阈值能区分出图像的前景和背景的主要区域所在，但在图像的细微处还没有很好的区分度。

但令人惊讶的是，对某些特定图象，微小数据的变化却会引起分割效果的巨大改变，两者的数据只是稍有变化，但分割效果却反差极大

经试验比较，对于直方图双峰明显，谷底较深的图像，迭代方法可以较快地获得满意结果。但是对于直方图双峰不明显，或图像目标和背景比例差异悬殊，迭代法所选取的阈值不如最大类间方差法。

4.2 最大类间方差的方法实现图像切割 由Otsu于1978年提出的最大类间方差法以其计算简单、稳定有效，一直广为使用。从模式识别的角度看，最佳阈值应当产生最佳的目标类与背景类的分离性能，此性能我们用类别方差来表征，为此引入类内方差、类间方差和总体方差。最大类间方差法计算简单、稳定有效，一直广为使用，是一种受到普遍欢迎的阈值选取方法。其基本思路是将直方图在某一阈值处分割成两组，当被分成的两组的方差为最大时，得到阈值。因为方差是灰度分布均匀性的一种量度，方差值越大，说明构成图像的两部分差别越大，当部分目标错分为背景或部分背景错分为目标都会导致两部分差别变小，因此使类间方差最大的分割意味着错分概率最小。 [3][6]

5 程序设计

5.1 程序简单介绍

程序除了有题目中要求的自编程序实现最大类间方差法外，还扩展了迭代法阈值分割，功能较多，为了使程序有良好的人机界面，主程序中提供了简单的菜单界面。同时为了增加程序的可读性，分模块编写，然后在主函数中调用。

5.2 程序代码

function main %主程序

clear;close all; %清除变量，关闭所有窗口

[name path]=uigetfile('C:\Documents and Settings\All Users\Documents\My Pictures\示例图片\Winter.jpg');%获取图片路径及文件名

I=imread([path name]);%读图

I=rgb2gray(I);%彩色转黑白

while 1 %循环执行

disp('0 原图 1 OTSU 2 全局阈值 3 迭代法 4 退出');%提示信息

select=input('请选择：'); %选择相应功能

switch select

case 0

yuantu(I); %输入0 显示原图

case 1 %输入1 OTSU法

fun1(I);

case 2 %输入2 全局阈值法

fun2(I);

case 3

fun3(I) %输入3 迭代法

case 4

return; %输入4 退出

otherwise %其他值不处理

end

end

% OTSU

function fun1(I) %阈值计算程序

Ni=imhist(I);%计算直方图数组

N=sum(Ni); %总像素点个数

delamax=0; %类间方差最大值

threshold=0; %阈值

for k=2:255

u=dot([0:255],Ni/N); %图像的总平均灰度级

w0=sum(Ni(1:k)/N); %C0类像素所占面积的比例

w1=1-w0; %C1类像素所占面积的比例

if w0==0|w0==1 %当w0为1或0时提前结束本次循环 continue

end

u0=dot([0:k-1],Ni(1:k)/N)/w0; %C0类像素的平均灰度

u1=dot([k:255],Ni(k+1:256)/N)/w1; %C1类像素的平均灰度 dela(k)=w0*(u-u0)^2+w1*(u-u1)^2; %类间方差公式

%求出类间方差的最大值，最大时的那个值对应的k值存入delamax if dela(k)>delamax

delamax=dela(k);

threshold=k-1;

end

end

BW1=im2bw(I,threshold/255); %阈值分割

figure

imshow(BW1),title('自编程序运行结果')%显示图片

disp('自编程序求的的阈值：')%显示提示信息

disp(threshold) %显示结果

%调用系统函数OTSU阈值分割

level=graythresh(I);%调用系统函数计算阈值

BW2=im2bw(I,level); %阈值分割

figure

imshow(BW2),title('调用库函数运行结果')

level=round(level*255);

disp('调用库函数求得的阈值：')

disp(level)

%全局阈值

function fun2(I) %阈值分割程序

[width height]=size(I);;%获取图片宽高

th=input('请输入阈值:');

for i=1:width

for j=1:height

if(I(i,j)

BW1(i,j)=0;

else % 灰度大于等于阈值时则为白色 BW1(i,j)=1;

end

end

end

figure

imshow(BW1),title('全局阈值')

%迭代法

function fun3(I) %迭代法求阈值

J=double(I);

T=(min(J(:))+max(J(:)))/2; %T的初始阈值

done=false; %done 的初始值为false

i=0;

while ~done %当done为false是则执行循环

r1=find(J<=t); %找出不大于t的所有像素点索引值="" r2="find(J">T); %找出大于T的所有像素点索引值 Tnew=(mean(J(r1))+mean(J(r2)))/2; %计算新的阈值 done=abs(Tnew-T)<1; %判断|t2-t1|是否小于t0="" t="">

i=i+1;

end

J(r1)=0; %小于等于阈值的为黑

J(r2)=1 ; %大于阈值的为白

figure

imshow(J),title('迭代法求阈值')

function yuantu(I) %显示原图

figure

imshow(I);title('原图')

6 结果与分析

图 6.1 原图，自编程结果图，调用库函数图

不同阈值的全局阈值图，迭代法求阈

值图的比较

结果分析：

调用系统函数与自编的程序求得的阈值一样，处理的结果完全一样。手动输入阈值分割的话输入的阈值不同处理的结果也不一样，阈值差别越大图像的区别越明显但如果输入的阈值和OTSU算出的阈值一样，处理结果就完全一样。另外采用迭代法，通过多次用不同的图片处理发现算出的阈值和最大类间方差法基本一样，但有时会有误差，处理结果基本一样。

结束语

通过本学期对数字图像处理的学习，在老师的指导下，对选取的图像进行迭代阈值法的处理得到处理结果，并对结果进行分析。这个过程分为三步：

第一步：熟悉所运用的仿真软件MATLAB，查阅相关资料了数字图像处理中的

迭代阈值法的原理。

第二步：设计相应的程序并选取合适的图像，运用仿真软件对图像进行仿真处理，

从而得到处理结果。

第三步：对得到的处理结果进行分析，比较不同的迭代法的优缺点，加深对这种

方法的理解。

通过这个课程设计的学习，使我对数字图像处理的理解更加深刻了。数字图像处理和我们的生活息息相关，理解它就是理解生活，同时我也感受到老师对我们学习科学知识的重要性，并且明白想要知道更多的东西就要在动脑的同时多动手进行实践，帮助我们去理解，即理论与实际相结合。

参考文献

[1] 王家文.MATLA.图形图像处理[M].北京，国防工业出版社.2001:12-17:30-55

[2]甲永红.数字图像处理[M] .武汉，武汉大学出版社.2011:71-84.:25-68

[3]李红梅.二值图像的阈值分割方法探讨[J] .上海，科技经济市

场.2007:38-42:68-79

[4] 刘刚. MATLAB数字图像处理[M].北京，机械工业出版社.2010：34-45.:128-205

[5]张强, 王正林.《精通MATLAB图像处理》. 电子工业出版社，1997:99-106

[6]陈怀琛.《MATLAB及其在课程中的应用指南》. 西安电子科技大学出版社，2000:256-306

[7]朱习军.《MATLAB在信号与系统与图象处理中的应用》. 电子工业出版社，2002:78-109

范文四：迭代阈值算法

(1) 求出图像的最大灰度值和最小灰度值，分别记为:ZMAX和ZMIN，令初始阈值

T0=(ZMAX+ZMIN)/2;

(2) 根据阈值T将图像分割为前景和背景，分别求出两者的平均灰度值ZO和ZB; (3) 求出新阈值T=(ZO+ZB)/2;

(4) 若两个平均灰度值ZO和ZB不再变化，即T为阈值，否则转2.迭代计算。 实现如下:

I=imread('gaoshanbei01.jpg'); >> ZMax=max(max(I));

>> ZMin=min(min(I));

>> TK=(ZMax+ZMin)/2;

>> bCal=1;

>> iSize=size(I);

>> while(bCal)

iForeground=0;

iBackground=0;

ForegroundSum=0;

BackgroundSum=0;

for i=1:iSize(1)

for j=1:iSize(2)

tmp=I(i,j);

if(tmp>=TK)

iForeground=iForeground+1;

ForegroundSum=ForegroundSum+double(tmp);

else

iBackground=iBackground+1;

BackgroundSum=BackgroundSum+double(tmp);

end

end

end

ZO=ForegroundSum/iForeground; ZB=BackgroundSum/iBackground; TKTmp=uint8((ZO+ZB)/2);

if(TKTmp==TK)

bCal=0;

else

TK=TKTmp;

end

end

disp(strcat('迭代后的阈值:',num2str(TK))); newI=im2bw(I,double(TK)/255); >> subplot(121),imshow(I)

>> subplot(122),imshow(newI)

范文五：基于回溯的迭代硬阈值算法

第37卷第3期2011年3月

自动化学报ACTAAUTOMATICASINICA

Vol.37,No.3March,2011

基于回溯的迭代硬阈值算法

杨海蓉1,2

方红3

张成1

韦穗1

摘要针对压缩传感(Compressedsensing,CS)理论中迭代硬阈值(Iterativehardthresholding,IHT)算法迭代次数多和时间长的问题,提出基于回溯的迭代硬阈值算法(Backtracking-basediterativehardthresholding,BIHT),该算法通过加入回溯的思想,优化了IHT算法迭代支撑的选择,减少支撑被反复选择的次数.模拟实验表明,在保证重建质量的前提下,相比较于IHT和正规化迭代硬阈值(NormalizedIHT,NIHT)算法,BIHT算法的重建时间降低了2个数量级.用本身稀疏的0-1随机信号的重建实验表明,若测量次数和稀疏度相同,BIHT算法的重建概率高于IHT算法.关键词DOI

压缩传感,迭代硬阈值,正规化迭代硬阈值,回溯,稀疏10.3724/SP.J.1004.2011.00276

IterativeHardThresholdingAlgorithmBasedonBacktracking

YANGHai-Rong1,2

FANGHong3

ZHANGCheng1

WEISui1

AbstractThebacktracking-basediterativehardthresholding(BIHT)algorithmisproposedtosolvetheproblemthatthenumberofiterationsistoolargeandtheiterationtimeistoolongwhentheiterativehardthresholding(IHT)algorithmisappliedtothecompressivesensing.TheBIHTalgorithmoptimizesthesub-optimalchoiceofsupportsforeachiterationandreducesthetimesofsomesupportsiteratedrepeatedlybyaddingtheideaofbacktracking.Thesimulationdemonstratesthatbacktracking-basedalgorithmensuresthereconstructionqualityanddecreasesthetimebytwoordersofmagnitudewhencomparedwithIHTandNormalizediterativehardthresholding(NHT)algorithmsforlownoiselevel.Simulationonthe0-1sparsesignaldemonstratesthatthereconstructionprobabilityofBIHTalgorithmishigherthanthatoftheIHTalgorithmifthemeasurementtimesandsparsityofthesignalarethesame.

KeywordsCompressivesensing,iterativehardthresholding(IHT),normalizediterativehardthresholding(NIHT),backtracking,sparsity

Nyquist-Shannon采样定理50多年来一直作为信号获取系统的基础.运用这个理论,信号以其带宽的两倍采样,而不用考虑可能是众所周知的额外信号结构.最近涌现的压缩传感理论(Compressedsensing,CS)[1],假设信号在某个变换域上(近)稀疏,运用线性变换(例如傅里叶变换、小波变换等)将信号投影到一个低维(比较于Nyquist比率)空间上,则少量随机的线性投影中包含了信号的大部分信息,通过非线性解码机制可以高概率精确重建原始信号.CS运用少量的随机投影代替以Nyquist比率采样,相比较于信号空间维数,稀疏限制大大

录用日期2010-08-20

ManuscriptreceivedMarch30,2010;acceptedAugust20,2010高等学校博士学科点专项科研基金(20070357003),“新一代宽带无线移动通信网”国家科技重大专项(2009ZX-03006-001-02)资助

SupportedbySpecializedResearchFundfortheDoctoralPro-gramofHigherEducation(20070357003)and“TheNextGen-erationBroadbandWirelessMobileCommunicationsNetwork”NationalScienceandTechnologyofMajorProjects(2009ZX-03006-001-02)

1.安徽大学计算智能与信号处理教育部重点实验室合肥2300392.合肥师范学院数学系合肥2300393.上海第二工业大学理学院上海201209

1.KeyLaboratoryofIntelligentComputingSignalProcess-ing,MinistryofEducation,AnhuiUniversity,Hefei2300392.DepartmentofMathematics,HefeiNormalUniversity,Hefei2300693.SchoolofScience,ShanghaiSecondPolytechnicUniversity,Shanghai201209收稿日期2010-03-30

减少了包含信号主要信息的空间尺度.Donoho[1]开创性工作的重要贡献之一是证明了只要采样矩阵Φ满足具有一个常数参数的限制等容条件(Restrictedisometryproperty,RIP),线性规划算法和贪婪算法等能高概率重建原始信号.压缩传感理论中,众多重建算法能有效地逼近最优解,例如传统的正交匹配追踪(Orthogonalmatchingpursuit,OMP)算法[2]、逐步正交匹配追踪(Stagewiseorthogonalmatchingpursuit,STOMP)[3]、子空间追踪(Sub-spacepursuit,SP)[4]以及有良好实现保证的迭代硬阈值(Iterativehardthresholding,IHT)[5]、正规化的迭代硬阈值(Normalizediterativehardthresh-olding,NIHT)[6]和压缩采样匹配追踪(Compres-sivesamplingmatchingpursuit,CoSaMP)算法[7]等.其中,逐步正交匹配追踪(STOMP)算法利用的是软阈值迭代,通过一个软阈值,也就是一个阈值参数确定迭代中保留的项,而硬阈值迭代通过一个固定的阈值确定的非线性算子进行迭代,每次迭代中保留的项数相等.软、硬阈值算子在压缩传感中的应用均已得到广泛研究,文献[8?9]中通过 1模确定了去噪的一个变分公式,从而诱导了一个简单的软阈值,由于其简单性引起了广泛的关注,文献[3]中STOMP算法运用软阈值对应的 1模求信号

的逼近解;而文献[5]中也确定了硬阈值与函数最小值间的关系,它通过满足向量非零元个数的约束,也就是 0模达到函数最小值的要求.国内关于阈值的理论及应用研究也有许多,例如文献[10]讨论了与阈值迭代下如何实现全变差,文献[11]讨论了CS在医学磁共振成像(Magneticresonanceimaging,MRI)当中如何实现软阈值迭代和算法的收敛性.本文针对IHT算法,通过在算法迭代过程中加入回溯思想[12],提出基于回溯的迭代硬阈值(Backtracking-basediterativehardthresholding,BIHT)算法,该算法在迭代过程中,将测量向量在迭代过程的前后两次支撑并集对应矩阵列张成的空间上重新做投影,再利用IHT算法中的非线性算子得到逼近解.该回溯过程优化了迭代过程中解向量支撑的选择.文献[6]中证明了NIHT算法在有IHT算法类似的实现保证基础上,算法稳定,并且收敛快于IHT算法.而本文的模拟实验和真实实验结果表明,BIHT算法在保证重建质量的前提下,迭代次数和时间明显少于IHT、NIHT算法.

1CS理论

本文中,x表示向量.supp(x)表示向量x=(x1,···,xN)的非零坐标,称为向量x的支撑.当|supp(x)|=K,K N时,称x是一个K-稀疏信号.通常,编码一个N维的K-稀疏离散信号,可通过计算一个包含向量x的M N个线性投影的测量向量y完成,由

y=Φx+e(1)

描述.e为噪声,这也是压缩传感有噪声情况下的模型,e=0为无噪声情况下的模型,其中Φ是M×N矩阵,因为不是一个方阵,故没有传统意义上的逆,所以从y=Φx中恢复信号x是一个病态求逆问题,由于信号的稀疏性,信号恢复可视为从所有可能的解中找出非零元最少的一种解形式,又根据 0模的概念, x 0表示向量x的非零元个数,即在一个 0模最小化内,寻找下列问题的一个解

min x 0

s.t.y=Φx(2)

不幸的是,解决上述 0最小化问题是一个NP-hard问题,因此非常困难.由此,Chen等在文献[13]中提出了基于最小化 1模的线性优化模型,也就是基追踪(Basispursuit,BP)方法:

min x 1

s.t.y=Φx(3)

不同于 0模, 1模具有凸性,并能很好地保持信号的稀疏性.模型(3)可以通过现有的优化算法获得高效求解.并且文献[13]中证明了如果支撑集

supp(x)满足条件:

|supp(x)|<>

2

(1+ν?1)

(4)

这里ν表示Φ与单位基的相干系数,那么式(3)的最小化解是唯一的且和式(2)的最小化解等价.同时,测量矩阵的限制等容条件(RIP)为从y=Φx中精确重建信号提供了理论保证.

定义1(限制等容条件).称一个测量矩阵Φ满足具有参数(K,δK),δK∈(0,1)的限制等容条件(RIP),如果对所有的K稀疏向量x,有

(1?δK) x 2≤ Φx 2≤(1+δK) x 2

(5)

成立.文献[4]中给出了截断的定义.

定义2(截断(Truncation)).设Φ∈RM×N,I?{1,···,N},矩阵ΦI表示Φ中指标i∈I对应的各列构成的矩阵,称为矩阵Φ的截断.通过ΦI的列张成的空间定义为张成空间(ΦI).

由上述定义,ΦI是一个秩为I的矩阵,由文献[14]中推论及定理8.2.2知,ΦI的伪逆唯一,且表示

为Φ?T?1I=(ΦIΦI)

ΦT

I.2迭代硬阈值算法

文献[5]为迭代硬阈值(IHT)算法用于压缩传感恢复问题提供了一系列的理论保证,证明算法成功运用较少的测量向量逼近最优的信号恢复,其迭代过程如下:假设x0=0,迭代

xn+1=HK(xn+ΦT(y?Φxn))

(6)

其中HK(α)是一个非线性算子,保留α值最大的K个分量,其余分量均设为零.如果按这样的机制所得非零支撑集合不唯一,则可随机选择其中一个集合或者预设要选择分量的支撑集合.在条件 Φ 2<>

min y?Φx2

x

2

s.t.

x 0≤K

(7)

即IHT算法对应解决的是被称为是NP-hard的 0模问题,直接考虑稀疏信号的非零元个数,寻找最能逼近稀疏信号的K项支撑.上述算法的更一般的形式是包含额外的步长尺度μ>0,即我们运用

xn+1=HK(xn+μΦT(y?Φxn))

(8)

进行迭代.并且,在文献[5]中,通过下列定理1和定理2保证了算法的收敛性以及实现.

定理1.如果测量矩阵Φ张成RM,以及 Φ 2

<1,则在x是k稀疏这个约束下,迭代硬阈值算法收敛于代价函数 y?φx="">

定理2.考虑一噪声测量y=Φx+e,其中x是一个任意向量.假设xK是xΦ满足具有参数δ√的最优K次逼近.如果3K<1>

x?xn n2≤2? xK 2+6?εK

(9)

其中

ε? x?x1

K=K 2+ x?xK 1+ e 2

而且,在最多

n?=log xK

2

2

ε?(10)

K

次迭代,IHT估计了x,具有精度

x?xn? 2≤7 x?xK 2+

1

K

x?xK 1+ e 2(11)更一般地,如果δ√

3K<1>

x?xn? 2≤c x?xK 2+

1

x?xK 1+ e 2(12)δ→1/√

3K时,c→∞.虽然IHT算法有良好的实现保证,但是由于限制等容条件以及迭代硬阈值对矩阵Φ的重新缩放是敏感的.当定理中的条件不满足,Φ的尺度变化时,IHT算法经常变得极不稳定,文献[6]中图1说明了这一点.从第5节的模拟实验中也可以看出,Φ取一般的随机高斯矩阵时,IHT算法对信号的重建概率甚至很低.Blumensath等在文献[6]中对此算法做了改进,在迭代中加入一个下降序列因子{μn},仍保证收敛于代价函数的一个局部最小值及保证算法实现,且取代IHT中通过限制 Φ 2以保证算法的收敛性,使得算法保持尺度独立,大大提高了算法的稳定性及对信号的重建概率,称之为正规化的迭代硬阈值(NIHT)算法,算法步骤见算法1.文献[6]中也给出了类似于定理1和定理2的实现保证,且相比较于IHT算法,NIHT算法虽然增加了下降序列因子{μn},但在迭代时间上不仅没有增加,在迭代次数较多时,迭代时间还少于IHT算法.本文基于IHT算法,提出了另一种改进,引入回溯思想,称之为基于回溯的迭代硬阈值

(BIHT)算法,虽然多了一步回溯的求伪逆过程,但是省略了步长序列的选择,计算复杂度和迭代次数及时间却远小于IHT、NIHT算法.

算法1.NIHT算法

1)初始化

迭代次数n=1,x1=0,Γ1=supp(HK(ΦTy))2)迭代

ggTn=ΦT(y?Φxn),μn=

ΓgΓn

gΓx

?n

ΦΓ

n

ΦΓngΓn

n+1=HK(xn+μngn),Γn+1=supp(x?n+1)若Γn+1=Γn,xn+1=x

?n+1若Γn+1=Γn

如果μ?n≤(1?c) xn?xn 2 Φ(?xn

?xn

) ,c为某个2

小的固定常数,则

xn+1=x

?n+1如果μn>(1?c) ?xn

?xn 22

Φ(?

x

n?xn) ,对某个

2

κ>1/(1?c),更新:

μn←μn/(κ(1?c)),x

?n+1=HK(xn+μngn)迭代直到μn≤(1?c) ?xn

?xn 22

nnΓ2

n+1=supp(x

?n+1),xn+1=x?n+13一种回溯的迭代硬阈值算法

受子空间追踪(SP)算法[4]的启发,将回溯的思想引入IHT算法,即BIHT算法的基本思想借用具有回溯(试探法)的顺序编码理论[12].在这个解码框架中,首先选择最可能张成编码空间的K个支撑集合,张成空间.如果得到的向量到这个空间的距离被认为很大,则算法根据它们的可靠值,递增地消除和增加新的基向量,直到确定一个充分靠近的候选支撑,但在增加和消除的过程中我们有一步回溯的思想,重新估计候选者的可靠性.BIHT算法使用一个类似的策略,但在每一个步骤,保留IHT算法的非线性算子原则,将相同数量的向量从候选名单中删除.因此,IHT、NIHT和BIHT重建策略之间的主要区别是前者根据一个非线性算子HK(α)依次产生一列候选者,没有回溯过程.在每一次迭代中,通过不断将信号残差投影到Φ列构成的原子库中最匹配它的K个向量,确定K个可信赖的候选,一旦坐标被认为是可信赖的,则被增加于名单,进行下一次迭代.这种迭代因为缺乏Φ列之间的正交性而不是最优的,向量α的某个支撑可能在之前的迭代中因为不在最大的K值之内已被设为零,但是在下一次的迭代中又因为加入残差投影,在最大的K个值之内而被重新保留下来,也就是说,α的某个支撑可能被反复选择和丢弃,所以导致的结果是IHT、NIHT算法需要较多的迭代次数才能高概率重建原始信号,增加了计算复杂度和重建时间.相比之下,BIHT算法在每一次迭代过程中增加了一步回溯的思想,将

前次迭代结果与非线性算子HK(α)产生的新向量支撑合并,再通过伪逆过程与非线性算子HK(α)重新估计.我们通过模拟实验发现,仅仅需要很少的迭代次数既能高概率恢复原始信号,又能大大减少精确重建所需的迭代时间,BIHT算法的主要步骤见算法2.相比较于IHT、NIHT算法,BIHT算法迭代过程是为了确定向量最大K个元的支撑集,然后利用伪逆过程将测量向量y在IHT算法中前后两次支撑并集对应矩阵列张成的空间中重新做一个投影,再利用非线性算子HK(α)选择最大的K项,得到信号的逼近解,所以下降因子对本算法没有影响,我们可以直接假设μ=1.这不同于IHT和NIHT算法,下降序列的选择影响了算法的稳定性.NIHT算法迭代次数可以选用具体的数值,当然也可以有别的停止条件.算法加入回溯的思想,从而优化了IHT、NIHT算法每次迭代支撑的选择,减少x中的某些支撑被反复迭代的情况.BIHT算法的主要贡献是在其每一次迭代中均取得重大的进步,所以经过少数几次的迭代就能达到良好的重建结果.

算法2.BIHT算法

1)初始化

x1=0,残差r0=y,迭代次数n=12)迭代

an+1=HK(xn+μΦTrn),Γn=supp(xn)∪supp(an+1)3)回溯

?n+1=Φ??n+1)xΓny,xn+1=HK(x

及BIHT算法对于真实图像重建的效果比较,结果

表明我们的算法在同样的迭代次数下,效果及运行速度均优于IHT和NIHT算法;其次通过对本身稀疏的0-1信号重建比较OMP、IHT、NIHT以及BIHT算法的重建概率与稀疏度及测量次数之间的关系关于信号稀疏度及所需测量数目之间的比较;最后考虑了有噪声情况下NIHT及BIHT算法的比较情况.

4.1真实图像重建

本节算法重建的真实图像是一个64像素×64像素大小的海螺黑白图(来自于计算机视觉图像库(http://sta?.science.uva.nl/?aloi),测量矩阵采

?,其用高斯随机矩阵,Φ∈RM×N:Φ(i,j)=1

ij

中?ij～N(0,1),从下面模拟实验中可知稀疏度、测量次数,模拟实验在一台CPU为IntelE5300(双核2.60GHz),内存为2.00GB的联想机上各运行100次,重建的平均结果如图1所示

.

4模拟实验

文献[5]中说明了测量矩阵Φ的敏感性,推导了对于不同的下降因子μ,Φ的要求也不同.原算法中,为了概念上的方便,固定IHT一般算法中的μ值为1.文献[6]中指出,相比较于IHT算法,NIHT算法不再通过限制 Φ 2保证算法收敛,对测量矩阵Φ的缩放不敏感,算法稳定的同时,保留了类似IHT算法的实现保证.但是我们也看到,NIHT算法增加了下降序列因子{μn},这是否会增加运算量和运算时间？本节做了三组实验,首先是IHT、NIHT以

表1

Table1

IHT

n10100200

t(s)133.411389.12531.1

e2.279E+003778.67696.82

t(s)132.781338.82245.0

图1IHT、NIHT、BIHT算法的重建结果((1)～(4)分别

为海螺原始图像、IHT、NIHT、BIHT算法经10次迭代的重建结果;(5)～(10)依次为IHT、NIHT、BIHT算法经

100次和200次迭代的重建结果)Fig.1

ImageconstructionsrespectivelybasedonIHT,

NIHT,andBIHTalgorithmswithdi?erentiterationtimes

(Originalimageandreconstructionresultswith10iterationtimesarefrom(1)to(4);reconstructionresultswith100,200iterationtimesarefrom(5)to(10).)

表1为对应于IHT、NIHT和BIHT算法分别经10,100,200次迭代各运行100次的重建结果的平均时间与误差数据.

算法时间与误差比较

Comparisonsoftimeanderrors

NIHT

e650.178.677E?0045.632E?005

t(s)13.70315.35917.525

BIHT

e3.590E?0123.140E?0122.740E?012

从图1和表1中看出,一方面,在算法运行时间上,虽然NIHT算法比IHT算法虽然多了一个下降序列因子,但是运行时间上却与之相当,而且在迭代次数增加时,运行时间少于IHT算法,我们的BIHT算法在任何迭代次数下的运行时间均少于前面两者,且在超过100次迭代时时间有两个数量级的下降;另一方面,在算法重建效果上,算法均取10次迭代时,IHT算法的效果较差,BIHT算法能以极小的误差重建海螺图像,正因为如此,我们取100次和200次迭代时,BIHT算法运行的时间和误差与10次迭代比较只是稍有增加,而NIHT算法达到100次以上时,能以较小的误差重建海螺图像,IHT算法由于其重建的不稳定性,重建效果却仍然不十分理想.由此可见,对于重建相同质量的海螺图像,采用O(NlnK)的测量次数,BIHT算法比IHT、NIHT算法在重建真实图像时的误差要小得多,同时在保证重建质量的情况下,所需的迭代次数相差几十倍,重建所需时间也降低了两个数量级.

4.2算法比较

为了比较IHT、NIHT、BIHT以及传统OMP算法关于稀疏度和测量次数的重建性能,我们采用本身稀疏的0～1随机信号,通过假设N=256,M=128,改变信号稀疏度,以及假设N=256,K=30,改变测量次数,分别来进行1000次模拟实验,所得结果如图2和图3所示

.

图2N=256,M=128时,OMP、IHT、NIHT、BIHT算法对应的重建概率的比较(X轴表示0～1信号的稀疏度,

Y轴表示1000次实验成功的百分比)Fig.2

Comparisonofrelativereconstructionprobabilities

basedonOMP,IHT,NIHTandBIHTalgorithms(X-axisrepresentsthesparsityof0～1signal.Y-axisrepresentsthepercentageofsuccessfulexperimentin1000times.)

从图2和图3可以看出,IHT算法极不稳定,导致了重建概率很低,NIHT,BIHT算法均优于传统的OMP算法,撇开信号的重建时间,BIHT算法在重建概率上也是高于IHT、NIHT算法.图2中,相

同测量次数的情况下,随着稀疏度的增加,BIHT算法的重建概率曲线下降速率慢于IHT、NIHT算法;同样地,图3中,保持稀疏度不变的情况下,随着测量次数的增加,BIHT算法在任意测量数量情况下的重建概率也高于NIHT算法

.

图3N=256,K=30时,OMP、IHT、NIHT、BIHT算法对应的重建概率的比较(X轴表示测量次数,Y轴表示

1000次实验成功的百分比)Fig.3

Comparisonofrelativereconstructionprobabilities

basedonOMP,IHT,NIHTandBIHTalgorithms(X-axisrepresentsthemeasurementtimes.Y-axisrepresentsthepercentageofsuccessfulexperimentin1000times.)

4.3有噪声测量下的NIHT、BIHT算法

由于IHT算法的不稳定性,本节主要讨论NIHT、BIHT算法对测量向量以及第4.1节中海螺灰度图加不同程度的噪声时图像的重建效果,本节选用的噪声为高斯噪声,σ2表示方差,SNR表示信噪比.因为从表1中看出,BIHT、NIHT算法100次以上迭代的重建误差和时间都在一个数量级上,故本节下面实验中的NIHT、BIHT算法中均取100次迭代.首先,采用无噪声海螺源图像,只对测量向量y加不同程度的高斯噪声,也就是式(1)的重建情况.为了便于比较,在下面的重建实验中,我们将测量值和高斯噪声均归一化,这不会影响重建结果,因为由重建式(1)知,归一化相当于对测量值y和信号x同时除以一个常数,我们将重建结果乘以常数,则仍重建出原来的信号.对不同程度噪声分别做100次模拟实验,所得的平均重建时间、重建误差及重建效果如表2和图4.

从表2和图4中看出,NIHT算法在有测量噪声情况下的重建时间与无噪声情况下相当,BIHT算法重建时间增加了一个数量级,但仍少于NIHT算法.重建效果方面,NIHT、BIHT算法在噪声方差较小时,重建误差均很小,NIHT算法略优于BIHT算法,但随着方差增大,重建误差增长较快,明显高于BIHT算法.其次,对图像加高斯噪声,经

过100次模拟试验,所得不同程度噪声所需的平均重建时间、重建误差及重建效果如表3和图

5.

图4IHT、NIHT、BIHT算法的重建结果(第1行分别为海螺原始图像、NIHT、BIHT算法在测量向量高斯噪声方差

为σ2

=0.05时对应的海螺重建图像;第2、3行依次为

NIHT、BIHT算法在σ2=0.10,σ2

=0.15以及σ2=0.20

时的重建结果)Fig.4

ImageconstructionsrespectivelybasedonIHT,NIHT,andBIHTalgorithmswiththedi?erentGaussian

noisesatmeasurements(OriginalimageandreconstructionresultswiththeGaussiannoisevarianceσ2

=0.05areinthe?rstrow;reconstructionresultswithσ2=0.10,σ2=0.15,andσ2=0.20areinthesecondand

thethird

rows.)

图5IHT、NIHT、BIHT算法的重建结果(第1行分别为

高斯噪声方差σ2

=5时的噪声海螺源图像及NIHT、BIHT算法对应的海螺重建图像;第2、3行依次为NIHT、BIHT

算法在σ2=10,σ2=15以及σ2

=20时的噪声海螺源图

像及NIHT、BIHT算法对应的海螺重建图像)Fig.5

ImageconstructionsrespectivelybasedonIHT,

NIHT,andBIHTalgorithmswiththedi?erentGaussiannoisesatimages(OriginalnoiseimageandreconstructionresultswiththeGaussiannoisevarianceσ2=5areinthe?rstrow;originalnoiseimagesandreconstructionresultswithσ2=10,σ2=15,andσ2=20areinthesecondand

thethirdrows.)

从表3和图5中看出,NIHT算法在有图像噪声情况下的重建时间与无噪声情况下相当,BIHT

算法重建时间增加了一个数量级,但仍少于NIHT

算法.重建效果方面,BIHT算法在方差较小时信噪比较高,明显高于BIHT算法.但随着方差增大,NIHT、BIHT算法对应信噪比均较差,BIHT算法的信噪比下降更快,希望下一步能在图像噪声变大时改善信噪比.

表2噪声测量下算法时间与误差

Table2Comparisonsoftimeanderrorswith

Gaussiannoiseatmeasurements

NIHT

BIHT

σ2

t(s)e(%)t(s)e(%)0.051374.51.8218325.232.64430.101336.23.5049330.144.56430.151418.823.898340.3610.1590.20

1441.2

40.717

346.42

28.738

表3噪声测量下算法时间与误差

Table3Comparisonsoftimeanderrorswith

Gaussiannoiseatimages

NIHT

BIHT

σ2t(s)SNRt(s)SNR51372.510.476297.5844.251101382.29.5033340.3815.795151426.14.3023355.253.014620

1455.4

2.3417

356.09

0.4494

5总结及展望

本文主要提出一种基于回溯的IHT算法,采用IHT的原子选择机制以及通过回溯的思想减少同一原子反复迭代的过程,从而减低了计算复杂度.从模拟实验中看出,相比较于IHT及NIHT算法,BIHT算法大大减少了迭代次数,从而使得重建时间也大大减少;同时,在相同的测量次数和信号稀疏度情况下,NIHT算法的重建概率也有所提高.压缩传感理论是当前信号处理的一个热点,国内也有许多此类研究,例如文献[15?16].而快速稳定的重建算法是将CS推向实用化的关键之一,也是CS的主要研究内容之一.目前关于CS的重建算法越来越多,我们也一直致力于进一步寻找更加实用有效的重建算法(例如文献[17]),在减少测量次数的同时还能降低对信号稀疏度的要求.

References

1DonohoDL.Compressedsensing.IEEETransactionson

InformationTheory,2006,52(4):1289?1306

2TroppJ,GilbertA.Signalrecoveryfromrandommeasure-mentsviaorthogonalmatchingpursuit.IEEETransactionsonInformationTheory,2007,53(12):4655?46663DonohoDL,TsaigY,DroriI,StarckJL.SparseSolutionofUnderdeterminedLinearEquationsbyStagewiseOrthogo-nalMatchingPursuit.TechnicalReportNo.2006-2,Depart-mentofStatistics,StanfordUniversity,USA,20064DaiW,MilenkovicO.Subspacepursuitforcompressivesensingsignalreconstruction.IEEETransactionsonInfor-mationTheory,2009,55(5):2230?22495BlumensathT,DaviesME.Iterativehardthresholdingforcompressedsensing.AppliedandComputationalHarmonicAnalysis,2009,27(3):265?2746BlumensathT,DaviesME.Normalizediterativehardthresholding:guaranteedstabilityandperformance.IEEEJournalofSelectedTopicsinSignalProcessing,2010,4(2):298?3097NeedellD,TroppJA.CoSaMP:iterativesignalrecoveryfromincompleteandinaccuratesamples.AppliedandCom-putationalHarmonicAnalysis,2009,26(3):301?3218ChambolleA,DeVoreRA,Nam-YongL,LucierBJ.Nonlinearwaveletimageprocessing:variationalproblems,compression,andnoiseremovalthroughwaveletshrink-age.IEEETransactionsonImageProcessing,1998,7(3):319?3359ZibulevskyM,EladM.L1-L2optimizationinsignalandimageprocessing.IEEESignalProcessingMagazine,2010,27(3):76?8810MaJW.Compressedsensingbyinversescalespaceand

curveletthresholding.AppliedMathematicsandComputa-tion,2008,206(2):980?98811QuXB,ZhangWR,GuoD,CaiCB,CaiSH,Chen

Z.IterativethresholdingcompressedsensingMRIbasedoncontourlettransform.InverseProblemsinScienceandEn-gineering,2010,18(6):737?75812HanYS,HartmannCRP,ChenCC.E?cientpriority-?rstsearchmaximum-likelihoodsoft-decisiondecodingoflinearblockcodes.IEEETransactionsonInformationThe-ory,1993,39(5):1514?152313ChenSS,DonohoDL,SaundersMA.Atomicdecomposi-tionbybasispursuit.SiamJournalonScienti?cComputing,1998,20(1):33?6114ShiRong-Chang,WeiFeng.MatrixAnalysis.Beijing:Bei-jingInstituteofTechnologyPress,2005.302?303

(史荣昌,魏丰.矩阵分析.北京:北京理工大学出版社.2005.302?303)

15ShiGuang-Ming,LiuDan-Hua,GaoDa-Hua,LiuZhe,Lin

Jie,WangLiang-Jun.Advanceintheoryandapplicationofcompressedsensing.ChineseJournalofElectronics,2009,37(5):1070?1081

(石光明,刘丹化,高大化,刘哲,林杰,王良君.压缩感知理论及其研究发展.电子学报,2009,37(5):1070?1081)16LiShu-Tao,WeiDan.Asurveyoncompressivesensing.Acta

AutomaticaSinica,2009,35(11):1369?1377

(李树涛,魏丹.压缩传感综述.自动化学报.2009,35(11):1369?1377)17FangHong,ZhangQuan-Bing,WeiSui.Imagereconstruc-tionbasedonimprovedbackwardoptimizedorthogonalmatchingpursuitalgorithm.JournalofSouthChinaUniver-sityofTechnology(NaturalScienceEdition),2008,36(8):23?27

(方红,章权兵,韦穗.改进的后退型最优正交匹配追踪图像重建方法.华南理工大学学报(自然科学版),2008,36(8):23?27)

杨海蓉安徽大学博士研究生.主要研究方向为信号处理.本文通信作者.E-mail:murong81@163.com(YANGHai-RongPh.D.candi-dateatAnhuiUniversity.Hermainresearchinterestissignalprocessing.Correspondingauthorofthispaper.)方红上海第二工业大学理学院副教授.主要研究方向为图像处理.E-mail:luckymars@gmail.com(FANGHongAssociateprofessorattheSchoolofScience,ShanghaiSec-ondPolytechnicUniversity.Hermainresearchinterestisimageprocessing.)张成安徽大学博士研究生.主要研究方向为信号处理.

E-mail:question1996@163.com

(ZHANGChengPh.D.candidateatAnhuiUniversity.Hismainresearchinterestissignalprocessing.)

韦穗安徽大学教授.主要研究方向为计算机视觉、图像处理和全息显示.E-mail:swei@ahu.edu.cn

(WEISuiProfessoratAnhuiUni-versity.Herresearchinterestcoverscomputervision,imageprocessing,andholographydisplay.)

转载请注明出处范文大全网 » 基于迭代方法的软阈值估计小波