曹玉林, 馬建萍

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

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

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

1 基本原理

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

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

圖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)都是廣義平穩(wěn)的,y(n)也就是廣義平穩(wě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)式的一個近似是

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

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

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

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

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

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

黃金分割法的搜索過程：

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

2) 按坐標點計算公式α1=b-λ(b-a),α2=a+λ(b-a),計算α1和α2,并計算其對應的函數(shù)值.比較函數(shù)值f1和f2,并縮短搜索區(qū)間;

(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 算法實現(xiàn)

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

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

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

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

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

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

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

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

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

%返回值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所采取的做法是將反卷積問題轉(zhuǎn)化為參數(shù)優(yōu)化問題.假定核函數(shù)為高斯函數(shù),則有

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

由此

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

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

達到最小.在這里使用一維搜索方法進行尋優(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的相對誤差即

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

實驗方法表明誤差品質(zhì)

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

3 實驗結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)語

利用XRD圖譜正性和儀器寬化函數(shù)類型已知等先驗知識,將迭代盲目反卷積算法用于X射線衍射圖譜,數(shù)值試驗表明迭代盲目反卷積算法能夠有效提高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]. 核電子學與探測技術(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] 鄒謀炎. 反卷積與信號復原[M]. 北京:國防工業(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] 楊懷棟,徐立,陳科新,等. 盲目反卷積光譜圖超分辨復原算法[J]. 光譜學與光譜分析,2007,27(7):1249-1253.

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

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

[11] 張瑩. 一維搜索的程序?qū)崿F(xiàn)[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在數(shù)學建模中的應用[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.