SLM & Wavefront Control2026 年 9 月 14 日 · 4 min read

基于 Zernike 多项式与共轭梯度法的 SLM 二维平顶光束整形与 λ/100 波前像差矫正手记

从菲涅尔衍射正反向建模出发，用 CUDA 加速的共轭梯度（CG）算法在 SLM 上生成 10×10 光镊阵列所需的二维平顶光束；再以 Zernike 多项式对系统像差做闭环拟合矫正，把波前残差 RMS 压到 λ/100 量级。含完整 PyTorch 实现与踩坑记录。

SLMCGZernikeTop-HatCUDA

0. 实验目标#

构建 10×10（并预留扩展到 20×20）规模均匀光镊阵列，核心瓶颈在两件事：

能量域：把高斯入射光整形为方形二维平顶（Top-Hat）分布，均匀性目标 $\sigma/\mu \le 3\%$ （ $\mu$ 为焦面能量均值， $\sigma$ 为标准差）；
相位域：消除光路系统像差，使波前残差 RMS 达到 $\lambda/100$ 量级 —— 否则每个光镊位点都会继承不同的相位倾斜，重排后阵列亮度立刻「花掉」。

为什么不用 Dammann 光栅

Dammann 光栅设计简单、对像差宽容，但衍射效率受二值相位限制（理论 ~80%，实测往往低于 70%），且均匀性对制造误差敏感。连续相位优化（CG / GS 变体）在效率与均匀性上都明显占优，代价是需要精确的波前模型 —— 这正是本文第二部分闭环矫正的意义。

1. 理论框架：菲涅尔衍射的正反向模型#

SLM 相位面到焦面的传播用单 FFT 菲涅尔积分描述。设 SLM 面场为 $U_0(x,y)=A_0(x,y)\,e^{i\phi(x,y)}$ ，焦面场为

U_f(u,v) = \frac{e^{ikf}}{i\lambda f}\, e^{\frac{i\pi}{\lambda f}(u^2+v^2)} \iint U_0(x,y)\, e^{\frac{i\pi}{\lambda f}(x^2+y^2)}\, e^{-\frac{i2\pi}{\lambda f}(xu+yv)} \, \mathrm{d}x\,\mathrm{d}y ,

其中 $f$ 为透镜焦距， $\lambda = 1064\,\mathrm{nm}$ 。注意两个二次相位因子：一个乘在 SLM 面上（ $e^{i\pi(x^2+y^2)/\lambda f}$ ），一个乘在焦面上 —— 单 FFT 形式的代价就是这两个 chirp，但可以把整个传播写成对角算符 + 傅里叶算符的复合：

U_f = Q_2 \cdot \mathcal{F}\{\, Q_1 \cdot U_0 \,\}, \qquad Q_1 = e^{\frac{i\pi}{\lambda f}(x^2+y^2)},\quad Q_2 = e^{\frac{i\pi}{\lambda f}(u^2+v^2)} .

伴随算符（反向传播）严格等于 $Q_1^{*} \cdot \mathcal{F}^{-1}\{ Q_2^{*} \cdot \bullet \}$ ，这是 CG 收敛的数学保证 —— 正反向必须是严格互伴随的，数值上差一个 chirp 采样误差，收敛曲线就会平台化。

2. 相位优化：从 Gerchberg–Saxton 到共轭梯度#

2.1 目标泛函#

GS 交替投影在平顶整形里收敛慢且容易陷入局部极小。我们把整形写成显式优化问题：

\phi^* = \arg\min_{\phi} \underbrace{\tfrac{1}{2}\,\big\|\, |U_f(\phi)|^2 - I_{\text{target}} \,\big\|_2^2}_{\text{焦面能量分布误差}} \; + \; \alpha\, \mathcal{R}(\phi),

其中 $I_{\text{target}}$ 是目标平顶强度（方形，边缘按 erf 过渡以放松数值刚度）， $\mathcal{R}$ 为正则项（我们用的是相位二阶差分，抑制高频「散斑化」相位）。

2.2 CUDA 加速的共轭梯度实现#

