基于眼底结构特征的彩色眼底图像视盘定位方法

肖志涛 邵一婷 张 芳* 温 佳 耿 磊吴 骏 尚丹丹 苏 龙 单春艳

1(天津工业大学电子与信息工程学院，天津 300387)2(天津医科大学第二医院，天津 300211)3(天津医科大学代谢病医院，天津 300070)4(天津市光电检测技术与系统重点实验室，天津 300387)

基于眼底结构特征的彩色眼底图像视盘定位方法

肖志涛1,4邵一婷1张 芳1*温 佳1耿 磊1,4吴 骏1尚丹丹1苏 龙2单春艳3

1(天津工业大学电子与信息工程学院，天津 300387)2(天津医科大学第二医院，天津 300211)3(天津医科大学代谢病医院，天津 300070)4(天津市光电检测技术与系统重点实验室，天津 300387)

眼底图像中视盘的大小和形状等参数是判断眼底病变的重要辅助参数，视盘的检测和定位对眼科疾病的计算机辅助诊断具有重要意义。提出一种基于眼底结构特征的彩色眼底图像视盘定位方法。首先采用基于低帽运算的方法，提取眼底图像中的静脉血管；然后基于静脉血管的结构特征，采用最小二乘抛物线拟合法初步定位视盘；最后通过滑动窗口灰度扫描的方法，精确定位视盘。在4个公开的眼底图像数据库(DRIVE、DIABETED0、STARE和MESSIDOR)中，对所提出的视盘定位方法进行测试，定位准确率分别为100%、98.6%、93.8%、99.75%，验证了该方法的准确性和通用性。

视盘定位；低帽变换；抛物线拟合；窗口扫描

引言

视盘是眼底图像的一个重要特征，它在眼底图像中是一个类圆形的黄色亮区域。在彩色眼底图像中，由于视盘是血管进入眼部的起始区域，因此主要血管在该区域汇合[1]，即任意两条主血管在视盘内部聚合。在眼底图像中，视盘的大小和形状等参数是判断糖尿病视网膜病变的重要辅助参数，所以视盘的检测和定位对眼底疾病的临床诊断具有重要意义。目前，视盘定位的方法主要有两类。

第1类方法基于视盘的自身特点，即视盘的亮度高，形状呈类圆形，视盘内部灰度差异大。Sinthanayothin等寻找灰度变化幅度最大的矩形区域的中心作为视盘位置[2]。Li等首先找到1%的亮度最大的像素点，再通过聚类分割出视盘候选区域，然后采用PCA分析找到真正的视盘[3]。邹北骥等在预处理中增加裁剪图像来去除或减少眼底图像边缘高亮环对视盘定位的影响，然后通过寻找亮度最大区域来定位视盘[4]。Haar等则依据视盘的形状特点，构建视盘区域结构模板，通过模板匹配方法找到视盘[5]。Ravishankar等利用视盘的类圆形特点，通过Hough变换圆形检测方法来检测视盘[6]。Tobin等则分析血管亮度、宽度和方向信息，定位视神经乳头[7]。郑绍华等采用基于定向对比度的方法，定位视盘[8]。这些方法都充分考虑了视盘的自身亮度和纹理特点。但是，在病变图像中，由于亮的病变区域干扰以及视盘本身的病变造成视盘不完整或亮度变化，使得这类方法很难正确检测视盘。

第2类方法基于血管和视盘的关系特性检测视盘。Foracchia等采用两条抛物线来描述血管在视盘左、右的走向，视盘则刚好位于两抛物线的公共顶点[9]。Hoover等计算血管的汇合度，将血管汇合度最大的点确定为视盘[10]。朱琳琳等通过视盘内部的血管分布和方向特点，构建匹配模板，将匹配值最大的点作为视盘定位点[11]。Youssif等利用眼底图像中的血管方向信息，确定视盘位置[12]。张东波等注意到视盘区域中的血管基本沿垂直方向延伸，由此根据血管的分布特性确定视盘的水平坐标，利用视盘外观特性确定垂直坐标，通过两次投影找到视盘中心[1]。这类算法充分利用血管的特性(如血管的走向、血管的结构等信息)，可以有效地解决由于病变导致视盘在自身特性不明显的情况下的检测问题，但是这些方法需要制定较为严格的几何模板，算法的复杂度高。而且模板制定依赖于血管特性的精确分割，但在病变图像中，血管特性的完整性容易被破坏。此外，Godse等结合视盘的亮度、形状、面积和密集度多种特征，定位视盘[13]。Harangi则提出结合现有多种方法进行视盘定位[14]。这类方法在一定程度上能提高算法的准确性。但在病变图像中，由于病变影响，所制定的标准很难具有普遍适应性，导致检测的鲁棒性较差。

