迭代盲目反卷积算法在X射线衍射谱中的应用

曹玉林, 马建萍

(1. 陕西师范大学 计算机科学学院, 陕西 西安 710062; 2. 青海师范大学 计算机学院, 青海 西宁 810008; 3. 青海师范大学 数学系, 青海 西宁 810008)

X射线衍射[1](XRD)图谱往往不是由一条条衍射线组成,而是由具有一定宽度的衍射峰组成,很多因素会导致衍射线变宽.首先,晶粒的细化及微观应变能够引起衍射线形增宽,即为物理宽化.再者,还存在着由于X射线源的几何尺寸、入射线发散及平板样品聚焦不良,以及接收狭缝的大小和衍射仪调整精度等原因而产生衍射线的宽化,即所谓的仪器宽化.在X射线衍射仪[2]上所观测到的样品的衍射线是上述两类宽化的合成.理论分析表明,这种合成并非是简单叠加,而是几个互不相关效应的卷积[3].通过对观测到的综合衍射线进行解卷处理,去除仪器引起的宽化就能够求得反映结构缺陷的物理宽化函数,从而对结构中各种形式的缺陷进行研究.微细晶粒的平均大小、粒度分布、微观应力、结构面的堆垛层错等信息,都能通过对物理宽化函数的分析得到.

为了去除物理宽化,提高分辨率,人们提出在Jade上使用标准样品制作一条半高宽补正曲线,再进行解卷处理[4-5].这种方法所使用的标准样品是完全退火态样品[6].该方法有两大缺陷:一方面半高宽补正曲线因仪器而异,往往随着仪器的使用环境、年限等变化,仪器宽化函数会不断的发生变化,导致此方法缺乏通用性;另一方面要保证标准样品的完全退火态和试验条件的一致性不容易.基于上述缺陷人们提出了另外一种解决方法,即盲目反卷积方法.常规的反卷积如Weiner滤波、Jansson迭代等方法必须要求卷积核函数(文中指仪器宽化函数)已知,而通过理论分析或测量等手段准确获得该函数并非易事.而盲目反卷积却能通过对观测信号本身的处理同时获得更接近真实的衍射谱和仪器宽化函数.文献[6-7]提出了迭代盲目反卷积算法,并通过对光谱图的处理,验证了该算法的可行性.而杨怀栋等[8]对光谱图迭代盲目反卷积算法的具体实现进行了论述.然而这些研究成果应用于X射线衍射图谱领域的文献屈指可数,文献[4,9]利用盲目反卷积算法对XRD图谱进行了研究,但也只是基于峰度的.本文尝试将迭代盲目反卷积算法应用于XRD图谱处理,以提高图谱分辨率,对有结构缺陷的晶面进行晶粒大小和微观应变等方面的研究有借鉴意义.

1 基本原理

1.1频域Weiner滤波器本文中所用的Weiner滤波器是离散非因果滤波器[10],其基本原理框图如图1所示.

给定观测序列y(n),它是一个非因果系统的

图1 Weiner 滤波器的基本原理框图Fig. 1 The basic principle diagram of Weiner filter

输出

(1)

式中,ξ(n)是零均值白噪声,h(n)是单位脉冲响应.希望找到一个非因果滤波器g(n),它用y(n)作输入,使其输出

满足

(2)

y(n-m)]=0, ∀m.

假定x(n)和ξ(n)都是广义平稳的,y(n)也就是广义平稳的,由此得到

这个式子的离散时间傅里叶变换给出

Sxy(ω)=Syy(ω)G(ω),

式中,G(ω)是g(n)的离散时间傅里叶变换(DTFT),Sxy(ω)和Syy(ω)分别是交叉功率谱和自功率谱.Weiner滤波器可以表达为

G(ω)=Sxy(ω)/Syy(ω).

另一方面,由(1)式可以证明

Syy(ω)=|H(ω)|2Sxx(ω)+Snn(ω),

