D O I :10.3969/j
i s s n .1001-5337.2024.1.053收稿日期:2022-07-05
基金项目：山东省自然科学基金(Z R 2022MA 081
)；枣庄学院博士科研启动基金.通信作者：孙敏，男,1980-，博士，副教授；研究方向：最优化理论;E -m a i l :z i y
o u x i a o d o u @163.c o m.基于P R P 共轭梯度的图像去噪算法
孙 敏①, 孙肖丽①, 田茂英②
(①枣庄学院数学与统计学院;②山东煤炭卫生学校,277160,
山东省枣庄市)  摘要：两阶段算法是去除椒盐噪声的一类非常有效的方法.该文设计了一类P R P 共轭梯度法来求解两阶段法中的无约束最优化问题.该方法具有两个显著特点：搜索方向具有充分下降性；结合一类定制的A r m i -
j
o 型非精确线搜索，其点列的梯度趋于零，而不是通常共轭梯度法的下极限趋于零.将所设计方法及一些经典的共轭梯度法应用到两阶段算法的最优化问题中，数值结果表明该文设计的方法在收敛速度与计算精度方面都具有优势.
关键词:P R P 共轭梯度法;
两阶段算法；椒盐噪声；全局收敛性中图分类号:O 221.1  文献标识码:A  文章编号:1001-5337(2024)01-0053-08
0 引 言
通过视觉感知图像是人类认识与了解现实世界的重要手段.数字图像利用一个数字序列表示图像，便
于使用计算机对其进行各种图像处理.在数字化采样及传输过程中，数字图像容易受到各种噪声的污染，降低了图像的可读性.在追求高品质视觉体验的时代，数字图像的降噪处理引起了业界的广泛关注，该问题已成为数字图像领域的一个重要研究课题.目前学术界提出了多种降低噪声的方法，包括利用线性拉普拉斯
滤波器的模糊图像处理；利用二维中值滤波器的椒盐噪声图像处理[1
]
；利用一阶原始-对偶方法的运动模糊图像处理[2]；利用共轭梯度法的模糊图像处理[3]；利用两阶段法的椒盐噪声图像处理[4]
；利用交替方向法处理含噪声及缺失像素的视频背景与前景分离问题[5]
等.
数值实验表明，即使椒盐噪声浓度达到90%，两阶段法也能得到比较高的峰值信噪比(P S N R /d B )
，这说明该方法修复的图像比较好地还原了原始图像.在第一阶段该方法利用自适应中值滤波器(AM F )出受椒盐噪声污染的像素坐标.在第二阶段该方法通过求解一个无约束最优化问题填补椒盐噪声处的像素值.通过先检测再填充的方式，该方法有效地提高了数字图像去噪处理的速度与质量.
由上面的分析可以看出，两阶段法的关键在于第二阶段的无约束最优化问题的求解.共轭梯度法是求
解无约束最优化问题的最优秀的算法之一.众多的共轭梯度法中,P R P 共轭梯度法又是数值表现最好的共
轭梯度法之一.P R P 共轭梯度法在理论上存在的一个缺点是其收敛性条件比较苛刻.对于一般的目标函数，即使采用精确线搜索，其也可能发散[6].本文设计了一类P R P 共轭梯度法,
该方法在较弱的条件下具有更强的全局收敛性，即其序列梯度趋于零，而不是通常共轭梯度法的序列梯度下极限趋于零.同时将该方法应用到两阶段法的无约束最优化问题.
1 修正P R P 共轭梯度法
考虑光滑无约束最优化问题
m i n f (
x ):x ɪℝn ,(1
) 第50卷 第1期2024年1月      曲阜师范大学学报J o u r n a l  o f  Q u f u  N o r m a l  U n i v e r s i t y
V o l .50 N o .1
J a n .2024
其中f (x ):ℝn ңℝ是可微函数，其梯度记为Ñf (
x )，简记为g (x )，进一步g (x k )简记为g k .由于计算量与存储量的优势，共轭梯度法是求解式(1)
最有效的数值方法之一，尤其是求解大规模无约束最优化问题.共轭梯度法属于一阶迭代算法，给定当前迭代点x k ɪℝn ,
其迭代格式为x k +1=x k +αk
d k ,k =0,1,2, ，其中αk ɪℝ是步长因子,
通常由非精确线搜索确定;d k ɪℝn 是搜索方向，由下式确定,d k =
-g k ,
k =0,-g k +βk
d k -1,k ȡ1,{
(2
)这里βk ɪℝ是区分不同共轭梯度法的参数，也就是说不同共轭梯度法的区别主要在于βk 的不同表达式.βk 经典的表达式包括
βF R
k
=g T k g k g T k -1g k -1,βC D k =g T k g k -g T k -1(g k -g k -1),βD Y
k =g T k g k d T
k -1(g k -g k -1),βP R P
k =g T k (g k -g k -1)g T k -1g k -1,βH S k =g T k (g k -g k -1)d T k -1(g k -g k -1)
,βL S k =g T
k (g k -g k -1)-g T
k -1d k -1.  在求解一般的无约束最优化问题时,
不同的共轭梯度法有着截然不同的理论结果与数值表现.一般而言，由βF R k ㊁βC D k ㊁βD Y k 确定的共轭梯度法与βP R P k ㊁βH S k ㊁β
L S
k 确定的共轭梯度法在理论分析与数值效果方面具有互补的表现，即由βF R k ㊁βC D k ㊁βD Y k 确定的共轭梯度法具有较强的收敛性，但数值表现差一些，而由βP R P k ㊁βH S k ㊁β
L S
k 确定的共轭梯度法具有更好的数值表现，但收敛性条件更强.基于共轭梯度法最近的一些研究成果，本文将设计一个修正的P R P 共轭参数,
使得由其确定的共轭梯度法在通常的条件下具有全局收敛性.下面给出该方法的设计思路.
首先，搜索方向d k 满足下降性g T
k d k <0，这是迭代下降类算法的基本要求.再进一步，如果d k 满足充分
下降性g T k d k ɤ-c  g k  2
(c >0)，则基于其设计的算法往往具有更好的理论分析结果与数值表现.为了产生具有充分下降性的搜索方向，有3位学者做出了突出的贡献:S h i 等[7]
设计了一个具有凸组合结构的下降方向，利用范数 ㊃ 2的凸性，该方向具有充分下降性；后来,H a g
e r 等[8]设计了一类共轭梯度法，其中共轭参数取为
β
H Z
k
=1
<d k -1,y k -1> y k -1-2d k -1 y k -1 2<d k -1,y k -1>
,g k  ，其中y k -1=g k -g k -1.
利用柯西-施瓦兹不等式可以证明该方法产生的方向满足g T
k d k ɤ-78正则化共轭梯度法
g k  2
.
如果采用精确线搜索，则有g T k d k -1=0，于是该方法退化为经典的H S 共轭梯度法，因此其是对H S 共轭梯度法的
改进；几乎与此同时,Z h a n g 等利用正交化的思想,
给出了一个修正的三项P R P 共轭梯度法，该方法产生的搜索方向满足g T k d k =- g k  2.
基于这三篇文献对共轭梯度法做出的开创性贡献，最近十多年学者们设计了许多修正的共轭梯度法.本文受文献[8]中方法的启发，设计了一个具有充分下降性的共轭梯度法，其共轭参数取为
β
N P R P
k =g T
k y
k -
1-μ
y k -1 2 g k -1 2
g T k d k -1 g k -1
2,(3
)其中μ>
14
.如果采用精确线搜索，则由式(3)可以得到新的βN P R P k 退化为经典的βP R P k
.将式(3)代入式(2)得到搜索方向
d k =
-g k ,
k =0,-g k +βN P R P
k
d k -1,k ȡ1.{
(4)该方法继承了P R P 共轭梯度法自动再启动的优点，即当迭代点前进量几乎为零时，该方法自动切换到最速下降法.同时其还具有P R P 共轭梯度法所不具备的优点,
即充分下降性.其次，为了得到收敛性分析中的一些关键不等式，修正P R P 共轭梯度法采用如下A r m i j
o 型线搜索：步45            曲阜师范大学学报(自然科学版)              2024年
长αk 是集合{υk ,υk ρ,υk ρ2,
}中满足下面不等式的最大的α[7]
:f (
x k +αd k )ɤf (x k )-σα2
d k  4,(5
)其中υk =-δg T
k
d k  d k  2
,δ>0,ρɪ(0,1),σ>0.容易看出，初始步长υk 的选取主要是受二次函数精确步长的启发.
最后给出修正P R P 共轭梯度法的算法流程.步0 选取常数μ>
1
4
,δ>0,ρɪ(0,1),σ>0.选取初始迭代点x 0ɪℝn ，计算g 0.令k :=0.步1 如果 g k  =0,
停止，否则转步2.步2 由A r m i j o 线搜索式(5)确定步长αk .
步3 由式(4)确定搜索方向d k .步4 令x k +1:=x k +αk d k ，计算g k +1.步5 令k :=k +1,
返回步1.2 修正P R P 共轭梯度法收敛性分析
为了分析修正P R P 共轭梯度法的收敛性，需要对式(1
)
施加如下假设:(H 1)水平集Ω={x ɪℝn |f (
x )ɤf (x 0)}有界;(H 2)在Ω的某个邻域中，梯度g (x )是L i p s c h i t z 连续的，即存在一个常数L >0，使得 g (x )-g (y ) ɤL  x -y  ,∀x ,y ɪℝn
.(6)  下面的引理说明由式(5
)生成的方向满足充分下降性.引理2.1 对于任意的k =0,1,2, ，由式(5)定义的搜索方向d k 满足
g T k d k ɤ-c  g k  2
,
(7
)其中c =1-1
4μ
.
证明 如果k =0，则有g T k d k =- g k  2
;
如果k ȡ1，则有g T k d k =- g k  2+βN P R P k g T k d k -1=- g k  2+
g T
k y k -
1-μ
y k -1 2 g k -1 2
g T k d k -1 g k -1
2
g T
k d k -1=- g k  2
+1 g k -1 2g T k y k -1()g T k d k -1()-μ y k -1 2 g k -1
2
g T k d k -1()2æèçöø÷=- g k  2
+1
g k -1 2 g k -1 2g T k y k -1()24μ y k -1 2
-μ y k -1 g T k d k -1 g k -1 - g k -1 g T k y k -12μ y k -1 æèçöø÷2æèçöø÷ɤ- g k  2
+1
g k -1 2 g k -1 2g T
k y k -1()24μ y k -1 2
ɤ-1-14μæèçöø
÷ g k  2，其中最后一个不等式使用了柯西-施瓦兹不等式.证毕.再由柯西-施瓦兹不等式与式(7)得c  g k  ɤ d k  ,k =0,1,2, ,(8)由此得
αk ɤ-δg T
k d k  d k  2ɤδ g k  d k  ɤδ
c .(9
)  下面分析如果x k 不是问题(
1)的稳定点，则A r m i j o 型线搜索式(5)会在有限步内终止.55第1期          孙敏，等：基于P R P 共轭梯度的图像去噪算法
引理2.2 当 g k  ʂ0时,
A r m i j o 型线搜索式(5)会在有限步内终止.证明 利用反证法.假设结论不成立，则对任意正整数j ，有
f (x k +υk ρj d k )>f (x k )-συ2k ρ2j
d k  4.由f (x )的光滑性及中值定理得，存在θj ɪ(0,1)，使得ρj g (x k +θj υk ρj d k )T d k >-συk ρ2j
d k  4，即g (x k +θj υk ρj d k )T d k >-συk ρj
d k  4.在此式两边令j ңɕ，由ρɪ(0,1)及g (x )的连续性得g T k d k ȡ0.结合式(7)得-c  g k  2
ȡ0.这说明 g k  =0.
与已知条件矛盾.因此假设不成立.证毕.引理2.3 假设(H 1)与(H 2)成立，则由式(5)确定的步长αk 满足
αk ȡm i n υk ,
ρ
c  g k  2L +σ
d k  2() d k  2
{}
.(10
)  证明 如果αk <υk ，则由式(5)的规则,α'k =αk /ρ不满足式(
5)式，即f (
x k +α'k d k )>f (x k )-σ(α'k )2
d k  4.由中值定理，存在θk ɪ(0,1)使得g (x k +θk α'k d k )T
d k >-σα'k
d k  4.于是-σα'k  d k  4<g (x k +θk α'k d k )-g k ()T d k +g T
k
d k ɤL  d k  2α'k +g T k d k ɤL  d k  2α'k -c  g k  2
,
整理得αk =ρα'k >ρc  g k  2
L +σ d k  2() d k  2
.此式结合αk =υk 的情况可得式(10)成立.证毕.下面的定理表明修正P R P 共轭梯度法具有全局收敛性.
定理2.1 如果修正P R P 共轭梯度法产生无穷迭代点列{x k }，则l i m k ңɕ
g k  =0.
证明 由引理2.1得σα2
k
d k  4
ɤf (x k )-f (x k +1).由(H 1)
得正项级数ðɕ
k =0
f (
x k
)-f (
x k +1)()的部分和是有界的.再结合正项级数的比较审敛法得级数σðɕ
k =0
α2
k  d k  4是收敛的,
这说明其通项满足l i m k ңɕ
αk
d k  2
=0.结合式(9)得l i m k ңɕ
αk
d k  ()2ɤδ
c
l i m k ңɕαk  d k  2=0.(11)因此再结合式(8
)得l i m k ңɕ
αk  g k  =0.
(12
)  下面利用反证法证明结论成立.假设结论不成立，则存在常数ε0>
0及无穷集合K ⊆N ，使得对于任意的k ɪK ，都有
g k -1 ȡε
0.(13
)由(H 2)及αk ɤδc 得 g k -g k -1 2ɤL
α2
k  d k -1 2ɤL δc
αk  d k -1 2ң0(k ңɕ).由上式及式(13)得，存在ε1>0，使得对于任意的k ɪK ，都有 g k  ȡε
1.(14)于是结合式(12
)得l i m k ңɕ,k ɪK
αk =
0.(15
)  另一方面，由搜索方向的定义式(4
)得 d k  ɤ g k  +β
N P R P
k  d k -1 ɤ g k  + y k -1 +μ
y k -1 2 g k -1
2 d k -1  g k -1
2
g k  d k -1 ɤ g k  +1ε20L αk -1 d k -1 2
+μL 2ε20α2k -1 d k -1 4æè
çöø÷ g k  .
(16
)由式(12)得，存在r >0，使得任意的k ɪK ，有αk  g k  ɤ
r .将此式代入式(16)得 d k  ɤ r  g k  ,
∀k ɪK ,65            曲阜师范大学学报(自然科学版)              2024年
其中 r =1+L r 2ε2
0+μL 2r 4
ε4
.将上式代入式(10)，再结合式(8)得αk ȡm i n -δg T
k d k  d k  2,ρc  g k  2L +σr 2 g k  2() r 2 g k
2{}ȡm i n δc  r 2,ρc L +σr 2ε21() r 2{}
>0,∀k ɪK .显然上式与式(15
)是矛盾的.这说明假设不成立.证毕.3 数值实验
下面将本文设计的P R P 共轭梯度法(记为N P R P C G )应用到两阶段图像处理模型中，并与F R 共轭梯度
法(记为F R C G )以及文献[11]中的P R P 共轭梯度法(记为C P R P C G )
进行数值对比.首先简单回顾一下两阶段法[4
].令X =[x i ,j ]
m ˑn 表示一个数字图像，共m ˑn 个像素.对于每个(i ,j )ɪA ={1,2, ,m }ˑ{1,2, ,n }，令V i ,j 表示(
i ,j )的4个最相邻点构成的集合，即V i ,j ={
(i ,j -1),(i ,j +1),(i -1,j ),(i +1,j )}.对于点(i ,j )，令y i ,j 表示该点的正是像素值x i ,j 经过椒盐噪声污染后的观测值.在第一阶段，两阶段法利用自适应中值滤波器检测椒盐噪声的位置.令N ⊆A 表示检测到的可能含有椒盐噪声的像素点构成的集合;
第二阶段，两阶段法通过求解下面的最优化问题填充椒盐噪声点处的像素值:
m i n G α(u )=
ð(i ,j )ɪN ð(m ,n )ɪV i ,j ɘN φα(u i ,j -y m ,n )+12ð(m ,n )ɪV i ,j ɘN
φα(u i ,j -y m ,n )éëêêù
ûúú，其中α是正则化参数,φ
α是保边函数,u =[u i ,j ](i ,j )ɪN 是一个列向量(按字典序排列).本实验中φα定义为H u b e r 函数[1
1]
:φα(t )=t 2/(2α),  如果t ɤα,
t -α/2, 如果t >α,
{
其中α=10.使用峰值信噪比(P e a ks i g
n a l -t o -n o i s e r a t i o ,P S N R )来衡量图像修复的质量,P S N R =10l g
255
2
1MN ði ,j
u *i ,j -x i ,j ()2
,
其中u *
i ,j 表示点(
i ,j )处经两阶段法修复后的像素值.显然,P S N R 越大表明修复的图像越接近真实的图像.测试算法的终止条件设置为
F α(u k )-F α(u k -1)F α(u k )ɤ10-4或u k -u k -1u k
ɤ
10-4.选取初始迭代点x 0取为图像观测值.本次实验采用如下测试图片L e n a (256ˑ256)㊁C a m e r a m a n (256ˑ
256)㊁B a r b a r a (512ˑ512)与B o a t (512ˑ512
)，见图1.图1 原图片L e n a ㊁C a m e r a m a n ㊁B a r b a r a 与B o a t
首先对N P R P C G 的参数μ做灵敏度分析.
理论上μ>1
4都可以保证算法的收敛性，因此取μ满足μɪ{
0.25,0.35,0.45, ,2.05}.同时与文献[11]类似，对该算法采用类似于MA T L A B 的t r y -c a t c h 的运行思想：首先取近似最优步长αk =
75第1期          孙敏，等：基于P R P 共轭梯度的图像去噪算法

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。