本研究结合较为稳定的基于血管和视盘关系的视盘定位方法并对其进行改进，提出基于静脉血管结构特征的视盘定位方法。此方法不需要制定复杂的几何模板，而且只需要静脉血管的主要轮廓特性，避免了必须精确提取整体血管或者主血管特征的复杂性和局限性。通过分析眼底的结构特征，利用视盘在其邻域所表现的自身特性，采用窗口扫描方法修正初定位结果，保证视盘定位的准确性。

1 方法

1.1 方法概述

图1所示为本方法的总体流程，分为视盘初定位和视盘精确定位两部分。

图1 所提出方法的流程Fig.1 The flowchart of the proposed method

笔者经过分析发现，视盘邻域的上下两条静脉血管连接成的形状类似于抛物线，如图2所示。图中用绿色抛物线描绘静脉血管，用蓝色十字叉的交点标识抛物线的顶点，可以发现此抛物线顶点在两条静脉血管的交点处，且落在视盘内部。依据此特点，本研究采用基于静脉血管结构特征的方法初定位视盘。首先通过预处理和形态学运算得到静脉血管的主要轮廓，然后采用最小二乘法进行抛物线拟合，接着对拟合结果进行判定并且初定位视盘。

图2 一幅彩色眼底图像Fig.2 A color fundus image

由于视盘在其邻域中表现出较为明显的高亮度特性，所以本研究依据这一特点采用灰度窗口扫描方法来精确定位视盘。

1.2 图像预处理

1.2.1 选取颜色通道

分别提取彩色眼底图像的RGB等3个颜色通道，如图3所示。可以看出：红色通道中的视盘区域最明显，但是整幅图像的亮度较大、对比度较低，且主血管不清晰；蓝色通道图像整体偏暗，图像对比度最差；而绿色通道图像中的血管和视盘最清晰，对比度较高，因此选取绿色通道图像进行后续处理。对于绿色通道图像中出现的各区域明暗不一致以及图像的渐晕现象，本研究采用文献[10]的方法解决。

图3 眼底图像的RGB通道图像。(a)红色通道；(b)绿色通道；(c)蓝色通道Fig.3 The RGB channel images. (a) Red channel; (b) Green channel; (c) Blue channel

1.2.2 提取静脉血管

视盘粗定位是在眼底静脉血管的基础上进行的，静脉血管具有如下特点：

1)宽度方面。静脉血管的主要轮廓属于主血管，其宽度大于周围的细小血管和静脉血管的分支。利用这点，采用以主血管宽度作为半径的圆盘模板，对预处理图像进行低帽运算；然后采用Otsu二值分割算法，提取静脉血管，得到分割结果BT1。

2)长度方面。静脉血管主要轮廓的长度较大，一般超过图像宽度的1/3，而且静脉血管亮度值较低，即病变区域和视盘周围血管的亮度高于视盘附近两条静脉血管段的亮度。依据此特点，对预处理图像进行多方向线性形态学低帽运算[8]，其定义为

BT2=max(f·Sednπ/12-f)

(1)

式中，Sed为线性模板，长度为L，每15°为一个方向，进而得到12个低帽变换结果。

将包含像素最多的结果作为多方向线性低帽运算结果，然后采用Otsu二值分割算法，提取静脉血管(见图4)，得到分割结果BT2。

为了消除病变和分支血管的干扰，提取到准确的静脉血管主要轮廓，将上述两个分割结果相与，并采用形态学算法提取其骨架线。另外，为了消除可能引入的边界的影响，对得到的二值图像进行掩膜处理，掩膜模板如图4(b)所示。最终得到的静脉血管轮廓如图4(c)所示。

图4 提取静脉血管。(a)彩色眼底图像；(b)掩膜图像；(c)静脉血管的主要轮廓；(d)抛物线拟合结果(染色线)Fig.4 Vein vessels extraction. (a) Retinal image; (b) Template image; (c) Main course of vein vessels; (d) The parabolic fitting result (red curve)