Sxy(ω)=H*(ω)Sxx(ω),

Sxx(ω)和Snn(ω)分别是输入信号和噪声的(自)功率谱.

于是Weiner滤波器及其给出的估计为

(3)

由此可见,为了能利用Weiner滤波器作信号复原,必须知道卷积核h(n)以及信号和噪声的功率谱.如果2个功率谱难以得到,对(3)式的一个近似是

式中,γ是观测序列的信噪比的倒数,它起到规整化的作用,用以消除核函数的频域奇异性造成的病态问题.在用此滤波器时,信噪比只能通过原始的XRD图谱数据,利用公式SNR=Es/En进行估计.一个经验方法是,填零扩展y(n)到尺寸m,计算离散傅里叶变换(DFT)获得其频谱|Y(ω)|.根据|Y(ω)|可以确定一个频率点kf,使得信号能量完全包含在0～kf的范围内.用下式估计信号和噪声能量

1.2一维搜索方法当采用数学规划法寻求多元函数的极值点时,一般要进行一系列如下格式的迭代计算

xk+1=xk+αksk,k=0,1,2,…

当方向sk给定,求最佳步长αk就是求一元函数

f(xk+1)=f(xk+αksk)=φ(αk)

的极值问题,这一过程被称为一维搜索[11].一维搜索的方法主要分为两大类:解析法和数值法.解析法要进行求导,这对于函数形式复杂的十分不便,尤其是表达式未知的离散序列就更无能为力了,如XRD图谱数据.而数值法的应用就没这样的限制,其基本思路是,确定所要求的极值点α*的搜索区间,在不断缩小区间,最终获得近似值.数值法也有很多种,在此只介绍所用的黄金分割法.

黄金分割法的搜索过程：

1) 给出初始搜索区间[a,b]及收敛精度ε,将λ赋以0.618;

2) 按坐标点计算公式α1=b-λ(b-a),α2=a+λ(b-a),计算α1和α2,并计算其对应的函数值.比较函数值f1和f2,并缩短搜索区间;

(i) 若y1≤y2,丢去[α2,b],取[a,α2].α2→b,α1→α2,取新点α1=b-λ(b-a);

(ii) 若y1>y2,丢去[a,α1],取[α1,b].α1→a,α2→α1,取新点α2=a+λ(b-a);

3) 判断迭代终止条件.

2 算法实现

本文中的迭代盲目反卷积算法并不要求核函数已知,而是迭代的估计真实信号和卷积核.一般盲目反卷积的解不是唯一的[12],但对于XRD图谱,有重要的先验限制能够使迭代盲目反卷积算法收敛到正确的解.第一,是正性;第二,是核函数接近于某一钟罩函数,并有一个可猜测的大致范围.

设y(n)为XRD图谱数据,h(n)为仪器宽化函数(即核函数),x(n)为真实XRD图谱.算法步骤如下：

步骤1从XRD图谱中选择最狭窄的衍射峰,并选择一钟罩函数对其进行拟合,作为一个仪器宽化函数的起始猜测h0(n);

步骤2运用Weiner滤波器对y(n)和h0(n)作反卷积处理,求估计的XRD图谱序列x0(n);

步骤3再对y(n)和x0(n)作反卷积,求改进的核函数h1(n);

步骤4重复步骤2和3以改进x(n)和h(n)的估计.

下面就算法的实现作如下几点说明：1) 仪器宽化函数的选择如下:

%该程序用于找到拟合离散度最小的钟罩函数

%h为衍射峰分布数据,m1、m2、m3分别为3个钟罩函数的参数数组

%返回值k=1时,选高斯；k=2时,选柯西；k=3时,选双柯西.

function k=dispersion(h,m1,m2,m3)

h=h’;

N=length(h);

n=1:N;

h1=m1(1).*exp(-m1(2).*(n-m1(3)).^2);

h2=m2(1)./(1+m2(2).*(n-m2(3)).^2);

