ISSN 
-
图像恢复是图像处理中最基本的问题之一，其目的是从得到的噪声图像或者噪声模糊图像中恢复出原图像，同时避免破坏图像的一些基本信息。目前大量的文献已对该问题进行了深入的研究，但主要以加性高斯白噪声为主。然而，在光量子计数成像系统中，如天文成像、医学成像、荧光共焦显微成像等，得到的图像往往受到量子噪声的污染。量子噪声服从泊松分布的统计法则，并非加性噪声，并且泊松噪声的噪声强度与方差具有信号依赖性。统计上，亮度大的像素受到更多的干扰，因此去除该类泊松噪声是一个困难的科学问题。
图像去噪方法中最著名的是文献[1]基于变分理论提出的全变分 (total variation, TV) 去噪模型。TV模型在去除加性噪声的同时能很好地保持图像的边缘，但是会产生阶梯效应。文献[2]针对泊松噪声恢复问题，提出了新的全变分模型来处理泊松噪声图像，新变分模型的保真项适合泊松噪声，该模型同样也会产生阶梯效应，此后学者们又进行了进一步的研究[3-4]。
众所周知，TGV正则项[5-6]可以有效地消除阶梯效应，但同时也会使图像的边缘以及纹理区域变得模糊。更严重的，图像的一些细节信息在恢复过程中会丢失。Shearlet变换作为多尺度几何分析中的一种，能够对图像进行稀疏表示产生最优逼近，并在图像去噪过程中可以有效地保持图像边缘、角落等信息结合多尺度方法提取图像的几何特征。上述正则化方法在文献[7]中首次提出，并应用于压缩感知中图像重建问题。本文把TGV和Shearlet变换相结合来处理乘性泊松噪声恢复问题。
-
定义1[5] 设 $\mathit{\Omega } \subset {R^d}$ 是一个开区域， $k \ge 1$ 且 $\alpha = ({\alpha _0},{\alpha _1}, \cdots ,{\alpha _{k - 1}}) > 0$ 为加权系数，则对任意的 $\mathit{\boldsymbol{u}} \in L_{{\rm{loc}}}^1(\mathit{\Omega })$ ，k阶TGV定义为：
$$\begin{array}{l} {\rm{TGV}}_\alpha ^k(\mathit{\boldsymbol{u}}) = \sup \left\{ {\int_\mathit{\Omega } {\mathit{\boldsymbol{u}}{\rm{di}}{{\rm{v}}^k}\nu {\rm{d}}x|\nu \in } } \right.C_c^k(\mathit{\Omega },{\rm{Sy}}{{\rm{m}}^k}({R^d})),\\ \quad \quad \quad \left. {{{\left\| {{\rm{di}}{{\rm{v}}^l}\nu } \right\|}_\infty } \le {\alpha _l},l = 0,1, \cdots ,k - 1} \right\}{\kern 1pt} \end{array}$$ (1) 式中，Symk(Rd) 为k阶对称张量空间。
通过运用Legendre-Fenchel对偶定理，式 (1) 也可转换成它的原始对偶形式：
$$\rm{TGV}_{\mathit{\alpha }}^{\mathit{k}}\left( \mathit{\boldsymbol{u}} \right)=\underset{\begin{smallmatrix} {{\mathit{\boldsymbol{u}}}_{l}}\in {{C}^{k-1}}\left( \mathit{\bar{\Omega },}\rm{Sy}{{\rm{m}}^{\mathit{l}}}\left( {{\mathit{R}}^{\mathit{d}}} \right) \right) \\ l=1,2,\cdots ,k-1,{{u}_{0}}=\mathit{\boldsymbol{u}},{{u}_{k}}=0 \end{smallmatrix}}{\mathop{\inf }}\,\sum\limits_{\mathit{l}=1}^{\mathit{k}}{{{\alpha }_{\mathit{k}-\mathit{l}}}}{{\left\| \varepsilon \left( {{\mathit{\boldsymbol{u}}}_{\mathit{l}-1}} \right)-{{\mathit{\boldsymbol{u}}}_{\mathit{l}}} \right\|}_{1}}$$ 式中， $\varepsilon ({\mathit{\boldsymbol{u}}_{l - 1}})$ 表示对称梯度算子：
$$\varepsilon ({\mathit{\boldsymbol{u}}_{l - 1}}) = \frac{{\nabla {\mathit{\boldsymbol{u}}_{l - 1}} + {{(\nabla {\mathit{\boldsymbol{u}}_{l - 1}})}^{\rm{T}}}}}{2}$$ 本文主要使用二阶的TGV，在文献[6]中二阶 ${\rm{TGV}}_\alpha ^2$ 的能量函数定义为：
$${\rm{TGV}}_\alpha ^2{\rm{ = }}\mathop {\min }\limits_{\mathit{\boldsymbol{u}} \in {\rm{BGV}}_\alpha ^2(\mathit{\Omega }),\mathit{\boldsymbol{p}} \in {\rm{BD}}(\mathit{\Omega })} {\alpha _1}\int_\mathit{\Omega } {\left| {\nabla \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{p}}} \right|} + {\alpha _0}\int_\mathit{\Omega } {\left| {\varepsilon (\mathit{\boldsymbol{p}})} \right|} $$ 式中，；D1、D2分别是x和y方向上的向前有限差分算子矩阵，详见文献[7]，有：
$$\mathit{\boldsymbol{\varepsilon }}(\mathit{\boldsymbol{p}}) = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{D}}_1}{p_1}} & {\frac{1}{2}({\mathit{\boldsymbol{D}}_2}{p_1} + \mathit{\boldsymbol{D}}{}_1{p_2})}\\ {\frac{1}{2}({\mathit{\boldsymbol{D}}_2}{p_1} + \mathit{\boldsymbol{D}}{}_1{p_2})} & {{\mathit{\boldsymbol{D}}_2}{p_2}} \end{array}} \right]$$ -
Shearlet变换[7-9]最初是文献[9]根据小波理论衍生而来的。连续Shearlet变换[7]定义为：
$${\rm{S}}{{\rm{H}}_\mathit{\psi }}(f)(a,s,t) = ＜ f,{\mathit{\psi }_{a,s,t}} > $$ 式中， ${\mathit{\psi }_{a,s,t}}$ 是Shearlet基函数； $\mathit{\psi } \in {L_2}({R^2})$ 定义为：
$$\begin{array}{l} {\mathit{\psi }_{a,s,t}}(x) = {a^{ - 3/4}}\mathit{\psi }(\mathit{\boldsymbol{A}}_a^{ - 1}\mathit{\boldsymbol{S}}_s^{ - 1}(x - t)){\rm{ = }}\\ \quad {a^{ - 3/4}}\mathit{\psi }\left( {\left( {\begin{array}{*{20}{c}} {1/a} & { - s/a}\\ 0 & {1/\sqrt a } \end{array}} \right)(x - t)} \right) \end{array}$$ 式中，Aa为各向异性膨胀矩阵；Ss是剪切矩阵:
$${\mathit{\boldsymbol{A}}_a} = \left( {\begin{array}{*{20}{c}} a & 0\\ 0 & {\sqrt a } \end{array}} \right),a \in {R^ + },{\mathit{\boldsymbol{S}}_s} = \left( {\begin{array}{*{20}{c}} 1 & s\\ 0 & 1 \end{array}} \right),{\kern 1pt} {\kern 1pt} s \in R$$ 在后面试验中，采用文献[8]中的方法进行相关计算。
-
设 $\mathit{\boldsymbol{u}}\in R_{+}^{N}$ 是原始图像， $\mathit{\boldsymbol{f}}\in {{R}^{N}}$ 是观察图像， $\mathit{\boldsymbol{K}}\in {{R}^{N\times N}}$ 是线性模糊算子，依据泊松分布的定义假定 $\mathit{\boldsymbol{u}}$ 是有界并且是正的，则该退化模型描述为：
$$\mathit{\boldsymbol{f}}=P(\mathit{\boldsymbol{Ku}})$$ 式中， $P$ 表示泊松分布。基于文献[2]，运用贝叶斯法则，有：
$$P(\mathit{\boldsymbol{u}}|\mathit{\boldsymbol{f}})=\frac{P(\mathit{\boldsymbol{f}}|\mathit{\boldsymbol{u}})P(\mathit{\boldsymbol{u}})}{P(\mathit{\boldsymbol{f}})}$$ 对于任意的 $\mathit{\boldsymbol{u}}\in \mathit{\Omega }$ ，可以得到：
$$P(\mathit{\boldsymbol{f}}|\mathit{\boldsymbol{Ku}})=\prod\limits_{i=1}^{N}{P({{f}_{i}}|{{(\mathit{\boldsymbol{Ku}})}_{i}})}=\prod\limits_{i=1}^{N}{\frac{{{\rm{e}}^{-{{(\mathit{\boldsymbol{Ku}})}_{i}}}}{{({{(\mathit{\boldsymbol{Ku}})}_{i}})}^{{{f}_{i}}}}}{{{f}_{i}}!}}$$ 这里假设先验分布 $P(\mathit{\boldsymbol{u}})$ 是TGV和Shearlet变换，那么得到的新模型的正则项为：
$$P(\mathit{\boldsymbol{u}}) = exp\left( { - \lambda \sum\limits_{j = 1}^M {{{\left\| {{\rm{S}}{{\rm{H}}_j}(\mathit{\boldsymbol{u}})} \right\|}_1}} - {\rm{TGV}}_\alpha ^2} \right)$$ 式中， $\lambda $ 是正则化参数。这样本文的泊松复原模型可表示为：
$$\begin{array}{l} \mathop {\min }\limits_\mathit{\boldsymbol{u}} \beta \int_\mathit{\Omega } {(\mathit{\boldsymbol{Ku}} - \mathit{\boldsymbol{f}}\log \mathit{\boldsymbol{Ku}}){\rm{d}}x + } \\ \lambda \sum\limits_{j = 1}^M {{{\left\| {{\rm{S}}{{\rm{H}}_j}(\mathit{\boldsymbol{u}})} \right\|}_1}} {\kern 1pt} + {\rm{TGV}}_\alpha ^2 + {l_{{R^ + }}}(\mathit{\boldsymbol{u}}) \end{array}$$ 式中， ${\rm{S}}{{\rm{H}}_j}(\mathit{\boldsymbol{u}}) \in {R^N}$ 表示第j变换子带 (subband)，总的子带数目M由分解尺度的个数决定，试验中选择的分解尺度为3；ls表示投影算子保证图像中每个元素都为正，它的形式表示为：
$$\begin{array}{l}\quad \quad {l_S}(\mathit{\boldsymbol{u}}) = \left\{ {\begin{array}{*{20}{c}} {0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{u}} \in S}\\ { + \infty {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{u}} \notin S} \end{array}} \right.\\ {\mathit{\boldsymbol{u}}_ + } = \max \left\{ {0,\mathit{\boldsymbol{u}}} \right\},\log (0) = - \infty ,{\kern 1pt} {\kern 1pt} {\kern 1pt} 0\log (0) = 0 \end{array}$$ 进一步极小化离散TGV，模型可表示为：
$$\begin{array}{l} \mathop {\min }\limits_\mathit{\boldsymbol{u}} \beta \int_\mathit{\Omega } {(\mathit{\boldsymbol{Ku}} - \mathit{\boldsymbol{f}}\log \mathit{\boldsymbol{Ku}}){\rm{d}}x} + \lambda \sum\limits_{j = 1}^M {{{\left\| {{\rm{S}}{{\rm{H}}_j}(\mathit{\boldsymbol{u}})} \right\|}_1}} + \\ \quad \quad {\alpha _1}{\left\| {\nabla \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{p}}} \right\|_2} + {\alpha _0}{\left\| {\varepsilon (\mathit{\boldsymbol{p}})} \right\|_F} + {l_{{R^ + }}}(\mathit{\boldsymbol{u}}) \end{array}$$ (2) -
交替方向乘子法 (ADMM)[10, 11]在求解约束的凸优化问题时有极快的收敛性，因此被广泛的应用于图像处理。下面讨论如何应用ADMM算法求解新模型。
首先，引入辅助变量x、y、z、w、v，问题式 (2) 转化为以下约束问题：
$$\begin{array}{l} \quad \quad \quad \quad \mathop {\min }\limits_\mathit{\boldsymbol{u}} \beta \int_\mathit{\Omega } {\mathit{\boldsymbol{w}} - \mathit{\boldsymbol{f}}\log \mathit{\boldsymbol{w}}} + \\ \quad \quad \lambda \sum\limits_{j = 1}^M {{{\left\| {{\mathit{\boldsymbol{x}}_j}} \right\|}_1}} + {\alpha _1}{\left\| \mathit{\boldsymbol{y}} \right\|_2} + {\alpha _0}{\left\| \mathit{\boldsymbol{z}} \right\|_F} + {l_{{R^ + }}}(\mathit{\boldsymbol{v}})\\ {\rm{s}}{\rm{.t}}{\rm{.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{Ku}} = w,{\rm{S}}{{\rm{H}}_j}(\mathit{\boldsymbol{u}}) = {x_j},\nabla \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{p}} = y,\varepsilon (\mathit{\boldsymbol{p}}) = \mathit{\boldsymbol{z}},\mathit{\boldsymbol{v}} = \mathit{\boldsymbol{u}} \end{array}$$ (3) 约束条件可以进一步转化为：
$$\begin{array}{l} {\mathit{\boldsymbol{A}}_1} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}} & 0\\ {{\rm{S}}{{\rm{H}}_j}} & 0\\ D & { - \mathit{\boldsymbol{I}}}\\ 0 & \varepsilon \\ \mathit{\boldsymbol{I}} & 0 \end{array}} \right],{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{A}}_2} = \left[ {\begin{array}{*{20}{c}} { - \mathit{\boldsymbol{I}}} & {} & {} & {} & {}\\ {} & { - \mathit{\boldsymbol{I}}} & {} & {} & {}\\ {} & {} & { - \mathit{\boldsymbol{I}}} & {} & {}\\ {} & {} & {} & { - \mathit{\boldsymbol{I}}} & {}\\ {} & {} & {} & {} & { - \mathit{\boldsymbol{I}}} \end{array}} \right]\\ \quad \quad \quad \quad {\mathit{\boldsymbol{x}}_1} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{u}}\\ \mathit{\boldsymbol{p}} \end{array}} \right],{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{x}}_2} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{w}}\\ \mathit{\boldsymbol{x}}\\ \mathit{\boldsymbol{y}}\\ \mathit{\boldsymbol{z}}\\ \mathit{\boldsymbol{v}} \end{array}} \right],{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{d}} = \left[ {\begin{array}{*{20}{c}} 0\\ 0\\ 0\\ 0\\ 0 \end{array}} \right] \end{array}$$ 从以上形式，可以清楚地看出问题式 (3) 满足ADMM算法条件，并且收敛性也可以得到保证[12]。
对约束优化问题式 (3)，它的增广拉格朗日乘子函数为：
$$\begin{array}{l} \quad \quad \quad \quad L(\mathit{\boldsymbol{w,x,y,z,v,}}{c_i}) = \\ \quad \quad \quad \beta \int_\mathit{\Omega } {\mathit{\boldsymbol{w}} - \mathit{\boldsymbol{f}}\log \mathit{\boldsymbol{w}}} + \lambda \sum\limits_{j = 1}^M {{{\left\| {{\mathit{\boldsymbol{x}}_j}} \right\|}_1}} + \\ \quad \quad {\alpha _1}{\left\| \mathit{\boldsymbol{y}} \right\|_2} + {\alpha _0}{\left\| \mathit{\boldsymbol{z}} \right\|_F} + {l_{{R^ + }}}(\mathit{\boldsymbol{v}}){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} + \frac{{{\mu _1}}}{2}\left\| {\mathit{\boldsymbol{Ku}} - \mathit{\boldsymbol{w}} - {c_1}} \right\|_2^2{\kern 1pt} + \\ \frac{{{\mu _2}}}{2}\sum\limits_{j = 1}^M {\left\| {{\rm{S}}{{\rm{H}}_j}(\mathit{\boldsymbol{u}}) - \mathit{\boldsymbol{x}} - {c_2}} \right\|_2^2} + \frac{{{\mu _3}}}{2}\left\| {\mathit{\boldsymbol{Du}} - \mathit{\boldsymbol{p}} - \mathit{\boldsymbol{y}} - {c_3}} \right\|_2^2{\kern 1pt} + \\ \quad \quad \quad \frac{{{\mu _4}}}{2}\left\| {\mathit{\boldsymbol{\varepsilon }}(\mathit{\boldsymbol{p}}) - \mathit{\boldsymbol{z}} - {c_4}} \right\|_2^2 + \frac{{{\mu _5}}}{2}\left\| {\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{v}} - {c_5}} \right\|_2^2 \end{array}$$ (4) 式中，ci是拉格朗日乘子，；i=1, 2, …, 5;μi是正的惩罚参数i=1, 2, …, 5。由ADMM算法，问题式 (3) 的迭代公式表示为：
$$\left\{ {\begin{array}{*{20}{l}} {({\mathit{\boldsymbol{w}}^{k + 1}},{\mathit{\boldsymbol{x}}^{k + 1}},{\mathit{\boldsymbol{y}}^{k + 1}},{\mathit{\boldsymbol{z}}^{k + 1}},{\mathit{\boldsymbol{v}}^{k + 1}},{\mathit{\boldsymbol{u}}^{k + 1}},{\mathit{\boldsymbol{p}}^{k + 1}}) = }\\ {\mathop {\arg \min }\limits_{\mathit{\boldsymbol{w,x,y,z,v,u,p}}} L(\mathit{\boldsymbol{w,x,y,z,v,u,p}})}\\ {c_1^{k + 1} = c_1^k + \gamma (\mathit{\boldsymbol{K}}{\mathit{\boldsymbol{u}}^{k + 1}} - {\mathit{\boldsymbol{w}}^{k + 1}})}\\ {c_2^{k + 1} = c_2^k + \gamma ({\rm{S}}{{\rm{H}}_j}({\mathit{\boldsymbol{u}}^{k + 1}}) - \mathit{\boldsymbol{x}}_j^{k + 1})}\\ {c_3^{k + 1} = c_3^k + \gamma (\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{u}}^{k + 1}} - {\mathit{\boldsymbol{p}}^{k + 1}} - {\mathit{\boldsymbol{y}}^{k + 1}})}\\ {c_4^{k + 1} = c_4^k + \gamma (\varepsilon ({\mathit{\boldsymbol{p}}^{k + 1}}) - {\mathit{\boldsymbol{z}}^{k + 1}})}\\ {c_5^{k + 1} = c_5^k + \gamma ({\mathit{\boldsymbol{u}}^{k + 1}} - {\mathit{\boldsymbol{v}}^{k + 1}})} \end{array}} \right.$$ (5) 式中，γ是迭代步长。由于 $c_1^k,c_2^k,c_3^k,c_4^k,c_5^k$ 很容易计算，因此主要求解式 (5) 中的极小化问题 $L(\mathit{\boldsymbol{w,x,y,z,v,u,p}})$ ，它可以分解为以下子问题：
$$\left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{w}}^{k + 1}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{w}} \beta \int_\mathit{\Omega } {\mathit{\boldsymbol{w}} - \mathit{\boldsymbol{f}}\log \mathit{\boldsymbol{w}}} + \frac{{{\mu _1}}}{2}\left\| {\mathit{\boldsymbol{K}}{\mathit{\boldsymbol{u}}^k} - \mathit{\boldsymbol{w}} - c_1^k} \right\|_2^2{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} }\\ {{\mathit{\boldsymbol{x}}^{k + 1}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{x}} \lambda \sum\limits_{j = 1}^M {{{\left\| {{\mathit{\boldsymbol{x}}_j}} \right\|}_1}} + \frac{{{\mu _2}}}{2}\sum\limits_{j = 1}^M {\left\| {{\rm{S}}{{\rm{H}}_j}({\mathit{\boldsymbol{u}}^k}) - {\mathit{\boldsymbol{x}}_j} - c_2^k} \right\|_2^2} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} }\\ {{\mathit{\boldsymbol{y}}^{k + 1}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{y}} {\alpha _1}{{\left\| \mathit{\boldsymbol{y}} \right\|}_2} + \frac{{{\mu _3}}}{2}\left\| {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{u}}^k} - {\mathit{\boldsymbol{p}}^k} - \mathit{\boldsymbol{y}} - c_3^k} \right\|_2^2{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} }\\ {{\mathit{\boldsymbol{z}}^{k + 1}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{z}} {\alpha _0}{{\left\| \mathit{\boldsymbol{z}} \right\|}_F} + \frac{{{\mu _4}}}{2}\left\| {\varepsilon ({\mathit{\boldsymbol{p}}^k}) - \mathit{\boldsymbol{z}} - c_4^k} \right\|_F^2{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} }\\ {{\mathit{\boldsymbol{v}}^{k + 1}} = \max \{ {\mathit{\boldsymbol{u}}^{k + 1}} + c_5^k,0\} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} }\\ \begin{array}{l} ({\mathit{\boldsymbol{u}}^{k + 1}},{\mathit{\boldsymbol{p}}^{k + 1}}) = \mathop {\arg \min }\limits_{\mathit{\boldsymbol{u,p}}} \frac{{\beta {\mu _1}}}{2}\left\| {\mathit{\boldsymbol{Ku}} - {\mathit{\boldsymbol{w}}^{k + 1}} - c_1^k} \right\|_2^2 + \\ \frac{{\lambda {\mu _2}}}{2}\sum\limits_{j = 1}^M {\left\| {{\rm{S}}{{\rm{H}}_j}(\mathit{\boldsymbol{u}}) - \mathit{\boldsymbol{x}}_j^{k + 1} - c_2^k} \right\|_2^2 + } \\ \frac{{{\alpha _1}{\mu _3}}}{2}\left\| {\mathit{\boldsymbol{Du}} - \mathit{\boldsymbol{p}} - {\mathit{\boldsymbol{y}}^{k + 1}} - c_3^k} \right\|_2^2 + \\ \frac{{{\alpha _0}{\mu _4}}}{2}\left\| {\varepsilon (\mathit{\boldsymbol{p}}) - {\mathit{\boldsymbol{z}}^{k + 1}} - c_4^k} \right\|_F^2{\kern 1pt} + \frac{{{\mu _5}}}{2}\left\| {\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{v}}^{k + 1}} - c_5^k} \right\|_2^2{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \end{array} \end{array}} \right.$$ (6) 式中，w子问题可直接求解，即对i=1, 2, …, N，有：
$$\begin{array}{l} \mathit{\boldsymbol{w}}_i^{k + 1} = \frac{1}{2}\left( {{{\left( {\mathit{\boldsymbol{K}}{\mathit{\boldsymbol{u}}^{k + 1}} + c_1^k - \frac{\beta }{{{\mu _1}}}} \right)}_i} + } \right.\\ \left. {\sqrt {\left( {\mathit{\boldsymbol{K}}{\mathit{\boldsymbol{u}}^{k + 1}} + c_1^k - \frac{\beta }{{{\mu _1}}}} \right)_i^2 + \frac{{4\beta {f_i}}}{{{\mu _1}}}} } \right) \end{array}$$ 子问题x、y、z的解，直接应用收缩阈值法：
1) 极小化问题x的解表示为：
$$\mathit{\boldsymbol{x}}_j^{k + 1} = {\rm{sgn}}({\rm{S}}{{\rm{H}}_j}({\mathit{\boldsymbol{u}}^k}) - c_2^k){\rm{max}}\left( {\left| {{\rm{S}}{{\rm{H}}_j}({\mathit{\boldsymbol{u}}^k}) - c_2^k} \right| - \frac{\lambda }{{{\mu _2}}},0} \right)$$ 2) y的解表示为：
$${\mathit{\boldsymbol{y}}^{k + 1}} = \max ({\left\| {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{u}}^k} - {\mathit{\boldsymbol{p}}^k} - c_3^k} \right\|_2} - \frac{{{\alpha _1}}}{{{\mu _3}}},0)\frac{{\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{u}}^k} - {\mathit{\boldsymbol{p}}^k} - c_3^k}}{{{{\left\| {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{u}}^k} - {\mathit{\boldsymbol{p}}^k} - c_3^k} \right\|}_2}}}$$ 3)z的解表示为：
$${\mathit{\boldsymbol{z}}^{k + 1}} = \max \left( {{{\left\| {\varepsilon ({\mathit{\boldsymbol{p}}^k}) - c_4^k} \right\|}_F} - \frac{{{\alpha _0}}}{{{\mu _4}}},0} \right)\frac{{\varepsilon ({\mathit{\boldsymbol{p}}^k}) - c_4^k}}{{{{\left\| {\varepsilon ({\mathit{\boldsymbol{p}}^k}) - c_4^k} \right\|}_F}}}$$ 另外，(u, p) 子问题是一个鞍点问题，可以分解为以下问题。
1) 对于u子问题，求解下面的极小化问题：
$$\begin{array}{l} {\mathit{\boldsymbol{u}}^{k + 1}} = \mathop {arg\min }\limits_\mathit{\boldsymbol{u}} \frac{{{\mu _1}\beta }}{2}\left\| {\mathit{\boldsymbol{Ku}} - {\mathit{\boldsymbol{w}}^{k + 1}} - c_1^k} \right\|_2^2 + \\ {\kern 1pt} \quad \quad \frac{{\lambda {\mu _2}}}{2}\sum\limits_{j = 1}^M {\left\| {{\rm{S}}{{\rm{H}}_j}(\mathit{\boldsymbol{u}}) - \mathit{\boldsymbol{x}}_j^{k + 1} - c_2^k} \right\|_2^2 + } \\ \frac{{{\alpha _1}{\mu _3}}}{2}\left\| {\mathit{\boldsymbol{Du}} - \mathit{\boldsymbol{p}} - {\mathit{\boldsymbol{y}}^{k + 1}} - c_3^k} \right\|_2^2 + {\kern 1pt} \frac{{{\mu _5}}}{2}\left\| {\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{v}}^{k + 1}} - c_5^k} \right\|_2^2 \end{array}$$ 通过求解下面的线性方程得到结果：
$$\begin{array}{l} \quad \quad \quad \quad \beta {\mu _1}{\mathit{\boldsymbol{K}}^{\rm{T}}}(\mathit{\boldsymbol{Ku}} - {\mathit{\boldsymbol{w}}^{k + 1}} - c_1^k) + \\ \quad \quad \quad \lambda {\mu _2}\sum\limits_{j = 1}^M {{\rm{SH}}_j^*} ({\rm{S}}{{\rm{H}}_j}(\mathit{\boldsymbol{u}}) - \mathit{\boldsymbol{x}}_j^{k + 1} - c_2^k) + \\ {\alpha _1}{\mu _3}\sum\limits_{j = 1}^2 {\mathit{\boldsymbol{D}}_j^{\rm{T}}({\mathit{\boldsymbol{D}}_j}\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{p}}_j} - \mathit{\boldsymbol{y}}_j^{k + 1} - c_{{3_j}}^k)} {\kern 1pt} {\kern 1pt} {\kern 1pt} + {\mu _5}(\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{v}}^{k + 1}} - c_5^k) = 0 \end{array}$$ 即u表示为：
$$\begin{array}{l} \quad \quad \quad {\mathit{\boldsymbol{u}}^{k + 1}} = (\beta {\mu _1}{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{K}} + \quad \quad \lambda {\mu _2}\sum\limits_{j = 1}^M {{\rm{SH}}_j^ * {\rm{S}}{{\rm{H}}_j}} + \\ \quad \quad {\mu _5}\mathit{\boldsymbol{I}} + {\alpha _1}{\mu _3}\sum\limits_{j = 1}^2 {\mathit{\boldsymbol{D}}_j^{\rm{T}}{\mathit{\boldsymbol{D}}_j}} {)^{ - 1}}(\beta {\mu _1}{\mathit{\boldsymbol{K}}^{\rm{T}}}({\mathit{\boldsymbol{w}}^{k + 1}} + c_1^k) + \\ \quad \quad \quad \quad \quad \lambda {\mu _2}\sum\limits_{j = 1}^M {{\rm{SH}}_j^*(\mathit{\boldsymbol{x}}_j^{k + 1} + c_2^k) + } \\ \quad \quad {\alpha _1}{\mu _3}\sum\limits_{j = 1}^2 {\mathit{\boldsymbol{D}}_j^{\rm{T}}} ({\mathit{\boldsymbol{p}}_j} + \mathit{\boldsymbol{y}}_j^{k + 1} + c_{{3_j}}^k)) + {\mu _5}({\mathit{\boldsymbol{v}}^{k + 1}} + c_5^k)) \end{array}$$ 2) 子问题p转化为以下极小化问题：
$$\begin{array}{l} {\mathit{\boldsymbol{p}}^{k + 1}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{p}} \frac{{{\alpha _1}{\mu _3}}}{2}\left\| {\mathit{\boldsymbol{Du}} - \mathit{\boldsymbol{p}} - {\mathit{\boldsymbol{y}}^{k + 1}} - c_3^k} \right\|_2^2 + \\ \quad \quad \quad \frac{{{\alpha _0}{\mu _4}}}{2}\left\| {\varepsilon (\mathit{\boldsymbol{p}}) - {\mathit{\boldsymbol{z}}^{k + 1}} - c_4^k} \right\|_F^2 \end{array}$$ 而子问题p又可分解成子问题p1、p2来求解：
① 子问题p1，可由下列线性方程获得：
$$\begin{array}{l} \quad \quad {\alpha _1}{\mu _3}({\mathit{\boldsymbol{p}}_1} - {\mathit{\boldsymbol{D}}_1}\mathit{\boldsymbol{u}} + \mathit{\boldsymbol{y}}_1^{k + 1} + c_{{3_1}}^{k + 1}) + \\ \quad \quad {\alpha _0}{\mu _4}\mathit{\boldsymbol{D}}_1^{\rm{T}}({\mathit{\boldsymbol{D}}_1}{\mathit{\boldsymbol{p}}_1} - \mathit{\boldsymbol{z}}_1^{k + 1} - c_{{4_1}}^{k + 1}) + \\ \frac{{{\alpha _0}{\mu _4}}}{2}\mathit{\boldsymbol{D}}_2^{\rm{T}}({\mathit{\boldsymbol{D}}_2}{\mathit{\boldsymbol{p}}_1} + {\mathit{\boldsymbol{D}}_1}{\mathit{\boldsymbol{p}}_2} - 2\mathit{\boldsymbol{z}}_3^{k + 1} - 2c_{{4_3}}^{k + 1}) = 0{\kern 1pt} {\kern 1pt} \end{array}$$ 因此 $\mathit{\boldsymbol{p}}_1^{k + 1}$ 的解可表示为：
$$\begin{array}{l} \mathit{\boldsymbol{p}}_1^{k + 1} = {({\alpha _1}{\mu _3} + {\alpha _0}{\mu _4}\mathit{\boldsymbol{D}}_1^{\rm{T}}{\mathit{\boldsymbol{D}}_1}{\kern 1pt} + \frac{{{\alpha _0}{\mu _4}}}{2}\mathit{\boldsymbol{D}}_2^{\rm{T}}{\mathit{\boldsymbol{D}}_2})^{ - 1}} \times \\ ({\alpha _1}{\mu _3}({\mathit{\boldsymbol{D}}_1}u - \mathit{\boldsymbol{y}}_1^{k + 1} - c_{{3_1}}^k) + {\alpha _0}{\mu _4}\mathit{\boldsymbol{D}}_1^T(\mathit{\boldsymbol{z}}_1^{k + 1} + c_{{4_1}}^k) + \\ \quad \frac{{{\alpha _0}{\mu _4}}}{2}\mathit{\boldsymbol{D}}_2^{\rm{T}}(2\mathit{\boldsymbol{z}}_3^{k + 1} + 2c_{{4_3}}^k - {\mathit{\boldsymbol{D}}_1}{\mathit{\boldsymbol{p}}_2} - {\mathit{\boldsymbol{D}}_2}{\mathit{\boldsymbol{p}}_1})) \end{array}$$ ② 类似的，子问题p2可表示为下列线性方程的解：
$$\begin{array}{l} \quad \quad \quad {\alpha _1}{\mu _3}({\mathit{\boldsymbol{p}}_2} - {\mathit{\boldsymbol{D}}_2}u + \mathit{\boldsymbol{y}}_2^{k + 1} + c_{{3_2}}^{k + 1}) + \\ \quad \quad \quad {\alpha _0}{\mu _4}\mathit{\boldsymbol{D}}_2^{\rm{T}}({\mathit{\boldsymbol{D}}_2}{\mathit{\boldsymbol{p}}_2} - \mathit{\boldsymbol{z}}_2^{k + 1} - c_{{4_2}}^{k + 1}) + \\ \frac{{{\alpha _0}{\mu _4}}}{2}\mathit{\boldsymbol{D}}_1^{\rm{T}}({\mathit{\boldsymbol{D}}_2}{\mathit{\boldsymbol{p}}_1} + {\mathit{\boldsymbol{D}}_1}{\mathit{\boldsymbol{p}}_2} - 2\mathit{\boldsymbol{z}}_3^{k + 1} - 2c_{{4_3}}^{k + 1}) = 0 \end{array}$$ 这样 $\mathit{\boldsymbol{p}}_2^{k + 1}$ 的解表示为：
$$\begin{array}{l} \mathit{\boldsymbol{p}}_2^{k + 1} = {({\alpha _1}{\mu _3} + {\alpha _0}{\mu _4}\mathit{\boldsymbol{D}}_2^{\rm{T}}{\mathit{\boldsymbol{D}}_2}{\kern 1pt} + \frac{{{\alpha _0}{\mu _4}}}{2}\mathit{\boldsymbol{D}}_1^{\rm{T}}{\mathit{\boldsymbol{D}}_1})^{ - 1}} \times \\ {\alpha _1}{\mu _3}({\mathit{\boldsymbol{D}}_2}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{y}}_2^{k + 1} - c_{{3_2}}^k) + {\alpha _0}{\mu _4}\mathit{\boldsymbol{D}}_2^{\rm{T}}(\mathit{\boldsymbol{z}}_2^{k + 1} + c_{{4_2}}^k){\kern 1pt} {\kern 1pt} + \\ \quad \frac{{{\alpha _0}{\mu _4}}}{2}\mathit{\boldsymbol{D}}_1^{\rm{T}}(2\mathit{\boldsymbol{z}}_3^{k + 1} + 2c_{{4_3}}^k - {\mathit{\boldsymbol{D}}_1}{\mathit{\boldsymbol{p}}_2} - {\mathit{\boldsymbol{D}}_2}{\mathit{\boldsymbol{p}}_1}) \end{array}$$ 这样通过ADMM交替求解式 (6)，并将其转化为以分块循环矩阵作为系数矩阵的线性系统，并通过傅里叶变换加速求解，最终获得复原图像。
-
为了验证新模型的有效性，分别对Lena和Peppers两幅图像进行实验，对Lena图像加高斯模糊，用文献[13]中的模糊函数psfGauss (5, 2)；Peppers图像加运动模糊，其中 $r = 6,\theta = 45$ ，然后再用Poissrnd函数加泊松噪声。实验结果与文献[3]中的方法PIDAL和文献[4]中的方法PID-Split进行比较，如图 1、图 2所示，其中图 1c~图 1e、图 2c~图 2e是用对应的算法获得的复原图像；图 2f~图 2h是局部放大图像。
图 1 Lena图像与其他复原方法的比较结果
图 2 Peppers图像与其他复原方法的比较结果
以上两个实验直观地展示了本文模型能够缓解甚至消除阶梯效应，并且很好地保持了图像的边缘细节信息。
为了更加客观地评价本文模型的优越性，用相对误差 (relative error, RelErr)、信噪比 (signal-to-noise ratio, SNR)、峰值信噪比 (peak signal-to-noise, PSNR)、结构相似度指数测量 (structural similarity index, SSIM)4种定量指标来评价本文模型的去噪效果。它们的计算公式分别为：
$$\begin{align} & \rm{RelErr=}\frac{{{\left\| \mathit{\boldsymbol{u}}-\mathit{\boldsymbol{\hat{u}}} \right\|}_{\rm{2}}}}{{{\left\| \mathit{\boldsymbol{u}} \right\|}_{\rm{2}}}},\rm{SNR=20}{{\log }_{\rm{10}}}\left( \frac{{{\left\| \mathit{\boldsymbol{u}} \right\|}_{\rm{2}}}}{{{\left\| \mathit{\boldsymbol{\hat{u}}}-\mathit{\boldsymbol{u}} \right\|}_{\rm{2}}}} \right) \\ & \quad \rm{PSNR=10}\lg (\rm{25}{{\rm{5}}^{\rm{2}}}\rm{/}\frac{\rm{1}}{\mathit{MN}}\ \sum\limits_{\mathit{i}\rm{=0}}^{\mathit{M}\rm{-1}}{\sum\limits_{\mathit{j}\rm{=0}}^{\mathit{N}\rm{-1}}{{{\left[ \mathit{\boldsymbol{u}}(\mathit{i,j})-\mathit{\boldsymbol{\hat{u}}}(\mathit{i,j}) \right]}^{\rm{2}}}}}) \\ & \quad \quad \rm{SSIM=}\frac{(\rm{2}{{\mu }_{\mathit{\boldsymbol{u}}}}{{\mu }_{{\mathit{\boldsymbol{\hat{u}}}}}}\rm{+}{{\mathit{C}}_{\rm{1}}})(\rm{2}{{\sigma }_{\mathit{\boldsymbol{u\hat{u}}}}}\rm{+}{{\mathit{C}}_{\rm{2}}})}{(\mu _{\mathit{\boldsymbol{u}}}^{\rm{2}}\rm{+}\mu _{{\mathit{\boldsymbol{\hat{u}}}}}^{\rm{2}}\rm{+}{{\mathit{C}}_{\rm{1}}})(\sigma _{\mathit{\boldsymbol{u}}}^{\rm{2}}\rm{+}\sigma _{{\mathit{\boldsymbol{\hat{u}}}}}^{\rm{2}}\rm{+}{{\mathit{C}}_{\rm{2}}})} \\ \end{align}$$ 式中，u和 ${\mathit{\boldsymbol{\hat u}}}$ 分别表示原始图像和恢复的图像； ${\mu _\mathit{\boldsymbol{u}}},{\mu _{\mathit{\boldsymbol{\hat u}}}}$ 、 ${\sigma _\mathit{\boldsymbol{u}}},{\sigma _{\mathit{\boldsymbol{\hat u}}}},{\sigma _{\mathit{\boldsymbol{u\hat u}}}}$ 分别是u、 ${\mathit{\boldsymbol{\hat u}}}$ 的均值、标准差；共变异数C1、C2为常数。一般来说，PSNR, SNR, SSIM值越大，RelErr值越小时图像的去噪效果越好。在表 1中，给出实验的数值结果，从中可以看出本文模型的优越性。
表 1 不同算法的SNR、PSNR、RelErr和MSE值的结果比较
图像 方法 SNR PSNR RelErr SSIM Lena PIDAL 24.01 29.76 0.063 0.866 PIDSplit 24.15 29.80 0.062 0.871 Proposed 24.54 30.16 0.059 0.894 Pepper PIDAL 25.02 30.37 0.056 0.936 PIDSplit 25.07 30.41 0.055 0.940 Proposed 25.53 30.87 0.052 0.942 -
本文基于TGV和Shearlet变换提出了一种新的泊松图像去噪模型。它在有效去除噪声的同时不仅可以消除阶梯效应，还能很好地保持图像的边缘细节信息。然而，模型引入了多个参数，如何自适应地选择参数将需要进一步的研究。
Research on Poisson Noise Image Restoration Problems Based on Shearlet Transform
-
摘要: 为了解决泊松噪声图像的复原问题，几种正则化方法已被提出，其中最著名的是全变差（TV）模型，但TV模型会引起阶梯效应。总广义变差（TGV）是全变差的推广，用TGV作为正则项来恢复泊松图像，可以消除阶梯效应，但图像的边缘细节信息不能很好地保持。为了克服这个缺点，基于TGV和Shearlet变换，该文提出了一种新的正则化模型，并用交替方向乘子法（ADMM）求解。数值结果有效地展示了该模型在保持图像边缘细节上的优越性。
-
关键词:
- 交替方向乘子法 /
- 泊松噪声 /
- Shearlet变换 /
- 阶梯效应 /
- 总广义变差
Abstract: Restoring Poisson noise images has been drawn a lot of attention in recent years. To solve this problem, several regularization methods have been put forward. One of the most famous methods is the Total variation (TV) model. However, the TV model will cause staircasing effects. The total generalized variation (TGV) is the extension of TV. Using TGV as a regularization term to recover the Poission image can eliminate staircase effects but the edge details will not preserved very well. In order to overcome this drawback, based on TGV and Shearlet transform, we propose a new regularization method. The proposed model is solved by the alternating direction method of multiplier (ADMM). The numerical results reflect the efficiency of the new model in dealing with Poisson noise image. -
表 1 不同算法的SNR、PSNR、RelErr和MSE值的结果比较
图像 方法 SNR PSNR RelErr SSIM Lena PIDAL 24.01 29.76 0.063 0.866 PIDSplit 24.15 29.80 0.062 0.871 Proposed 24.54 30.16 0.059 0.894 Pepper PIDAL 25.02 30.37 0.056 0.936 PIDSplit 25.07 30.41 0.055 0.940 Proposed 25.53 30.87 0.052 0.942 
下载: 导出CSV
-
[1] RUDIN L I, OSHER S, FATERMI E. Nonlinear total variation based noise removal algorithms[J]. Physica D: Nonlinear Phenomena, 1992, 60(1): 259-268. [2] LE T, CHARTRAND R, ASAKI T J. A variational approach to reconstructing images corrupted by Poisson noise[J]. Journal of Mathematical Imaging and Vision, 2007, 27(3): 257-263. doi: 10.1007/s10851-007-0652-y [3] FIGUEIREDO M A T, BIOUCAS-DIAS J M. Restoration of Poissonian images using alternating direction optimization[J]. IEEE Transactions on Image Processing, 2010, 19(12): 3133-3145. doi: 10.1109/TIP.2010.2053941 [4] LIU X, HUANG L. Total bounded variation‐based Poissonian images recovery by split Bregman iteration[J]. Mathematical Methods in the Applied Sciences, 2012, 35(5): 520-529. doi: 10.1002/mma.v35.5 [5] BREDIES K, KUNISCH K, POCK T. Total generalized variation[J]. SIAM Journal on Imaging Sciences, 2010, 3(3): 492-526. doi: 10.1137/090769521 [6] BREDIES K, VALKONEN T. Inverse problems with second-order total generalized variation constraints[C]//Proceedings of SampTA 2011 – 9th International Conference on Sampling Theory and Applications. Singapore: IEEE Press, 2011: 201-205. [7] GUO W, QIN J, YIN W. A new detail-preserving regularization scheme[J]. SIAM Journal on Imaging Sciences, 2014, 7(2): 1309-1334. doi: 10.1137/120904263 [8] ETTER S, GROHS P, OBERMEIER A. FFRT: a fast finite ridgelet transform for radiative transport[J]. SIAM Journal on Multiscale Modeling & Simulation, 2015, 13(1): 1-42. https://www.researchgate.net/publication/276235323_FFRT_A_Fast_Finite_Ridgelet_Transform_for_Radiative_Transport [9] LABATE D, KUTYNIOK G. Sparse multidimensional representation using shearlets[J]. Proceedings of SPIE-the International Society for Optical Engineering, 2005, 5914(1): 254-262. https://www.researchgate.net/profile/Demetrio_Labate2/publication/220132077_Optimally_Sparse_Multidimensional_Representation_Using_Shearlets/links/00b7d5395d06dc356f000000.pdf?origin=publication_detail [10] GABAY D, MERCIER B. A dual algorithm for the solution of nonlinear variational problems via finite element approximation[J]. Computers & Mathematics with Applications, 1976, 2(1): 17-40. https://www.researchgate.net/publication/222625249_A_dual_algorithm_for_the_solution_of_nonlinear_variational_problems_via_finite_element_approximation [11] GLOWINSKI R, MARROCO A. Sur l'approximation, par éléments finis d'ordre un, et la résolution, par pénalisation-dualité d'une classe de problèmes de Dirichlet non linéaires[J]. ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique, 1975, 9(R2): 41-76. http://citeseerx.ist.psu.edu/showciting?cid=1045342 [12] BERTSEKAS D P, TSITSIKLIS J N. Parallel and distributed computation: numerical methods[M]. Englewood Cliffs, New Jersey: Prentice Hall Inc, 1989: 243-260. [13] NAJY J G, PALMER K, PERRONE L. Iterative methods for image deblurring: a Matlab object-oriented approach[J]. Numerical Algorithms, 2004, 36(1): 73-93. doi: 10.1023/B:NUMA.0000027762.08431.64 -
点击查看大图
图(2) / 表(1)
计量
- 文章访问数: 3947
- HTML全文浏览量: 1451
- PDF下载量: 189
- 被引次数: 0
