一种半自动X射线晶面指数标定方法和晶胞常数计算方法

具体实施方式

本发明提供了一种用于LaNi₅基贮氢材料的基于MATLAB的半自动X射线晶面指数标定方法和晶胞常数计算方法，依据LaNi₅基贮氢材料具有与LaNi₅在X射线谱上相同的峰形的基本原理编制了计算机半自动标定程序实现晶面标定，同时利用该晶面标定结果采用最小二乘法编制了晶胞常数精确计算程序，包括如下步骤：

1)将样品制粉并过20～30um细筛，然后在X射线衍射仪上，使用CuKα辐射(波长λ＝1.54056)对样品进行连续扫描，扫描速度：0.12～2(秒/步)，扫描步长：0.02～0.04(度/步)，2θ的范围：20-90(度)得到样品的X射线衍射谱。

2)根据最小二乘法数学原理推导出针对LaNi₅基贮氢材料的用于晶胞常数计算的正则方程组：

$>\begin{array}{}>>>\Sigma >>\alpha >i>>>sin>2>>>\theta >i>>=>Asup>>\Sigma \alpha >i>2sup>>+>C\Sigma >>\alpha >i>>>\gamma >i>>+>D>>\Sigma \alpha >i>>>\delta >i>>>>>>>\Sigma \gamma >i>>>sin>2>>>\theta >i>>=>A\Sigma >>\alpha >i>>>\gamma >i>>+>Csup>>\Sigma \gamma >i>2sup>>+>>D\Sigma \gamma >i>>>\delta >i>>>>>>\Sigma >>\delta >i>>>sin>2>>>\theta >i>>=>A>>\Sigma \alpha >i>>>\delta >i>>+>C>>\Sigma \gamma >i>>>\delta >i>>+>D\Sigma sup>>\delta >i>2sup>>>>>>>\end{array}$

其中 $>>A>=>>>\lambda >2>sup>>>3>a>>0>2sup>>>,>>>\alpha ＝h2+hk+k2，$ >>C>=>>>\lambda >2>sup>>>4>c>>0>2sup>>>,>>>\gamma ＝l2，$ >>\delta >=>>sin>2>>2>\theta >>(>>1>>sin>\theta >>>+>>1>\theta >>)>>,>>>\lambda 为X射线波长，\theta 为衍射角，h、k、l为晶面指数，a$$$₀、c₀为晶胞常数。

3)依据该方程组编制MATLAB函数cohenh，其输入变量为样品的衍射角和对应晶面指数和衍射线的波长，输出变量为样品的晶胞常数a₀、c₀、V₀。

4)由于LaNi₅基贮氢材料具有与LaNi₅相类似的晶胞结构，所以表现在X射线谱上LaNi₅基贮氢材料具有与LaNi₅相同的峰形，只是其峰位发生了改变，依据该原理编制MATLAB函数indexing，其输入变量为样品的X射线衍射谱的衍射角和对应的衍射强度，以及LaNi₅标准X射线衍射峰的衍射角与对应晶面指数，输出变量为样品的衍射角和对应晶面指数。

5)编制MATLAB函数cal_cell_para，完成实验数据的分配，显示样品和LaNi₅合金X射线谱以及LaNi₅标准峰位分布，以及分别调用indexing函数和cohenh函数，完成标定及计算后将数据显示在窗口中。indexing函数用于标定样品的X射线衍射谱，输入变量为样品的X射线衍射谱的衍射角和对应的衍射强度，以及LaNi₅标准X射线衍射峰的衍射角与对应晶面指数，输出变量为样品的衍射角和对应晶面指数；cohenh函数用于计算样品的晶胞常数，输入变量为样品的衍射角和对应晶面指数和衍射线的波长，输出变量为样品的晶胞常数a₀、c₀、V₀，所用数学计算方法为最小二乘法。

下面分别就建立正则方程组的数学推导过程和晶面指数标定的程序实现思想进行详细描述。

最小二乘法的目的是求出被测量数据的最佳值，最佳值的条件是它与测量值之差的平方和为最小，即：∑(最佳值-测量值)2＝最小值

