CN102012466A - 数字x射线成像系统的噪声测量方法 - Google Patents

Abstract

本发明属于生物医学工程及计算机领域，涉及一种数字X射线成像系统的噪声测量方法：包括下列步骤：设置数字放射成像系统的曝光参数，放置脂肪仿体，采集图像；设置噪声功率谱研究的感兴趣区域，分别对各幅采集图像的感兴趣区域进行二维经验模态分解，提取各自对应的背景趋势，而后将各自的对应趋势去除，得到各幅图像的去趋势图像；计算所有去趋势图像感兴趣区域的平均图像，依次从各幅去趋势图像中去除该平均图像，获取各幅相应的噪声图像；根据空间分辨率要求进行图像分块，计算所有子块噪声功率谱并进行叠加平均，分别获取不同分辨率下数字X射线成像系统的二维噪声功率谱；选取二维噪声功率谱频率轴上下各7条谱线，叠加平均以获取不同分辨率下数字X射线成像系统的一维噪声功率谱。本发明为正确识别系统随机噪声特性，有效评估放射成像系统性能提供有力支持。

Description

数字X射线成像系统的噪声测量方法

技术领域

本发明属于生物医学工程及计算机领域，涉及一种数字X射线成像系统的噪声测量方法。

背景技术

噪声功率谱是成像系统性能评估的关键技术指标，已被广泛应用于放射影像系统成像质量定量评价的临床实践和研究，精确测量噪声功率谱是实现成像系统准确评估的关键。近年来随着数字影像设备的快速发展，针对数字X射线成像系统的噪声评估已经成为国内外放射影像学研究的重点内容之一。实现噪声功率谱测量的重点和关键在于通过有限的图像样本获取准确光滑的噪声功率谱，并且最大程度的保证其频率分辨率。

到目前为止，针对数字X射线成像系统实现噪声功率谱测量的技术难点主要有以下两方面。第一，理论上要获得真实的噪声功率谱必须通过测量无限数量的成像图像样本，但是实际进行噪声功率谱测量时只能获取有限数量的图像数据。因此，为了通过有限的样本数据来得到当前条件下最准确的噪声功率谱估计，必须根据实际情况采用恰当的测量及评估手段。比如对于二维X射线图像数据，在实际评估噪声功率谱时，必须处理好图像尺寸以及图像分块数之间的关系，兼顾噪声功率谱准确度与分辨率方面的要求，Dobbins、Mark等学者在近年来的研究中对此作了大量的工作。第二，成像系统评估针对的是系统随机噪声，但是实际噪声图像还包括了一些静态伪迹。Mark、Dobbins以及Samei等专家在研究中发现这些静态伪迹作用于傅立叶功率谱的低频部分，却不反映实际需要测量的系统随机噪声，最终导致噪声功率谱的低频成份产生系统性的上升，不利于真实噪声功率谱的准确评估。

目前，针对有限图像样本与噪声功率谱理论定义之间的矛盾，临床实践和研究已经展开了大量的研究并取得了可喜的成果。而对于第二个技术难点，即噪声中的静态伪迹问题，研究进展缓慢。如何有效去除噪声功率谱低频成份的系统性上升，准确实现数字X射线成像系统的噪声测量，实现噪声精确评估，已经成为成像系统性能评估的研究重点，这对于成像系统的发展应用具有十分重要的意义。

发明内容

本发明的主旨是克服现有技术的上述不足，提出一种数字X射线成像系统的噪声测量方法，以实现数字X射线成像系统的噪声测量，为进一步全面评估放射成像系统性能提供有力条件。本发明的技术方案如下：

一种数字X射线成像系统的噪声测量方法：包括下列步骤：

①设置数字放射成像系统的曝光参数，放置脂肪仿体，采集图像；

②设置噪声功率谱研究的感兴趣区域，分别对各幅采集图像的感兴趣区域进行二维经验模态分解，提取各自对应的背景趋势，而后将各自的对应趋势去除，得到各幅图像的去趋势图像；

③计算所有去趋势图像感兴趣区域的平均图像，依次从各幅去趋势图像中去除该平均图像，获取各幅相应的噪声图像；

