CN101135722A - 基于重构信号替代频谱数据的信号去噪方法 - Google Patents

Abstract

本发明涉及一种基于重构信号替代频谱数据的信号去噪方法，包括从观测信号g(x)中提取奇异点和奇异值、根据得到的奇异点和奇异值构造相应的奇异函数、根据系统预设的频谱置零阈值构造相应的去噪函数和置换函数、根据上述的结果利用重构信号替代法得到去噪后的信号。采用该种基于重构信号替代频谱数据的信号去噪方法，有效地避免了信号除噪过程中的高频信号丢失或信号失真，从而在保证图像的高分辨率和高精确度条件下，能够有效去除图像噪声，提高了信噪比；而且该方法高效实用，工作性能稳定可靠、适用范围较为广泛，给信号采集和处理带来很大的便利，并且也为医学图像处理技术的进一步发展和大范围普及应用奠定了坚实的理论和实践基础。

Description

基于重构信号替代频谱数据的信号去噪方法

技术领域

本发明涉及信号及医学图像处理技术领域，特别涉及磁共振成像保真信号除噪声领域，具体是指一种基于重构信号替代频谱数据的信号去噪方法。

背景技术

随着现代医学技术的不断发展，核磁共振成像(MRI)技术已经成为医学成像检测领域中不可或缺的手段，其中，磁共振信号空间(原始数据空间)称为K空间，即为傅里叶变换空间，K空间采样到信号经过傅里叶反变换后再取模，即得到核磁共振(MR)图像。

通常的图像信号中均包含有各种噪声，而噪声可分为加性噪声和乘性噪声，乘性噪声大小与其染污的信号大小成比例，加性噪声大小与其染污的信号大小无关。含加性噪声的观测信号f(x)数学模型可表述为：

f(x)＝g(x)+ns(x)，x＝0，1，...，N-1......(0)

其中g(x)和ns(x)分别表示为无噪声真实信号序列和噪声信号序列。大多情况下，g(x)序列是一个非平稳信号，因而观测信号f(x)一般也是一个非平稳信号。

现有技术中的除噪声方法大致有以下三类：

第一类：多重观测信号平均法。其主要思想是根据ns(x)序列的各个元素可认作相互独立的、具有同分布的、零均值的、平稳的随机变量的假设。这样将多次采集到的观测信号f(x)序列迭加平均时，随机噪声ns(x)会相互抵消弱，从而达到去噪目的。这种方法是目前公认的保真去噪方法，在医学设备中广泛使用。但是这种方法有如下缺陷：

(1)要重复采集同一信源的观测信号是很困难的，有时甚至是不可能的。例如，当信源有随机运动(如心脏跳动)、信源是时变系统(心率失常病人的心电信号)等情况(请参阅文献：He Wei，Xie Zhengxian，Studying Waveform Distortion of Electricardiac Signals IntroducedBy Signal-Averaged，CHINESE J MED PHYS Vol.16No.1，Jan.1999，pp26-27)；

(2)重复采集观测信号费时费设备占用时间，降低设备使用效。

为此，人们发展出了以下的单一观测信号的除噪方法。

第二类：单一观测信号邻域估计法。这类方法的基本思想是基于f(x)在局部小邻域内也可认作近似相互独立的、具有同分布的、零均值的、平稳的随机变量的假设。这样就可以用局部空间邻域估计值(如均值，中值，拟合值)对f(x)进行去噪。但是在绝大多数情况下，此类方法常常会使信号细节的丢失，从而失真。为此，人们又提出了提高信号保真的信号空间域方法(请参阅文献：Charles，D.；Davies，E.R.，Distance-weighted median filters and theirapplication to colour images，Visual Information Engineering，2003.VIE2003.InternationalConference on7-9July2003Page(s)：117-120；Nai-Xiang Lian；Zagorodnov，V.；Yap-Peng Tan，Edge-preserving image denoising via optimal color space proj ection，Image Processing，IEEETransactions on，Volume15，Issue9，Sept.2006Page(s)：2575-2587；Balster，E.J.；Zheng，Y.F.；Ewing，R.L.，Combined spatial and temporal domain wavelet shrinkage algorithm for videodenoising，Circuits and Systems for Video Technology，IEEE Transactions on Volume16，Issue2，Feb.2006Page(s)：220-230；Rosiles，J.G.；Smith，M.J.T.，Image denoising using directional filterbanks，Image Processing，2000.Proceedings.2000Intemational Conference on Volume3，10-13Sept.2000Page(s)：292-295；和Han Liu；Yong Guo；Gang Zheng，Image Denoising Based onLeast Squares Support Vector Machines，Intelligent Control and Automation，2006.WCICA2006.The Sixth World Congress on，Volume1，21-23June2006Page(s)：4180-4184)。但是这类方法的缺陷是空间去噪的信号失真问题无法根本解决。

