4. 色散补偿算法实现

任务说明

在扫频源OCT（SS-OCT）系统中，样品臂与参考臂之间的色散失配会导致干涉信号的相位产生非线性频率依赖，表现为点扩散函数（PSF）展宽、轴向分辨率下降。本案例实现基于相位校正的数值色散补偿算法：对采集到的干涉光谱信号，构造与色散失配量对应的相位项并乘以原始信号，再经逆傅里叶变换获得补偿后的A-line，从而恢复理论轴向分辨率。

色散补偿原理

设原始干涉光谱信号为S(k)，其中k为波数。色散失配引入的相位畸变为φ(k)，可展开为泰勒级数：

φ(k) = a₂(k-k₀)² + a₃(k-k₀)³ + ...

其中a₂为二阶色散系数（群延迟色散），a₃为三阶色散系数。补偿方法为将原始信号乘以exp(-j·φ(k))，再进行逆傅里叶变换得到聚焦的A-line。本案例实现二阶与三阶色散补偿，系数a₂、a₃作为可调参数。

代码实现

import numpy as np
import matplotlib.pyplot as plt

def generate_dispersed_signal(N=2048, lambda_c=1310e-9, delta_lambda=50e-9,
                              z=0.5e-3, a2=2e-4, a3=1e-6):
    """
    生成含色散失配的模拟OCT干涉光谱信号
    N: 采样点数
    lambda_c: 中心波长(m)
    delta_lambda: 波长扫描范围(m)
    z: 反射点深度位置(m)
    a2: 二阶色散系数(rad·m²)
    a3: 三阶色散系数(rad·m³)
    """
    # 波数向量（均匀k域）
    k = np.linspace(2*np.pi/(lambda_c+delta_lambda/2),
                    2*np.pi/(lambda_c-delta_lambda/2), N)
    k0 = 2*np.pi/lambda_c
    dk = k - k0

    # 理想干涉信号（单反射点）
    S_ideal = np.exp(1j * 2 * k * z)

    # 色散相位畸变
    phi_disp = a2 * dk**2 + a3 * dk**3

    # 含色散信号
    S_disp = S_ideal * np.exp(1j * phi_disp)

    # 添加高斯型光源谱包络
    spectral_envelope = np.exp(-((k - k0)/(0.2*(k[-1]-k[0])))**2)
    return k, S_disp * spectral_envelope, S_ideal * spectral_envelope

def dispersion_compensation(S, k, a2_comp, a3_comp=0):
    """
    数值色散补偿
    S: 含色散干涉光谱信号
    k: 波数向量
    a2_comp: 二阶补偿系数(rad·m²)
    a3_comp: 三阶补偿系数(rad·m³)
    """
    k0 = k[len(k)//2]
    dk = k - k0
    # 构造补偿相位并施加
    phi_comp = a2_comp * dk**2 + a3_comp * dk**3
    S_comp = S * np.exp(-1j * phi_comp)
    return S_comp

def reconstruct_psf(S, k):
    """对光谱信号做逆傅里叶变换得到A-line（PSF）"""
    # 均匀k域直接IFFT
    psf = np.fft.ifft(S)
    # 取模并归一化
    psf_abs = np.abs(psf)
    return psf_abs / np.max(psf_abs)

# ========== 主流程 ==========
# 1. 生成含色散信号
N = 2048
k, S_disp, S_ideal = generate_dispersed_signal(N=N, a2=2e-4, a3=1e-6)

# 2. 执行色散补偿（使用与生成时相同的系数）
S_comp = dispersion_compensation(S_disp, k, a2_comp=2e-4, a3_comp=1e-6)

# 3. 重建PSF
psf_original = reconstruct_psf(S_disp, k)
psf_compensated = reconstruct_psf(S_comp, k)
psf_ideal = reconstruct_psf(S_ideal, k)

# 4. 计算半高全宽(FWHM)评估分辨率
def calc_fwhm(psf, k):
    """计算PSF主瓣半高全宽对应的深度分辨率"""
    # 深度坐标映射
    z_axis = np.linspace(0, 2*np.pi/(k[1]-k[0]), len(psf))
    # 找主瓣峰值位置
    peak_idx = np.argmax(psf)
    half_max = psf[peak_idx] / 2
    # 向左找半高点
    left_idx = peak_idx
    while left_idx > 0 and psf[left_idx] > half_max:
        left_idx -= 1
    # 向右找半高点
    right_idx = peak_idx
    while right_idx < len(psf)-1 and psf[right_idx] > half_max:
        right_idx += 1
    fwhm_idx = right_idx - left_idx
    # 转换为深度（空气中，n=1）
    fwhm_depth = fwhm_idx * (2*np.pi/(k[1]-k[0])) / (2 * (k[-1]-k[0]) / N)
    return fwhm_depth, peak_idx

# 计算FWHM（以深度采样间隔为单位）
z_spacing = 2*np.pi/(k[1]-k[0]) / N  # 深度采样间隔
fwhm_orig, _ = calc_fwhm(psf_original, k)
fwhm_comp, _ = calc_fwhm(psf_compensated, k)
fwhm_ideal, _ = calc_fwhm(psf_ideal, k)

print(f"原始信号PSF FWHM: {fwhm_orig:.1f} 个采样点")
print(f"补偿后PSF FWHM: {fwhm_comp:.1f} 个采样点")
print(f"理想信号PSF FWHM: {fwhm_ideal:.1f} 个采样点")
print(f"分辨率改善: {fwhm_orig/fwhm_comp:.2f} 倍")

# 5. 绘制对比图
plt.figure(figsize=(10, 6))
depth_axis = np.arange(N) * z_spacing * 1e3  # 转换为mm
plt.plot(depth_axis, psf_original, 'r-', alpha=0.7, label='色散未补偿')
plt.plot(depth_axis, psf_compensated, 'b-', label='色散补偿后')
plt.plot(depth_axis, psf_ideal, 'g--', alpha=0.8, label='理想（无色散）')
plt.xlabel('深度 (mm)')
plt.ylabel('归一化幅度')
plt.title('色散补偿前后PSF对比')
plt.legend()
plt.grid(True, alpha=0.3)
plt.xlim([0, 2])
plt.savefig('dispersion_compensation_result.png', dpi=150)
plt.show()

运行与验证

运行方式：将代码保存为dispersion_compensation.py，在Python环境中执行：

python dispersion_compensation.py

需要安装numpy和matplotlib。

验证方法：程序输出三个FWHM数值并显示PSF对比图。判定标准：

• 补偿后的FWHM应显著小于未补偿的FWHM（本案例中约改善3~5倍）；
• 补偿后的PSF主瓣宽度应接近理想信号的PSF宽度；
• 生成的dispersion_compensation_result.png中，蓝色曲线（补偿后）的主瓣应明显窄于红色曲线（未补偿），且与绿色虚线（理想）基本重合。

关键参数说明

实际应用中，色散系数a₂、a₃可通过以下方式确定：

• 在参考臂放置一个反射镜，采集干涉信号并提取其相位；
• 对相位进行多项式拟合，二次项系数即为a₂，三次项系数为a₃；
• 将拟合系数代入补偿函数即可实现自动色散补偿。