h3=m3(1)./((1+m3(2).*(n-m3(3)).^2).^2);

d1=sum((h1-h).^2)/N;

d2=sum((h2-h).^2)/N;

d3=sum((h3-h).^2)/N;

if d1

k=1;

elseif d2

k=2;

elseif d3

k=3;

end

2) 在编写Weiner滤波器程序时,要对计算出的图谱施以正性限制.

3) 步骤3所采取的做法是将反卷积问题转化为参数优化问题.假定核函数为高斯函数,则有

式中,a为归一化系数,它使得

由此

于是反卷积问题可以转化为单个参数b的寻优,使得误差函数

E=‖y(n)-h(n)x(n)‖2

达到最小.在这里使用一维搜索方法进行寻优,其程序实现如下:

function c=search(dir1,dir2,a,b,aqr)

f1=fopen(dir1,‘r’);

fseek(f1,25,0);

p=importdata(dir1);

f2=fopen(dir2,‘r’);

fseek(f2,25,0);

q=importdata(dir2);

y=p(:,2);

x=q(:,2);

namta=0.618;

t1=b-namta*(b-a);k1=Error_cauchy(y,x,t1);

t2=a+namta*(b-a);k2=Error_cauchy(y,x,t2);

while abs((b-a)/b)>aqr||abs((k2-k1)/k2)>aqr

if k1>=k2

a=t1;t1=t2;k1=k2;

t2=a+namta*(b-a);k2=Error_cauchy(y,x,t2);

else

b=t2;t2=t1;k2=k1;

t1=b-namta*(b-a);k1=Error_cauchy(y,x,t1);

end

end

c=(a+b)/2;

4) 该算法的迭代终止条件为核函数的特征参数b的相对误差即

|(b(k+1)-b(k))/b(k)|.

实验方法表明误差品质

的曲线拐点也可以用来指示迭代的终结,其中,hinit是核函数h的起始猜测.故可以用误差品质曲线的拐点来验证上述迭代终止条件的合理性,亦可验证该算法所得结果的正确性.

3 实验结果与分析

用Y-4Q型X射线衍射仪分别对SiO2粉末和藏药银矿石粉末进行了测试.测试条件为：步宽0.01°,管压30 kV,管流20 mA ,铜靶的Kα辐射,滤波片为镍,扫描速度为0.03(°)/s,连续扫描,时间常数0.5 s,扫描范围10°～90°.

实验得到的石英石和银矿石XRD图谱不能直接用于盲目反卷积算法的数值试验.因为实验得到的衍射峰是Kα1和Kα2的双重线,在低角度区域它们严重重叠,而仅在高角度区域才能分离,Kα1、Kα2的重叠会妨碍求算单一波长的剖面.故在进行数值实验前,需对得到的图谱应用Jade 5.0软件作扣除背景和分离Kα2的处理.

首先对石英石和银矿石图谱中最细锐的峰进行作非线性拟合,并进行拟合离散度的计算,如图2和图3所示.

然后,用MATLAB编写的基于迭代盲目反卷积的程序对两图谱进行解卷积处理.图4和图5为石英石和银矿石原始图谱及部分迭代估计图谱.

现进行如下分析：

1) 在数值实验中分别用高斯函数、柯西函数、双柯西函数分别对石英石和银矿石原始图谱中最细锐谱线进行拟合,如图2和图3所示,并作拟合离散度计算.对于石英石各函数拟合离散度分别为1 884.4、6 670.1和2 763.5,对于银矿石分别为333.402 8、715.988 1和441.419 2.因此将石英石和银矿石所对应的仪器宽化函数类型定为高斯函数型.

2) 分离出Kα2的石英石和银矿石图谱经迭代盲目反卷积算法处理后,在一些衍射峰附近分离出一些细小的峰,而一些峰被清楚地分开,如石英石图谱中位于2θ=26.530°的衍射峰,及银矿石图谱中位于2θ=24.753°的衍射峰,如图4和图5所示.这是因为经过算法处理后,衍射峰剖面变窄,使得其附近不明显的重叠峰凸显出来,即分辨率提高了.