④根据空间分辨率要求进行图像分块，计算所有子块噪声功率谱并进行叠加平均，分别获取不同分辨率下数字X射线成像系统的二维噪声功率谱；

⑤选取二维噪声功率谱频率轴上下各7条谱线，叠加平均以获取不同分辨率下数字X射线成像系统的一维噪声功率谱。

本发明解决了数字X射线成像系统的噪声功率谱估计中低频静态伪迹问题，从而为正确识别系统随机噪声特性，有效评估放射成像系统性能提供有力支持。采用本发明方法显著抑制了噪声评估结果中噪声功率谱低频成份的系统性上升，大大降低了功率谱估计低频成份的不确定性(即功率谱方差)，因此能有效抑制背景静态伪迹对系统噪声评估的影响，从而实现了数字X射线成像系统的噪声特性精确测量。本发明的应用，将为放射影像医学的临床实践和研究提供有效地技术指导，并为实际影像诊断时成像设备参数设置以及影像设备设计提供技术支持，因此该发明的实现具有广泛的应用前景和重要的社会意义。

附图说明

图1.(a)去除背景趋势前，15副图像ROI的二位噪声功率谱(b)背景趋势的功率谱密度曲线。

图2.背景趋势去除前后，噪声功率谱曲线及其不确定度比较。谱估计分析时采用15副图像ROI区域(256×256)。(a)噪声功率谱对比(b)噪声功率谱方差对比。

图3.在图像ROI区域分成2×2子块的情况下，背景趋势去除前后的噪声功率谱曲线及其不确定度比较。谱估计分析时采用60个图像子块(128×128)。(a)噪声功率谱对比(b)噪声功率谱方差对比。

图4.在图像ROI区域分成4×4子块的情况下，背景趋势去除前后的噪声功率谱曲线及其不确定度比较。谱估计分析时采用240个图像子块(64×64)。(a)噪声功率谱对比(b)噪声功率谱方差对比。

具体实施方式

本发明针对数字X射线成像系统性能评估中噪声功率谱估计的低频静态伪迹问题，以及目前相关技术不能依据不同成像图像实现自适应伪迹去除，也没有关于伪迹去除前后噪声功率谱定量性能评估。为此，本发明提出应用扩展的二维经验模态分解方法提自适应取放射曝光图像的低频背景趋势，以去除低频静态伪迹的图像作为噪声功率谱测量估计对象，从而正确识别系统随机噪声特性，评估系统成像性能。最后得到的技术方案如下：①设置数字放射成像系统的曝光参数，放置脂肪仿体，采集图像；②设置噪声功率谱研究的感兴趣区域，分别对各幅采集图像的感兴趣区域进行二维经验模态分解，提取各自对应的背景趋势，而后将对应趋势去除；③计算所有去趋势图像感兴趣区域的平均图像，依次从各幅去趋势图像中去除该平均图像，获取各幅相应的噪声图像；④根据空间分辨率要求进行图像分块，计算所有子块噪声功率谱并捷星叠加平均，分别获取不同分辨率下数字X射线成像系统的二维噪声功率谱。⑤选取二维噪声功率谱频率轴上下各7条谱线，叠加平均以获取不同分辨率下数字X射线成像系统的一维噪声功率谱，以此进行数字X射线成像系统的噪声特性评估。下面从几个方面对本发明的测量方法进行详细说明。

1.噪声功率谱估计

噪声功率谱被定义为空间频率域上随机信号在每一个频率位置的方差。根据维纳-辛钦理论，通过计算自相关函数的傅立叶变换就可以获得噪声功率谱。但在实际计算中，噪声功率谱通常通过对噪声密度图像的二维傅立叶变换取模来直接获取，具体计算公式如下：

$\mathrm{NPS}(u,v)=\underset{X,Y\→\∞}{\lim }\frac{1}{X\·Y}{\left|{\∫}_{-X/2}^{X/2}{\∫}_{-Y/2}^{Y/2}f(x,y)\exp \right[-2\mathrm{\πi}(\mathrm{ux}+\mathrm{vy})\left]\mathrm{dxdy}\right|}^2=\underset{X,Y\→\∞}{\lim }\frac{1}{X\·Y}{\left|\mathrm{FT}\right[f(x,y)\left]\right|}^2---\left(1\right)$