第三类：单一观测信号变换域系数分离法。这类方法的基本假定是：噪声污染信号可以在变换域里区分为信号变换域系数和噪声变换域系数，可以将噪声变换域系数置零，然后用反变换法重构出无噪声信号，达到去噪声目的。常见办法有：Fourier变换，小波变换(请参阅文献：Yunyi Yan；Baolong Guo；Wei Ni，Image Denoising：An Approach Based on WaveletNeural Network and Improved Median Filtering，Intelligent Control and Automation，2006.WCICA2006.The Sixth World Congress on Volume2，21-23June2006Page(s)：10063-10067；Hui Cheng；Qiuze Yu；Jinwen Tian；Jian Liu，Image denoising using wavelet and supportvector regression，Image and Graphics，2004.Proceedings.Third Intemational Conference on18-20Dec.2004Page(s)：43-46；Wink，A.M.；Roerdink，J.B.T.M.，Denoising functional MR images；acomparison of wavelet denoising and Gaussian smoothing，Medical Imaging，IEEETransactions onVolume23，Issue3，March2004Page(s)：374-387；Nai-Xiang Lian；Zagorodnov，V.；Yap-Peng Tan，Color image denoising using wavelets and minimum cut analysis，Signal ProcessingLetters，IEEE Volume12，Issue11，Nov.2005Page(s)：741-744；Chen，G.Y.；Bui，T.D.；Krzyzak，A.，Image denoising using neighbouring wavelet coefficients，Acoustics，Speech，andSignal Processing，2004.Proceedings.(ICASSP′04).IEEE International Conference on Volume2，17-21May2004Page(s)：917-202；Zhang，S.；Salari，E.，Image denoising using a neural networkbased non-linear filter in wavelet domain，Acoustics，Speech，and Signal Processing，2005.Proceedings.(ICASSP′05).IEEE International Conference on Volume2，18-23March2005Page(s)：989-992)，稀疏变换(Sparse Tranform)(请参阅文献：Li Shang；Deshuang Huang，Image denoising using non-negative sparse coding shrinkage algorithm，Computer Vision andPattern Recognition，2005.CVPR2005.IEEE Computer Society Conference on，Volume1，20-25June2005Page(s)：1017-1022；Guleryuz，O.G.，Nonlinearapproximation based imagerecovery using adaptive sparse reconstructions and iterated denoising-part I：theory，ImageProcessing，IEEE Transactions on，Volume15，Issue3，March2006Page(s)：539-554；Elad，M.；Aharon，M.，Image Denoising Via Sparse and Redundant Representations Over LeamedDictionaries，Image Processing，IEEE Transactions on，Volume15，Issue12，Dec.2006Page(s)：3736-3745)，Hilbert-Huang变换(请参阅文献：Zhuo-Fu Liu；Zhen-Peng Liao；En-FangSang，Speech enhancement based on Hilbert-Huang transform，Machine Learning and Cybernetics，2005.Proceedings of2005InternationalConferenceon，Volume8，18-21Aug.2005Page(s)：4908-4912；Xiaojie Zou；Xueyao Li；Rubo Zhang，Speech Enhancement Based onHilbert-Huang Transform Theory，Computer and Computational Sciences，2006.IMSCCS′06.First International Multi-Symposiums on，Volume1，20-24June2006Page(s)：208-213；WuWang；Xueyao Li；Rubo Zhang，Speech Detection Based on Hilbert-Huang Transform，Computerand Computational Sciences，2006.IMSCCS′06.First International Multi-Symposiums on，Volurne1，20-24June2006Page(s)：290-293)等等。

事实上到目前为止，没有一种变换可以将观测信号中的真实信号和噪声的变换系数严格区分开来，这样难免会有或多或少的信号丢失，总之，单一信号去噪声核心问题——信号失真问题依然没有得到较好的解决。

发明内容

本发明的目的是克服了上述现有技术中的缺点，提供一种能够有效避免信号除噪过程中的高频信号丢失或信号失真、同时兼有其它除噪方法的优越性、信号采集容易、有效去除图像噪声、精确显示原磁共振图像、高效实用、工作性能稳定可靠、适用范围较为广泛的基于重构信号替代频谱数据的信号去噪方法。

