迭代盲目反卷積算法在X射線衍射譜中的應(yīng)用

曹玉林, 馬建萍

(1. 陜西師范大學(xué) 計(jì)算機(jī)科學(xué)學(xué)院, 陜西 西安 710062; 2. 青海師范大學(xué) 計(jì)算機(jī)學(xué)院, 青海 西寧 810008; 3. 青海師范大學(xué) 數(shù)學(xué)系, 青海 西寧 810008)

X射線衍射[1](XRD)圖譜往往不是由一條條衍射線組成,而是由具有一定寬度的衍射峰組成,很多因素會(huì)導(dǎo)致衍射線變寬.首先,晶粒的細(xì)化及微觀應(yīng)變能夠引起衍射線形增寬,即為物理寬化.再者,還存在著由于X射線源的幾何尺寸、入射線發(fā)散及平板樣品聚焦不良,以及接收狹縫的大小和衍射儀調(diào)整精度等原因而產(chǎn)生衍射線的寬化,即所謂的儀器寬化.在X射線衍射儀[2]上所觀測(cè)到的樣品的衍射線是上述兩類寬化的合成.理論分析表明,這種合成并非是簡(jiǎn)單疊加,而是幾個(gè)互不相關(guān)效應(yīng)的卷積[3].通過(guò)對(duì)觀測(cè)到的綜合衍射線進(jìn)行解卷處理,去除儀器引起的寬化就能夠求得反映結(jié)構(gòu)缺陷的物理寬化函數(shù),從而對(duì)結(jié)構(gòu)中各種形式的缺陷進(jìn)行研究.微細(xì)晶粒的平均大小、粒度分布、微觀應(yīng)力、結(jié)構(gòu)面的堆垛層錯(cuò)等信息,都能通過(guò)對(duì)物理寬化函數(shù)的分析得到.

為了去除物理寬化,提高分辨率,人們提出在Jade上使用標(biāo)準(zhǔn)樣品制作一條半高寬補(bǔ)正曲線,再進(jìn)行解卷處理[4-5].這種方法所使用的標(biāo)準(zhǔn)樣品是完全退火態(tài)樣品[6].該方法有兩大缺陷:一方面半高寬補(bǔ)正曲線因儀器而異,往往隨著儀器的使用環(huán)境、年限等變化,儀器寬化函數(shù)會(huì)不斷的發(fā)生變化,導(dǎo)致此方法缺乏通用性;另一方面要保證標(biāo)準(zhǔn)樣品的完全退火態(tài)和試驗(yàn)條件的一致性不容易.基于上述缺陷人們提出了另外一種解決方法,即盲目反卷積方法.常規(guī)的反卷積如Weiner濾波、Jansson迭代等方法必須要求卷積核函數(shù)(文中指儀器寬化函數(shù))已知,而通過(guò)理論分析或測(cè)量等手段準(zhǔn)確獲得該函數(shù)并非易事.而盲目反卷積卻能通過(guò)對(duì)觀測(cè)信號(hào)本身的處理同時(shí)獲得更接近真實(shí)的衍射譜和儀器寬化函數(shù).文獻(xiàn)[6-7]提出了迭代盲目反卷積算法,并通過(guò)對(duì)光譜圖的處理,驗(yàn)證了該算法的可行性.而楊懷棟等[8]對(duì)光譜圖迭代盲目反卷積算法的具體實(shí)現(xiàn)進(jìn)行了論述.然而這些研究成果應(yīng)用于X射線衍射圖譜領(lǐng)域的文獻(xiàn)屈指可數(shù),文獻(xiàn)[4,9]利用盲目反卷積算法對(duì)XRD圖譜進(jìn)行了研究,但也只是基于峰度的.本文嘗試將迭代盲目反卷積算法應(yīng)用于XRD圖譜處理,以提高圖譜分辨率,對(duì)有結(jié)構(gòu)缺陷的晶面進(jìn)行晶粒大小和微觀應(yīng)變等方面的研究有借鑒意義.

1 基本原理

1.1頻域Weiner濾波器本文中所用的Weiner濾波器是離散非因果濾波器[10],其基本原理框圖如圖1所示.

給定觀測(cè)序列y(n),它是一個(gè)非因果系統(tǒng)的

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

輸出

(1)

式中,ξ(n)是零均值白噪聲,h(n)是單位脈沖響應(yīng).希望找到一個(gè)非因果濾波器g(n),它用y(n)作輸入,使其輸出