这里FT[f(x，y)]是噪声图像f(x，y)的傅立叶变换，X和Y分别是图像在x方向和y方向的最大尺寸。

在数字成像系统中，上述的噪声功率谱的连续表达形式通常需要转换为离散形式进行计算。因此，需要重新计算连续噪声图像函数f(x，y)在一系列离散位置点的像素值，这些位置点是x＝mΔx，y＝nΔy，m＝0ΛM，n＝0ΛN，这里X＝MΔx，Y＝NΔy。这样就可以把连续图像表示为离散形式f(m，n)。相应的，空间频率变量u和v的函数也需要在离散频率位置重新计算，这些位置包括u＝kΔu，v＝lΔv，k＝0，±1，±2Λ±K，l＝0，±1，±2Λ±L。最大空间频率即奈奎斯特采样频率在两个频率轴方向分别为U＝1/2Δx及V＝1/2Δy。

为了估计数字成像系统地噪声功率谱，上述计算功率谱的连续表达式也需要转换为离散形式。根据上述定义，可以将公式(1)改写为：

${\mathrm{NPS}}_d(u,v)=\mathrm{NPS}(\mathrm{k\Δu},\mathrm{l\Δv})=\underset{\mathrm{M\Δx},\mathrm{N\Δy}\→\∞}{\lim }\frac{\mathrm{\Δx\Δy}}{\mathrm{MN}}{\left|F{T}_d\right[f(m,n)\left]\right|}^2---\left(2\right)$

这里FT_d[f(m，n)]是数字图像f(m，n)的二维傅立叶变换。

2.基于经验模态分解方法的静态伪迹去除

数字X射线系统成像结果中的静态伪迹可能表现为成像过程足跟效应或者是曝光量平方反比作用而导致的背景伪影，也可能表现为探测器结构因素造成的固定趋势。这些伪迹成份通常叠加在由量子数决定的随机噪声上，表现为低频背景趋势。一些学者正在开展低频背景趋势去除的研究工作，目前主要发展的技术有如下几种：(1)去除数字图像的低通滤波成份；(2)去除数字图像的二维一阶平面拟合；(3)去除数字图像的二维二阶曲面拟合。这几类方法在各类场合下都所应用，但是不同方法的选择对系统噪声最终评估影响很大，而到目前为止，尚无未见伪迹去除前后噪声功率谱定量评估以及如何根据成像条件选择伪迹去除方法的相关报道。而在实际临床研究与实践，如何恰当的去除静态伪迹，直接影响到图像质量的评估，以及进一步的图像恢复。因此，发展一套具有自适应特性的伪迹去除方法，针对不同成像条件实现放射成像系统的噪声特性准确估计，可以为当前影像医学的发展提供有力支持。

在传统的Fourier分析中，频率被定义为整个分析数据长度中具有一定幅度的正、余弦函数。受这种固有观念的影响，人们在认识和接受瞬时频率的意义和概念时，总是从正、余弦函数的有关角度来分析。这样当人们定义局部频率值时就需要多于一个周期的正、余弦波动，基于这个逻辑，少于一个周期长度的信号将无法给出其频率的定义。而对于非线性和非平稳信号来说，其主要特征频率是时变的，即仅仅是在某一局部时间内存在或曾在某一时刻出现过，在描述频率随时间的变化关系上，Fourier变换显然已经无能为力。为了弥补Fourier变换对时变信号分析的不足，人们对原始信号加窗，认为在某个“窄带”内的信号是平稳的或近似平稳的，然后再对窗内的信号进行分析，如短时傅立叶变换、小波分析等，这些方法不同程度上对非线性和非平稳信号的时变性进行了描述，大大改进了Fourier变换的不足。但是由于受Heisenberg不确定原理的制约，在时间和频率上的分辨率不能同时达到最小，因此，所得的结果是窗内信号的平均结果，同样也摆脱不了Fourier变换的局限性。

