偏振频域光学相干层析成像光谱校准方法与流程

本发明涉及频域光学相干层析成像(fourier-domainopticalcoherencetomography，简称fd-oct)，特别是一种偏振频域光学相干层析成像光谱校准方法，更确切地说是一种自动寻找波数k的最优校准系数进行偏振频域光学相干层析成像光谱校准的方法。

背景技术：

光学相干层析成像(opticalcoherencetomography，oct)是一种非侵入、基于低相干光干涉的生物医学成像技术。早期的时域oct需要移动参考臂来探测样品的不同深度位置信号。随后发展起来的频域oct(fd-oct)无需移动参考臂，只需一次探测即可获得样品内部的深度信息(a-line信号)，再进行横向扫描即可获得样品的纵切面图像(b-scan图像)，具有高速高分辨的特点。偏振频域oct(fourierdomainpolarization-sensitiveoct，fd-ps-oct)(参见在先技术[1]，michaelr.heeanddavidhuang，“polarization-sensitivelow-coherencereflectometerforbirefringencecharacterizationandranging”，j.opt.soc.am.b，vol.9，no.6，903-908，1992)是在此基础上发展起来的一种功能oct技术，在得到双折射样品强度图像的同时可获得样品的延迟图像和快轴图像，能够额外提供样品的双折射信息，已被广泛应用在皮肤、牙科和眼科等领域。fd-oct通过傅里叶变换对样品不同深度位置一次成像，傅里叶变换中深度位置z对应波数k(k＝2π/λ，λ为波长)。但在fd-oct中，ccd等光电探测器探测的是波长等间隔的数据。要得到准确的b-scan图像，需要对波数k进行等间隔数据处理，这就需要进行准确的波长校准，即准确映射光谱仪探测器的像素与波长之间的对应关系。错误的波长映射会使点扩展函数展宽，产生类似色散对点扩展函数的影响。因此，传统的fd-oct需要对波长进行校准。而fd-ps-oct除获得强度图像外，还要定量描述样品的延迟和快轴等偏振信息，波长校准不准确会引入干涉相位误差，导致偏振参数计算精度下降，因此精确的波长校准更为重要。另外，fd-ps-oct为了获得样品偏振参数，至少需要同时采集正交的水平偏振干涉谱h和垂直偏振干涉谱v两路干涉信号。而两路信号由于各自波长范围或光谱分辨率不同等原因，对于同一深度位置，会出现两路a-line信号错位现象，使相位计算误差增大。因此除上述波长校准外，还需要对h和v进行对准(参见在先技术[2]，m.wojtkowski,r.leitgeb,a.kowalczyk,t.bajraszewski,anda.f.fercher，“invivohumanretinalimagingbyfourierdomainopticalcoherencetomography”，j.biomed.opt.vol.7，no.3，457–463，2002)，使两路信号保持相同的k分布，提高计算精度。

wojtkowski等人(参见在先技术[2])第一次指出fd-oct中光谱校准的重要性。park等人指出错误的波长校准会产生双折射计算误差(参见在先技术[3]，b.h.park,m.c.pierce,b.cense,s.h.yun,m.mujat,g.j.tearney,b.e.bouma,andj.f.deboer,“real-timefiber-basedmulti-functionalspectral-domainopticalcoherencetomographyat1.3μm”，opt.express，13(11),3931–3944，2005)。在典型的偏振频域oct中，干涉信号分为h和v两路，即使很小的波长错位也会产生双折射相关参数的计算误差。这样光经过一段传播距离后，误差会累积增大，可能导致无法与样品本身产生的双折射参数区分。mujat等人提出采用h和v两路a-line信号的峰值位置来求出波数k错位的比例关系，然后以一路为基准，校准另一路信号(参见在先技术[4]，mujatm,parkbh,censeb,etal,“autocalibrationofspectral-domainopticalcoherencetomographyspectrometersforinvivoquantitativeretinalnervefiberlayerbirefringencedetermination,”journalofbiomedicaloptics，12(4),041205-041205-6，2007)。该方法适合高斯谱型较好的光源，因为a-line信号具有较好的单一峰值，可以准确求出k的比例系数。但对于非高斯谱型光源，a-line有多个峰值位置或者最高峰位置不在中心位置处，会导致校准误差较大。

