基于多对比度磁共振颈动脉图像的配准算法研究

吴玉霞 徐肖攀 张 曦 刘 洋 张国鹏 陈慧军 卢虹冰 李宝娟*

1(第四军医大学生物医学工程系，西安 710032)2(清华大学生物医学影像研究中心，北京 100084)

基于多对比度磁共振颈动脉图像的配准算法研究

吴玉霞1徐肖攀1张 曦1刘 洋1张国鹏1陈慧军2卢虹冰1李宝娟1*

1(第四军医大学生物医学工程系，西安 710032)2(清华大学生物医学影像研究中心，北京 100084)

基于在血管壁增强显示方面的强大潜力，近年来多对比度核磁共振成像已成为斑块分析研究的有力工具，但其效能受到多序列图像血管不匹配的影响。为实现多序列图像的准确血管配准，在管腔分割的基础上提出一种由粗到精的两步配准策略：先采用迭代最近点实现多对比度图像中心线的刚性配准，再采用薄板样条实现基于血管边界点的非刚性配准。在第二步配准中，为准确寻找不同序列血管边界的匹配点，创新使用形状上下文描述子对边界点进行筛选，并应用确定性退火技术进行全局优化。采用新型三维多对比度磁共振血管成像序列，对提出算法的有效性进行定量评价。结果表明，配准后不同序列的边界重合度均达到95%以上，平均表面距离0.12 mm，可有效提高配准精度，为后续斑块成分的分析奠定基础。

颈动脉粥样硬化；易损斑块；磁共振图像；配准

引言

大量临床研究表明，颈动脉粥样硬化斑块成分及易损性与脑卒中等心脑血管疾病的发生密切相关[1]，对颈动脉斑块及其成分的诊断分析成为目前影像学技术研究的热点。其中，核磁共振成像具有无创、无辐射、软组织分辨率高等优势，特别是，多对比度MRI成像既能够显示管腔狭窄，又能显示斑块形态及成分信息，与组织病理学具有高度的一致性[2]，从而为颈动脉斑块的诊断和评估提供了一种强有力的检测手段。但是，多对比度成像存在扫描时间长、空间分辨率低等局限[3]，加上血管自身拓扑复杂，存在直径、曲率、分支等改变，造成多序列间血管的不完全匹配，这对于像素级别的颈动脉粥样硬化斑块成分识别是致命的。因此，多序列图像准确配准是此类研究的关键基础。

根据采用相似性测度的不同，图像配准可分为基于灰度和基于特征的两大类算法。其中，基于灰度信息(如互信息[4]、归一化互信息[5]等)算法的核心是像素的联合概率分布模型。对于多序列图像的血管配准，由于血管及斑块在不同序列中的灰度分布不同，且体积很小，当空间信息被忽略时，配准的准确性有待提升。基于特征的配准主要利用图像的几何结构特征(如点、边缘、轮廓等)的对应关系作为相似性测度，在血管等复杂结构的配准中应用广泛。Can等采用血管结构中的分支和交叉点作为特征点进行配准[6]。Stewart等则加入血管中心线特征，包括中心线的位置、切线方向和宽度，改进配准的精度[7]。Guo等与笔者类似，基于管腔中心线和管腔边界点，采用改进的迭代最近点(iterative closest point，ICP)和薄板样条(thin plate spline, TPS)函数进行配准[8]。由于采用均匀采样法对不同序列的匹配点进行筛选，对于具有复杂形变的血管，其配准精度有待提高。

为实现高精度的血管配准，本研究在管腔分割的基础上，提出了基于特征点的两步自动配准方法。由于中心线能够很好地表达血管的拓扑结构，因此首先选取血管中心线为特征点，采用迭代最近点算法进行三维刚性粗配准；在此基础上，以管腔内边界点作为特征点，创新使用形状上下文(shape context)方法进行匹配点筛选，并采用薄板样条函数进行二维非刚性精配准，从而实现不同序列图像中颈动脉血管的更准确匹配。

1 方法

1.1 影像数据的采集