为了把信号的频谱分析精确到每一个时间点、频率点上，美国国家航空航天局，美国工程院院士HuangN E等人提出了经验模态分解(Emprical Mode Decomposition，EMD)法。1996年，Huang在一次国际学术会议上首次提出这种适于非平稳信号分析的新方法的设想-基于经验的模式分解方法。Huang认为，对于瞬态与非平稳现象，频率与能量一般都是时间的函数，因此需要给出瞬时频率和瞬时能量的定义。目前信号瞬时能量的概念已被广泛接受，但是瞬时频率的概念和意义却一直存有争议。当可以使数据解析化的Hilbert变换出现之后，人们根据Hilbert变换提供的能够完全表达原始数据幅度和相位的函数，给出了瞬时频率的统一定义，从定义可以看出瞬时频率是时间的单值函数，即在任意时刻只存在一个振荡模式。所以在使用瞬时频率这一概念时，对应的数据受到了一定的限制。这主要因为任何一个时刻，数据中可能包含多个振荡模式，此时Hilbert变换不能给出该信号完全的频率内容，所得到的结果只是多个振荡模式的平均效果，从而瞬时频率的意义变得模糊。为了从复杂信号中得到有意义的瞬时频率，Huang根据瞬时频率物理意义上的必要条件，提出把含有多个振荡模式的数据分解为满足一定条件的多个单一振荡模式分量的线性叠加，每个单一振荡模式分量又叫做一个基本模式分量，并提出了一种基于经验的模式分解方法。每一个单一模式分量都满足Hilbert变换的必要条件，使得用Hilbert变换求解信号的瞬时频率成为可能。

经验模态分解法的意义在于：在信号分析中信号频率的定义是基于波形的局部特征和瞬时特征，它能在信号数据的每个时间点上，从点与点之间的变化特征来给出瞬时频率值，而不是需要多个振荡周期的波形才能给出一个频率值。经验模态分解法分析中若存在一个频率，仅仅表示该频率对应的信息在某一局部时间内存在或者曾在某一时刻出现过。这样无论从概念上还是从信号分析本质上来看，这种分析方法打破了传统频率思想，给出了一个全新的频率概念。经验模态分解法对信号分析具有重大的意义，同时也为非平稳信号处理提供了新思路，开辟了新途径。

为了能把一般数据分解成固有特征振荡模式，Huang N E提出了经验模态分解的方法。和以前几乎所有的分解方法不同，该新方法是直观的、直接的、后验的和自适应的，其分解的基函数立足于数据并且来源于数据本身。

经验模态分解方法是建立在以下的假设之上的：(1)信号至少有两个极值点，一个最大值和一个最小值；(2)特征时间尺度是通过两个极值点之间的时间间隔来定义；(3)若数据缺乏极值点但有变形点，则可通过数据微分一次或几次获得极值点，然后再通过积分来获得分解结果。

这种方法的本质是通过数据的经验特征时间尺度来获得其内在的振荡模式，然后据此分解数据。根据Drazin的经验，数据分析的第一步是人工观察数据，有两种能直接区分不同尺度振荡模式的方法：观察依次交替出现的极大、极小值点间的时间间隔；和观察依次出现的过零点的时间间隔。交替的局部极值点和过零点形成了复杂的数据：一个波动骑在另一个波动上，同时它们又骑在其他的波动上，依此类推，每个波动都定义了数据的一个特征尺度，这个特征尺度是内在的。采取依次出现的极值点的时间间隔作为振荡模式的时间尺度，因为这个方法对振荡模式不但有更高的分辨率，而且能应用于非零均值的数据，例如没有过零点的，全部数据点是正的或负的数据。为了把各种振荡模式依次从数据中提取出来，使用一种系统的方法，即经验模态分解方法，或形象地称之为“筛选”的过程。对实信号s(t)进行经验模态分解的步骤为：

1)确定s(t)的所有极大值和极小值；

2)根据极大值和极小值作三次样条差值来构造s(t)的上下包络线；

3)根据上下包络线，计算出s(t)的局部均值(上下包络的平均值)m₁(t)，以及s(t)和m₁(t)的差值h₁(t)＝s(t)-m₁(t)；

