做定量相位成像这几年最让我头疼的不是光学平台而是重建算法。单波长下跑跑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_{l1}^{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}}) 是三对角矩阵我直接写了一个稀疏矩阵来配合pcgfunction 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_z0.1)(\rho_v0.05)(\beta0.01)最大迭代300轮支持域取一个直径80像素的圆。4.2 重建质量对比我把四种方法跑了同样的数据经典GS、HIO、不带光谱约束的ADMM、带光谱近邻算子的ADMM。评价指标用峰值信噪比、结构相似性和相位均方根误差。算法迭代数PSNR(dB)SSIM相位RMSE(rad)GS50024.60.840.42HIO50027.30.910.28ADMM无光谱约束30028.90.930.23ADMM光谱近邻算子30032.10.970.12带光谱近邻算子的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稳但相位恢复本质上还是非凸问题初始化不好一样会掉进坏局部最小。我的经验是三步走对每个波长单独用支持域约束的Gerchberg-Saxton跑30轮得到一个中等质量的初始估计。对初始估计做相邻波长平均抹掉一部分独立恢复带来的光谱噪声。把这个平均值作为ADMM的X起点Z起点设为它的傅里叶变换对偶变量全部置零。这个流程比随机初始化稳定得多。我试过纯随机相位初始化ADMM大概有30%的几率收敛到伪影严重的解用GS预热后失败率降到5%以下。代价只是GS的30轮迭代很划算。还有一种更省事的初始化直接用各波长平均强度的平方根乘以随机相位。这个方案在支持域约束很紧时也能用但如果支持域不够紧建议还是走GS预热。5.2 rho和beta的调整策略(\rho) 的取值直接影响收敛速度。我最初固定 (\rho_z1)结果原始残差降得很慢。后来改成动态调整策略每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) 调过头了回去把阈值放宽一个量级再跑一轮通常就对了。