LaNi₅和LaNi₅基贮氢材料的晶胞均属于六方晶系结构，在建立六方晶系的正则方程组之前必须建立与之相应的精确误差函数。

误差函数的表示：

用X射线衍射法测定点阵参数的依据是衍射线的位置，即2θ角。在衍射花样已经指数化的基础上，可通过Bragg方程和晶面间距公式计算点阵参数。根据布拉格方程：

$>>>sin>2>>\theta >=>>>\lambda >2>>>>4>d>>2>>>·>·>·>>(>1>)>>>>$

两边取对数：

$>>>>In>sin>>2>>\theta >=>In>>>\lambda >2>>4>>->2>Ind>·>·>·>>(>2>)>>>>$

两边微分(为避免与晶面间距d混淆，微分符号用Δ)：

$>>>>>>\Delta >sin>>2>>\theta >>>>sin>2>>\theta >>>=>2>>(>>\Delta \lambda >\lambda >>->>\Delta d>d>>)>>·>·>·>>(>3>)>>>>$

从而得到：

$>>>\Delta d>d>>=>>\Delta \lambda >\lambda >>->cot>\theta >·>\Delta \theta >·>·>·>>(>4>)>>>>$

可见d值的误差主要来源于衍射角误差Δθ与波长误差Δλ，而X射线波长的精度已达10^-6且对任何测量此误差均相等，故忽略不计，即Δλ＝0，所以：

$>>>>>>\Delta >sin>>2>>\theta >>>>sin>2>>\theta >>>=>->2>>\Delta d>d>>·>·>·>>(>5>)>>>>$

对于晶面间距的相对误差，Nelson和Riley以及Taylor和Sinelair等研究了不同吸收性质的试样，得出更为严格的误差表达式：

$>>>\Delta d>d>>=>k>·>>cos>2>>\theta >>(>>1>>sin>\theta >>>+>>1>\theta >>)>>·>·>·>>(>6>)>>>>$

将(5)式代入(6)式得(令D＝-(k/2))：

$>>>>\Delta >sin>>2>>\theta >=>->>k>2>>>sin>2>>\theta >>(>>1>sin>>+>>1>\theta >>)>>=>D>>sin>2>>2>\theta >>(>>1>>sin>\theta >>>+>>1>\theta >>)>>·>·>·>>(>7>)>>>>$

又由于

Δsin²θ＝sin²θ(具有一定系统误差的最佳值)-

sin²θ₀(真实值)                           ，根据六方晶系的衍射方程：

$>>>sin>2>>>\theta >0>>=>>>\lambda >2>sup>>>3>a>>0>2sup>>>>(>>h>2>>+>hk>+>>k>2>>)>>+>>>\lambda >2>sup>>>4>c>>0>2sup>>>·>>1>2>>·>·>·>>(>8>)>>>>$

所以由(7)式和(8)式得：

$>>+>D>>sin>2>>2>\theta >>(>>1>>sin>\theta >>>+>>1>\theta >>)>>·>·>·>>(>9>)>>>>$

这就是六方晶系sin²θ误差函数的表达式。

正则方程组的确定：

在sin²θ误差函数的基础上，利用最小二乘法来确定目标方程。

            sin²θ(最佳值)＝Aα+Cγ+Dδ                 (10)

式中 $>>A>=>>>\lambda >2>sup>>>3>a>>0>2sup>>>,>>>\alpha ＝h2+hk+k2，$ >>C>=>>>\lambda >2>sup>>>4>c>>0>2sup>>>,>>>\gamma ＝l2，$ >>\delta >=>>sin>2>>2>\theta >>(>>1>>sin>\theta >>>+>>1>\theta >>)>>.>>>$$$

对于衍射谱中的n个衍射峰，可分别求出它们的sin²θ，α，γ，δ代入(10)式，可以得到n个这样的方程，将每个方程乘以未知数A的系数α并相加得：

       ∑α_isin²θ_i＝A∑α_iu+C∑α_iγ_i+D∑α_iδ_i                 (11)

