4. OCT光谱信号色散补偿与相位校正

色散补偿与相位校正：从PSF展宽到分辨率恢复

在完成波数域线性化后，OCT光谱干涉信号仍可能因系统色散不平衡而出现相位非线性。理想情况下，参考臂与样品臂的色散应完全匹配，干涉信号在波数域中表现为等间隔余弦条纹，其FFT结果是一个接近理论极限的窄点扩散函数（PSF）。但实际系统中，光纤长度差异、透镜材料色散、分光棱镜色散等因素会引入随波数变化的附加相位 φ_disp(k)，使条纹相位偏离线性，导致PSF展宽、旁瓣升高、轴向分辨率下降。

色散补偿的核心思路是：从点目标（或强反射面）的干涉信号中提取相位残差，拟合出非线性色散系数，构造共轭相位进行校正。常用模型为泰勒展开：

φ_disp(k) ≈ a₂(k - k₀)² + a₃(k - k₀)³ + a₄(k - k₀)⁴ + ...

其中 k₀ 为中心波数，a₂ 为二阶色散系数（单位 rad·mm²），a₃ 为三阶色散系数，更高阶项在多数系统中可忽略。补偿时，将光谱信号乘以 exp(-j·φ_disp(k))，再执行FFT，即可得到压缩后的PSF。

下面给出一个可运行的Python示例，使用模拟的色散展宽光谱信号，演示相位提取、多项式拟合、补偿和PSF评估的完整过程。示例参数为教学假设值，用于说明方法逻辑，不代表真实系统实测结果。

import numpy as np
from scipy.signal import hilbert
from scipy.optimize import curve_fit

# ========== 1. 模拟参数 ==========
k0 = 7.5e6          # 中心波数 (rad/m)，对应约 838 nm
dk = 0.5e6          # 波数范围半宽 (rad/m)
N = 2048            # 采样点数
k = np.linspace(k0 - dk, k0 + dk, N)  # 波数轴 (rad/m)

# 真实色散系数（教学假设值）
a2_true = 1.2e-12   # 二阶色散 (rad·m²)
a3_true = 3.0e-19   # 三阶色散 (rad·m³)

# 理想干涉信号（单点目标，深度 z0）
z0 = 0.5e-3         # 目标深度 0.5 mm
S_ideal = np.cos(2 * z0 * k)  # 理想余弦条纹

# 加入色散相位
phi_disp = a2_true * (k - k0)**2 + a3_true * (k - k0)**3
S_disp = np.cos(2 * z0 * k + phi_disp)

# 加入少量噪声
np.random.seed(42)
S_disp += 0.02 * np.random.randn(N)

# ========== 2. 提取相位并拟合 ==========
# 使用Hilbert变换获取解析信号相位
analytic = hilbert(S_disp)
phi_meas = np.unwrap(np.angle(analytic))

# 去除线性项（对应深度z0），保留非线性残差
# 线性拟合
p_lin = np.polyfit(k - k0, phi_meas, 1)
phi_linear = p_lin[0] * (k - k0) + p_lin[1]
phi_residual = phi_meas - phi_linear

# 对残差做多项式拟合（2阶和3阶）
def poly_model(x, a2, a3):
    return a2 * x**2 + a3 * x**3

popt, _ = curve_fit(poly_model, k - k0, phi_residual, p0=[0, 0])
a2_est, a3_est = popt
print(f"估计色散系数: a2 = {a2_est:.3e}, a3 = {a3_est:.3e}")
print(f"真实色散系数: a2 = {a2_true:.3e}, a3 = {a3_true:.3e}")

# ========== 3. 补偿 ==========
phi_comp = a2_est * (k - k0)**2 + a3_est * (k - k0)**3
S_comp = S_disp * np.exp(-1j * phi_comp)

# ========== 4. FFT与PSF评估 ==========
def compute_psf(signal):
    """对光谱信号做FFT，返回深度轴和归一化幅度"""
    spectrum = np.fft.fft(signal)
    spectrum = np.fft.fftshift(spectrum)
    depth = np.fft.fftshift(np.fft.fftfreq(len(signal), d=(k[1]-k[0])))
    mag = np.abs(spectrum)
    mag = mag / mag.max()
    return depth, mag

depth_ideal, psf_ideal = compute_psf(S_ideal)
depth_disp, psf_disp = compute_psf(S_disp)
depth_comp, psf_comp = compute_psf(S_comp)

# 计算FWHM（以理想PSF为参考）
def fwhm(depth, psf):
    idx = np.argmax(psf)
    half = 0.5 * psf[idx]
    left = idx
    while left > 0 and psf[left] > half:
        left -= 1
    right = idx
    while right < len(psf)-1 and psf[right] > half:
        right += 1
    return abs(depth[right] - depth[left])

print(f"理想PSF FWHM: {fwhm(depth_ideal, psf_ideal)*1e6:.2f} μm")
print(f"色散PSF FWHM: {fwhm(depth_disp, psf_disp)*1e6:.2f} μm")
print(f"补偿后PSF FWHM: {fwhm(depth_comp, psf_comp)*1e6:.2f} μm")

运行上述代码，预期可观察到：色散使PSF明显展宽，FWHM增大；补偿后PSF重新压缩，FWHM接近理想值。实际系统中，色散系数需从实测点目标信号中提取，且可能随温度、光纤弯曲等因素缓慢变化，需定期校准。

补偿效果评估与判断依据

补偿是否有效，不能仅凭视觉观察，需要定量评估。最直接的指标是PSF的FWHM（半高全宽），它直接对应轴向分辨率。理论上，OCT轴向分辨率约为 δz = 2ln2·λ₀²/(π·Δλ)，其中 λ₀ 为中心波长，Δλ 为光源带宽。补偿后FWHM应接近该理论值。

若补偿后FWHM仍显著大于理论值，可能原因包括：色散系数估计不准确、存在高阶色散未补偿、波数域线性化不完善、或光源本身相干长度不足。此时应检查相位拟合残差是否呈随机分布，若仍有系统性弯曲，需增加多项式阶数或重新标定波数轴。