滿足

(2)

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

假定x(n)和ξ(n)都是廣義平穩(wěn)的,y(n)也就是廣義平穩(wěn)的,由此得到

這個(gè)式子的離散時(shí)間傅里葉變換給出

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

式中,G(ω)是g(n)的離散時(shí)間傅里葉變換(DTFT),Sxy(ω)和Syy(ω)分別是交叉功率譜和自功率譜.Weiner濾波器可以表達(dá)為

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

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

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

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

Sxx(ω)和Snn(ω)分別是輸入信號(hào)和噪聲的(自)功率譜.

于是Weiner濾波器及其給出的估計(jì)為

(3)

由此可見(jiàn),為了能利用Weiner濾波器作信號(hào)復(fù)原,必須知道卷積核h(n)以及信號(hào)和噪聲的功率譜.如果2個(gè)功率譜難以得到,對(duì)(3)式的一個(gè)近似是

式中,γ是觀測(cè)序列的信噪比的倒數(shù),它起到規(guī)整化的作用,用以消除核函數(shù)的頻域奇異性造成的病態(tài)問(wèn)題.在用此濾波器時(shí),信噪比只能通過(guò)原始的XRD圖譜數(shù)據(jù),利用公式SNR=Es/En進(jìn)行估計(jì).一個(gè)經(jīng)驗(yàn)方法是,填零擴(kuò)展y(n)到尺寸m,計(jì)算離散傅里葉變換(DFT)獲得其頻譜|Y(ω)|.根據(jù)|Y(ω)|可以確定一個(gè)頻率點(diǎn)kf,使得信號(hào)能量完全包含在0～kf的范圍內(nèi).用下式估計(jì)信號(hào)和噪聲能量

1.2一維搜索方法當(dāng)采用數(shù)學(xué)規(guī)劃法尋求多元函數(shù)的極值點(diǎn)時(shí),一般要進(jìn)行一系列如下格式的迭代計(jì)算

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

當(dāng)方向sk給定,求最佳步長(zhǎng)αk就是求一元函數(shù)

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

的極值問(wèn)題,這一過(guò)程被稱為一維搜索[11].一維搜索的方法主要分為兩大類:解析法和數(shù)值法.解析法要進(jìn)行求導(dǎo),這對(duì)于函數(shù)形式復(fù)雜的十分不便,尤其是表達(dá)式未知的離散序列就更無(wú)能為力了,如XRD圖譜數(shù)據(jù).而數(shù)值法的應(yīng)用就沒(méi)這樣的限制,其基本思路是,確定所要求的極值點(diǎn)α*的搜索區(qū)間,在不斷縮小區(qū)間,最終獲得近似值.數(shù)值法也有很多種,在此只介紹所用的黃金分割法.

黃金分割法的搜索過(guò)程：

1) 給出初始搜索區(qū)間[a,b]及收斂精度ε,將λ賦以0.618;

2) 按坐標(biāo)點(diǎn)計(jì)算公式α1=b-λ(b-a),α2=a+λ(b-a),計(jì)算α1和α2,并計(jì)算其對(duì)應(yīng)的函數(shù)值.比較函數(shù)值f1和f2,并縮短搜索區(qū)間;

(i) 若y1≤y2,丟去[α2,b],取[a,α2].α2→b,α1→α2,取新點(diǎn)α1=b-λ(b-a);

(ii) 若y1>y2,丟去[a,α1],取[α1,b].α1→a,α2→α1,取新點(diǎn)α2=a+λ(b-a);

3) 判斷迭代終止條件.

2 算法實(shí)現(xiàn)

本文中的迭代盲目反卷積算法并不要求核函數(shù)已知,而是迭代的估計(jì)真實(shí)信號(hào)和卷積核.一般盲目反卷積的解不是唯一的[12],但對(duì)于XRD圖譜,有重要的先驗(yàn)限制能夠使迭代盲目反卷積算法收斂到正確的解.第一,是正性;第二,是核函數(shù)接近于某一鐘罩函數(shù),并有一個(gè)可猜測(cè)的大致范圍.

設(shè)y(n)為XRD圖譜數(shù)據(jù),h(n)為儀器寬化函數(shù)(即核函數(shù)),x(n)為真實(shí)XRD圖譜.算法步驟如下：