将每个方程乘以C的系数γ并相加得：

       ∑γ_is_in²θ_i＝A∑α_iγ_i+C∑γ_i²+D∑γ_iδ_i                  (12)

将每个方程乘以D的系数δ并相加得：

      ∑δ_is_in²θ_i＝A∑α_iδ_i+C∑γ_iδ_i+D∑δ_i²                   (13)

相应的正则方程组为：

$>\begin{array}{}>>>\Sigma >>\alpha >i>>>sin>2>>>\theta >i>>=>A\Sigma >>>\alpha >i>>2>>+>C>\Sigma >>\alpha >i>>>\gamma >i>>+>D\Sigma >>\alpha >i>>>\delta >i>>>>>>\Sigma >>\gamma >i>>>sin>2>>>\theta >i>>=>A\Sigma >>\alpha >i>>>\gamma >i>>+>C\Sigma >>>\gamma >i>>2>>>+>D\Sigma >>\gamma >i>>>\delta >i>>>>>>\Sigma >>\delta >i>>>sin>2>>>\theta >i>>=>A\Sigma >>\alpha >i>>>\delta >i>>+>C\Sigma >>\gamma >i>>>\delta >i>>+>D\Sigma >>>\delta >i>>2>>>>>>>>\end{array}$

解出正则方程组得到A和C，进而获得六方晶胞的点阵参数a₀和c₀。

指数标定的程序实现思想：

由于LaNi₅基贮氢材料具有与LaNi₅相类似的晶胞结构，所以表现在X射线谱上LaNi₅基贮氢材料具有与LaNi₅相同的峰形，只是其峰位发生了改变。所以，首先通过作出LaNi₅标准xrd峰值图和具有同样晶型LaNi₅基贮氢材料的xrd谱。依据LaNi₅xrd谱和标准xrd峰值图用户可以很快得确定样品xrd谱的峰位和晶面指数的对应关系。该方法很好得体现了人机协作的概念：由计算机确定样品xrd谱的峰位然后由人判定峰位对应的晶面与LaNi₅标准xrd峰值图中对应晶面的相似关系，从而最终由计算机确定峰位与晶面指数一一对应关系。

如图2a-c所示MATLAB程序流程图，整个程序包含3个函数cal_cell_para函数、cohenh函数、indexing函数；其中indexing函数用于标定样品的X射线衍射谱，输入变量为样品的X射线衍射谱的衍射角和对应的衍射强度，以及LaNi₅标准X射线衍射峰的衍射角与对应晶面指数，输出变量为样品的衍射角和对应晶面指数；cohenh函数用于计算样品的晶胞常数，输入变量为样品的衍射角和对应晶面指数和衍射线的波长，输出变量为样品的晶胞常数a₀、c₀、V₀，以下是整个标定及计算过程的MATLAB程序；cal_cell_para函数用于完成实验数据的分配，显示样品和LaNi₅合金X射线谱以及LaNi₅标准峰位分布，以及分别调用indexing函数和cohenh函数，完成标定及计算后将数据显示在窗口中。

该部分程序用于标定样品X射线衍射谱的晶面指数：

function[out1，out2，out3，out4]＝indexing(thet2a，intensity，standxrd)

％thet2a，intensity--分别为样品X射线衍射谱的2theta，强度

％standxrd为LaNi5标准衍射数据

％LaNi5标准衍射数据来源于“粉末衍射标准联合委员会(JCPDS)”的PDF卡片

％获得所需峰的2theta值及对应晶面指数

％copyright@alex12cheng，Cheng Honghui

％------------------------

％load the standard-xrd of LaNi5

％get the standard thet2a and h，k，l

intensity＝intensity．/max(intensity).*100；％normlizing

thet2a_standard＝standxrd(：，1)；

h_standard＝standxrd(：，3)；

k_standard＝standxrd(：，4)；

l_standard＝standxrd(：，5)；

％plot the peak of the standxrd and the requested xrd

plot(standxrd(：，1)，standxrd(：，2)，′R+′，thet2a，intensity，′black-′)；

set(gcf，′menubar′，′none′，′Position′，[1 39 1024 698])；％make the interface window fully open