为了实现上述的目的，本发明的基于重构信号替代频谱数据的信号去噪方法如下：

该基于重构信号替代频谱数据的信号去噪方法，其主要特点是，所述的方法包括以下步骤：

(1)从观测信号g(x)中提取奇异点和奇异值；

(2)根据得到的奇异点和奇异值构造相应的奇异函数；

(3)根据系统预设的频谱置零阈值，构造相应的去噪函数和置换函数；

(4)根据上述的观测信号g(x)、去噪函数、置换函数和奇异函数，利用重构信号替代法得到去噪后的信号。

该基于重构信号替代频谱数据的信号去噪方法中的从观测信号g(x)中提取奇异点和奇异值，包括以下步骤：

(11)对观测信号g(x)进行差分；

(12)在上述差分中取差分绝对值大于系统预设的噪声阈值的那些点作为奇异点b₁，b₂，...，b_Q，其中Q为奇异点的数量；

(13)取所述的奇异点b₁，b₂，...，b_Q的差分值作为相应的奇异值a₁，a₂，...，a_Q。

该基于重构信号替代频谱数据的信号去噪方法中的根据得到的奇异点和奇异值构造相应的奇异函数，具体为：

根据所述的奇异点b₁，b₂，...，b_Q和相应的奇异值a₁，a₂，...，a_Q，通过以下公式得到奇异函数：

$\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}{a}_q{w}_{b_q}\left(x\right);$