4)以h₁(t)代替原始信号s(t)，重复以上三步骤k次，直到所得的平均包络趋于零为止(h_1，(k-1)(t)与h_1，k(t)之间的方差小于设定值)，即认为h_1，k(t)是一个IMF分量，记c₁(t)＝h_1，k(t)，r₁(t)＝s(t)-c₁(t)，s(t)＝r₁(t)；第一个IMF分量代表原始数据中最高频率的分量。将原始数据序列s(t)减去第一个IMF分量c₁(t)，可以得到一个去掉高频分量的差值数据序列r₁(t)。

5)对r₁(t)重复以上四步平稳化处理过程，可以得到第二个IMF分量c₂(t)，如此重复下去直到最后一个差值序列r_n(t)不可分解为止(r_n(t)小于一设定值，或者变成一个单调函数时)，原始信号的经验模态分解结束，得到s(t)的分解式如下：

$s\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_i\left(t\right)+r_n\left(t\right)---\left(3\right)$

其中r_n(t)为残余函数，代表原始数据的趋势或均值，

由于每一个IMF分量都是代表一组特征尺度的数据序列，因此这个平稳化的处理过程实际上是将原始数据序列分解为不同特征波动的叠加。需要说明的是，每一个IMF分量既可以是线性的也可以是非线性的。

经验模态算法实际是一个筛选的过程，首先将信号中频率最高的成份筛选出来，而后从原信号中将该成份去除，再从新的信号中选出频率最高的成份，依此类推，直到信号不可分解为止。这种过程可以看作为一系列的滤波器族，从算法的以上描述可知，信号可以通过有限的筛选步骤完成分解。每次筛选时，新信号的最大、最小值的个数都在减少。为了减少提取IMF的筛选步骤，定了了参数SD：

$\mathrm{SD}=\underset{0}{\overset{T}{\mathrm{\Σ}}}\frac{{|h_{1,(k-1)}\left(t\right)-h_{1,k}\left(t\right)|}^2}{h_{1,(k-1)}{\left(t\right)}^2}---\left(4\right)$

当SD小于某一常数时停止筛选，一般SD的值在0.2至0.3之间。另外在筛选过程中，由于该算法采用的是三次样条插值，所以当信号的极大值或者极小值的个数小于2时，停止筛选。这些有限的IMF经过希尔伯特变换后产生了以时间为参数的具有实际意义的瞬时能量和瞬时频率。这样可以在时域和频域同时对信号进行有效的分析，这种特性是以往以傅立叶变换为基础的信号分析方法所不具有的。

由于经验模态方法是依据数据本身的时域信息进行的时域分解，得到的IMF通常个数是有限的和平稳的，而且是具有实际意义的窄带信号，基于这些IMF分量进行的Hilbert变换其结果反映了真实的物理信息，而且基于经验模态的Hilbert变换得到的每个IMF的振幅和频率是随时间变化的，消除了经典频谱分析法中为反映非线性、非平稳过程而引入的多余无物理意义的简谐波。因此，其Hilbert谱也能够准确反映出信号能量、频率在空间或时间尺度上的分布。经验模态方法是基于信号的局部特征时间尺度实现分解的。与小波分析法相比，经验模态-Hilbert方法具有小波分析的全部优点，并克服了小波变换的非自适应性，因此这种基于经验模态的Hilbert频谱分析方法，在非线性和非平稳过程的分析中具有很高的应用价值。

从本质上说，经验模态是对数据在进行Hilbert变换之前作的一个预处理。通过经验模态，数据被分解为若干固有模态函数(IMF)的集合，每一个IMF刻画了信号的一个简单振荡模式。从表达形式上来看，一个IMF类似于Fourier分解表达式中的一个简谐振动，但是，它比简谐振动更加一般化。通过经验模态分解，可以准确获取信号的趋势项，该趋势项有效反映了信号的低频静态伪迹，基线漂移等。

为了实现基于经验模态分解方法的二维静态伪迹消除，我们需要把上述的一维经验模态分解方法扩展到二维。对于二维的信号如图像，二维经验模态分解可以将原图像分解成频率从高到低的有限个二维的固有模态函数和趋势图像。固有模态函数应该满足两个约束条件：首先是对称于原二维信号的局部均值并且自身均值为0；其次是它的极大值点都是正的，极小值点都是负的(没有骑波)。二维经验模态分解的算法步骤如下：