set(gca，′position′，[0.03，0.03，0.96，0.96])％make the axis more big

axis([20，90，0，102])；

hold on

thet2a peak＝

p＝；％permutation vector

n＝0；

but＝1；％when push the ESC key，the choosing peak manipulation terminate

while but～＝27％ESC key

    [xi，yi，but]＝ginput(1)；

          if but～＝27

              ％Tb improve the accurancy,the peak will be acquired near the

              ％locafion of cursor

              upi＝xi+0.2；

              downi＝xi-0.2；

Irange＝find(thet2a＞downi&thet2a＜upi)；

[int_peaki，Ipeak]＝max(intensity(Irange))；

thet2a_peaki＝thet2a(Irange(Ipeak))；

plot(thet2a_peaki，int_peaki，′m^′)；

n＝n+1；

thet2a_peak(n，1)＝thet2a_peaki；

ifbut＝＝1％click the left button of mouse if the standard peak is close

         ％to the left side ofthe sample xrd peak

    Ip＝find(thet2_astandard＞(thet2a_peaki-5)&thet2a_standard＜thet2a_peaki)；

    p(n)＝Ip(end)；

    hi＝h_standard(Ip(end))；

    ki＝k_standard(Ip(end))；

    li＝l_standard(Ip(end))；

elseifbut＝＝97％click a key on keyboard if the standard peak is not

                ％so close to the left side ofthe sample xrd peak

    Ip＝find(thet2a standard＞(thet2a_peaki-5)&thet2a_standard＜thet2a_peaki)；

    p(n)＝Ip(end-1)；

    hi＝h_standard(Ip(end-1))；

    ki＝k_standard(Ip(end-1))；

    li＝l_standard(Ip(end-1))；

elseif but＝＝32％click spacebar if the manipulator can not figure out the location

                ％relation between the standardxrd peak and the sample xrd peak

    Ip＝find(thet2a_standard＞(thet2a_peaki-0.5)&thet2a_standard＜(thet2a_peaki+0.5))；

    [C，I]＝min(abs(thet2a_standard(Ip)-thet2a_peaki))；

    p(n)＝Ip(I)；

    hi＝h_standard(Ip(I))；

    ki＝k_standard(Ip(I))；

    li＝l_standard(Ip(I))；

elseif but＝＝3％click right button of mouse if the standard peak is close to the

                     ％right side of the sample xrd peak

            Ip＝find(thet2a_standard＞thet2a_peaki&thet2a_standard＜(thet2a_peaki+5))；

            p(n)＝IP(1)；

            hi＝h_standard(Ip1))；

            ki＝k_standard(Ip(1))；

            li＝l_standard(Ip(1))；

        elseifb but＝＝100％click d key on keyboard if the standard peak is not so

                          ％close to the right side of the sample xrd peak

            Ip＝find(thet2a_standard＞thet2a_peaki&thet2a_standard＜(thet2a_peaki+5))；

            p(n)＝Ip(2)；

            hi＝h_standard(Ip(2))；

            ki＝k_standard(Ip(2))；

            li＝l_standard(Ip(2))；

        end

        string＝[num2str(thet2a_peaki)，′(′，num2st(hi)，′′，num2str(ki)，′′，num2str(li)，′)′]；

        text(thet2a_peaki，int_peaki，string，...

        ′VerticalAlignment′，′bottom′，′HorizontalAlignment′，′center′)；

        ％display the result of indexing

    end

end

holdo off

h＝h_standard(p)；

k＝k_standard(p)；

l＝l_standard(p)；

hkl＝[h，k，l]；

thet2a_peakhkl＝[thet2a_peak，h，k，l]；

if nargout＝＝1

      out1＝thet2a_peakhkl；

elseif nargout＝＝2

    out1＝thet2a_peak；

    out2＝hkl；

elseif nargout＝＝4

    out1＝thet2a_peak；

    out2＝h；

    out3＝k；

    out4＝1；

end

该部分程序用于样品的晶胞常数计算：

function[a0，c0，V0]＝cohenh(h，k，l，thet2a，lumda)