技术实现要素：

本发明提供一种偏振频域光学相干层析成像光谱校准方法，更确切地说，是一种自动寻找波数k的最优校准系数进行偏振频域光学相干层析成像光谱校准的方法，对采集到的干涉谱，先粗略校准再进行精确校准，精确校准以一路信号为基准，通过自动寻找k的比例系数来对另一路信号进行校准，并采用两路信号的互相关结果作为自动寻找的标准评判最优校准系数。本发明克服了上述在先技术的不足，采用相关度来评价结果，对光源的谱型没有要求，对高斯型、非高斯型光谱均适用。

本发明的技术解决方案如下：

一种自动寻找波数k的最优校准系数进行偏振频域光学相干层析成像光谱校准的方法，该方法先对各个光谱进行光电探测器像素与波长对应关系的粗略校准，然后以一路光谱(v)的波数为基准，将另一路光谱(h)的波数乘以最优系数实现精确校准；最优系数通过评判光谱v和校准后光谱h的a-line信号相关度来自动确定，互相关最大时，认为校准效果最好，此时对应的系数即为最优校准系数；样品测量时，用得到的最优校准系数对光谱h的波数k校准，然后用校准后的k对h进行重采样，结合评判光谱v，实现样品偏振参数的计算。

本发明的自动寻找波数k的最优校准系数进行偏振频域光学相干层析成像光谱校准的方法，具体步骤如下：

①偏振频域oct系统探测的两路正交干涉信号，分别为水平偏振干涉谱ihr和垂直偏振干涉谱ivr；对于每一路干涉谱，分别进行光电探测器像素与波长λ对应关系的粗略校准；

②采用一简单样品，光电探测器探测该样品的水平干涉光谱ih和垂直干涉光谱iv：

其中，下标h和v分别表示水平干涉光谱h和垂直干涉光谱v；ih0表示水平干涉光谱的直流项和自相关项；iv0表示垂直干涉光谱的直流项和自相关项；s(k)表示光源光谱；rr表示参考镜的反射率；rsn表示样品第n层的反射率；δzn表示样品臂各层与参考臂之间的光程差；δn表示从样品表面到样品第n层的延迟量；θn表示样品第n层的快轴方向；

③判断光电探测器探测到的数据ih和iv是否是以波数等间隔分布的，如果数据是波数等间隔分布的，则不用对波数k进行插值预处理，直接转到步骤④；如果数据不是波数等间隔分布的，则对波数k进行插值预处理，得到等间隔的水平干涉光谱的波数kh和等间隔的垂直干涉光谱的波数kv；