1)初始化，设R(m，n)为待分解的图，设定当前每次迭代筛选前后的图像之间的标准偏差值SD的阈值；

2)如果待分解图像R(m，n)单调或达到图像的分解层数，则算法停止，此时得到的图像R(m，n)即为所求得的背景趋势图像；否则，令H(m，n)＝R(m，n)；

3)对图像H(m，n)进行极值点求解，找出区域极大值点集和区域极小值点集；

4)分别对区域极大值点集和区域极小值点集进行平面插值，得出图像的上、下包络面，根据上、下包络面求出图像H(m，n)的均值M(m，n)；

5)令H(m，n)＝H(m，n)-M(m，n)，判断图像H(m，n)与迭代前的图像之间的标准偏差值小于所设定的阈值，不满足则转步骤3)，否则当前迭代过程停止；

6)令D(m，n)＝H(m，n)，迭代得到一个二维固有模态函数D(m，n)；

7)令R(m，n)＝R(m，n)-D(m，n)，转步骤2。

在上述算法中，极值点求解、平面插值和筛选的停止条件是算法的核心。经过J层二维经验模态分解后，最终的分解过程可以表示为

$I(m,n)=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{D}_j(m,n)+R_J(m,n),J\∈N---\left(5\right)$

其中，D_j是第j个二维固有模态函数，R_J是经过J层分解后的趋势图像。

求解二维图像信号极值点是找出区域极值点而不是局部极值点。局部极值点是指该点在其周围邻居点中是极大值或极小值，而区域极值点是指该点在其周围邻居及其连通的部分是极大值或极小值。可以很容易看出区域极值点是局部极值点，而反过来则不一定。例如，一个点周围是具有相同值的平面，这一点是局部极大值点，但这个平面上可能会有一个邻居点的值高于这个平面，这样，这个点就不是区域极大值点。求区域极值可以采用Vincent提出的形态学重构方法，这种方法在图像中查找区域的极大值或极小值点。

在平面插值时，为了防止二维经验模态分解的端点效应问题，本发明采用极值点对称延拓方法。首先将图像的4条边周围较近的极值点以图像的边为对称中心向外延拓。延拓的目的是为了分解数据的完整性，从而抑制图像的4条边周围的误差向内传播。

在进行平面插值时基于对极值点集进行Delaunay三角划分的方法，然后在划分后的三角形内采用三次插值的方法。基于三角划分和三次插值的方法有如下优点：首先，被插值的极值点(极大值或极小值)不一定局限于方形的网格，而图像中的极值点往往是离散的、不规则的；其次，该方法计算出来的插值平面是很光滑的，是真正意义上的二维平面插值。

因为二维经验模态分解过零点的数目是无法统计的，所以分解时可以用固有模态分量的约束条件作为筛选过程的停止条件，也可以用Cauchy-type收敛条件]和能量差跟踪的方法作为筛选过程的停止条件。本发明选用Cauchy-type收敛条件，并且取SD＝0.3。

3.数字X射线成像系统

实验成像系统为Pixarray 100小动物数字放射成像系统，由美国BIOPTICS公司制造。该系统的探测器为1024×1024的CCD阵列，像素大小为50μm×50μm，14级灰度。横向及纵向的空间分辨率均为每毫米20像素。X射线管的焦斑尺寸为50μm。实验中，X射线源的工作电压为33kVp，工作电流为0.5mA，X射线源到探测器的距离保持830mm。曝光时间设置为10秒，成像物体选取1.4cm厚的脂肪仿体。在上述系统设置条件下，连续对脂肪仿体曝光成像15幅，用于进一步的噪声功率谱计算及分析。

4.数字X射线成像系统的噪声评估算法

本发明的数字X射线成像系统的噪声精确评估方法的流程描述如下：

1)设置数字放射成像系统的曝光参数，放置脂肪仿体，连续采集15幅图像，并设置图像中央256×256区域作为噪声功率谱研究的感兴趣区域(Region of interest，ROI)；

2)分别对15幅脂肪仿体图像的ROI进行二维经验模态分解，获取各自对应的背景趋势，而后将背景趋势从15幅脂肪仿体曝光图像ROI中去除；

