ADMM与光谱近邻算子在定量相位成像中的应用及Matlab实现

做定量相位成像这几年，最让我头疼的不是光学平台，而是重建算法。单波长下跑跑Gerchberg-Saxton或者HIO还能糊弄过去，可一旦把照明换成高光谱宽带光源，同时采集多个波长的衍射强度，问题立刻变得棘手：每个波长既共享样本信息，又各自带着色散和噪声，独立恢复的结果互相打架，相位图连起来看完全是乱的。后来我把主力算法换成ADMM（交替方向乘子法），在波长维上引入光谱近邻算子做约束，整套流程用Matlab实现后，重建稳定性和相位精度都有明显提升。这篇就把我从建模、推导到调参的完整过程写出来，给同样在做高光谱宽带相位恢复和定量相位成像的朋友一条可参考的技术路线，也顺带把Matlab实现里的关键代码和踩坑记录都留在下面。

文章不搬理论证明，重点说清楚三件事：ADMM在高光谱宽带相位恢复里到底怎么拆问题，光谱近邻算子的来龙去脉和代码怎么写，以及定量相位成像里如何从恢复结果换算出你真正关心的物理量。适合正在做衍射成像、无透镜成像、定量相位成像的科研人员和工程师参考，也适合刚接触相位恢复但想跳过漫长试错阶段的研究生。

1. 高光谱宽带相位恢复到底难在哪

1.1 相位恢复的基础测量模型

先捋一下基本设定。假设样本的复透射场是 (x)，照明波长为 (\lambda)，探测器只能记录强度：

[ y = |\mathcal{F}(x)|^2 ]

这里的 (\mathcal{F}) 在具体系统里可能是傅里叶变换，也可能是角谱传播算子。相位恢复的目标，就是从 (y) 反解出 (x) 的幅度和相位。早期方法基本都是迭代投影思路：不断在空间域约束和傅里叶域约束之间来回投影，直到两个约束同时满足。

单波长情形下，这套思路已经有不少坑。最典型的是孪生像问题——恢复结果里经常出现物体翻转叠加的伪影；其次是对初始猜测极其敏感，随机相位起步，十次里有三次收敛到明显错误的解。我当时做单波长衍射成像时，为了压制这些伪影，各种初始化策略、支持域约束、频域过采样全都试过，最后还是靠多次随机重启取最优才勉强稳定。

1.2 单波长算法的短板与多波长信息互补

到了高光谱宽带场景，单个波长的恢复本身就不稳，还要同时处理多个波长，问题复杂度直接乘上通道数。这里有三个现实难点：

• 每个波长独立恢复，得到的相位图互相不一致。因为色散关系，同一物理位置的样本在不同波长下相位延迟本来就不同，但如果完全各自为政，恢复出来的相位差和理论色散曲线对不上，后续定量换算就没法做。
• 波长之间的信息冗余没有被利用。相邻波长的复振幅场在绝大多数光学样本上高度相关，这种相关性是天然的约束条件，经典交替投影法完全没考虑。
• 宽带照明下衍射图案的噪声结构更复杂。不同波长的量子效率、光源功率都不一样，有些波长信噪比高，有些低，独立处理时低信噪比波长会把整次重建拖垮。

多波长数据其实给了你一条额外的路：既然样本是同一个，光谱维上的解不应该剧烈跳变。只要用光谱近邻算子在波长维上施加平滑约束，就能让强波长帮弱波长、高信噪比通道带动低信噪比通道，获得比单波长恢复更稳的结果。这也是我最终转向ADMM的核心原因——它天然适合把"每个波长的数据保真"和"光谱维上的结构约束"拆成两个子问题交替求解。

1.3 光谱近邻算子的作用位置

光谱近邻算子，简单说就是把相邻波长之间的差分作为正则项，通过近邻映射（proximal mapping）在每轮迭代里对变量做一次软阈值操作。它的作用位置在ADMM的辅助变量更新里，负责执行"让相邻波长看起来像同一个样本"这一约束。相比直接把约束写进目标函数，近邻算子的好处是计算代价低，对复数变量也成立，且不需要调整主迭代的收敛特性。

2. ADMM把问题拆成了哪三个子问题

2.1 目标函数和变量引入

我的高光谱宽带相位恢复目标函数写成如下形式：

[ \min_{X,Z,V} ; \sum_{l=1}^{L} \frac{1}{2}|Z_l - b_l|^2 + \beta |V|_1 ]

