阅读说明：本技术 基于矩形窗函数优化的光学相干层析成像信号处理方法 (Optical coherence tomography signal processing method based on rectangular window function optimization ) 是由 张茜 李中梁 王向朝 南楠 何东航 欧阳君怡 杨晨铭 刘腾 于 2020-08-14 设计创作，主要内容包括：在利用光学相干层析成像系统对被测样品进行成像时,样品臂与参考臂的色散失配会影响系统分辨率,因此进行深度分辨的色散补偿是非常重要的。此时需要将样品不同深度处的干涉信号从A-line信号中截取再做傅里叶变换(FFT),通过对截取所用的矩形窗函数进行优化,可有效减小信号加窗FFT中加窗产生的相位误差,准确提取被测样品在不同深度处的色散系数,提高色散补偿精度。(When an optical coherence tomography system is used for imaging a measured sample, the dispersion mismatch between a sample arm and a reference arm affects the system resolution, so that depth-resolved dispersion compensation is very important. Interference signals at different depths of the sample are intercepted from the A-line signal and then subjected to Fourier transform (FFT), and a rectangular window function used for intercepting is optimized, so that phase errors generated by windowing in the signal windowing FFT can be effectively reduced, dispersion coefficients of the sample to be measured at different depths are accurately extracted, and dispersion compensation precision is improved.)

基于矩形窗函数优化的光学相干层析成像信号处理方法

技术领域

本发明涉及频域光学相干层析成像(Fourier Domain Optical CoherenceTomography，简称FD-OCT)技术，尤其是一种基于矩形窗函数优化的光学相干层析成像信号处理方法。

背景技术

光学相干层析成像(OCT)通过探测样品背向散射光的强度来获取样品的结构信息，是一种高分辨、非侵入、可在体检测组织内部微观结构的光学断层成像技术。频域OCT技术通过对频域干涉谱信号的逆傅里叶变换(IFT)获得被测物体的层析图，具有微米或亚微米量级的空间分辨率。但是在高分辨率系统中由于采用的超宽光谱光源带宽超过100nm甚至达到几百nm，会导致色散效应加剧。色散会导致OCT相干信号的展宽与畸变，使系统实际分辨率小于理论值，因此色散补偿是实现高分辨率OCT的关键技术之一。

色散的算法补偿是通过对OCT获取的数据进行后处理来消除色散展宽，具有灵活、便捷等优势。主要方法包括解卷积算法、迭代算法、自聚焦算法和全深度色散补偿方法等。由于大部分被测组织存在深度变化的色散，张仙玲、黄炳杰、A.F.Fercher等人提出不同种用于频域光学相干层析成像的深度分辨色散补偿方法，通过动态滤波，分离出被测样品各个深度的干涉谱，进而提取出各层深度的色散系数进行分别补偿，可以在全深度范围内获得较好的色散补偿效果。这类方法需要提取样品不同深度的信号进行FFT分析相位信息，得到不同深度的二阶及以上的色散系数，所以通常需要选用窗函数对样品结构信号进行截取得到某一深度的信号，再进行后续信号处理。在进行深度分辨色散补偿时，需要分离出不同深度处的干涉信号，加窗后仅要求精确读出主瓣频率，不考虑幅值精度；而矩形窗主瓣窄、旁瓣大、频率识别精度最高，幅值识别精度最低，因此使用矩形窗。但是在加窗截取信号过程中，如果窗口宽度过小，可能会使某些结构层的信息丢失，如果窗口宽度太宽则可能削弱后处理深度色散补偿算法的效果。并且为了更精确分析不同深度处频谱信息，加窗截取局部深度处的干涉信号时，在原则上需要尽可能增大窗口宽度，以提高频谱的频率分辨率。因此在对不同深度处干涉信号的相位提取时合适的窗口选择是非常重要的。

发明内容