3)计算15幅去趋势ROI图像的平均ROI图像；

4)将15幅经过去除背景趋势的脂肪仿体图像分别减去步骤3获得的平均图像，以此获取15幅噪声图像；

5)按照噪声图像不分块，分成4块以及分成16块这三种情况，分别对所有子块计算噪声功率谱，将上述三种情况下的噪声功率谱叠加平均，分别获取不同分辨率下数字X射线成像系统的二维噪声功率谱。

6)从二维噪声功率谱图中心沿横轴向右至终点，取横轴上下各7行功率谱值，而后将这14行功率谱值叠加平均，计算不同分辨率下数字X射线成像系统的一维噪声功率谱。

下面结合具体实施例和附图对本发明做进一步介绍。

本发明采用美国BIOPTICS公司生产的Pixarray 100小动物数字放射成像系统。根据国际电子协会(International Electrotechnical Commission，IEC)标准中的建议，本发明选取曝光图像(1024×1024)中央区域的256×256区域作为图像感兴趣区域(ROI)来计算噪声功率谱。图1(a)给出了在背景趋势去除之前，利用15幅ROI图像获取相应的15幅噪声图像，而后分别计算噪声功率谱图，最后通过叠加平均所获取的二维噪声功率谱。这里，为了更清楚的显示结果，我们对二维功率谱进行了取对数操作，计算公式为log₁₀(1+NPS_d(u，v))。从图中可以看出，沿横轴及纵轴方向的噪声功率谱基本一致，因此，在下面的分析中，我们以沿横轴方向的一维噪声功率谱为例说明本发明方法的有益效果。在本发明中，我们从二维噪声功率谱图中心沿横轴向右至终点，取横轴上下各7行功率谱值，而后将这14行功率谱值叠加平均，计算不同分辨率下数字X射线成像系统的一维噪声功率谱。

我们对15幅脂肪仿体曝光图像ROI区域(256×256)分别进行二维经验模态分解，获取对应的背景趋势。图1(b)给出了15幅图像ROI区域(256×256)背景趋势的一维功率谱(横轴方向)，这里为了更好的显示结果，同样进行了取对数操作。

在获取15幅图像ROI区域的背景趋势之后，我们将其分别从对应的脂肪仿体曝光图像ROI区域(256×256)中去除，而后对背景趋势去除前后的噪声功率谱进行评估和比较。图2(a)给出了背景趋势取出前后，采用15幅图像ROI区域(256×256)最终获取的一维噪声功率谱曲线。从图中可以看出，经过背景趋势去除之后，噪声功率谱低频成份的系统性上升受到了显著压制。为了评估背景趋势去除前后噪声功率谱的准确性，我们计算了15条噪声功率谱估计曲线的方差衡量功率谱估计的“不确定度(uncertainty)”，计算结果如图2(b)所示。从图2(b)可以看出，经过背景趋势去除后，噪声功率谱低频成份谱估计的方差(不确定度)有了显著地降低，即去除背景趋势可以有效改善噪声功率谱估计的准确性。

我们将上述15幅噪声图像ROI区域(256×256)分割成互不重叠的2×2子块(子块大小为128×128)，而后分别对每幅噪声图像的4个子块计算噪声功率谱。将每幅噪声图像下4个子块对应的1维噪声功率谱曲线进行叠加平均，叠加平均后的噪声功率谱曲线即为当前噪声图像ROI(256×256)的噪声功率谱估计。依照上述方法，分别计算15幅噪声图像的噪声功率谱估计，而后将15个噪声功率谱估计曲线进行叠加平均，以获得在当前空间分辨率(128×128)下的噪声功率谱曲线。图3(a)与图3(b)分别给出了去背景趋势前后，采用60个噪声图像子块(128×128)计算获得的一维噪声功率谱以及相应的功率谱方差曲线。从图中可以看出，在背景趋势去除之前，通过将噪声图像分成2×2子块(128×128)，噪声功率谱低频成份的系统性上升以及所有频段噪声功率谱成份的方差(包括低频功率谱成份)都得到了降低。经过背景趋势去除之后，噪声功率谱低频成份的系统性上升以及低频功率谱成份的方差相对背景趋势去除前有了进一步降低。由此可知，在通过图像分块求解噪声功率谱的情况下，去除背景趋势可以进一步有效改善噪声功率谱估计的准确性。