约束条件：

[ Z_l = \mathcal{F}(X_l), \quad V = D_{\text{spec}} X, \quad X_l \in \mathcal{C} ]

其中 (X) 是尺寸为 ([N_y, N_x, L]) 的复振幅场张量，(b_l = \sqrt{y_l}) 是第 (l) 个波长的实测幅度，(D_{\text{spec}}) 是作用在波长维上的差分算子（对相邻波长的复振幅做差），(\mathcal{C}) 是支持域约束（样本在视场内有限区域）。(\beta) 控制光谱平滑强度。

增广拉格朗日写出来就是：

[ \begin{aligned} \mathcal{L} = &\sum_l \frac{1}{2}|Z_l - b_l|^2 + \beta|V|1 \ &+ \frac{\rho_z}{2}|\mathcal{F}(X) - Z + U_z|^2 \ &+ \frac{\rho_v}{2}|D{\text{spec}}X - V + U_v|^2 \end{aligned} ]

这里 (U_z)、(U_v) 是缩放对偶变量。拆开看就是三个子问题：(X) 子问题、(Z) 子问题、(V) 子问题。下面逐个说。

2.2 X子问题：带支持域约束的最小二乘

固定 (Z)、(V) 和两个对偶变量，(X) 子问题变成：

[ \min_X ; \frac{\rho_z}{2}|\mathcal{F}(X) - Z + U_z|^2 + \frac{\rho_v}{2}|D_{\text{spec}}X - V + U_v|^2 ]

因为 (\mathcal{F}) 是傅里叶变换，(\mathcal{F}^H\mathcal{F} = I)，这个子问题可以写成如下正规方程：

[ \left(I + \gamma D_{\text{spec}}^T D_{\text{spec}}\right) X = \mathcal{F}^{-1}(Z - U_z) + \gamma D_{\text{spec}}^T (V - U_v) ]

其中 (\gamma = \rho_v / \rho_z)。(D_{\text{spec}}^T D_{\text{spec}}) 在光谱维上是一个三对角矩阵，系数矩阵是稀疏的，不需要直接求逆。我在Matlab里用预条件共轭梯度法（pcg）解，通常十几次迭代就收敛，因为每轮ADMM里 (X) 的变化本来就不大，以上一轮解作为pcg初值，实际计算量很低。

解完之后再乘上支持域掩膜 ( \text{support} )，把视场外的值清零。这一步是对样本有限尺寸的先验，也是从编码孔径或光瞳约束里继承来的。

2.3 Z与V子问题：傅里叶幅度投影与光谱软阈值

(Z) 子问题是数据保真项，目标函数：

[ \min_Z ; \frac{1}{2}|Z - b|^2 + \frac{\rho_z}{2}|\mathcal{F}(X) - Z + U_z|^2 ]

这是一个逐点二次函数，闭式解是：

[ Z = \frac{b + \rho_z |\mathcal{F}(X) + U_z|}{1 + \rho_z} \cdot e^{i \angle(\mathcal{F}(X) + U_z)} ]

也就是把傅里叶域当前估计的幅度朝实测幅度方向拉，但保留相位。这就是经典的傅里叶域幅度投影，不过多了一个 (\rho_z) 平衡项。(\rho_z) 越大，越信任当前迭代的傅里叶值；越小，越贴近实测幅度。

(V) 子问题则完全是光谱近邻算子的主场：

[ \min_V ; \beta|V|1 + \frac{\rho_v}{2}|D{\text{spec}}X - V + U_v|^2 ]

它的闭式解是软阈值：

[ V = \operatorname{soft}\left(D_{\text{spec}}X + U_v, \frac{\beta}{\rho_v}\right) ]

这一步就是文章标题里"光谱近邻算子"的实际样貌：在波长维上对相邻波长的差分做收缩。差分的幅度大于阈值的部分保留，小于阈值的部分直接清零，让相邻波长的解趋于一致，但又不会把真实的色散差异全部抹掉。

剩下的对偶更新也很机械：

Uz = Uz + (FX - Z); Uv = Uv + (DspecX - V);

整个ADMM循环就是：更新X，更新Z，更新V，更新对偶变量，重复直到残差收敛。从结构上看，数据保真和光谱约束被安排到了不同类型的子问题里，互不干扰，这是ADMM相比单层投影法最大的优势。

2.4 停止准则与残差监控