本发明的目的在于提供一种基于矩形窗函数优化的光学相干层析成像信号处理方法。在利用光学相干层析成像系统对被测样品进行成像时，样品臂与参考臂的色散失配会影响系统分辨率，因此进行深度分辨色散补偿是非常重要的。此时需要将样品不同深度处的干涉信号从A-line信号中截取后再做FFT，通过对截取所用的矩形窗函数进行优化，即将初级矩形窗得到的干涉信号FFT后提取相位，根据存在色散时相位与波数的函数关系将得到的相位用最小二乘法拟合二次多项式，并以拟合二次多项式的标准误差为判定条件，不断改变矩形窗口长度和中心位置，其中最小标准误差即对应最优化的矩形窗，该矩形窗下截取的干涉信号FFT后可得到逼近真实的相位。减小矩形窗的引入所带来的相位误差，更准确提取被测样品在不同深度处的色散系数，来提高色散补偿精度。

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

一种基于矩形窗函数优化的光学相干层析成像信号处理方法，其特征在于该方法包括如下步骤：

①首先利用频域光学相干层析成像系统对样品进行扫描，该系统的光电探测阵列记录样品原始的干涉信号并传送计算机；

②对原始干涉信号去背景后做逆傅里叶变换，得到A-line信号；

③对A-line信号中任一深度的信号加初级矩形窗进行信号截取，并对截取到的信号进行傅里叶变换后求相位角得到该深度信号对应的初级相位；

④根据色散时相位与波数的函数关系，以初级相位为目标函数，用最小二乘法拟合得到初级相位与波数的二次多项式，并计算拟合结果的标准误差；

⑤改变初级矩形窗的窗口长度和中心位置，得到当前矩形窗；对A-line信号加当前矩形窗截取该深度范围内的信号，并对截取到的信号进行FFT后求相位角得到该深度信号对应的相位，以得到的相位为目标函数，用最小二乘法拟合得到相位与波数的二次多项式，并计算拟合结果的标准误差，即该矩形窗对应的相位的标准误差；然后继续改变矩形窗的窗口长度和中心位置，截取该深度范围内的信号，在矩形窗遍历所有设定情况后得到不同矩形窗对应的相位的标准误差；

⑥步骤⑤中得到的标准误差最小值所对应的矩形窗即为最优矩形窗，利用该最优矩形窗即求出该深度信号对应的最终相位；

⑦由该深度信号的最终相位计算对应的二阶色散系数；

⑧对A-line信号中不同深度的信号重复步骤③④⑤⑥⑦，得到不同深度信号对应的二阶色散系数，并计算由色散引起的相位偏差；从该深度对应的相位中减去偏差量，得到色散补偿后的相位；

⑨利用各个深度补偿后的相位重建出色散补偿后的频域干涉信号。

所述步骤③中初级矩形窗的窗口长度和中心位置选取方法为：

无样品时，A-line信号的平均幅值为噪声的均值；

设置样品时，A-line信号中某一深度的信号的极大值，以该极大值所在的深度为初级矩形窗的中心位置，并向两侧扩展窗口，直到两侧的信号强度等于噪声的均值，此时初级矩形窗的窗口长度为该深度信号的窗口长度。

所述步骤⑤中改变矩形窗的窗口长度和中心位置方法为：

以A-line信号的单位长度为一个步长，将初级矩形窗起始位置向右移动1个步长，窗口末端位置依次向左和向右从1个步长移到N个步长，共得到2N个矩形窗，截取每个窗口下的干涉信号进行傅里叶变换提取相位，重复上述操作直到窗口起始位置向右移动N个步长；同样的，将初级矩形窗起始位置向左移动1个步长，窗口末端位置依次向左和向右从1个步长移到N个步长，共得到2N个矩形窗，截取每个窗口下的干涉信号进行傅里叶变换提取相位，重复到窗口起始位置向左移动N个步长结束，得到所有不同矩形窗截取信号时该深度信号对应的相位，N取15-30。

所述步骤⑦中计算对应深度的二阶色散系数，具体是：

步骤7.1对该深度处的干涉信号进行傅里叶变换得到最终的相位信息