采用的多对比度MRI数据来自清华大学生物医学影像研究中心，包含3个新型3D序列，即MERGE (motion sensitized driven equilibrium prepared rapid gradient echo)[9]、SNAP (simultaneousnoncontrast angiography and intraplaque hemorrhage)[10]和T2W-VISTA (T2 weighted volume isotropic TSE acquisition)序列[11]。由于3个序列的管壁管腔信号互补，其中MERGE的管壁信号较强，有利于斑块负荷的测量，SNAP管腔信号较强，又可以清晰地显示出血病灶，T2W-VISTA和SNAP可用于检测高危斑块成分，因此结合分析，可为临床提供更客观、准确的判断信息[12]。

实验中采用Philips 3T Achieva TX扫描仪和清华大学自制神经血管36通道线圈，对具有不同程度粥样硬化斑块的7名患者进行扫描，各序列扫描参数如表1所示。

表1 多对比度MRI序列扫描参数

Tab.1 Scanning parameters of multi-contrast MRI sequences

MERGESNAPT2W-VISTA(TR/TE)/ms9.2/4.310/4.82500/268层厚/mm0.80.80.8层间距/mm0.40.40.4像素间距/mm0.34720.39060.4883采集方向冠状位冠状位冠状位采集矩阵/像素720×720×100640×640×100512×512×124

1.2 算法流程

本研究提出的配准算法流程如图1所示。首先采用交互式图割方法，从待配准颈动脉MRI图像中分割出血管内壁，进而提取血管三维中心线；然后采用迭代最近点算法，实现基于中心线的粗配准；再用形状上下文描述子，筛选出相匹配的边界点；最后对获得的匹配点进行薄板样条变换，从而对图像进行精配准，进一步提高配准精度。

图1 算法流程Fig.1 The diagram of the proposed method

1.3 颈动脉血管分割及中心线的提取

在本研究中，血管壁的准确分割是多对比度血管图像配准的基础。为实现血管内壁的准确提取，笔者采用交互式图割方法[13]。该方法对于处理血管壁等模糊边界具有较强的优势，能在较复杂的颈动脉图像中通过少量交互工作得到满意的分割效果。

该方法将一幅图像定义为G=<γ，ε>，其中γ是图像中的所有节点，ε是连接节点的每一条边。对图像中的任意点i，定义：xi=1，i属于图像的前景；xi=0，i属于图像的背景。在用户手动勾勒图像前景(血管区域)与背景区域的基础上，通过定义能量函数对图像中的像素属于前景还是背景进行判断，有

(1)

式中，E1是衡量将一个节点i定义为前景或者背景所需要消耗的能量，E2是当两个相邻的像素标记不一样时付出的代价，λ是平衡这两个能量函数的一个参数。有关算法的具体描述及各项的定义计算详见文献[13]。

获得颈动脉内壁边界点后，血管中心线则由每一层边界的几何中心点确定。

1.4 颈动脉血管配准

考虑到血管拓扑结构及变化的复杂性，本研究基于提取的管腔内边界和血管中心线，提出由粗到细的两步配准法，实现多对比度血管图像的准确配准。首先采用迭代最近点算法，对不同序列的血管中心线进行三维粗配准，获得较为准确的初始位置。由于中心线难以反映血管壁的各种变化，因此将血管边界点作为特征点，进一步进行精配准。考虑到不同序列图像的灰度差异及形变，其血管边界点的数目和位置有所不同，为提高匹配的准确度，提出基于形状上下文的匹配点筛选方法。通过比较两个形状的相似性，对边界点进行筛选，以充分利用边界点的相对位置信息。在此基础上，对获得的匹配点进行薄板样条变换，并用确定性退火算法进行全局优化，这样不仅增强了算法的鲁棒性，而且大大减小了陷入局部极小值的可能性，达到更好的配准效果。

1.4.1 基于ICP算法的粗配准

为获得不同序列颈动脉中心线间的刚性变换矩阵，首先采用ICP算法进行三维几何粗配准。

设P和Q分别为参考图像和浮动图像颈动脉中心线三维点集，点数分别为Np和Nq，对Q中每个点qk，在P中寻找与该点欧式距离最近的点pk，设这两组点集间的旋转矩阵和平移向量分别为R和T，则求解以下最优化问题就能获得R和T，有

(2)

本研究基于改进的ICP方法，利用k-d树寻找邻近点加速迭代过程[14]，再用四元数奇异值分解求得R和T，并用立方插值法获得配准后图像。

1.4.2 基于TPS算法的精配准