ADMM的停止准则用原始残差和对偶残差双指标。原始残差衡量约束满足程度：

[ r_{\text{prim}} = |\mathcal{F}(X) - Z|F + |D{\text{spec}}X - V|_F ]

对偶残差衡量最优性条件满足程度：

[ s_{\text{dual}} = \rho_z |\mathcal{F}(X^k - X^{k-1})|F + \rho_v |D{\text{spec}}(X^k - X^{k-1})|_F ]

实际跑的时候，我习惯每轮都画一下这两个残差的对数曲线。如果原始残差下降但幅度很小，多半是 (\rho) 选大了；对偶残差振荡大，多半是 (\rho) 选小了。这个经验比任何理论收敛条件都直观。

3. 光谱近邻算子的Matlab实现细节

3.1 近邻算子的数学定义与软阈值写法

严格来说，一个凸函数 (g) 的近邻算子定义为：

[ \operatorname{prox}_g(v) = \arg\min_u \frac{1}{2}|u - v|_2^2 + g(u) ]

当 (g(u) = \tau |u|_1) 时，近邻算子就是软阈值：

[ \operatorname{prox}_{\tau|\cdot|_1}(v) = \operatorname{sign}(v) \cdot \max(|v| - \tau, 0) ]

对复数变量，这里的 (\operatorname{sign}) 要改成相位项，也就是 (v / |v|)。这一点最容易出错，很多人直接套实数软阈值，把复数幅度信息弄丢了。

在ADMM的 (V) 子问题里，变量是一个 ([N_y, N_x, L-1]) 的差分张量，软阈值作用在每一个像素点、每一对相邻波长上。形式上就是在光谱维上"瘦身"，所以叫光谱近邻算子。

3.2 核心函数代码

下面给出我实际在Matlab里用的核心实现。第一是复数软阈值：

function V = soft_threshold(A, tau) amp = abs(A); scale = max(amp - tau, 0) ./ max(amp, eps); V = scale .* A; end

这段代码保留了复数的相位，只对幅度做收缩。( \text{eps} ) 是为了防止零幅度处除零。

第二是 (X) 子问题的求解。因为 (D_{\text{spec}}) 是光谱差分算子，(D_{\text{spec}}^T D_{\text{spec}}) 是三对角矩阵，我直接写了一个稀疏矩阵来配合pcg：

function X = update_X(Z, Uz, V, Uv, support, rho_z, rho_v) gamma = rho_v / rho_z; RHS = ifft2(Z - Uz) + gamma * adjoint_spectral_diff(V - Uv); A = @(x) x + gamma * spectral_diff_adjoint(spectral_diff_forward(x)); [X, ~] = pcg(A, RHS, 1e-6, 20, [], [], X_prev); X(~support) = 0; end

其中：

function d = spectral_diff_forward(X) d = diff(X, 1, 3); end function x = spectral_diff_adjoint(d) x = cat(3, -d(:, :, 1), -diff(d, 1, 3), d(:, :, end)); end

这两个函数一正一伴，配合使用。注意Matlab的diff(X, 1, 3)会丢掉最后一层，邻接算子要补回来边界项，否则正规方程不对称，pcg会直接发散。我第一次写的时候就忘了补边界，结果残差死活降不下去。

第三是 (Z) 子问题的更新：

function Z = update_Z(FX, Uz, b, rho_z) temp = FX + Uz; amp = (b + rho_z * abs(temp)) / (1 + rho_z); Z = amp .* exp(1i * angle(temp)); end

主循环：

for k = 1:opts.maxIter X_prev = X; X = update_X(Z, Uz, V, Uv, support, rho_z, rho_v); FX = fft2(X); Z = update_Z(FX, Uz, b, rho_z); Uz = Uz + (FX - Z); Dx = spectral_diff_forward(X); V = soft_threshold(Dx + Uv, beta / rho_v); Uv = Uv + (Dx - V); r_prim(k) = norm(FX - Z, 'fro') + norm(Dx - V, 'fro'); s_dual(k) = rho_z * norm(fft2(X - X_prev), 'fro') ... + rho_v * norm(spectral_diff_forward(X - X_prev), 'fro'); if r_prim(k) < opts.tol && s_dual(k) < opts.tol break; end end

这套代码跑下来的稳定性和速度都还不错，128x128像素、3个波长的数据，在普通台式机上300轮迭代大概十几秒。

3.3 为什么要用L1范数而不是L2