其中，w_bq(x)为奇异点b_q的奇异函数， $w_{b_q}\left(x\right)=w(x-x_{b_q})=\left\{\begin{array}{c}1,x\&GreaterEqual;x_{b_q}\\ 0,x<x_{b_q}\end{array}.\right.$

该基于重构信号替代频谱数据的信号去噪方法中的根据系统预设的频谱置零阈值构造去噪函数和置换函数，包括以下步骤：

(31)根据以下公式得到观测信号g(x)的频谱函数G(k)：

G(k)＝DFT[g(x)]；

其中DFT[·]为离散付里叶变换算子；

(32)将观测信号频谱幅度值|G(k)|与系统预设的频谱置零阈值T进行比较，并根据以下条件得到高信噪比频谱区域Ω和其余频谱区域Ω

如果|G(k)|≤T，则k∈Ω；

如果|G(k)|<T，则k∈Ω；

(33)根据高信噪比频谱区域Ω和其余频谱区域Ω通过以下公式得到指示函数R(k)和R(k)：

$R\left(k\right)=\left\{\begin{array}{c}1,k\&Element;\mathrm{\&Omega;}\\ 0,k\&NotElement;\mathrm{\&Omega;}\end{array}\right.;$

$\stackrel{\&OverBar;}{R}\left(k\right)\left\{\begin{array}{c}1,k\&Element;\stackrel{\&OverBar;}{\mathrm{\&Omega;}}\\ 0,k\&NotElement;\stackrel{\&OverBar;}{\mathrm{\&Omega;}}\end{array}\right.;$

(34)根据以下公式得到去噪函数r(x)和置换函数r(x)：

r(x)＝IDFT[R(k)]；

r(x)＝IDFT[R(k)]；

其中，IDFT[·]为离散付里叶反变换算子。

该基于重构信号替代频谱数据的信号去噪方法中的利用重构信号替代法得到去噪后的信号，具体为：

根据观测信号g(x)、去噪函数r(x)、置换函数r(x)和奇异函数

通过以下公式得到去噪后的信号g_r(x)：

$g_r\left(x\right)=g\left(x\right)*r\left(x\right)+\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_q{w}_{b_q}\left(x\right)*\stackrel{\&OverBar;}{r}\left(x\right);$

其中，“*”为卷积积分运算符。

采用了该发明的基于重构信号替代频谱数据的信号去噪方法，用于实际磁共振图像去噪，其方法步骤是：首先从对观测信号提取奇异点和奇异值，构造去噪函数和置换函数，最后由观测信号与去噪函数卷积积分、置换函数与奇异函数卷积积分的加权和共同得到去噪声后的信号。因而这种方法有效地避免了信号除噪过程中的高频信号丢失或信号失真，同时具有其它信号去噪方法法的优越之处，并巧妙克服了单一观测信号去噪法的信号失真问题，从而在保证图像所具有的高分辨率和高精确度条件下，而且能够有效去除图像噪声，提高了信噪比，为医学核磁共振成像检测提供了高质量的可靠图像信息；同时，本发明的方法高效实用，工作性能稳定可靠、适用范围较为广泛，给信号采集和处理带来了很大的便利，并且也为医学图像处理技术的进一步发展和大范围普及应用奠定了坚实的理论和实践基础。

附图说明

图1a、1c分别为仿真图像去噪声实验中的原始标准无噪图像和标准差为5、零均值的高斯白噪声图像示意图。

图1b为将图1c的高斯白噪声加入图1a的原始标准无噪图像后的用于算法测试的图像示意图。

图1d、1e、1f分别为对图1b采用中值滤波法去噪、小波滤波法去噪和本发明的基于重构信号替代频谱数据的信号去噪方法去噪后得到的图像比较示意图。

图2a、2b、2c分别为图1d、1e、1f相对于图1a的误差图像比较示意图。

图2d、2e、2f分别为图1d、1e、1f相对于图1b的误差图像比较示意图。

图2g、2h、2i分别为图2d、2e、2f进行付里叶变换后的频谱图像比较示意图。

图3为图1b、2d、2e、2f相对于图1a的误差直方图统计示意图。

图4为对于加入标准差为1～9的零均值高斯白噪声的测试图像分别采用中值滤波法去噪、小波滤波法去噪和本发明的基于重构信号替代频谱数据的信号去噪方法去噪后得到的图像与原始标准图像的误差标准差比较示意图。

图5a、5b分别为实际水模磁共振成像试验中采用本发明的图像去噪方法前后的图像比较示意图。

图5c为图5a、5b的误差图像示意图。

图5d为图5a、5b的频谱幅度误差图像示意图。

图5e、5f分别为图5a、5b在第216列的线图比较示意图。

图6a、6b分别为实际人头磁共振成像试验中采用本发明的图像去噪方法前后的图像比较示意图。

图6c为图6a、6b的误差图像示意图。

图6d为图6a、6b的频谱幅度误差图像示意图。

图7为本发明的基于重构信号替代频谱数据的信号去噪方法的算法流程示意图。

具体实施方式

为了能够更清楚地理解本发明的技术内容，特举以下实施例详细说明。

在阐述本发明的整体工作过程及工作原理之前，为了更加明确其技术含义，首先需要对奇异谱模型做如下介绍：

根据奇异谱分析理论(请参阅文献：Luo Jianhua and Yuemin Zhu，MR image ReconstructionFrom truncated k-space using a layer singular point extraction technique.IEEE Transactions onNuclear Science，2004.1，Vol.51，No.1；157：169)，任意信号f(x)都可以表示为奇异函数加权和，即有：

$f\left(x\right)=\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_q{w}_{b_q}\left(x\right),x=\mathrm{0,1}...,M-1......\left(1\right)$

其中：

$w(x-x_{b_1})=\left\{\begin{array}{c}1,x\&GreaterEqual;x_{b_1}\\ 0,x<x_{b_1}\end{array}\right.......\left(2\right)$

为奇异函数，{b₁，b₂，...，b_Q}为f(x)上的全体奇异点，数量为Q，即f(x)差分不为零的点，a₁，a₂，...，a_Q是该Q个奇异点{b₁，b₂，...，b_Q}上的奇异值，即差分值，{w_b1(x)，w_b2(x)，...，w_bQ(x)}分别为关于{b₁，b₂，...，b_Q}为奇异点的Q个奇异函数。记DFT[·]为离散付里叶变换算子，F(k)＝DFT[f(x)]。对上述公式(1)两边取付里叶变换得，任意付里谱函数F(k)可以表示奇异谱函数加权和，即有：

$F\left(k\right)\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_q{W}_{b_q}\left(k\right),k=\mathrm{0,1}...,M-1......\left(3\right)$

其中：

$W_{b_q}\left(k\right)=\mathrm{DFT}\left[w_{b_q}\left(x\right)\right]......\left(4\right)$

称为W_bq(k)奇异谱函数。上述的公式(1)和(3)构成奇异谱分析模型。根据此式，就可以利用奇异点和奇异值来构造信号f(x)的付里叶变换频谱数据。

在进行本发明的重构信号置换去噪方法时，首先设观测信号g(x)为无噪信号f(x)及加性高斯白噪声n(x)之和，即：

g(x)＝f(x)+n(x)            ......(5)

分别记为G(k)＝DFT[g(x)]、F(k)＝DFT[f(x)]和N(k)＝DFT[n(x)]，对上述公式(5)两边同取付里叶变换，得到：

G(k)＝F(k)+N(k)               ......(6)

由于高斯白噪声的功率谱密度|N(k)|为常量，因而在观测信号频谱|G(k)|幅度大的频段信噪比高，观测信号频谱幅度小的频段信噪比低。因此可以将全体高信噪比的频谱区域用Ω表示，而其余频谱区域Ω表示。这样就可以将高信噪的Ω区域频谱数据保留，而将其余的频谱区域Ω数据用关于无噪信号f(x)的奇异谱分析模型重构的数据置换，从而达到去噪目的。为方便获取Ω和Ω中频谱数据，分别定义指示函数如下为：

$R\left(k\right)=\left\{\begin{array}{c}1,k\&Element;\mathrm{\&Omega;}\\ 0,k\&NotElement;\mathrm{\&Omega;}\end{array}\right.......\left(7\right)$

和：

$\stackrel{\&OverBar;}{R}\left(k\right)=\left\{\begin{array}{c}1,k\&Element;\stackrel{\&OverBar;}{\mathrm{\&Omega;}}\\ 0,k\&NotElement;\stackrel{\&OverBar;}{\mathrm{\&Omega;}}\end{array}\right.......\left(8\right)$

这样总有：

$R\left(k\right)+\stackrel{\&OverBar;}{R}\left(k\right)=1......\left(9\right)$

考虑到， $G\left(k\right)=\underset{i=1}{\overset{Q_g}{\mathrm{\&Sigma;}}}{a}_i{W}_{b_i}\left(k\right)$ 观测信号可以表示为：

g(x)＝IDFT[G(k)]

＝IDFT[(R(k)+R(k))G(k)]      ......(10)

＝IDFT[R(k)G(k)+R(k)G(k)]

由于R(k)G(k)部分频谱已经舍去，因此可以用奇异谱

替换，则重构信号置换去噪可以由下式表示：

$g_r\left(x\right)=\mathrm{IDFT}[R\left(k\right)G\left(k\right)+\stackrel{\&OverBar;}{R}\left(k\right)\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_l{W}_{b_l}\left(k\right)]......\left(11\right)$

上述公式(11)右边第一项表示体现为去噪声作用，第二项表示对丢失信号的置换。重构信号置换去噪法可以写成以下卷积形式：

$\hat{f}\left(x\right)=g\left(x\right)*r\left(x\right)\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_{i}w(x-x_i)*\stackrel{\&OverBar;}{r}\left(x\right)......\left(12\right)$

其中r(x)＝IDFT[R(k)]和r(x)＝IDFT[R(k)]分别称为去噪函数和置换函数，“*”表示卷积。

在正确获得无噪声信号f(x)的奇异点b₁，b₂，...，b_Q和奇异值a₁，a₂，...，a_Q的条件下，那些低信噪比频谱中的噪声便完全除去。但事实上无法得到无噪声真实信号f(x)的奇异点和奇异值，只有靠观测信号g(x)进行估计。其方法是：对g(x)进行差分，并取绝对差分值大于噪声阈值的那些点作为f(x)的奇异点，其差分值作为奇异值。

这样做不可避免地给奇异值和奇异点引入误差。如果这种误差产生的负面影响大于噪声产生的负面影响，那么这种取代是不可取的。为此，只是差分绝对值小的奇异点被略去，因为奇异值绝对值越小，它对

的贡献也就越小。另外，这种替代仅在较高频段进行，这种误差不会被放大，因为绝大多数情况下高频段总有较低的信号—噪声谱能比，这种替代将会是大有帮助的。

据此，请参阅图7所示，该基于重构信号替代频谱数据的信号去噪方法，包括以下步骤：

(1)从观测信号g(x)中提取奇异点和奇异值；包括以下步骤：

(a)对观测信号g(x)进行差分；

(b)在上述差分中取差分绝对值大于系统预设的噪声阈值的那些点作为奇异点b₁，b₂，...，b_Q，其中Q为奇异点的数量；

(c)取所述的奇异点b₁，b₂，...，b_Q的差分值作为相应的奇异值a₁，a₂，...，a_Q；

(2)根据得到的奇异点和奇异值构造相应的奇异函数，具体为：

根据所述的奇异点b₁，b₂，...，b_Q和相应的奇异值a₁，a₂，...，a_Q，通过以下公式得到奇异函数：

$\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_q{w}_{b_q}\left(x\right);$

其中，w_bq(x)为奇异点b_q的奇异函数， $w_{b_q}\left(x\right)=w(x-x_{b_q})=\left\{\begin{array}{c}1,x\&GreaterEqual;x_{b_q}\\ 0,x<x_{b_q}\end{array}\right.;;$

(3)根据系统预设的频谱置零阈值，构造相应的去噪函数和置换函数，包括以下步骤：

(31)根据以下公式得到观测信号g(x)的频谱函数G(k)：

G(k)＝DFT[g(x)]；

其中DFT[·]为离散付里叶变换算子；

(32)将观测信号频谱幅度值|G(k)|与系统预设的频谱置零阈值T进行比较，并根据以下条件得到高信噪比频谱区域Ω和其余频谱区域Ω

如果|G(k)|≤T，则k∈Ω；

如果|G(k)|<T，则k∈Ω；

(33)根据高信噪比频谱区域Ω和其余频谱区域Ω，通过以下公式得到指示函数R(k)和R(k)：

$R\left(k\right)=\left\{\begin{array}{c}1,k\&Element;\mathrm{\&Omega;}\\ 0,k\&NotElement;\mathrm{\&Omega;}\end{array}\right.;$

$\stackrel{\&OverBar;}{R}\left(k\right)=\left\{\begin{array}{c}1,k\&Element;\stackrel{\&OverBar;}{\mathrm{\&Omega;}}\\ 0,k\&NotElement;\stackrel{\&OverBar;}{\mathrm{\&Omega;}}\end{array}\right.;$

(34)根据以下公式得到去噪函数r(x)和置换函数r(x)：

r(x)＝IDFT[R(k)]；

r(x)＝IDFT[R(k)]；

其中，IDFT[·]为离散付里叶反变换算子；

(4)根据上述的观测信号g(x)、去噪函数、置换函数和奇异函数，利用重构信号替代法得到去噪后的信号，具体为：

根据观测信号g(x)、去噪函数r(x)、置换函数r(x)和奇异函数通过以下公式得到去噪后的信号gr(x)：

$g_r\left(x\right)=g\left(x\right)*r\left(x\right)+\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_q{w}_{b_q}\left(x\right)*\stackrel{\&OverBar;}{r}\left(x\right);$

其中，“*”为卷积积分运算符。

在实际应用当中，为了测试信号置换的去噪声方法的效果，本发明中对仿真数据和实际磁共振图像数据进行了的实验。在仿真数据实验中，还与小波滤波和中值滤波去噪声方法进行了比较。实验中，将这方法用于图像去噪声，只要将图像按行展开，便成为一维信号，进行去噪声处理。实验结果比较手法有：图像比较，图像误差比较等。现分别下面将实验结果讨论分析如下：

仿直图像去噪声实验分析

实验中，标准无噪图像为256X256图像(如图1a)。对标准图像加上标准差为5的、零均值的、高斯白噪声作为用于测试算法的图像(如图1.b)，其噪声图像如图1c所示。中值滤波使用3x3滤波窗口，小波滤波使用小小波函数为sym8，取全局软阀值，有关参数由函数ddencmp(′den′，′wv′，X)计算获得。本实验中的信号置换的阀值T按重构像素概率为98％计算获得，即取一个T，使得98％频谱数据用重构频谱数据来取代。中值滤波法、小波滤波法和信号置换法去噪声后图像分别如图1d、1e、1f所示。从图1a、1b、1d、1e和1f上因图像太小几乎看不出什么差。为此，我们作各去噪声后图像与标准图像分别进行减运算，得到误差图像，结果如图2a、2b和2c。从图2a中发现，中值滤波的图像边缘误差大，但数量不多，范围在(-200，200)之间，对边缘有一定的损坏。从图2b中发现，小波分析法的图像也在边缘误差较大，且点数也多，误差范围在(-80，80)之间，边缘损坏严重。而从图2c中发现，信号置换法去噪声的图像，在边缘处也有一定的误差，但幅度仅为(-25，25)，在加入噪声强度范围之内，没有因去噪声对图像边缘造成损坏。我们把各种方法去噪声后的图像也去噪声前图像图1b也进行了减运算，得出误差图像如图2d、2e和2f所示，得和图2a、2b和2c相似的结果。对图2d、2e和2f进行付里叶变换，并用log(1+|E(k_x，k_y)|)频谱幅度图像来观察测频谱幅度的变化，如图2g、2h和2i图所示。从图2g、2h频谱幅度图上发现中值滤波和小波分析去噪的频谱图明有显高频亮区，说明高频分量损失较大，而i图上发现频谱重构法去噪的频谱幅度图除中心区域低频保持少误差的特征外，符合随机噪声频谱特征，对高频没有特别的损失。

为了进一步分析误差的分布，本发明对误差图的误差进行了直方图统计，结果如图3，其中横坐标为误差值，纵坐标为误差值的频率。误差直方图的曲线成为0点位置上的脉冲时表示两图像的误差为零，直方图曲线越平坦，表示离开零误差越远，即误差越大。从误差直方图上可以清晰看到，本发明的重构信号替代频谱数据的信号去噪方法的直方曲线2最陡，误差最小。其次为中值滤波的曲线3，小波分析的曲线4最差，它比去噪声前的图像曲线1还要平坦，因而质量还要差，因为它比被去噪声图像质量还要差。

另外，为了测试算法对噪声大小的敏感性，本发明对加入标准差为1～9的零均值高斯白噪声图像，用中值滤波法，小波滤波法和本发明的信号置换去噪声方法进行了测试，结果如图4所示，横坐标为加入的高斯白噪声图像的标准差，纵坐标为各个去噪后的图像与无噪标准图像的标准差。

其中：

·“○”表示标准差为1～9的零均值高斯白噪声值

·“△”表示中值滤波法的误差标准差值

·“+”表示小波滤波法的误差标准差值

·“*”表示本发明的信号置换去噪声方法的误差标准差值

从图中可以发现信号置换去噪声方法能适应各种噪声大小情况，都有很好的去噪效果。而中值和小波的误差标准差都大于加入噪声的标准差。

通过仿真实验分析表明信号置换去噪声方法能很好适应种强度噪声污染图象的去噪声情况，并有很好的保图像边缘的特征。这是说明本方法和现有单一观测信号去噪声法完全不同，有不可比拟的优点。

实际图像去噪声实验

1、实际水模磁共振图像

实验数据用FLASH方法扫描，图像大小为512×512，Te＝0.172ms，Tr＝400ms，平均次数为2。该图像多用于测试成像系统的成像分辨率，其特点是图像内有一光栅，可以用来观测的图像空间分辨率情况。

图5a为去噪声前原始图像，图5b为频谱重构法去噪声后图像，图5c为去噪声前后图像的误差，即图5a和图5b的差，图5d为去噪声前后误差的频谱幅度图，即图5c的频谱图。从图5b中发现，噪声已经基本去掉了，而且比图5a清晰，图5a和5b的信噪比分别是24.05和48.54，说明本方法有很强的去噪能力。图5c实际上反映了去噪声后算法从原始图像中抽取的噪声，从这一图像中不难发现，它近于均匀噪声，没有边缘信号混杂其中的痕迹，即图像边缘没有显现。图5d就是原始图像的频谱与去噪后图像频谱的差，反映了去噪过程去掉的频谱。图5e和5f分别是图5a和图5b的第216列线图比较。从图中可以发现图5f线图的噪声比图5e上的少很多，在光栅处又保持完好。也说明了本方法具有去噪声而保护图像细节的能力。

2、实际人头磁共振图像

采用永磁0.4特斯拉的开放成像系统的T2加权图像，尺寸256×256×12。图像除噪声前后，整卷图像的信噪比增强了2.16倍左右。现将第1片图像分析进行分析，如图6所示。图6a和6b分别为去噪声前后的图像，6c为去前后图6a和6b的误差图像，图6d为6c的频谱幅图像。从图6中可以明显看出去后图6b比去噪声前图6a更加清晰，它们的差图像说明这种频谱重构去法有相当好的去噪声效果并没有伤及图像细节，图6d是差图像在频域上显示，说明了高频中的噪声分量大大减弱。

采用了上述的基于重构信号替代频谱数据的信号去噪方法，用于实际磁共振图像去噪，其方法步骤是：首先从对观测信号提取奇异点和奇异值，构造去噪函数和置换函数，最后由观测信号与去噪函数卷积积分、置换函数与奇异函数卷积积分的加权和共同得到去噪声后的信号。因而这种方法有效地避免了信号除噪过程中的高频信号丢失或信号失真，同时具有其它信号去噪方法法的优越之处，并巧妙克服了单一观测信号去噪法的信号失真问题，从而在保证图像所具有的高分辨率和高精确度条件下，而且能够有效去除图像噪声，提高了信噪比，为医学核磁共振成像检测提供了高质量的可靠图像信息；同时，本发明的方法高效实用，工作性能稳定可靠、适用范围较为广泛，给信号采集和处理带来了很大的便利，并且也为医学图像处理技术的进一步发展和大范围普及应用奠定了坚实的理论和实践基础。

在此说明书中，本发明已参照其特定的实施例作了描述。但是，很显然仍可以作出各种修改和变换而不背离本发明的精神和范围。因此，说明书和附图应被认为是说明性的而非限制性的。

Claims (5)

1.一种基于重构信号替代频谱数据的信号去噪方法，其特征在于，所述的方法包括以下步骤：

(1)从观测信号g(x)中提取奇异点和奇异值；

(2)根据得到的奇异点和奇异值构造相应的奇异函数；

(3)根据系统预设的频谱置零阈值，构造相应的去噪函数和置换函数；

(4)根据上述的观测信号g(x)、去噪函数、置换函数和奇异函数，利用重构信号替代法得到去噪后的信号。

2.根据权利要求1所述的基于重构信号替代频谱数据的信号去噪方法，其特征在于，所述的从观测信号g(x)中提取奇异点和奇异值，包括以下步骤：

(11)对观测信号g(x)进行差分；

(12)在上述差分中取差分绝对值大于系统预设的噪声阈值的那些点作为奇异点b₁，b₂，...，b_Q，其中Q为奇异点的数量；

(13)取所述的奇异点b₁，b₂，...，b_Q的差分值作为相应的奇异值a₁，a₂，...，a_Q。

3.根据权利要求2所述的基于重构信号替代频谱数据的信号去噪方法，其特征在于，所述的根据得到的奇异点和奇异值构造相应的奇异函数，具体为：

根据所述的奇异点b₁，b₂，...，b_Q和相应的奇异值a₁，a₂，...，a_Q，通过以下公式得到奇异函数：

$\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_q{w}_{b_q}\left(x\right);$

其中，w_bq(x)为奇异点b_q的奇异函数， $w_{b_q}\left(x\right)=w(x-x_{b_q})=\left\{\begin{array}{cc}1,& x{\&GreaterEqual;x}_{b_q}\\ 0,& x<x_{b_q}\end{array}\right..$

4.根据权利要求3所述的基于重构信号替代频谱数据的信号去噪方法，其特征在于，所述的根据系统预设的频谱置零阈值构造去噪函数和置换函数，包括以下步骤：

(31)根据以下公式得到观测信号g(x)的频谱函数G(k)：

G(k)＝DFT[g(x)]；

其中DFT[·]为离散付里叶变换算子；

(32)将观测信号频谱幅度值|G(k)|与系统预设的频谱置零阈值T进行比较，并根据以下条件得到高信噪比频谱区域Ω和其余频谱区域

如果|G(k)|≤T，则k∈Ω；

如果|G(k)|＜T，则 $k\&Element;\stackrel{\&OverBar;}{\mathrm{\&Omega;}};$

(33)根据高信噪比频谱区域Ω和其余频谱区域

通过以下公式得到指示函数R(k)和

$R\left(k\right)=\left\{\begin{array}{cc}1,& k\&Element;\mathrm{\&Omega;}\\ 0,& k\&NotElement;\mathrm{\&Omega;}\end{array}\right.;$

$\stackrel{\&OverBar;}{R}\left(k\right)=\left\{\begin{array}{cc}1,& k\&Element;\stackrel{\&OverBar;}{\mathrm{\&Omega;}}\\ 0,& k\&NotElement;\stackrel{\&OverBar;}{\mathrm{\&Omega;}}\end{array}\right.;$

(34)根据以下公式得到去噪函数r(x)和置换函数

r(x)＝IDFT[R(k)]；

$\stackrel{\&OverBar;}{r}\left(x\right)=\mathrm{IDFT}\left[\stackrel{\&OverBar;}{R}\left(k\right)\right];$

其中，IDFT[·]为离散付里叶反变换算子。

5.根据权利要求4所述的基于重构信号替代频谱数据的信号去噪方法，其特征在于，所述的利用重构信号替代法得到去噪后的信号，具体为：

根据观测信号g(x)、去噪函数r(x)、置换函数

和奇异函数 $\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_q{w}_{b_q}\left(x\right),$ 通过以下公式得到去噪后的信号g_r(x)：

$g_r\left(x\right)=g\left(x\right)*r\left(x\right)+\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_q{w}_{b_q}\left(x\right)*\stackrel{\&OverBar;}{r}\left(x\right);$

其中，“*”为卷积积分运算符。

Images