1.3 视盘定位

1.3.1 拟合静脉血管

图5 坐标系结构示意图。(a)坐标系建立；(b)左眼眼底静脉血管骨架；(c)右眼眼底静脉血管骨架Fig.5 The schematic diagram of rectangular coordinate system. (a) Rectangular coordinate system; (b) Main course of vein vessels (left eye); (c) Main course of vein vessels(right eye)

在眼底图中，以图像左上角为原点，水平方向作为Y轴的正方向，竖直向下作为X轴的正方向，建立坐标系，如图5(a)所示。这样静脉血管骨架线上点的行标号映射为X坐标，列标号映射为Y坐标，得到映射后的数据点集(xi,yi)(1≤i≤N)，N为数据点个数。根据眼底图像的成像特点，在左眼的眼底图像中，上下两条静脉血管连接而成的弧形类似于一条开口向上的抛物线，如图5(b)所示；对于右眼的眼底图像，上下两条静脉血管曲线类似于开口向下的抛物线，如图5(c)所示。

抛物线的一般式为一元二次形式，有

f(x)=ax2+bx+c

(2)

因此，本研究采用式(2)拟合静脉血管骨架线映射的数据点(xi,yi)(1≤i≤N)。

为了使拟合出的近似曲线能尽量反映拟合数据的变化趋势，采用最小二乘拟合法，使得所有数据点残差的误差平方和S(a,b,c)最小，有

(3)

当式(3)取最小值时，S对于(a,b,c)的一阶导数为零，由此可得

(4)

求解式(4)即可得到所拟合抛物线方程的参数(a,b,c)。

1.3.2 视盘的初定位

依据眼底血管走向特征，以及所建立的坐标系的结构特点，制定静脉血管拟合评价准则，判定抛物线拟合方法是否有效，并确立视盘初定位方法。

由所建立的坐标系模型以及眼底图像血管结构走向的特点得知，对于拟合的抛物线的顶点，理论上应落在坐标轴的第一象限，因此抛物线顶点的两个坐标值应当均为非负值。于是，将静脉血管骨架线映射的待拟合数据点(xi,yi)(1≤i≤N)以及拟合结果采用两个评估标准进行评价。

1) 将待拟合数据代入方程组式(4)，方程有解；

若两个评价标准不能同时满足，则认为采用抛物线拟合的定位方法失效，拟合失败的主要原因是血管病变。由于血管病变和视盘自身病变同时存在的情况极少，此时依据文献[4]，采用基于血管自身特性的方法初定位视盘。

1.3.3 视盘的精确定位

在正常的彩色眼底图像中，视盘区域亮度最高且形状呈类圆形。在病变图像中，渗出物的亮度以及形状都可能类似视盘，此时仅采用视盘的灰度或形状信息进行视盘定位则准确性较低。但是，在视盘的邻域，视盘仍然呈现两个明显特点：一是视盘区域整体亮度较大，这在LAB颜色空间的A通道图像中较明显，如图6(a)所示；二是视盘与其邻域灰度差异大，由于HSI颜色空间的I分量对图像的轮廓信息反映较好，该特点在HSI颜色空间的I通道中最明显，如图6(b)所示。

图6 两个颜色通道图像。(a)A通道图像；(b)I通道图像Fig.6 Two channel images. (a) A channel image of LAB color space; (b) I channel image of HSI color space

依据上述视盘邻域特点，对初定位的结果进行修正。首先，在原始彩色图像中截取以初定位的位置为中心、长度为两个视盘直径的方形区域，作为预处理图像B；然后，定义一个直径为视盘直径的圆形滑动窗口，在B中从左向右、从上到下、以步长1个像素进行扫描。具体步骤为：

步骤1，在每个滑动窗口内，计算其灰度响应值，有

(5)

式中，(x,y)为滑动窗口的中心坐标，I(x,y)为此窗口图像所对应的LAB颜色空间中A通道图像的灰度均值，δ(x,y)为此窗口图像所对应的HSI颜色空间中I通道图像的灰度方差值。

步骤2，计算响应值最大的滑动窗口所对应的中心坐标(x0,y0)，即为最终的视盘位置。

当视盘定位点位于视盘边界内时，认为定位正确。

1.4 样本及参数设定