经过粗配准后，参考图像和浮动图像中血管的位置已比较接近，分叉处在同一层，因此精配准阶段可以在二维空间中逐层进行。在对两边界点集之间变换关系求解时，需要解决两个问题：匹配点的筛选和匹配点间变换关系的求解。

1.4.2.1 匹配点的筛选

人们已研究和实现了多种不同的形状匹配算法，其中最常用的是均匀采样法，通过对两组边界点进行插值，从而分别均匀地选取相同个数的样本点，组成一一对应的匹配点。由于血管存在非刚性形变，为提高匹配点的准确度，本研究提出采用基于形状上下文的算法对边界点进行优化筛选[15]。

作为近年来提出的最成功的形状描述子之一，形状上下文的定义如图2所示。设血管边界点坐标集合P={p1,p2,…，pN}，对其进行极坐标变换；对于每一点pi，与其余的(N-1)个点可构成(N-1)个方向向量，这(N-1)个向量提供了丰富的信息，决定了目标的形状。

图2 内壁上点pi到其余点的向量Fig.2 The vectors from pi to other points on the inner wall

构造如图3所示的对数极坐标系：它将整个平面空间平均划分出m个方向，同时在半径上从(N-1)个向量长度的最小值到最大值按对数等分为n份，整个平面空间就被划分为K=m×n个区域。对于边界上的任意一点pi，可以用以pi点为中心的极坐标系中处于每一个区域bin(k)的像素点个数hi(k) 对它的属性进行描述，即得到该点的形状上下文，有

蚀花肉红石髓珠5粒。蚀花肉红石髓珠又称“光玉髓”，常常被古代人制作成小颗的串珠，作为随身佩戴的饰物。此类蚀花工艺最早出现于西亚和南亚一带。

hi(k)={q≠pi∶(q-pi)∈bin(k)}

(3)

k={1,2,…，K}

图3 对数极坐标系Fig.3 The log-polar coordinates

对于参考图像上任一点pi与浮动图像上任一点qj，定义这两个点的匹配代价为

(4)

式中，hi(k)和hj(k)分别代表pi和qj的归一化直方图。

如果两边界点的数目不相等，则为较小点集增加虚拟节点以达到数目相等，每个虚拟节点与另一边界上所有点之间的匹配代价值为εd。匹配代价值越小，表示两个点的形状上下文越接近，匹配的相似度越高，反之，匹配相似度越低。要达到完全匹配，必须使两个点集总的匹配代价最小。这是一个典型的平方分配问题(或称“加权二分图匹配”)，使用匈牙利算法可有效解决[16]。

1.4.2.2 匹配点间变换关系的求解

获取匹配点后，采用TPS算法[17]，求出其映射关系，从而实现血管图像的非刚性配准。为了避免陷入局部极小值，采用确定性退火技术作为全局优化算法[18]，根据退火过程，将求解优化问题的最优点转化为求一系列随温度变化的物理系统的自由能函数的极小值，最终避开局部极小值而得到全局极小值。

1.5 实验设置

为评测算法的有效性，在实验中将本算法与基于均匀采样法的匹配点筛选方法进行比较。为保证公平对比，两算法中相同意义的参数设置相同。其中，图像分割过程中，根据文献[13]及分割效果，本研究将参数设定为：迭代次数为6，前景点的聚类数为5，背景点的聚类数为3，平衡参数λ=50。在粗配准中，根据收敛情况设置迭代次数为35。在精配准阶段，根据文献报道[19-20]及初步实验结果，本研究将参数设定为：对数极坐标系的方向数m=12，半径划分n=5份，虚拟节点匹配代价值εd=0.15，退火速率r=0.95，迭代次数t=5。

1.6 评价准则

为了客观评价提出算法的配准效果，本研究采用边界重合度(overlap ratio，OR)、平均表面距离(mean surface distance,MSD)和最大表面距离(maximum surface distance,MAXD)作为定量评价指标，其计算方法如下：

(5)

式中，UF和UM分别代表参考图像和浮动图像的管腔边界内区域，OR范围是0 ～1，其值越大代表重合度越高，配准效果越好。

式中，LM表示浮动图像管腔内边界点集，pi表示参考图像管腔内边界上的某一点，n表示参考图像管腔内边界点的个数，q表示浮动图像管腔内边界上的点。