步驟1從XRD圖譜中選擇最狹窄的衍射峰,并選擇一鐘罩函數(shù)對(duì)其進(jìn)行擬合,作為一個(gè)儀器寬化函數(shù)的起始猜測(cè)h0(n);

步驟2運(yùn)用Weiner濾波器對(duì)y(n)和h0(n)作反卷積處理,求估計(jì)的XRD圖譜序列x0(n);

步驟3再對(duì)y(n)和x0(n)作反卷積,求改進(jìn)的核函數(shù)h1(n);

步驟4重復(fù)步驟2和3以改進(jìn)x(n)和h(n)的估計(jì).

下面就算法的實(shí)現(xiàn)作如下幾點(diǎn)說(shuō)明：1) 儀器寬化函數(shù)的選擇如下:

%該程序用于找到擬合離散度最小的鐘罩函數(shù)

%h為衍射峰分布數(shù)據(jù),m1、m2、m3分別為3個(gè)鐘罩函數(shù)的參數(shù)數(shù)組

%返回值k=1時(shí),選高斯；k=2時(shí),選柯西；k=3時(shí),選雙柯西.

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濾波器程序時(shí),要對(duì)計(jì)算出的圖譜施以正性限制.

3) 步驟3所采取的做法是將反卷積問(wèn)題轉(zhuǎn)化為參數(shù)優(yōu)化問(wèn)題.假定核函數(shù)為高斯函數(shù),則有

式中,a為歸一化系數(shù),它使得

由此

于是反卷積問(wèn)題可以轉(zhuǎn)化為單個(gè)參數(shù)b的尋優(yōu),使得誤差函數(shù)

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

達(dá)到最小.在這里使用一維搜索方法進(jìn)行尋優(yōu),其程序?qū)崿F(xiàn)如下:

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) 該算法的迭代終止條件為核函數(shù)的特征參數(shù)b的相對(duì)誤差即

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

實(shí)驗(yàn)方法表明誤差品質(zhì)

的曲線拐點(diǎn)也可以用來(lái)指示迭代的終結(jié),其中,hinit是核函數(shù)h的起始猜測(cè).故可以用誤差品質(zhì)曲線的拐點(diǎn)來(lái)驗(yàn)證上述迭代終止條件的合理性,亦可驗(yàn)證該算法所得結(jié)果的正確性.

3 實(shí)驗(yàn)結(jié)果與分析

用Y-4Q型X射線衍射儀分別對(duì)SiO2粉末和藏藥銀礦石粉末進(jìn)行了測(cè)試.測(cè)試條件為：步寬0.01°,管壓30 kV,管流20 mA ,銅靶的Kα輻射,濾波片為鎳,掃描速度為0.03(°)/s,連續(xù)掃描,時(shí)間常數(shù)0.5 s,掃描范圍10°～90°.

實(shí)驗(yàn)得到的石英石和銀礦石XRD圖譜不能直接用于盲目反卷積算法的數(shù)值試驗(yàn).因?yàn)閷?shí)驗(yàn)得到的衍射峰是Kα1和Kα2的雙重線,在低角度區(qū)域它們嚴(yán)重重疊,而僅在高角度區(qū)域才能分離,Kα1、Kα2的重疊會(huì)妨礙求算單一波長(zhǎng)的剖面.故在進(jìn)行數(shù)值實(shí)驗(yàn)前,需對(duì)得到的圖譜應(yīng)用Jade 5.0軟件作扣除背景和分離Kα2的處理.

首先對(duì)石英石和銀礦石圖譜中最細(xì)銳的峰進(jìn)行作非線性擬合,并進(jìn)行擬合離散度的計(jì)算,如圖2和圖3所示.

然后,用MATLAB編寫的基于迭代盲目反卷積的程序?qū)蓤D譜進(jìn)行解卷積處理.圖4和圖5為石英石和銀礦石原始圖譜及部分迭代估計(jì)圖譜.

現(xiàn)進(jìn)行如下分析：

1) 在數(shù)值實(shí)驗(yàn)中分別用高斯函數(shù)、柯西函數(shù)、雙柯西函數(shù)分別對(duì)石英石和銀礦石原始圖譜中最細(xì)銳譜線進(jìn)行擬合,如圖2和圖3所示,并作擬合離散度計(jì)算.對(duì)于石英石各函數(shù)擬合離散度分別為1 884.4、6 670.1和2 763.5,對(duì)于銀礦石分別為333.402 8、715.988 1和441.419 2.因此將石英石和銀礦石所對(duì)應(yīng)的儀器寬化函數(shù)類型定為高斯函數(shù)型.