从DRIVE[15]、STARE[16]、DIABETED0[17](DB0)和MESSIDOR[18]图库中选取图像样本作为实验材料来进行方法测试。DRIVE图库共40幅图像，包括33幅健康眼底图像和7幅糖尿病病变眼底图像。STARE图库共400幅图像，选取其中81幅图像用作本次实验样本，这些图像是视盘定位算法通用的测试样本，包括31幅健康眼底图像和50幅病变眼底图像。DIARETDB0(DB0)是公开的用于糖尿病视网膜病变检查的标准图库，共130幅眼底图像，包括20幅健康眼底图像和110幅糖尿病视网膜病变眼底图像。MESSIDOR图库共1 200幅彩色眼底图像，包括540幅健康眼底图像和660幅病变眼底图像。图像样本的具体选取情况如表1所示。

表1 实验图像集Tab.1 Experimental images

在所提的方法中，有几个重要的参数，这些参数的合理选取对实验结果起关键作用。一是多方向线性模板中的模板长度L，由于眼底图像的成像特点，以及静脉血管自身的特性，在眼底图像中，静脉血管长度一般超过整幅图像宽度的1/3，因此根据图像大小自动设定L为图像宽度的1/3。二是决定低帽变换模板大小的主血管宽度，以及精确定位时决定扫描窗口大小的视盘直径D。在本研究中，这些参数根据彩色眼底图像中独特的比例关系[19]来设定，视盘的直径同眼底图像宽度大小之间的比例为1/8～1/5，而主血管宽度同视盘直径的比例为1/7～1/6。例如，DRIVE数据库中的图像大小为768像素×584像素，L自动取为768/3，取整得256；视盘的直径D自动取为768/7，而主血管宽度取D的1/7。由此，本方法依据彩色眼底图像的大小自动获取参数的取值。

2 结果

本方法在上述4个图库中的定位准确率分别为100%、93.8%、98.6%和99.75%。对包含不同对比度图像的DRIVE数据库、包含严重病变图像的STARE数据库、包含各种糖尿病视网膜病变图像的DB0数据库、包含不同眼底病变图像的MESSIDOR数据库，本方法均具有较好的视盘定位效果。例如，选取2幅典型的健康眼底图像和6幅典型病变眼底图像进行视盘定位结果展示(见图7)，结果表明本方法对亮度和病变均具有较强的适应性。

图7 不同类型眼底图像视盘定位结果。(a)较暗的健康眼底；(b)较亮的健康眼底；(c)存在渗出物病变的眼底；(d)存在渗出物和出血点病变的眼底；(e)存在大量出血的眼底；(f)存在激光斑的眼底；(g)存在大面积亮斑的眼底；(h)视盘发生病变的眼底Fig.7 OD localization results with different fundus images. (a) Dark healthy; (b) Bright healthy fundus; (c) Exudation lesions; (d) Exudation lesions and hemorrhage; (e) Retinal hemorrhage; (f) Laser spots; (g) Bright speck; (h) Optic neuropathy

表2给出了不同方法在STARE数据库中的测试结果，其中对于时间复杂度的测试均在Windows XP系统、CPU主频为2.5 GHz、内存为2 GB条件下进行，编程环境为Matlab 2008。文献[2,5,8]分别利用眼底图像中视盘亮度大、灰度差异大、呈类圆形特征定位视盘，但由于STARE样本中包含较为严重的病变图像，而视盘的自身特性受病变影响，因此这些定位的准确率均不高。文献[1,7,10]利用血管和视盘的关系定位视盘，由于血管特征稳定，所以这类方法的视盘定位准确率有所提高，但定位精度也低于本方法。文献[9,12]的共同特点是充分利用了血管的整体特征，检测准确率高达98.8%，但这都需要制定较为严格的几何模板，而且模板制定依赖于对血管特性的精确分割，算法的复杂度高，耗时长，文献[12]对单幅图像定位视盘用时达3.5 min，文献[9]为2 min。本方法在STARE库中的检测率为93.8%，用时却仅为3 s，耗时远远低于其他方法。另外，由表3可知，相比其他方法，本方法在DB0和MESSIDOR数据库具有更高的定位正确率。

表2 不同方法在STARE数据库中的测试结果比较Tab.2 The comparison of localization results in STARE dataset

表3 不同方法在DB0和MESSIDOR数据库中的视盘定位准确率比较

Tab.3 The comparison of OD localization success rate in DB0 and MESSIDOR datasets