3) 由图4和图5可以看出,两图谱的第一次迭代估计和原始图谱相比在线形上有很大的变化,但第一次迭代估计和其余迭代估计相比,表面上看没有太大的区别,但实际上衍射峰的半高宽却变窄了,如表1所示.各次迭代估计围绕着最终迭代估计来回摆动并最终收敛于最终迭代即真实图谱.

4) 若仪器宽化函数的起始猜测和真实的仪器宽化函数很相近,收敛将进行很快；若猜测的差距很大,则迭代次数也会随之增大.如对于石英石XRD图谱,迭代进行了16次,而银矿石只有8次.图6和图7为核函数(仪器宽化函数)的迭代估计情况,石英石和银矿石图谱卷积核函数最终的迭代估计很接近,但仍有差别,这是因为2次实验毕竟不能保证完全一样,必然会有一定的误差存在.

5) 从表1可以看出迭代盲目反卷积算法能很大程度上改善XRD图谱的分辨率,从而去除仪器所引起的宽化.

表 1 应用盲目反卷积算法处理前后半高宽的比较及分辨率变化

4 结语

利用XRD图谱正性和仪器宽化函数类型已知等先验知识,将迭代盲目反卷积算法用于X射线衍射图谱,数值试验表明迭代盲目反卷积算法能够有效提高XRD图谱的分辨率,去除宽化,得到真实的XRD图谱.

致谢衷心感谢审稿专家及编辑为本文修改提出了宝贵建议,谨致谢意.

[1] Takehira R, Momose Y, Yamamura S. Quantitative analysis of crystalline pharmaceuticals in tablets bypattern-fitting procedure using X-ray diffraction pattern[J]. Int J Pharm,2010,398(1/2):33-38.

[2] 黄继武. MDI Jade使用手册:X射线衍射实验操作指导[M]. 长沙:中南大学出版社,2006:31-35.

[3] 中国国家标准化管理委员会. GB/T 23413-2009,纳米材料晶粒尺寸及微观应变的测定X射线衍射线宽化法[S]. 北京:中国标准出版社,2009.

[4] 胡自强,袁景和,徐宝龙,等. X射线粉末衍射谱的峰度盲目反卷积[J]. 核电子学与探测技术,2007,27(6):1156-1158.

[5] Fiori S. Geodesic-based and projection-based neural blind deconvolution algorithms[J]. Signal Processing,2008,88(3):521-538.

[6] 邹谋炎. 反卷积与信号复原[M]. 北京:国防工业出版社,2001:93-95,151-155.

[7] Vural C, Sethares W A. Blind image deconvolution via dispersion minimization[J]. Digital Signal Processing,2006,16(2):137-148.

[8] 杨怀栋,徐立,陈科新,等. 盲目反卷积光谱图超分辨复原算法[J]. 光谱学与光谱分析,2007,27(7):1249-1253.

[9] 李玲,胡学刚,蒋伟. 一种基于LIP的全变分图像去噪新模型[J]. 四川师范大学学报:自然科学版,2011,34(2):134-138.

[10] 杨丽娟,张白桦,叶旭桢. 快速傅里叶变换FFT及其应用[J]. 光电工程,2004,31(S1):1-3.

[11] 张莹. 一维搜索的程序实现[J]. 沈阳教育学院学报,2002,4(4):104-106.

[12] Vural C, Sethares W A. Blind image deconvolution via dispersion minimization[J]. Digital Signal Processing,2006,16(2):145-148.

[13] 卓金武,魏永生,秦健,等. MATLAB在数学建模中的应用[M]. 北京:北京航空航天大学出版社,2011:5-15.

[14] Branicki M, Gershgorin B, Majda A J. Filtering skill for turbulent signals for a suite of nonlinear and linear extended Kalman filters[J]. J Comput Phys,2012,231(4):1462-1498.