为了进一步分析本发明的有益效果，我们将上述15幅噪声图像ROI区域(256×256)继续分割为互不重叠的4×4子块(子块大小为64×64)。依次对每幅噪声图像的16个子块进行噪声功率谱计算，而后将每幅图像下所有子块对应的1维噪声功率谱曲线进行叠加平均，叠加平均获得的噪声功率谱曲线即为当前噪声图像ROI区域的噪声功率谱估计结果。将上述获得的15个噪声功率谱估计曲线再进行叠加平均，从而获得当前空间分辨率(64×64)下的噪声功率谱曲线。图4(a)与图4(b)分别给出了去背景趋势前后，采用240个噪声图像子块(子块大小为64×64)计算获得的一维噪声功率谱以及相应的功率谱方差曲线。从图中可以看出，当没有进行背景趋势去除时，噪声功率谱低频成份的系统性上升以及对应方差在进一步图像分块后基本没有改善。而在背景趋势去除情况下，噪声功率谱低频成份的系统性上升以及低频功率谱成份的方差得到了进一步显著改善，并且相对于图3所示的2×2分块情况，改善结果更为显著。

最终结果表明，在进行数字X射线成像系统的噪声评估时，通过二维验模态分解方法去除曝光图像背景趋势，而后根据系统评估所需的空间分辨率，进行恰当的图像分块，可以获得相对背景趋势去除前更为准确的噪声功率谱估计。该方法可以有效降低由于低频静态伪迹引起的噪声功率谱低频成份系统性上升，改善噪声功率谱估计的不确定度，从而实现了数字X射线成像系统的噪声功率谱准确估计。该方法的应用，将为正确识别系统随机噪声特性，有效评估放射成像系统性能，深入开展放射影像学临床实践和研究提供有力支持。

Claims (2)

1.一种数字X射线成像系统的噪声测量方法，包括下列步骤：

第一步：设置数字放射成像系统的曝光参数，放置脂肪仿体，采集图像；

第二步：设置噪声功率谱研究的感兴趣区域，分别对各幅采集图像的感兴趣区域进行二维经验模态分解，提取各自对应的趋势图像，而后将各自对应的趋势去除，得到各幅图像的去趋势图像；

第三步：计算所有去趋势图像感兴趣区域的平均图像，依次从各幅去趋势图像中去除该平均图像，获取各幅相应的噪声图像；

第四步：根据空间分辨率要求进行图像分块，计算所有子块噪声功率谱并进行叠加平均，分别获取不同分辨率下数字X射线成像系统的二维噪声功率谱；

第五步：选取二维噪声功率谱频率轴上下各7条谱线，叠加平均以获取不同分辨率下数字X射线成像系统的一维噪声功率谱。

2.跟根据权利要求1所述的数字X射线成像系统的噪声测量方法，其特征在于，其中的第二步采用下面的方法得到图像的背景趋势图像：

1)设R(m，n)为待分解的图，设定当前每次迭代筛选前后的图像之间的标准偏差值SD的阈值；

2)如果待分解图像R(m，n)单调或达到图像的分解层数，则算法停止，此时得到的图像R(m，n)即为所求得的背景趋势图像；否则，令H(m，n)＝R(m，n)；

3)对图像H(m，n)进行极值点求解，找出区域极大值点集和区域极小值点集；

4)分别对区域极大值点集和区域极小值点集进行平面插值，得出图像的上、下包络面，根据上、下包络面求出图像H(m，n)的均值M(m，n)；

5)令H(m，n)＝H(m，n)-M(m，n)，判断图像H(m，n)与迭代前的图像之间的标准偏差值小于所设定的阈值，不满足则转步骤3)，否则当前迭代过程停止；

6)令D(m，n)＝H(m，n)，迭代得到一个二维固有模态函数D(m，n)；

7)令R(m，n)＝R(m，n)-D(m，n)，转步骤2。

Images