MSD、MAXD的值越小，表示配准效果越好。

2 结果

某患者的MERGE和SNAP序列中心线粗配准结果如图4所示。经过粗配准，血管中心线的整体位置以及两个分叉点位置均比配准前重合度高，平均误差为0.517 3 mm，为精配准阶段提供良好的初始位置。

经过粗配准和精配准，3个序列的最终配准结果如图5所示。其中，上图是血管分叉部分同一层的灰度图像，下图是血管主干部分同一层的灰度图像，图中的边界线是SNAP序列图像的血管内壁分割结果。可以看出，配准后另外两个序列中血管的位置和形状与SNAP序列基本一致。

图4 中心线粗配准结果。(a)配准前两序列中心线的位置；(b)配准后两序列中心线的位置Fig.4 Coarse registration of centerlines. (a) The position of two centerlines before registration; (b) The position of two centerlines after registration

图5 MERGE 和 T2W-VISTA序列图像分叉部分(上行)和主干部分(下行)基于边界的配准结果。(a) SNAP序列; (b) MERGE序列配准后图像; (c) T2W-VISTA序列配准后图像Fig.5 Boundary-based registration results for the bifurcations (top row)and trunks (bottom row) of MERGE and T2W-VISTA images.(a) Images of SNAP; (b) Images of MERGE after registration; (C) Images of T2W-VISTA after registration.

为了验证匹配点筛选中形状上下文的优势，本研究对同样的数据采用文献[8]的均匀采样法进行匹配点筛选，精配准前后7例数据的平均重合度如图6所示，其中每组数据依次表示粗配准后、形状上下文精配准后、均匀采样精配准后的平均重合度。可以看出，只采用粗配准时重合度较低，进一步精配准后重合度得到极大提高。在本方法中，MERGE序列和SNAP序列精配准后平均重合度从66.90%提高到了95.51%，T2W-VISTA序列和SNAP序列精配准后平均重合度从59.25%提高到了95.08%。而采用均匀采样方法后，两序列与SNAP序列精配准后平均重合度分别为93.54%和93.94%。两种方法MSD和MAXD的结果如表2所示。

图6 不同方法精配准前后平均重合度测量Fig.6 The mean value of overlap ratios before and after the fine registration with different methods

表2 配准后的表面距离测量

Tab.2 MSD andMAXDafter the registration using the proposed method

配准序列形状上下文方法均匀采样方法MSD/mmMAXD/mmMSD/mmMAXD/mmSNAP-MERGE0.11±0.001.35±0.510.15±0.001.51±0.52SNAP-T2W-VIS-TA0.12±0.001.32±0.050.14±0.001.62±0.53

从图6和表2可以看出，采用形状上下文描述子进行匹配点筛选，序列边界的重合度均在95%以上，平均表面距离≤0.12 mm，配准误差更低，充分表明所提出算法的有效性。

3 讨论

从患者多序列图像配准结果可以看出，本研究提出的两步配准方法可有效提高配准精度。如图6所示，只采用粗配准时，由于只利用了血管中心线的拓扑结构信息及刚性配准算法，难以反映患者呼吸、脉动、体位改变等引起的非刚性形变和血管形状改变，血管重合度较低，因此需要进一步基于血管边界点进行非刚性精配准。

文献[21]采用归一化互信息作为相似性测度、三维B样条模型作为变换函数，对多对比度MRI颈动脉图像进行配准，得到的相似系数(类似于边界重合度)小于90%，平均表面距离为0.288±0.128 mm。与形状相似性不同，互信息方法统计的是两个图像间联合直方图的分布，反映的是信息量的变化，配准效果会受图像噪声的影响。由于空间信息被忽略，难以反映像素位置的改变。从图4可以看出，本研究采用的3个对比度图像，相比传统MRI成像，血管边界的清晰度更高，更有利于血管边界的准确分割，进一步改进了基于形状特征的血管配准效果。考虑到不同序列图像灰度分布的差异，今后研究中可考虑将灰度信息和空间信息的配准相结合，进一步提高弱边界处的配准精度。