方法准确率/%DB0MESSIDOR基于投影[1]95.5—基于亮度[4]—98.7基于多特征[13]96.92—基于圆形检测[6]97.7—基于投票[14]—98.33本方法98.699.75

3 讨论

本研究提出的视盘定位方法充分利用了彩色眼底图像中的眼底结构特征。根据眼底图像中静脉血管的结构特点，提出了基于静脉血管的视盘定位方法。由于本方法只需粗略提取静脉血管的骨架结构，并不需要精确地提取血管特征，所以降低了算法的复杂度。与依赖精确血管分割结果的传统的视盘定位方法不同，本方法只要依据局部血管的主要轮廓信息，就可以克服传统的基于血管的视盘定位方法中必须精确提取整体血管或者主血管特征的局限性，避免血管分割结果在一定程度上影响视盘检测正确率的问题。如图8所示，由于病变影响，血管的特征模糊，传统方法需要依赖精确的血管分割结果，所以针对这类图像定位失败，如文献[1]基于血管在视盘部分的分散度特征进行视盘定位，结果定位失败(见图8(a))；但由于图8中血管的主要骨架轮廓特征仍然保留，所以对此类图像本方法仍能定位成功(见图8(b))。此外，本研究依据眼底血管走向特征以及所建立的坐标系的结构特点，制定静脉血管拟合评价准则，并根据评价结果融合视盘自身特性以及血管特性初定位视盘，保证了算法的鲁棒性和适应性。

图8 视盘定位结果比较。(a)文献[1]定位结果；(b)本算法定位结果Fig.8 The comparison of OD localization results. (a) OD localization result of Ref[1]; OD localization result of the proposed method

虽然本研究提出的方法相比其他方法在4个实验数据库中整体定位效果较好，但仍存在一定的局限性。针对静脉血管信息缺失、主要轮廓信息不完整的情况(如图9所示，眼底图像仅出现上半部分主血管)，通过最小二乘抛物线拟合法所拟合出的抛物线顶点不落在视盘邻域，导致视盘定位出错。因此，需要深入研究抗干扰性更强、稳定性更高的抛物线拟合算法，以进一步提升算法的鲁棒性。

图9 本算法定位失败的结果示例Fig.9 Example of incorrect result

4 结论

本研究提出一种基于眼底图像结构特征的视盘定位方法。首先，分析了眼底静脉血管结构的特点，通过预处理和形态学运算得到静脉血管主要轮廓；然后，采用最小二乘法，对静脉血管进行抛物线拟合；接着，依据眼底血管走向特征以及的建立的坐标系的结构特点，制定静脉血管拟合评价准则，并根据评价结果辅助视盘定位；最后，基于视盘在其邻域中的高亮度特性，通过灰度窗口扫描方法修正初定位结果，得到最终的定位结果。本方法在分析眼底结构特点的基础上，充分利用视盘的自身特性以及视盘与血管的位置关系来定位视盘。定位精度高，运算复杂度低。

[1] 张东波, 易瑶, 赵圆圆. 基于投影的视网膜眼底图像视盘检测方法[J]. 中国生物医学工程学报, 2013, 32(4): 477-483.

[2] Sinthanayothin C, Boyce JF, Cook HL, et al. Automated localisation of the optic disc, fovea, and retinal blood vessels from digital colour fundus images[J]. British Journal of Ophthalmology, 1999, 83(8): 902-910.

[3] Li Huiqi, Chutatape O. Automatic location of optic disk in retinal images[C] // IEEE International Conference on Image Processing. Thessaloniki: IEEE, 2001, 2: 837-840.

[4] 邹北骥, 张思剑, 朱承璋. 彩色眼底图像视盘自动定位与分割[J]. 光学精密工程, 2015, 23(4): 1187-1195.

[5] Haar FT. Automatic localization of the optic disc in digital colour images of the human retina[D]. Utrecht: Utrecht University, 2005.

[6] Ravishankar S, Jain A, Mittal A. Automated feature extraction for early detection of diabetic retinopathy in fundus images[C]// IEEE Conference on Computer Vision and Pattern Recognition. Miami: IEEE, 2009: 210-217.

[7] Tobin KW, Chaum E, Govindasamy VP, et al. Detection of anatomic structures in human retinal imagery[J]. IEEE Transactions on Medical Imaging, 2007, 26(12): 1729-1739.