2) 分離出Kα2的石英石和銀礦石圖譜經(jīng)迭代盲目反卷積算法處理后,在一些衍射峰附近分離出一些細(xì)小的峰,而一些峰被清楚地分開(kāi),如石英石圖譜中位于2θ=26.530°的衍射峰,及銀礦石圖譜中位于2θ=24.753°的衍射峰,如圖4和圖5所示.這是因?yàn)榻?jīng)過(guò)算法處理后,衍射峰剖面變窄,使得其附近不明顯的重疊峰凸顯出來(lái),即分辨率提高了.

3) 由圖4和圖5可以看出,兩圖譜的第一次迭代估計(jì)和原始圖譜相比在線形上有很大的變化,但第一次迭代估計(jì)和其余迭代估計(jì)相比,表面上看沒(méi)有太大的區(qū)別,但實(shí)際上衍射峰的半高寬卻變窄了,如表1所示.各次迭代估計(jì)圍繞著最終迭代估計(jì)來(lái)回?cái)[動(dòng)并最終收斂于最終迭代即真實(shí)圖譜.

4) 若儀器寬化函數(shù)的起始猜測(cè)和真實(shí)的儀器寬化函數(shù)很相近,收斂將進(jìn)行很快；若猜測(cè)的差距很大,則迭代次數(shù)也會(huì)隨之增大.如對(duì)于石英石XRD圖譜,迭代進(jìn)行了16次,而銀礦石只有8次.圖6和圖7為核函數(shù)(儀器寬化函數(shù))的迭代估計(jì)情況,石英石和銀礦石圖譜卷積核函數(shù)最終的迭代估計(jì)很接近,但仍有差別,這是因?yàn)?次實(shí)驗(yàn)畢竟不能保證完全一樣,必然會(huì)有一定的誤差存在.

5) 從表1可以看出迭代盲目反卷積算法能很大程度上改善XRD圖譜的分辨率,從而去除儀器所引起的寬化.

表 1 應(yīng)用盲目反卷積算法處理前后半高寬的比較及分辨率變化

4 結(jié)語(yǔ)

利用XRD圖譜正性和儀器寬化函數(shù)類型已知等先驗(yàn)知識(shí),將迭代盲目反卷積算法用于X射線衍射圖譜,數(shù)值試驗(yàn)表明迭代盲目反卷積算法能夠有效提高XRD圖譜的分辨率,去除寬化,得到真實(shí)的XRD圖譜.

致謝衷心感謝審稿專家及編輯為本文修改提出了寶貴建議,謹(jǐn)致謝意.

[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使用手冊(cè):X射線衍射實(shí)驗(yàn)操作指導(dǎo)[M]. 長(zhǎng)沙:中南大學(xué)出版社,2006:31-35.

[3] 中國(guó)國(guó)家標(biāo)準(zhǔn)化管理委員會(huì). GB/T 23413-2009,納米材料晶粒尺寸及微觀應(yīng)變的測(cè)定X射線衍射線寬化法[S]. 北京:中國(guó)標(biāo)準(zhǔn)出版社,2009.

[4] 胡自強(qiáng),袁景和,徐寶龍,等. X射線粉末衍射譜的峰度盲目反卷積[J]. 核電子學(xué)與探測(cè)技術(shù),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] 鄒謀炎. 反卷積與信號(hào)復(fù)原[M]. 北京:國(guó)防工業(yè)出版社,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] 楊懷棟,徐立,陳科新,等. 盲目反卷積光譜圖超分辨復(fù)原算法[J]. 光譜學(xué)與光譜分析,2007,27(7):1249-1253.

[9] 李玲,胡學(xué)剛,蔣偉. 一種基于LIP的全變分圖像去噪新模型[J]. 四川師范大學(xué)學(xué)報(bào):自然科學(xué)版,2011,34(2):134-138.

[10] 楊麗娟,張白樺,葉旭楨. 快速傅里葉變換FFT及其應(yīng)用[J]. 光電工程,2004,31(S1):1-3.

[11] 張瑩. 一維搜索的程序?qū)崿F(xiàn)[J]. 沈陽(yáng)教育學(xué)院學(xué)報(bào),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在數(shù)學(xué)建模中的應(yīng)用[M]. 北京:北京航空航天大學(xué)出版社,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.