式中k表示波数，Δz_n为样品第n层相对于参考臂反射镜的光程差，k₀为光源中心波长对应的波数，n_n是样品第n层处的有效折射率，n_g,n是样品第n层处的有效群折射率，a₂是样品第n层处的二阶色散补偿系数。

步骤7.2以(k-k₀)为自变量，对进行多次项数值拟合，得到的二次相位项即为该深度处信号对应的二阶色散系数。

本发明基于矩形窗函数优化的光学相干层析成像信号处理方法的特点是针对被测样品在进行深度分辨色散补偿时，对样品不同深度处干涉信号做加窗FFT，对其A-line信号截取所用的单一阈值确定的初级矩形窗函数的窗口长度和中心位置进行优化，不断逼近相位的真实值，减小信号加窗FFT中窗口引入的相位误差，以此得到样品在不同成像深度处的最优色散系数，达到更好的色散补偿效果。

本发明的技术解决方案原理：

目前对数字干涉信号时-频域相位提取的方法主要为傅里叶变换法，但是傅里叶变换是一种全局变换，较适用分析平稳信号；而在对非平稳信号进行分析时，会产生不同局域的频谱混叠，以致于不能精确地提取基频。加窗傅里叶变换能够有效地使频谱局域化，改善不同级次频谱的叠加现象，提高测量的精确度，但是如果加窗傅里叶变换中窗口大小固定不变，也会在频谱信息的分析中引入误差。

理想的FFT要求时域信号是无限长的，而在实际计算测量中，只能对有限长的信号进行FFT，这样就相当于对无限长的信号进行截断。如果输入信号是无限长的，那么FFT所得到的频谱就是完全正确的频谱，但是如果输入信号是截断的有限长记录样本，并且所截取的信号又不是整周期时，离散谱线向两边展宽就会产生频谱泄露。由于矩形函数的傅里叶变换谱在两个方向都无限扩展，会产生振铃现象，则信号加矩形窗后用FFT算法计算频谱会产生旁瓣效应。所以加窗截取产生的频谱泄漏会使谱线模糊，谱分辨率降低，而且泄露的程度随截取长度的增加而降低，截取长度越长，泄露越小，谱分辨率越高。而当截取长度较短时，由于矩形窗函数谱存在许多旁瓣，从而使窗函数频谱与主信号频谱卷积以后也形成了许多旁瓣，这些旁瓣起谱间干扰作用，同样也会降低谱分辨率。

在对样品光学相干层析成像进行深度分辨色散补偿时，需要选取适当的矩形窗的窗口长度和中心位置来截取不同深度处的干涉信号。如果窗口宽度过小，可能会使某些结构层的信息丢失，如果太宽则可能导致窗口宽度太宽削弱深度色散补偿效果，又因为减小频谱泄漏需要截取的长度越长才能提高谱分辨率，所以选择合适的矩形窗长度和中心位置是至关重要的，这样才能还原更接近真实的频谱。

在频域OCT中通过探测样品背向散射光进行成像，光源发出的低相干光通过参考臂与样品臂分别照射到反射镜和样品上，从反射镜返回的参考光与从样品不同深度处返回的样品光发生干涉，去除直流项，样品的干涉信号为：

式中k表示波数，Re表示取复数的实部，I_n(k)表示样品第n层散射回的光强，I_r(k)表示反射镜返回的光强，Δz_n为样品第n层相对于参考臂反射镜的光程差，是样品第n层散射光相对于参考光的相位差，包括高阶色散相位Φ(k,Δz_n)。高阶色散相位的引入是导致干涉信号包络展宽和畸变、降低系统分辨率的主要原因，色散补偿的目的就是消除高阶色散相位。

第n层散射光相对于参考光的相位差

可以表示为：

式中β_n(k)是样品第n层处的有效传播系数，对于参考臂和样品臂中各光学元件已经给定的OCT系统来说，β_n(k)可以在光源中心波长对应的波数k₀附近做泰勒级数展开，于是可得：