％[a0，c0，V0]＝cohenh(h，k，l，thet2a，lumda)

％a0，c0，V0为所要获得的点阵参数

％cohenh为函数名

％h，k，l为晶面指数

％thet2a为对应的衍射角，单位为度

％lumda为X射线波长，单位为埃

theta＝thet2a./2；

alpha＝h.^2+h.*k+k.^2；

gama＝1.^2；

deta＝sin(thet2a.*pi./180).^2.*(sin(theta.*pi./180)+theta.*pi./180)./(theta.*pi./180)./sin(theta.*pi./180)；

a(1，1)＝sum(alpha.^2)；

a(1，2)＝sum(alpha.*gama)；

a(1，3)＝sum(alpha.*deta)；

a(2，2)＝sum(gama.^2)；

a(2，3)＝sum(gama.*deta)；

a(3，3)＝sum(deta.^2)；

for i＝2:3

    for j＝1:2

        ifi～＝j

            a(i，j)＝a(j，i)；

        end

    end

end

b(1)＝sum(alpha.*sin(theta.*pi./180).^2)；

b(2)＝sum(gama.*sin(theta.*pi./180).^2)；

b(3)＝sum(deta.*sin(theta.*pi./180).^2)；

x＝a\b′；

A＝x(1)；

C＝x(2)；

a0＝sqrt(lumda^2./3./A)；

c0＝sqrt(lumda^2./4./C)；

V0＝a0^2.*c0.*sin(pi/3)；

该部分程序完成实验数据的分配，显示样品和LaNi₅合金X射线谱以及LaNi₅标准峰位分布，以及分别调用indexing函数和cohenh函数，完成标定及计算后将数据显示在窗口中：

function cal_cell_para(Sample_xrd，lumda，standxrd_JCPDS，LaNi5xrd_reference)

plot(LaNi5xrd_reference(：，1)，LaNi5xrd_reference(：，2)，′b-′)；

text(60，40，′蓝色为LaNi5-x射线谱′，′color′，′blue′，′fontsize′，15)；

hold  on

thet2a＝Sample_xrd(：，1)；

intensity＝Sample_xrd(：，2)；

[thet2a，h，k，l]＝indexing(thet2a，intensity，standxrd_JCPDS)；

[a0，c0，V0]＝cohenh(h，k，l，thet2a，lumda)

s＝sprintf(′a0＝％0.4f＼nc0＝％0.4f＼n V0＝％0.4f，a0，c0，V0)；

text(60，80，s，′color′，′red′，′fontsize′，15)；

以纯度(at％)为La 99.3，Ni 99.9，Al 99.7的原料按设计成分34∶61.05∶4.95配重。将配好的料置于水冷铜模中，在氩气气氛保护下的电弧炉中进行熔炼。熔炼过程中将合金翻面重熔五次，同时进行电磁搅拌，以保证合金的均匀性。熔炼后将合金锭封于真空石英管中，真空度＜10^-2a，放入热处理炉随炉升温至1100℃，保温8个小时后淬火制得样品。研磨样品过20～30um细筛，使用日本理学(Rigaku)D/max-2500pc X射线衍射仪和CuKα辐射(λ＝1.54056)对粉末进行连续扫描，扫描速度：0.12(秒/步)，扫描步长：0.04(度/步)，2θ的范围：20-90(度)。整个实验共得到1750个数据点。将晶体参数及指数标定文件夹内文件名为LaNi5-CuKa-standard-xrd的数据文件导入到MATLAB的工作空间(workspace)中，然后将样品X射线衍射谱的衍射数据替代Sample_xrd已有数据即可(如图1所示)。然后在MATLAB的命令窗口(command window)中键入cal_ cell_para(Sample_xrd，lumda，standxrd_JCPDS，LaNi5xrd_reference)命令后屏幕上会出现一个人机交互式窗口。进行晶面指数标定及晶胞常数计算，得到LaNi_4.25Al0.75合金的晶胞常数a₀＝5.0525；c₀＝4.0535；V₀＝89.6124³(V₀为晶胞体积)