基于特征点的配准算法的一个局限是往往需要选取角点、高曲率点等作为对应特征点，以更好地反映形状的改变。例如，文献[6]就选取分叉点作为特征点，有些变化大的区域甚至需要手工选取对应点。本研究的创新性在于基于提取的血管边界，自动使用形状上下文对边界点进行特征描述，无需选取特殊标记点或关键点就能得到边界点的准确对应关系，具有更广泛的适用性。

为了验证形状上下文在边界匹配中的优势，笔者将其与文献[8]中的均匀采样方法进行了比较。由于血管存在非刚性形变，通过均匀采样得到的匹配点实际上难以一一对应。形状上下文能够对边界点的形状进行全局描述，充分利用边界点的相对位置信息，并且具有尺度不变性和旋转不变性，有效提高了配准的准确度。

耗时是非刚性配准算法的另一个瓶颈，为了避免陷入局部极小值，本研究采用确定性退火算法来求解最优的全局变换，一般只需迭代5次就能达到较好的配准效果。在配置为3.60 GHz CPU、8 GB RAM的工作站上，使用Matlab 2012b，以一套大小为(84×83×144)像素的三维数据为例，粗配准耗时1.2 s，精配准耗时1 min 50 s。可在今后工作中结合并行运算和算法优化，进一步提高配准速度。

本研究尚存在一些不足之处，希望在今后的研究中加以改进：

1)血管图像分割及配准中所用参数(如前背景点聚类数目、迭代次数等)对结果有一定影响。在实验中，根据文献中的参数值以及初步试验结果，选取了效果较好的参数，下一步研究希望能采用参数优化算法自动选取参数。

2)目前精配准过程是在二维空间中完成，对粗配准的层定位要求高。在后续工作中，拟将形状描述子扩展到三维，以获得边界点的三维空间变换。

3)本研究用于试验验证的样本例数较少，但试验时考虑到对患者的纳入具有一定的随机性，且每个患者均有3个不同的序列，相互间的配准效果已初步证实了算法的有效性。在后续研究中，会采集更多的图像数据，在验证算法有效性的同时，对斑块进行定量分析。

4 结论

本研究基于多对比度颈动脉MRI成像，在血管分割的基础上，提出了基于颈动脉血管中心线和管腔边界的两步配准算法，并在精配准过程中创新使用形状上下文描述子对匹配点进行筛选，提高了颈动脉配准的准确度。基于患者的三维多对比度MRI数据进行实验，结果表明：与常用的均匀采样法相比，采用形状上下文描述子进行匹配点筛选更有优势，提出的方法配准重合度更高。

[1] Naghavi M, Libby P, Falk E, et al. From vulnerable plaque to vulnerable patient: a call for new definitions and risk assessment strategies: part I [J]. Circulation, 2003, 108(14):1664-1672.

[2] CaiJianming, Hatsukami TS, Ferguson MS, et al. Classification of human carotid atherosclerotic lesions with in vivo multicontrast magnetic resonance imaging [J]. Circulation, 2002, 106(11):1368-1373.

[3] 李睿，陈慧军，苑纯. 磁共振颈动脉斑块成像技术的现状与进展[J]. 中国卒中杂志, 2014, 9(2):129-134.

[4] Viola P, Wells III WM. Alignment by maximization of mutual information [J]. International Journal of Computer Vision, 1995, 24(2):16-23.

[5] Studholme C, Hill D, Hawkes DJ. An overlap invariant entropy measure of 3D medical image alignment [J]. Pattern Recognition. 1999, 32(1):71-86.