其中n_n是样品第n层处的有效折射率，n_g,n是样品第n层处的有效群折射率，β”_n、β”'_n分别是样品第n层处的二阶有效色散系数、三阶有效色散系数，a₂、a₃分别称为二阶色散补偿系数、三阶色散补偿系数。

样品不同深度处有不同的色散系数β”_n和β”'_n。为了实现深度分辨色散补偿，需要先得到z_n处的相位，将频域OCT干涉谱信号作IFT变换得到样品的A-line信号，然后通过矩形窗函数截取得到不同深度z_n处的结构层信号包络，将这个包络信号沿深度z轴作FFT得到不同深度z_n处的干涉谱信号，提取得到相位。将z_n处的相位信息进行数值拟合，即得到二阶、三阶有效色散系数，计算色散引入的相位并从总的相位中减去即可得到补偿后的干涉谱信号。再对补偿后的干涉谱信号进行逆傅里叶变换得到样品无色散的A-line信号。

通常在实际光学相干层析成像系统中三阶色散以及更高阶色散的影响很小，与二阶色散项带来的影响相比可以忽略不计，所以色散补偿主要是消除干涉信号中存在的二阶色散相位，则上式得到的相位差可以表示为：

可以看出理论上相位为波数的二次多项式。传统做法中利用单一阈值确定的窗口对A-line信号进行截取，这种情况下某些深度处截取的信号会由于不合适的窗口长度和中心位置造成提取相位存在误差，影响后续算法计算精度。所以为了更精确提取样品在不同深度的相位，先利用初级矩形窗截取A-line信号得到某一深度相位，根据存在色散时相位与波数的函数关系将得到的相位用最小二乘法拟合二次多项式，以拟合二次多项式的标准误差为判定条件，不断改变矩形窗长度和中心位置，当标准误差最小时即对应最优化的矩形窗，该矩形窗截取干涉信号后进行FFT可得到逼近真实的相位。在这里初级矩形窗的确定方法为找到该深度的A-line信号高于噪声最大值的极大值，以极大值为中心向两侧扩展窗口，直到两侧的信号强度等于噪声的均值，此时的窗口宽度为该极大值处的初级窗口宽度，中心位置为极大值处。

本发明具有以下优点：

在光学相干层析成像的深度色散补偿中，通过优化矩形窗函数得到更准确的相位信息，可有效减小干涉信号加窗引入的相位误差。

此方法也适用于在光学相干层析成像干涉信号相位分析提取中其他窗函数的优化。

附图说明

图1为本发明模拟样品结构的A-line图，其中(a)为无色散A-line图，(b)为有色散A-line图。

图2为模拟样品无色散时不同深度的相位误差图，其中(a)为初级矩形窗截取信号并处理后得到的相位误差图，(b)为优化窗口截取信号并处理后得到的相位误差图。

图3为模拟样品有色散时利用初级矩形窗截取信号并处理后得到的不同深度的相位误差图。

图4为模拟样品有色散时随深度变化的相位图，其中(a)为初级窗口截取信号并处理后得到的相位图，(b)为优化窗口截取信号并处理后得到的相位图。

图5为模拟样品色散补偿后的A-line图，其中(a)为初级窗口截取信号并处理后得到二阶色散系数进行色散补偿效果图，(b)为优化窗口截取信号并处理后得到二阶色散系数进行色散补偿效果图。

具体实施方式

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