对上千次随机初始化做批量优化时，PyTorch + CUDA 的优势是碾压级的。核心代码骨架：

import torch

class FresnelPropagator(torch.nn.Module):
    """单 FFT 菲涅尔传播（对角算符 + FFT），batch 维并行。"""
    def __init__(self, n_pix: int, dxy: float, focal: float, wavelength: float,
                 device: str = "cuda"):
        super().__init__()
        k = 2 * torch.pi / wavelength
        x = (torch.arange(n_pix, device=device) - n_pix // 2) * dxy
        xx, yy = torch.meshgrid(x, x, indexing="ij")
        self.Q1 = torch.exp(1j * torch.pi / (wavelength * focal) * (xx**2 + yy**2))
        self.Q2 = torch.exp(1j * torch.pi / (wavelength * focal) * (xx**2 + yy**2))
        self.fx = torch.fft.fftfreq(n_pix, d=dxy)
        self.scale = 1j * wavelength * focal  # 频率→焦面坐标

    def forward(self, u0: torch.Tensor) -> torch.Tensor:
        # u0: (B, N, N) complex64
        v = torch.fft.fft2(self.Q1 * u0)
        fu = (torch.fft.fftshift(v, dim=(-2, -1))
              * self.Q2)                       # 焦面复振幅（已 shift 到中心）
        return fu

    def adjoint(self, fu: torch.Tensor) -> torch.Tensor:
        v = torch.fft.ifftshift(fu, dim=(-2, -1)) / self.Q2
        return torch.conj(self.Q1) * torch.fft.ifft2(v)


def conjugate_gradient_phase(target_amp: torch.Tensor, propagator: FresnelPropagator,
                             n_iter: int = 200, restarts: int = 256) -> torch.Tensor:
    """批量 CG：B 个随机初相并行优化，取损失最小者。"""
    B = restarts
    phi = (torch.rand(B, *target_amp.shape[-2:], device=target_amp.device) * 2 * torch.pi)
    phi.requires_grad_(True)

    # 共轭梯度作用在实值相位参数上，梯度由 autograd 沿复振幅链路自动获得
    opt = torch.optim.LBFGS([phi], max_iter=n_iter, history_size=24,
                            line_search_fn="strong_wolfe")

    def closure():
        opt.zero_grad()
        u0 = torch.exp(1j * phi)
        amp = propagator(u0).abs()
        loss = ((amp - target_amp) ** 2).flatten(1).sum(1).mean()
        loss.backward()
        return loss

    opt.step(closure)

    with torch.no_grad():
        u0 = torch.exp(1j * phi)
        losses = ((propagator(u0).abs() - target_amp) ** 2).flatten(1).sum(1)
        best = losses.argmin()
    return phi[best]

光路校准注意：chirp 采样与 FFT shift

三个高频踩坑：(1) $Q_1$ 、 $Q_2$ 的网格坐标必须用 fftfreq/fftfreq 对应的同一坐标系，混用 fftshift 前后的坐标会让正向「看起来对」但伴随算符错位；(2) SLM 像素间距 $d_{xy}$ 与焦面采样 $\Delta u = \lambda f/(N d_{xy})$ 的换算要做一次解析验证（点源打点测试），否则平顶的物理尺寸直接错；(3) batch 维用 complex64，显存随 $B \times N^2$ 增长，1080p SLM、 $B=256$ 时约 4 GB，够用但别再翻倍。

2.3 数值结果#

256 组随机初始化并行优化 200 步（RTX 4090，单轮 ~11 s），最优解的焦面能量分布：

指标	GS 基线（2000 步）	CG / LBFGS（200 步）
平顶均匀性 $\sigma/\mu$	8.9%	2.6%
衍射效率（平顶区内能量占比）	71.4%	84.1%
优化耗时（含 256 restarts）	~40 min（CPU）	~11 s（GPU）

实测在 SLM 上回放相位后，CCOS 相机直接测得 $\sigma/\mu \approx 3.4\%$ （模拟 2.6% → 实测 3.4% 的差值主要来自系统像差，这正是下一节要解决的）。

TOP-HAT (CG) — 图 1 · SLM 相位回放前后的焦面光束轮廓对比（模拟数据示意，拖动滑杆查看）：左为含像差的高斯光斑，右为 CG 优化相位后的二维平顶分布。

RAW GAUSSIAN — 图 1 · SLM 相位回放前后的焦面光束轮廓对比（模拟数据示意，拖动滑杆查看）：左为含像差的高斯光斑，右为 CG 优化相位后的二维平顶分布。

3. 波前像差闭环矫正：Zernike 模式空间#

3.1 像差模型#

把系统波前像差展开为 Zernike 多项式：

W(\rho, \theta) = \sum_{j=3}^{J} a_j Z_j(\rho, \theta), \qquad W_{\mathrm{RMS}} = \sqrt{\sum_j a_j^2},

$Z_j$ 采用 Noll 索引（ $j=1,2$ 的 piston 与 tilt 不影响焦面强度）。矫正目标即把 $W_{\mathrm{RMS}}$ 压到 $\lambda/100$ 。

3.2 闭环流程#

我们的闭环不依赖夏克-哈特曼传感器，而是以焦面形态为观测量的间接闭环：

在 SLM 上叠加试探性 Zernike 相位 $\delta\phi_j = \beta_j Z_j$ （幅度 ±0.1 rad 阶梯扫描）；
每步拍摄焦面光斑，计算锐度度量 $S = \sum |I - I_{\text{blur}}|^2$ 与 Strehl 比代理量；
信赖域内拟合 $S(\{a_j\})$ 的二次响应面，更新 $\{a_j\}$ ；
迭代直至波前残差判据 $\hat{W}_{\mathrm{RMS}}/\lambda \le 0.01$ 。

判据验证

闭环收敛后，用「相位恢复双焦面法」交叉验证：在两个离焦距离（±2% f）分别采集光斑，用Transport-of-Intensity方程反演波前。两种方法给出的 $W_{\mathrm{RMS}}$ 分别为 0.011λ 和 0.009λ —— 平均 ~λ/100，达标。

3.3 主要像差分量（矫正前）#

Zernike 模式	物理来源	矫正前幅度 (λ)
$Z_4$ Defocus	透镜-焦面离焦失配	0.18
$Z_5, Z_6$ Astigmatism	反射镜夹角应力	0.11 / 0.07
$Z_7, Z_8$ Coma	SLM 与透镜光轴横向偏心	0.06 / 0.04
$Z_{11}$ Spherical	透镜球差	0.09

闭环 6 轮后各分量均压缩到 ≤ 0.004λ，总残差 $\hat{W}_{\mathrm{RMS}} \approx \lambda/105$ 。

4. 阵列验证与下一步#

把平顶相位 + Zernike 补偿相位叠加后经 Dammann 型分裂相位（×10×10 复制）生成 100 个光镊位点，实测位点间功率分布：

阵列均匀性 $\sigma/\mu$ ：3.1%（未加像差矫正时为 9.8%）
平均单镊功率：350 mW 输入下 ~2.6 mW/site
稳定性：1 小时漂移 < 1.5%（受限于 SLM 相位 flicker，见下）

SLM 相位闪烁与 2π 回绕

(1) 液晶 SLM 的相位 flicker（~kHz 级、幅度 ~0.02λ）是均匀性长时漂移的主导项，锁相触发相机采样可以部分平均掉；(2) 优化相位回放前必须确认 2π 回绕映射与 SLM 的标定 LUT 一致 —— 我们第一版直接把连续相位 mod 2π 上载，均匀性凭空差了 2%，排查了一周才发现是 LUT 在 2π 附近的非单调段。

下一步：把该相位管线接入 AOM / 光纤耦合排查日志里的功率稳定链路，并探索三维阵列的轴向分裂相位设计。里德堡激发光路（420/480 nm）的整形需求与 1064 nm 差异较大，届时另开一篇记录。

相关阅读：魔术态蒸馏与表面码解码侧的理论笔记见 15-to-1 Bravyi-Kitaev 魔术态蒸馏协议笔记。