光谱维上的平滑约束如果换成L2范数，(V) 子问题的解就变成普通收缩，不做阈值截断，效果差异很大。L1正则允许少数相邻波长之间存在较大差异，不会把所有真实色散细节全都抹平；L2则倾向于把所有差异均匀摊薄，结果就是相位曲线的光谱细节被过度平滑。我用模拟数据对比过L1和L2在强色散样本上的表现，L1的相位RMSE比L2低约40%。所以光谱近邻算子里的阈值操作不是一个实现选择，而是一个建模选择。

3.4 复数域软阈值的坑

实数阈值 (\operatorname{sign}(v)\max(|v|-\tau,0)) 里的符号在复数域要替换成归一化相位。如果直接对实部和虚部分别做软阈值，会把幅度和相位耦合到一起，导致恢复结果出现明显的方格状伪影。这个坑我在一开始踩过，后来把所有处理都改成"幅度软阈值、相位保留"之后，伪影立刻消失。另外阈值 (\tau = \beta/\rho_v) 的单位是"复数场的幅度差"，如果你的数据没有做归一化，阈值量纲对不上，也会出现约束过强或过弱的情况。我的习惯是把各波长的衍射强度先归一化到总能量一致，再调 (\beta)。

4. 定量相位成像的仿真验证与物理量换算

4.1 仿真实验设置

为了验证这套方法在定量相位成像里的表现，我做了三层模拟实验。成像系统设定为透射式衍射成像，探测器距离样本约100个波长距离，像素数128x128，波长取488nm、561nm、640nm三个通道，模拟高光谱宽带照明的离散采样。

样本用的是模拟的双高斯相位球，等效厚度约2微米，折射率差0.05，模拟活细胞的量级。各波长衍射强度加入泊松噪声和5%高斯噪声混合，信噪比约20dB。算法参数：(\rho_z=0.1)，(\rho_v=0.05)，(\beta=0.01)，最大迭代300轮，支持域取一个直径80像素的圆。

4.2 重建质量对比

我把四种方法跑了同样的数据：经典GS、HIO、不带光谱约束的ADMM、带光谱近邻算子的ADMM。评价指标用峰值信噪比、结构相似性和相位均方根误差。

带光谱近邻算子的ADMM在相位误差上几乎是HIO的1/3，这也在意料之中——因为多波长数据里共享的结构信息被显式利用了。值得一提的还有收敛速度：GS和HIO在500轮时已经基本不下降，ADMM在300轮内就达到更高精度，而且没有出现GS常见的振荡。

从恢复图像上看，最明显的区别在样本边缘。GS和HIO恢复的边缘容易出现高频伪影环绕，而ADMM加光谱约束后边缘干净很多，这应该归功于光谱维上的差分正则承担了一部分高频噪声的抑制。

4.3 从相位恢复结果解耦厚度和色散

定量相位成像最终要回答的问题是：样本的厚度是多少、折射率分布如何。单波长相位图只能给出光程差，没法同时解出厚度和折射率；多波长的优势就在这里。

每个波长恢复出的相位满足：

[ \varphi_l(x,y) = \frac{2\pi}{\lambda_l} \left(n(\lambda_l) - n_m\right) d(x,y) ]

其中 (n_m) 是介质折射率，(d) 是样本厚度。把折射率色散用Cauchy形式近似：

[ n(\lambda) = A + \frac{B}{\lambda^2} + \frac{C}{\lambda^4} ]

那么对每个像素，把三个波长的相位带入，整理成一个线性方程组：

[ \frac{\varphi_l \lambda_l}{2\pi} = d \cdot \left(A - n_m + \frac{B}{\lambda_l^2} + \frac{C}{\lambda_l^4}\right) ]

三个波长三个未知数（(A-n_m)、(B)、(C)），加上 (d) 的耦合实际上是一个变量分离的拟合问题——因为 (d) 乘在括号外，需要联立求解。我实际是用一个小的最小二乘迭代，先固定色散系数求厚度，再固定厚度求色散，两轮就收敛。三个波长刚好够用，四个波长会更稳。

这一步里相位解包裹必须先做。ADMM恢复出来的是缠绕相位，直接带入公式会得到跳变的厚度图。我用的质量图引导路径跟随法，以各波长幅度投影的置信度为权重，质量高的像素先解，最后处理低信噪比区域。

5. 调参实战与踩坑记录

5.1 初始化决定成败