[6] Can A, Stewart CV, Roysam B, et al. A feature-based, robust, hierarchical algorithm for registering pairs of images of the curved human retina [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2002, 24(3):347-364.

[7] Stewart CV, Tsai CL, Roysam B. The dual-bootstrap iterative closest point algorithm with application to retinal image registration [J]. IEEE Transactions on Medical Imaging, 2003, 22(11):1379-1394.

[8] GuoHengkai, Wang Guilin, Huang Lingyun, et al. A robust and accurate two-step auto-labeling conditional iterative closest points (TACIP) algorithm for three-dimensional multi-modal carotid image registration [J]. PLoS ONE, 2016, 11(2):1-22.

[9] Balu N, Yarnykh VL, Chu B, et al. Carotid plaque assessment using fast 3D isotropic resolution black-blood MRI [J]. Magnetic Resonance in Medicine, 2011, 65(3):627-637.

[10] Wang Jinnan, Börnert P, Zhao Huilin, et al. Simultaneous noncontrast angiography and intraplaque hemorrhage (SNAP) imaging for carotid atherosclerotic disease evaluation [J]. Magnetic Resonance in Medicine, 2013, 69(2):337-345.

[11] Mihai G, Chung Y, Merchant A, et al. T1-weighted-SPACE dark blood whole body magnetic resonance angiography (DB-WBMRA): Initial experience [J]. Journal of Magnetic Resonance Imaging, 2010, 31(2):502-509.

[12] Zhao Huilin, Wang Jinnan, Liu Xiaosheng, et al. Assessment of carotid artery atherosclerotic disease by using three-dimensional fast black-blood MR imaging: comparison with DSA [J]. Radiology, 2015, 274(2): 508-516.

[13] Boykov Y, Veksler O, Zabih R. Fast approximate energy minimization via graph cuts [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2001, 23(11):1222-1239.

[14] 刘江,张旭,朱继文. 一种基于K-D树优化的ICP三维点云配准方法[J].测绘工程, 2016, 25(6):15-18.

[15] Belongie S, Malik J, Puzicha J. Shape matching and object recognition using shape contexts [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2002, 4(24):509-522.

[16] Papadimitriou CH. Combinatorial optimization: algorithms and complexity [M]//New York: Docer Pub, 1982: 248-254.

[17] Bookstein FL. Principal warps: thin-plate splines and the decomposition of deformations [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1989, 11(6):567-585.

[18] Rose K, Gurewitz E, Fox CG. Statistical mechanics and phase transitions in clustering [J]. Physical Review Letters, 1990, 65 (8):945-948.

[19] Belongie S, Malik J, and Puzicha J. Shape context: A new descriptor for shape matching and object recognition [C] // Advaces in Neural Information Processing Systems13. Barcelona: MIT Press, 2001: 831-837.

[20] Chui H, Rangarajan A. A new point matching algorithm for non-rigid registration [J]. Computer Vision and Image Understanding, 2003, 89(2-3):114-141.

[21] van’tKlooster R, Staring M, Klein S, et al. Automated registration of multispectral MR vessel wall images of the carotid artery [J]. Medical Physics, 2013, 40(12): 1-12.

Multi-Contrast Magnetic Resonance Imaging Registration of Carotid Arteries

Wu Yuxia1Xu Xiaopan1Zhang Xi1Liu Yang1Zhang Guopeng1Chen Huijun2Lu Hongbing1Li Baojuan1*

1(DepartmentofBiomedicalEngineering,FourthMilitaryMedicalUniversity,Xi′an710032,China)2(CenterforBiomedicalImagingResearch,TsinghuaUniversity,Beijing100084,China)

Recently multi-contrast magnetic resonance imaging (MRI) has become a powerful tool for plaque analyzing due to its strong potential in enhanced demonstration of carotid wall, but its performance is hampered by the misalignment of difference imaging sequences. To achieve accurate vessel registration of multi-sequence images, a two-step coarse-to-fine registration strategy was proposed in this study. First, the iterative closest point was applied to realize the rigid registration of the centerlines of the multi-contrast images, and then the thin plate spline was used to realize the non-rigid registration based on the boundaries of carotid artery. In the second step, to find the corresponding points of the boundaries in different sequences, a shape context descriptor was innovatively introduced to screen the boundary points. In addition, the deterministic annealing technique was employed to find a globally optimized solution. The effectiveness of the proposed algorithm was quantitatively evaluated by novel three-dimensional (3D) multi-contrast vessel wall MRI sequences. The results indicated that after registration, the overlap of two boundaries from different sequences was more than 95%, and their mean surface distance was 0.12 mm, which improved the accuracy of registration effectively and laid the foundation for further component analysis.

carotid atherosclerosis; vulnerable plaque; magnetic resonance imaging; registration

10.3969/j.issn.0258-8021. 2017. 02.001

2016-08-22， 录用日期:2016-11-22

国家自然科学基金青年基金(81301199)；国家自然科学基金重点项目(81230035);陕西省社会发展科技攻关项目(2016SF302)

R318

A

0258-8021(2017) 02-0129-07

*通信作者(Corresponding author)，E-mail: libjuan@163.com