利用计算机模拟分析本发明光学相干层析成像信号处理中加窗FFT的矩形窗函数优化方法的可行性。图(1)为一个有16层散射层的样品结构A-line图，样品中不同层玻璃都具有色散性质，并且每一层的厚度有所差异，参阅图(a)无色散时的样品A-line图，样品各结构层信号清晰，可以看到在样品中厚度薄的层两个信号峰相离非常近，厚度相对比较厚的层两个信号峰则相离较远。图(b)为引入色散时的样品A-line图，可以看出由于色散效应导致样品各结构层信号展宽，某些深度的结构层已经混叠，无法分辨开。模拟该样品无色散时，对干涉信号加初级矩形窗截取信号并处理后得到不同深度相位误差参阅图(2)中(a)所示，其中不同线型表示样品在不同深度的相位误差，可以看出加窗FFT在各个深度频谱分析时都会引入误差，误差为水平波动在5rad范围内。这是由于频谱泄漏和旁瓣效应是不可避免的。其中实线和点划线表示的两个深度处由于初级矩形窗的窗口长度和中心位置的不合适选择引起的相位误差很大，其最大误差随波数增加到接近35rad。经过窗函数优化方法提取信号，参阅图(2)中(b)所示，这两个深度处相位误差大大减小，误差分别降至0.3rad和2rad附近，波动在正常范围5rad范围内。再模拟引入色散情况下，根据存在色散时相位与波数的函数关系。理论上相位误差与波数呈二次曲线，频谱泄漏和旁瓣效应会引入一些波动，参阅图(3)可以看出由于色散的引入相位误差会随波数增加而变大，在合适的矩形窗下相位误差会随波数增加到26-42rad，但是初级矩形窗的窗口长度和中心位置的不合适选择会增大误差，如图中实线、虚线和点划线所示，这些深度的相位与实际相位偏差较大。对初级窗口优化并提取二阶色散系数具体步骤如下：

(1)对原始干涉信号去背景后做逆傅里叶变换，得到A-line信号；

(2)对A-line信号中任一深度的信号加初级矩形窗进行信号截取，并对截取到的信号进行傅里叶变换后求相位角得到该深度信号对应的初级相位；

(3)根据色散时相位与波数的函数关系，以初级相位为目标函数，用最小二乘法拟合得到初级相位与波数的二次多项式，并计算拟合结果的标准误差；

(4)改变初级矩形窗的窗口长度和中心位置，得到当前矩形窗；对A-line信号加当前矩形窗截取该深度范围内的信号，并对截取到的信号进行FFT后求相位角得到该深度信号对应的相位，以得到的相位为目标函数，用最小二乘法拟合得到相位与波数的二次多项式，并计算拟合结果的标准误差，即该矩形窗对应的相位的标准误差；然后继续改变矩形窗的窗口长度和中心位置，截取该深度范围内的信号，在矩形窗遍历所有设定情况后得到不同矩形窗对应的相位的标准误差；

(5)步骤(4)中得到的标准误差最小值所对应的矩形窗即为最优矩形窗，利用该最优矩形窗即求出该深度信号对应的最终相位；

(6)由该深度信号的最终相位计算对应的二阶色散系数；

(7)对A-line信号中不同深度的信号重复步骤(2)(3)(4)(5)(6)，得到不同深度信号对应的二阶色散系数，并计算由色散引起的相位偏差；从该深度对应的相位中减去偏差量，得到色散补偿后的相位；

(8)利用各个深度补偿后的相位重建出色散补偿后的频域干涉信号。

在引入色散情况下，将初级窗口和优化窗口截取信号进行FFT计算得到随深度变化的相位，理论上由可知相位与深度为线性关系，参阅图(4)可知初级窗口得到的相位在Z1、Z2、Z3和Z4深度处相位出现误差，对比图(a)和图(b)发现优化窗口提取的相位更准确。将初级窗口和优化窗口提取到的相位计算二阶色散系数进行色散补偿，参阅图(5)可以发现由初级窗口截取信号得到的各个深度的二阶色散系数在对应的这四个层并不准确，所以色散补偿效果较差甚至没有效果，而经过窗函数优化方法截取信号所得到的相位误差小，对该样品每一深度的相位都能准确还原，所以提取的随深度变化的二阶色散系数更准确，样品在各个结构层的色散补偿效果好，每一结构层都能分辨出。

以上所述只是本发明的一个具体实施例，该实施例仅用以说明本发明的技术方案而非对本发明的限制。凡本领域技术人员依本发明的构思通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在本发明的保护范围之内。