ADMM虽然比GS稳，但相位恢复本质上还是非凸问题，初始化不好一样会掉进坏局部最小。我的经验是三步走：

1. 对每个波长单独用支持域约束的Gerchberg-Saxton跑30轮，得到一个中等质量的初始估计。
2. 对初始估计做相邻波长平均，抹掉一部分独立恢复带来的光谱噪声。
3. 把这个平均值作为ADMM的X起点，Z起点设为它的傅里叶变换，对偶变量全部置零。

这个流程比随机初始化稳定得多。我试过纯随机相位初始化，ADMM大概有30%的几率收敛到伪影严重的解；用GS预热后，失败率降到5%以下。代价只是GS的30轮迭代，很划算。

还有一种更省事的初始化：直接用各波长平均强度的平方根乘以随机相位。这个方案在支持域约束很紧时也能用，但如果支持域不够紧，建议还是走GS预热。

5.2 rho和beta的调整策略

(\rho) 的取值直接影响收敛速度。我最初固定 (\rho_z=1)，结果原始残差降得很慢。后来改成动态调整策略：每50轮比较原始残差和对偶残差，如果原始残差偏大就增大 (\rho)，对偶残差偏大就减小 (\rho)。这里有个经验公式：

[ \rho \leftarrow \rho \cdot \min\left(2, \max\left(0.5, \frac{|r_{\text{prim}}|}{|s_{\text{dual}}|}\right)\right) ]

实际运行中，这个自适应策略能把迭代轮数减少40%左右。(\beta) 的选择也很有讲究。(\beta) 太小，光谱近邻算子基本不起作用，恢复结果接近独立ADMM；(\beta) 太大，相邻波长被强行拉成一样，真实的色散信息被破坏。我的标定方法是：先用无光谱约束的ADMM跑一遍，统计相邻波长恢复幅度差的平均绝对值，把这个值的5%到15%作为 (\beta) 的合理区间。这样标出来的 (\beta) 通常很稳。

5.3 相位解包裹与低信噪比区域处理

定量相位成像里，恢复出的相位分布超过 (2\pi) 就必然遇到解包裹。ADMM本身的输出并不会自动解决这个问题，它只负责给出最可能的缠绕相位。解包裹时最头疼的是低信噪比区域——样本边缘和视场外背景噪声会让质量图迅速恶化。我最后总结的流程是：

• 先对幅值图做一个简单的分割，把背景像素标记为不可信。
• 解包裹时只对样本区域内做路径跟随，背景区域用样条插值填充。
• 解完包裹再做一次中值滤波，去掉孤立的 (2\pi) 跳变点。

这三个步骤做完，厚度图基本没有明显的解包裹伪影。如果还有孤立坏点，多半是初始相位恢复时相位跳变本身搞错了，需要回到ADMM参数上找原因，而不是在解包裹阶段硬修。

5.4 低信噪比波长的权重处理

高光谱宽带数据里经常有一个波长特别弱。比如640nm在大多数探测器上量子效率偏低，衍射图案噪声很大。如果所有波长在目标函数里权重一样，弱波长就会拖累整体光谱平滑约束。我的做法是在数据保真项里按波长加权重系数 (w_l)，正比于该波长实测强度的对数均值。弱波长的权重降下来之后，它主要依靠光谱近邻算子从相邻强波长那里获得信息，而不是用自己的噪声强行主导重建。

这一点在实际成像里比仿真更容易被忽视。仿真里所有波长信噪比接近时，不加权重没什么感觉；一旦拿到真实系统数据，弱波长通道的高频噪声几乎可以让整个光谱维约束失效。所以如果你在真实系统上做，建议一开始就把权重项写进代码，省得后面回过头改。

整套方法跑通之后，我心里最深的感受是：高光谱宽带相位恢复真正困难的地方，不在于某个波长的相位恢复本身，而在于如何让多个波长的解既保持各自的物理特性，又共享合理的光谱结构。ADMM加光谱近邻算子恰好提供了一个干净利落的框架，把这两种诉求拆到不同子问题里交替解决。最后再分享一个小习惯：每次重建完，我都会把相邻波长的相位差画出来和理论Cauchy色散曲线叠加对比一下。如果相位差曲线和色散曲线趋势一致，说明光谱约束没有越界；如果出现明显背离，九成是 (\beta) 调过头了，回去把阈值放宽一个量级再跑一轮，通常就对了。