④对等间隔的水平干涉光谱的波数kh进行缩放：选取一组校准系数p＝[p1，p2，…，pi，…，pn-1，pn]，其中，i，n为自然数，表示校准系数p的序号，以等间隔的垂直干涉光谱的波数kv为基准，通过公式k'h＝kh+pi·(kh-kh0)对水平干涉光谱的波数kh进行缩放，得到一组新的波数k'h，其中kh0是kh的中心波数，“﹒”表示相乘，用新的k'h对水平干涉光谱ih进行重采样，ih变为关于k'h的干涉信号ih(k'h)：

对重采样后的ih(k'h)信号进行傅里叶变换，得到z域中重采样后的水平干涉光谱的a-line信号；

同样，判断光电探测器探测到的数据iv是以波数等间隔分布的，是，则直接对垂直干涉光谱iv进行傅里叶变换等处理，得到z域中垂直干涉光谱的a-line信号；光电探测器探测到的数据iv不是以波数等间隔分布的，则先用步骤③得到的垂直干涉光谱的波数kv对垂直干涉光谱iv进行重采样后再进行傅里叶变换等处理，得到z域中垂直干涉光谱的a-line信号；

⑤最后对重采样后的水平干涉光谱的a-line信号和垂直干涉光谱的a-line信号进行互相关计算，得到校准系数pi对应的互相关值c(pi)；重复上述过程，计算得到p中每一个系数pi对应的互相关值，当两路a-line信号互相关最大时认为校准效果最好，即max(c(pi),pi∈p)＝c(po)，此时所用的校准系数即为最优系数po，按下列光束计算最优系数po校准后的水平干涉光谱的波数k'hpo：

k'hpo＝kh+po·(kh-kh0).(3)

⑥经过上述校准后，z域中的水平干涉光谱的a-line信号和垂直干涉光谱的a-line信号可以完全重叠，但还需要对k域中干涉信号的谱型进行进一步校准：用步骤⑤得到的k'hpo对水平干涉光谱ih进行重采样，然后进行傅里叶变换，得到校准后的水平干涉光谱对应的a-line信号，用滤波器滤出与简单样品相关的信号；同时用滤波器滤出步骤④中得到的垂直干涉光谱的a-line信号中与样品相关的信号；对滤出的这两路与样品相关的信号分别进行傅里叶逆变换，分别得到水平干涉光谱h和垂直干涉光谱v在k域中的信号轮廓图；对这两幅信号轮廓图作互相关计算，互相关最大时，水平干涉光谱h的信号轮廓图移动的量即为k'hpo需要移动的像素数nshift(kv保持不变)，按公式(4)将移动的像素数k'hpo转化成移动量kshift：

kshift＝nshift·[(k'hpomax-k'hpomin)/n]，(4)

其中，k'hpomax和k'hpomin分别表示k'hpo的最大值和最小值，n表示k'hpo的总数据量，即光电探测器的像素数；

最终校准后的kh通过下式(5)计算出并记为k″hpo：

k″hpo＝k'hpo+kshift.(5)

⑦测量任意待测样品，获得样品的两路正交干涉信号后，用步骤⑥得到的k″hpo对样品的h路信号进行重采样，得到样品校准后的h光谱；对样品校准后的h光谱和样品的v光谱进行常规的傅里叶变换等步骤对偏振频域oct数据处理，得到样品的结构信息、延迟信息和快轴方向信息。

所述的粗略校准，可以采用光谱特征曲线法或干涉谱比较法等方法。光谱特征曲线法是通过测量得到几个已知的特征谱线，并找出与其对应的光电探测器像素，通过这种对应关系直接将多项式系数解出，初步确定各自的波长范围。干涉谱比较法采用商业光谱仪和系统中需要校准的光谱仪分别探测一样品反射镜对应的干涉谱，通过比较两个干涉谱的波峰和波谷得到多组波长和像素的对应关系，然后分析得到多项式系数，实现光谱校准。粗略校准不限定于上述两种方法，可以是任何其他能够获得光电探测器像素与波长对应关系的校准方法。

所述的光电探测器，是ccd或cmos阵列或其他的具有光电信号转换功能的探测器阵列。

所述的简单样品，可以是平面反射镜或其他类似样品。

所述的对波数k进行插值预处理，为提高校准精度，可以对k进行4倍插值处理或其他倍数的插值处理，对相关的干涉谱也要进行相应的补零插值处理。

所述的校准系数p，是预先给出的系数范围，与系统中光电探测器的光谱响应特性有关。当水平干涉光谱的波数kh和垂直干涉光谱的波数kv相差较大时，增大p的取值范围以保证可以找到两路的波数最佳对应关系，反之可以减小p的取值范围。

所述的以垂直干涉光谱的波数kv为基准，不限定于以kv为基准，也可以是以水平干涉光谱的波数kh基准，步骤中涉及到的内容都相应的把h和v对调。

与在先技术相比，本发明具有以下优点：

1.本发明采用相关度来评价校准结果，对光源谱型没有要求，对高斯型和非高斯型光谱均适用，拓展了应用范围，更具有普适性。

2.本发明采用互相关最大值作为评价函数，避免了寻找峰值位置时的误差，可以达到相同或者更高的校准精度。

附图说明

图1为本发明所采用的光谱校准的流程图。

图2为本发明实施例1中的非高斯光源光谱。

图3为本发明实施例1中两路偏振通道h和v的a-line信号。

图4为本发明实施例1中的校准前后的h和v的a-line信号，其中，图(a)为采用不同方法校准前后的h和v的a-line信号，图(b)为图(a)中虚线方框b部分的放大图，图(c)为图(a)中虚线方框c部分的放大图，图(d)为图(a)中虚线方框d部分的放大图。

图5为发明实施例1中延迟量和快轴方位角的计算误差绝对值，其中图(a)为延迟量的计算误差绝对值，图(b)为快轴方位角的计算误差绝对值。

图6为本发明实施例1中延迟量和快轴方位角的误差绝对值的平均值和标准差，其中图(a)为延迟量的误差绝对值的平均值和标准差，图(b)为快轴方位角的误差绝对值的平均值和标准差。

图7位本发明实施例2中的高斯光源光谱。

图8为本发明实施例2中延迟量和快轴方位角的误差绝对值的平均值和标准差，其中(a)为延迟量的误差绝对值的平均值和标准差，图(b)为快轴方位角的误差绝对值的平均值和标准差。

具体实施方式

下面结合实施例和附图对本发明作进一步说明，但不应以此实施例限制本发明的保护范围。

图1是本发明所采用的光谱校准的流程图。由图可见，本发明自动寻找波数k的最优校准系数进行偏振频域光学相干层析成像光谱校准的方法，该方法包括如下步骤：

①偏振频域oct系统探测的两路正交干涉信号，分别为水平偏振干涉谱ihr和垂直偏振干涉谱ivr；对于每一路干涉谱，分别进行光电探测器像素与波长λ对应关系的粗略校准。

②采用一简单样品，光电探测器探测得到该样品的h和v两路正交干涉信号，分别为水平干涉光谱ih和垂直干涉光谱iv，

其中，下标h和v分别表示水平干涉光谱h和垂直干涉光谱v；ih0表示水平干涉光谱的直流项和自相关项；iv0表示垂直干涉光谱的直流项和自相关项；s(k)表示光源光谱；rr表示参考镜的反射率；rsn表示样品第n层的反射率；δzn表示样品臂各层与参考臂之间的光程差；δn表示从样品表面到样品第n层的延迟量；θn表示样品第n层的快轴方向。

③判断光电探测器探测到的数据ih和iv是否是以波数等间隔分布的：如果数据是波数等间隔分布的，则不用对波数k进行插值预处理，直接转到步骤④；如果数据不是波数等间隔分布的，则对波数k进行插值预处理，得到等间隔的水平干涉光谱的波数kh和等间隔的垂直干涉光谱的波数kv；

④对等间隔的水平干涉光谱的波数kh进行缩放：选取一组校准系数p＝[p1，p2，…，pi，…，pn-1，pn]，其中，i，n为自然数，表示校准系数p的序号；以等间隔的垂直干涉光谱的波数kv为基准，通过公式k'h＝kh+pi·(kh-kh0)对水平干涉光谱的波数kh进行缩放，得到一组新的波数k'h，其中kh0是kh的中心波数，“﹒”表示相乘；用新的k'h对水平干涉光谱ih进行重采样，ih变为关于k'h的干涉信号ih(k'h)：

对重采样后的ih(k'h)信号进行傅里叶变换，得到z域中重采样后的水平干涉光谱的a-line信号；

同样，对垂直干涉光谱iv进行处理，得到z域中垂直干涉光谱的a-line信号；

⑤最后对重采样后的水平干涉光谱的a-line信号和垂直干涉光谱的a-line信号进行互相关计算，得到校准系数pi对应的互相关值c(pi)；重复上述过程，计算得到p中每一个系数pi对应的互相关值。当两路a-line信号互相关最大时认为校准效果最好，即max(c(pi),pi∈p)＝c(po)。此时所用的校准系数即是要求出的最优系数po。通过po得到校准后的水平干涉光谱的波数k'hpo：

k'hpo＝kh+po·(kh-kh0).(8)；

⑥经过上述校准后，z域中的水平干涉光谱的a-line信号和垂直干涉光谱的a-line信号可以完全重叠，但还需要对k域中干涉信号的谱型进行校准。用步骤⑤得到的k'hpo对水平干涉光谱ih进行重采样，然后进行傅里叶变换，得到校准后的水平干涉光谱对应的a-line信号，用滤波器滤出与简单样品相关的信号；同时用滤波器滤出步骤④中得到的垂直干涉光谱的a-line信号中与样品相关的信号；对滤出的这两路与样品相关的信号分别进行傅里叶逆变换，分别得到水平干涉光谱h和垂直干涉光谱v在k域中的信号轮廓图，对这两幅信号轮廓图作互相关计算，互相关最大时，水平干涉光谱h的信号轮廓图移动的量即为k'hpo需要移动的像素数nshift(kv保持不变)。移动的像素数转化成k'hpo的移动量kshift：

kshift＝nshift·[(k'hpomax-k'hpomin)/n]，(9)

其中，k'hpomax和k'hpomin分别表示k'hpo的最大值和最小值，n表示k'hpo的总数据量，即光电探测器的像素数。

最终校准后的kh通过下式计算并记为k″hpo：

k″hpo＝k'hpo+kshift.(10)；

⑦测量任意待测样品，获得样品的两路正交干涉信号后，用步骤⑥得到的的k″hpo对样品的h路信号进行重采样，得到样品校准后的h光谱；对样品校准后的h光谱和样品的v光谱进行常规的傅里叶变换等偏振频域oct数据处理，得到样品的结构信息、延迟信息和快轴方向信息。

实施例中需要计算简单样品的强度图像r(z)(结构信息)，延迟图像δ(z)(延迟信息)和快轴图像θ(z)(快轴方向信息)，可以通过下边公式获得：

其中，ah,v(z)表示校准后的水平干涉光谱h和垂直干涉光谱v各自对应的a-line信号的振幅，φh,v(z)表示校准后的水平干涉光谱h和垂直干涉光谱v各自对应的a-line信号的相位，下标分别表示h和v两通道。

所述的粗略校准，可以采用光谱特征曲线法或干涉谱比较法等方法。光谱特征曲线法是通过测量得到几个已知的特征谱线，并找出与其对应的光电探测器像素，通过这种对应关系直接将多项式系数解出，初步确定各自的波长范围。干涉谱比较法采用商业光谱仪和系统中需要校准的光谱仪分别探测一样品反射镜对应的干涉谱，通过比较两个干涉谱的波峰和波谷得到多组波长和像素的对应关系，然后分析得到多项式系数，实现光谱校准。粗略校准不限定于上述两种方法，可以是任何其他能够获得光电探测器像素与波长对应关系的校准方法。

所述的光电探测器，是ccd或cmos阵列或其他的具有光电信号转换功能的探测器阵列。

所述的简单样品，可以是平面反射镜或其他类似样品。

所述的对波数k进行插值预处理，为提高校准精度，可以对k进行4倍插值处理或其他倍数的插值处理，对相关的干涉谱也要进行相应的补零插值处理。

所述的校准系数p，是预先给出的系数范围，与系统中光电探测器的光谱响应特性有关。当水平干涉光谱的波数kh和垂直干涉光谱的波数kv相差较大时，增大p的取值范围以保证可以找到两路的波数最佳对应关系，反之可以减小p的取值范围。

所述的以垂直干涉光谱的波数kv为基准，不限定于以kv为基准，也可以是以水平干涉光谱的波数kh基准，步骤中涉及到的内容都相应的把h和v对调。

本实施例1，光源光谱为非高斯光谱，如图2所示。图中横坐标表示波长，纵坐标表示相对光强。本实施例中，光源中心波长为835nm，半高全宽为45nm。样品为3层，代表样品与参考臂的不同光程差位置，即实际深度位置，每一层与参考臂的光程差分别为200μm,700μm,1200μm。每层对应的延迟量δ分别为0.4969rad，0.8065rad，1.1282rad，对应的快轴方位角θ分别为1.2226rad，0.7083rad，0.4501rad。光电探测器采用ccd。理论上对应同一深度位置，h和v两路的a-line信号应重叠在一起。但由于两路通道的光谱分辨率不同、衍射角度不同等原因，h和v光谱出现错位，如图3所示。错位导致偏振参数计算误差增大，降低成像质量。图中，虚线表示h光谱的a-line信号，实线表示v光谱的a-line信号。

按照本发明上述校准方法的步骤对实施例1中的h和v光谱进行校准。校准后的a-line信号结果如图4所示。图4中“h”和“v”表示校准前的h和v两路的a-line信号，“peaks”表示采用峰值法(在先技术[4])校准后得到的h光谱的校准结果，“cor”表示本发明提供的最优系数法校准得到的h光谱校准结果，“△”表示实际深度位置。图(b)、(c)和(d)分别是图(a)中对应虚线方框部分的放大图。从图中可以看出，经过两种方法校准后，h的a-line信号中心都基本与v的a-line信号中心重合，但采用本发明方法校准后，两路信号重叠得更好，h的信号形状更接近于v的信号，在图(c)中尤为明显。图(c)中在峰值位置处，采用本发明方法校准后，变化趋势更倾向于v的信号，校准效果更好。

校准后计算实际深度位置处的延迟量δ和快轴方位角θ，各自的计算误差绝对值如图5所示。这里只需要得到计算值与真值之间的偏移量，因此本实施例中的误差计算都是针对误差的绝对值进行计算。从图中可以看出，在接近零光程差位置处，本发明校准方法的延迟量误差值小于峰值法，快轴方位角的计算误差与其接近。随着光程差增大，在第二层位置处，峰值法的延迟量和快轴方位角误差分别为0.025rad和0.033rad，而本发明分别为0.015rad和0.005rad，优势明显，能够更精确地校准h和v光谱。在非高斯光谱中，傅里叶变换后获得的a-line信号的峰值位置可能偏离实际深度位置，利用峰值方法校准可能会有较大的误差。在图5(b)的第三层位置处，可以明显看到，利用峰值法校准后的快轴方位角的计算误差较大，误差绝对值达到了0.091rad。而采用本发明方法校准后，快轴方位角的计算误差绝对值降到了0.029rad，精度能够提高50％以上。

图6是实验重复进行20次后获得的数据，利用该数据计算得到的延迟量和快轴方位角的误差绝对值的平均值和标准差。“peaks”是峰值方法计算得到的误差结果，“cor”是用本发明提供的方法计算得到的误差结果。图6(a)是延迟量δ的误差情况，峰值法和本发明提供的方法分别计算得到的误差平均值都在10^-2级次上，但是本发明方法误差值明显小于峰值方法。在第3层位置处，峰值方法的标准差平均值为0.00041rad，而本发明为0.00026rad，减少了36.6％，具有更高的测量准确性和稳定性。图6(b)是快轴方位角θ的误差情况，在接近零光程差位置处，本发明提供的校准方法的误差平均值和标准差与峰值法接近。随着光程差增大，两种方法的误差值都在增大。在第3层位置处，峰值法的误差平均值为0.092rad，本发明方法为0.030rad，远小于峰值方法，能够更精确地计算样品的快轴方位角，具有更高的校准精度。重复实验验证了本发明适用于光源为非高斯谱型的情况，能够提高校准精度。

实施例2，光源为高斯光谱，其他参量与实施例1一样，保持不变。实验重复进行20次。图7是实施例2所用的光源光谱，为高斯型。图8是计算的延迟量和快轴方位角的误差绝对值的平均值和误差结果。图8(a)是延迟量的误差绝对值的平均值和标准差，误差平均值都在10^-3量级上，误差很小，从图中可明显看出本发明优于峰值法。图8(b)为快轴方位角的计算误差，随着光程差增大，误差平均值和标准差都变大，与实施例1中的非高斯谱型类似。两种方法的标准差都在10^-4量级上，峰值法的误差平均值最大约为0.041rad，本发明方法约为0.026rad，比峰值法精度提高了36.6％。

实施例1和实施例2验证了本发明提供的方法的可行性。本发明不仅适用于高斯光谱，在非高斯光谱中也有较好表现，计算精度比峰值法有明显提高。