[8] 郑绍华, 陈健, 潘林, 等. 基于定向局部对比度的眼底图像视盘检测方法[J]. 中国生物医学工程学报, 2014, 33(3): 289-296.

[9] Foracchia M, Grisan E, Ruggeri A. Detection of optic disc in retinal images by means of a geometrical model of vessel structure[J]. IEEE Transactions on Medical Imaging, 2004, 23(10): 1189-1195.

[10] Hoover A, Goldbaum M. Locating the optic nerve in a retinal image using the fuzzy convergence of the blood vessels [J]. IEEE Transactions on Medical Imaging, 2003, 22(8): 951-958.

[11] 朱琳琳, 唐延东. 基于眼底特征的视盘自动检测[J]. 仪器仪表学报, 2010, 31(8): 118-122.

[12] Youssif A, Ghalwash A, Ghoneim A. Optic disc detection from normalized digital fundus images by means of a vessels’ direction matched filter [J]. IEEE Transactions on Medical Imaging, 2008, 27(1): 11-18.

[13] Godse DA, Bormane DS. Automated localization of optic disc in retinal images[J]. International Journal of Advanced Computer Science and Applications, 2013, 4(2): 65-71.

[14] Harangi B, Hajdu A. Detection of the optic disc in fundus images by combining probability models[J]. Computers in Biology and Medicine, 2015, 65:10-24.

[15] Max V. DRIVE: Digital retinal images for vessel extraction[EB/OL]. http://www.isi.uu.nl/Research/Databases/DRIVE, 2014-03-12/2015-03-18.

[16] Hoover A. Structured analysis of the retina[EB/OL]. http://www.ces.clemson.edu/～ahoover/stare, 2015-05-18/2015-06-09.

[17] Kauppi T, Kalesnykiene V, Kamarainen JK, et al. DIARETDB0-Standard Diabetic Retinopathy Database Calibration level 0[EB/OL]. http://www2.it.lut.fi/project/imageret/diaretdb0, 2014-09-06/2015-03-07.

[18] Klein J C. Methods to evaluate segmentation and indexing techniques in the field of retinal ophthalmology (MESSIDOR) [EB/OL]. http://www.adcis.net/en/Download-Third-Party/Messidor.html, 2015-01-11/2015-06-20.

[19] Tasman W, Jaeger EA. Duane’s Ophthalmology [M]. (15thEdition). Philadelphia: Lippincott Williams & Wilkins, 2009.

Optic Disc Localization in Color Fundus Images Based on Fundus Structure Feature

Xiao Zhitao1,4Shao Yiting1Zhang Fang1*Wen Jia1Geng Lei1,4Wu Jun1Shang Dandan1Su Long2Shan Chunyan3

1(SchoolofElectronicsandInformationEngineering,TianjinPolytechnicUniversity,Tianjin300387,China)2(TheSecondHospitalofTianjinMedicalUniversity,Tianjin300211,China)3(MetabolicDiseaseHospitalofTianjinMedicalUniversity,Tianjin300070,China)4(TianjinKeyLaboratoryofOptoelectronicDetectionTechnologyandSystem,Tianjin300387,China)

The size and the shape of optic disc (OD) in fundus images are very important for diagnosing the fundus diseases, therefore the detection and localization of OD have important significance for computer-aided diagnosis of the ophthalmology diseases. In this paper, the OD localization method in color fundus images based on the characteristics of fundus structure was proposed. Firstly, vein vessels in fundus image were extracted based on the bot-hat transformation. Then, parabolic fitting algorithm based on the least square method was applied to get the preliminary localization of OD. Finally, the accurate localization of OD was determined through window-scanning. Comparative experimental results were obtained from four publicly available datasets (DRIVE, DIABETED0, STARE and MESSIDOR). Testing results were 100%, 98.6%, 93.8%, 99.75% respectively, which proved the accuracy and robust of the proposed method.

optic disc localization; bot-hat transformation; parabolic fitting; window scanning

10.3969/j.issn.0258-8021. 2016. 03.001

2015-11-12， 录用日期:2016-01-05

国家自然科学基金(61401439)；天津市科技支撑计划重点项目(13ZCZDGX02100)；天津市应用基础与前沿技术研究计划一般项目(15JCYBJC16600)

R318

A

0258-8021(2016) 03-0257-07

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