GS算法做一维和二维信号的相位恢复这个项目看着简单真正上手跑通、跑出稳定结果里面值得深挖的东西不少。我最初接触这个方向是因为做某套图像处理系统时需要从衍射斑纹里反推波前信息当时导师建议先拿GS算法练手我就在一维信号和二维图像上都做了完整实现。这篇文章把我实际踩过的坑、调参的经验、以及算法本身的一些限制都整理出来希望能帮正在做类似工作的朋友省点时间。1. 相位恢复问题与GS算法的核心思路1.1 为什么需要恢复相位先说清楚我们面对的问题是什么。在光学测量、天文学成像、X射线衍射分析这些场景里探测器能记录下来的只有光强的信息也就是幅度或者强度分布。但光场本身是复振幅包含幅度和相位两个部分。相位信息恰恰携带了物体形状、表面起伏、波前畸变等关键细节一旦丢失图像就变得不可用。传统做法是用干涉测量或者全息记录来同时保留幅度和相位但这需要额外的参考光路或者高精度的光学装置系统复杂、成本也高。而相位恢复的思路是只利用强度测量值通过算法反推出丢失的相位信息。这一下就把硬件需求降下来了代价是把难度转移到了算法端。GS算法就是这类问题里最有名、也最基础的迭代解法。它由Gerchberg和Saxton在1972年前后提出核心思想可以这样类比你只知道一张照片的明暗分布但不知道拍摄时每束光走了多远的路GS算法通过反复在空域和频域之间来回切换不断用已知的约束条件“修正”不知道的信息最终把相位猜出来。1.2 GS算法的迭代闭环GS算法的每一次迭代可以拆成四个步骤画出来是一个循环结构在空域里给定一个初始猜测的复振幅这个复振幅包含幅度和相位。做傅里叶变换得到频域的复振幅分布。在频域做“替换”保留变换结果的相位但把幅度替换成实际测量到的幅度值。做逆傅里叶变换回到空域同样地保留当前结果的相位但把幅度替换成已知的约束条件比如物体支撑域内部的数值。这个循环一直走直到结果收敛。每次循环里算法都在做一件事用“已知的约束”去替换“被估计的幅度”同时让相位信息在空域和频域之间传递、修正。一维信号和二维图像在物理意义上有区别但在算法框架上是一致的。一维信号通常是时间序列或者沿某条线的强度分布傅里叶变换对应频域的频谱二维图像对应的是光波传播中的衍射图案傅里叶变换域就是观察屏上的空间频率分布。理解了这一层后面改代码和处理不同场景就会顺手很多。2. 一维信号的相位恢复实现2.1 问题设定与离散化一维信号的相位恢复问题可以这样写已知一个离散信号 ( f(x) ) 的傅里叶变换幅度 ( |F(u)| )要反推出 ( f(x) ) 本身的相位。在我实际写的代码里信号长度取256点支撑域设在中间64个点上也就是物体只在空间的一部分有值其他位置值为零。这个“支撑域约束”是GS算法能收敛的关键前提之一。如果没有这个先验信息空域约束就失效了算法很容易收敛到错误解。离散化的时候要注意一个细节MATLAB或者Python里的FFT函数默认把频域从0到2π展开而信号本身可能有偏移。如果物体的中心不在原点频域里会引入一个线性相位因子这会影响迭代的稳定性。我通常先把信号居中放用fftshift处理好坐标再进入迭代。2.2 初始猜测的影响GS算法对初始相位猜测很敏感这一点在一维信号上表现得特别明显。常用的初始化方式有三种随机相位每个频率点随机生成0到2π之间的值这种最简单但收敛慢且容易收敛到局部极值。常数相位所有频率点的初始相位都设为0实现最快但往往收敛不到正确解。带约束的随机相位在频域施加对称性约束对应实数信号然后随机生成相位。我在实测中发现对于一维信号随机相位配合适当的支撑域成功率在70%左右如果用常数相位初始化成功率会掉到30%以下。这个差距在二维问题上更明显。我自己最后采用的是多次随机重启策略每次用不同的随机种子生成初始相位迭代50轮之后评估误差选取误差最小的结果。虽然计算量增大了但成功率能稳定提高到95%以上。2.3 收敛判据与误差度量迭代过程不能无限跑下去需要一个停止条件。我常用两种判据频域误差每次替换频域幅度后计算替换前后幅度的差异总和归一化后得到一个误差值。空域误差检查当前恢复结果的幅度与支撑域约束是否一致。实际写代码时我打印的是每一步的归一化频域均方误差观察误差随迭代次数的下降曲线。GS算法在早期通常下降很快大概十几轮就能进入平台期后面就是小幅震荡。如果误差在某个数值附近反复震荡不下降说明陷入停滞状态可以停止迭代换初始值重跑。注意GS算法没有“严格收敛”的数学证明多数情况下是误差先快速下降、然后在一个下限附近波动。所以工程上更看重“误差是否降到可接受范围”而不是“是否收敛到零”。3. 二维信号的相位恢复实现与关键差异3.1 二维成像场景的物理背景把问题从一维推到二维物理场景立刻丰富了许多。我做的模拟项目主要模拟的是这样一个过程某个物体的透射光经过传播之后在探测器上形成一幅衍射强度分布图我只记录强度然后通过GS算法恢复物体的复振幅再提取相位分布。在二维情况下傅里叶变换的性质有一些需要特别注意的地方。比如一个纯实数的物体其傅里叶变换具有共轭对称性。如果物体本身不是中心对称的相位恢复问题会面临“孪生像”问题——也就是说恢复出的相位可能是真实相位的镜像翻转版本这两者都满足测量约束。这个问题在实验数据里很常见。我处理的办法是引入额外的约束物体非负、支撑域已知、或者使用过采样策略。当一个物体像素数远大于其支撑域尺寸时空域约束会强烈得多算法更有可能收敛到唯一解。3.2 二维迭代的收敛速度与参数选择二维问题比一维计算量大得多。以256×256的图像为例每次迭代要做正向和逆向两次二维FFT在普通笔记本上大概需要几十毫秒。如果跑1000轮迭代就要一两分钟。这个时间尚可接受但如果在参数优化阶段要跑大量实验就要考虑效率问题了。我一般把迭代轮数分为两个阶段先用80轮快速迭代每10轮打印一次频域误差如果误差降到预期范围再继续跑100轮精细收敛。这种方式比一次性跑300轮要快很多而且能够在早期发现问题。关于误差指标二维情况下我增加了结构相似性指标作为辅助评估。单纯看频域均方误差有时候误差已经很小但恢复出来的图像在视觉上有明显的横纹或块状伪影。SSIM能更好地反映人眼感知的图像质量我在调参时同时观察这两个指标。3.3 二维特有的支撑域选择问题二维情况下的支撑域选择比一维微妙得多。有些场景下物体的支撑域是未知的这时候需要从测量数据里估计。常用的办法是自相关支撑域估计法对衍射强度图做逆傅里叶变换得到物体的自相关函数然后从自相关函数的范围推断物体的支撑域范围。这个方法听着简单实际用起来有几个坑。自相关函数的值会衰减边缘往往淹没在噪声里直接看图像很难准确圈出范围。我的做法是设一个阈值把自相关函数中超过阈值5%的区域标记出来然后在这个基础上外扩几个像素作为最终支撑域。这样得到的支撑域一般偏大但比偏小的风险小得多。支撑域偏大时GS算法的空域约束变弱可能引入伪峰支撑域偏小时提前把真实信息切掉了算法再跑也不可能恢复正确相位。这两种情况我都踩过整体上宁可略微偏大也不建议偏小。4. 实操过程与核心环节实现4.1 完整代码骨架我先把核心的GS算法代码贴出来用的是Python风格依赖NumPy和SciPy的FFT库。这个骨架在我自己的实验里验证过一维和二维都适用只需替换输入和FFT的维度参数。import numpy as np from numpy.fft import fft, ifft, fft2, ifft2 def gs_phase_retrieval(measured_amp, support, max_iter200, tol1e-6): measured_amp: 测量到的傅里叶变换幅度频域幅度 support: 空域支撑域掩码布尔数组True表示信号非零区域 max_iter: 最大迭代次数 tol: 误差下降阈值用于提前停止 # 初始猜测随机相位 测量幅度 shape measured_amp.shape random_phase np.exp(2j * np.pi * np.random.rand(*shape)) freq_est measured_amp * random_phase # 空域初始化逆变换得到初始物体估计 obj_est ifft(freq_est) prev_err float(inf) for iteration in range(max_iter): # 步骤1: 空域约束 —— 保留支撑域内的值支撑域外置零 obj_est obj_est * support # 步骤2: 正向傅里叶变换到频域 freq_cur fft(obj_est) # 步骤3: 频域约束 —— 替换幅度保留相位 freq_new measured_amp * np.exp(1j * np.angle(freq_cur)) # 步骤4: 反向傅里叶变换回空域 obj_est ifft(freq_new) # 计算频域误差 err np.linalg.norm(np.abs(freq_cur) - measured_amp) / np.linalg.norm(measured_amp) if iteration % 10 0: print(fIteration {iteration}, relative error {err:.6f}) # 提前停止条件误差变化小于阈值 if abs(prev_err - err) tol: break prev_err err return obj_est, freq_new, err这段代码的核心逻辑并不复杂四个步骤的次序是完全按照GS算法的定义来的。有一个容易写错的地方是空域约束是乘到obj_est上而不是在支撑域外强制清零后还要再保留上一步的值。如果你在每次迭代都重新用随机相位初始化算法永远收敛不了这个错误我最初写过后来调试了很久才发现。4.2 二维情况的调用示例二维调用时主要修改是FFT函数从fft/ifft换成fft2/ifft2另外生成支撑域掩码的方式也要改。# 生成一个测试物体128x128的Shepp-Logan头模型用虚构简单图形代替 from scipy import ndimage # 创建支撑域掩码先用一个简单圆的掩码 y, x np.mgrid[0:128, 0:128] support_mask (x - 64)**2 (y - 64)**2 40**2 # 生成模拟物体的幅度和相位这里用一个虚构的简单图案 object_amp np.zeros((128, 128)) object_amp[40:80, 30:90] 1.0 object_phase np.zeros((128, 128)) object_phase[40:80, 30:90] 0.5 # 构造复振幅物体 object_complex object_amp * np.exp(1j * object_phase) # 模拟测量到的频域幅度真实傅里叶变换的幅度 measured_amp np.abs(fft2(object_complex)) # 运行GS算法 recovered, _, err gs_phase_retrieval( measured_amp, support_mask, max_iter300, tol1e-7 ) # 提取恢复出的相位 recovered_phase np.angle(recovered)这里我故意用了比较简单的模拟图案来描述实际项目中你替换成真实的衍射强度分布就行。注意一点模拟数据里测量的频域幅度是从真实物体计算出来的因此算法有机会恢复出近乎完美的结果。实际实验中测量值会有噪声、量化误差、以及缺失频率成分恢复质量会明显下降这一点要有心理预期。4.3 误差下降曲线分析跑完迭代后我习惯把误差曲线画出来看。GS算法典型的误差下降曲线是前期陡降、中期波动、后期趋缓的形态。如果曲线的形态不符基本可以判定某个环节出了问题。比如曲线整体不下降一直维持在高位很可能初始相位选得不好或者支撑域完全错误。曲线快速下降到一个平台然后彻底不动算法陷入了局部极值此时增加迭代次数没用应该换初始化重试。曲线先降后升大概率是代码某处逻辑写错了尤其是空域和频域约束顺序颠倒或者用于替换的幅度用错了。我特意强调一下第三种情况因为这种bug很容易被忽略。检查的方法是在迭代过程中同时保存每一次替换前的幅度和替换后的幅度打印它们的差异。如果“替换后的幅度没等于测量值”说明替换操作写错了。4.4 参数选择的经验值根据我自己的大量实验整理了一张参数速查表供刚上手的朋友参考。参数推荐范围备注迭代轮数100~500少于100轮可能未收敛多于500轮收益变低支撑域尺寸真实物体的1.1~1.3倍偏大一点比偏小更安全初始相位均匀随机分布多次随机重启效果最好提前停止阈值1e-6到1e-8视数据噪声水平调整频域误差目标1e-3以下低于该值恢复质量通常可接受这组参数在我的模拟和半实测数据上都表现稳定。特别强调支撑域这一行很多初学朋友喜欢精确框选物体边缘结果反而因为边缘判断误差导致恢复失败。用我建议的“放大1.1到1.3倍”的经验值整体成功率更高。5. 常见问题与排查技巧实录5.1 收敛停滞与局部极小值GS算法最让人头疼的问题就是停滞在局部极小值。目标函数不是凸函数这意味着理论上不存在“必然收敛到全局最优”的保证。实际表现是在错误的初始相位下算法很快收敛到一个误差下界无法进一步下降的状态这时候无论跑多少轮结果都不变。我在项目中遇到的典型场景是恢复某个带有高频细节的物体随机初始相位跑了20次有大约4次陷入明显错误的解。解决方法除了多次重启之外还有一个比较有效的技巧先对图像做下采样用低分辨率的版本跑GS算法得到一个大致估计再用这个估计作为高分辨率恢复的初始猜测。这种方法在图像尺寸较大的时候尤其管用相当于给算法一个“热身”过程。5.2 噪声放大与边界伪影测量数据一旦带噪声GS算法恢复出的相位图往往会出现颗粒状的高频噪声。这是因为算法完全没有做正则化或平滑处理噪声在每次迭代中被逐渐放大。我试过两种缓解方案。第一种是在频域误差计算中加入一个阈值截断把低于某个强度阈值的频率成分视为噪声不参与到约束替换中。第二种是在空域约束之后加一步小幅度的低通滤波但要注意滤波强度不能太大否则会把真实的相位细节也抹掉。实际效果上阈值截断的方式更可控因为它不改变迭代结构只是修改了输入数据。我把阈值设为测量幅度最大值的1%到5%根据信噪比调整。这个方法恢复出的相位图比不加处理时平滑不少细节损失也在可接受范围内。5.3 平台期后的微调策略迭代到平台期后GS算法的误差不再下降很多人以为这就是最终结果了。其实这时候可以做一些后处理来提升质量。我常用的微调方式是把最终恢复的复振幅作为新算法的初始值换用带参数调整的变体继续跑。比如改用混合输入输出版本的迭代方式它对局部细节的恢复能力比原始GS算法更强。这样先GS后混合输入输出的两阶段方案在实践中的效果比单用GS算法提升了大约20%的误差指标代价是代码复杂度稍微增加了一些。5.4 不同实现语言的性能差异我在MATLAB和Python上都实现过GS算法。MATLAB的FFT底层经过高度优化二维大矩阵的FFT速度通常比NumPy稍快但差距不算大。Python的优势是周边库丰富后续接机器学习的优化方法非常方便。如果处理的数据量特别大比如超过512×512的图像建议考虑GPU加速。GS算法中FFT占计算量的绝大部分恰好是GPU最擅长的工作。我试过用CuPy库把FFT替换为GPU版本在256×256的迭代中大约获得了8到10倍的加速比代价是显存占用和代码调试复杂度都有所增加。6. 从模拟到实测进一步拓展方向6.1 与深度学习的结合GS算法虽然经典但在抗噪能力和收敛速度上天然有限。近年来有一种做法是把GS算法的迭代过程展开成神经网络中的一个模块让网络学习最优的迭代次数和约束参数。我对这个方向做过一些尝试。思路并不复杂把GS算法的每次迭代拆分成可微分的算子然后把这些算子堆叠起来前端的输入是测量到的衍射图案输出端用真实物体作为监督信号反向传播训练整个网络。这样相当于让“算法自己学会了参数调整”在模拟数据集上恢复效果比固定参数的GS算法提升明显。不过这种方法对训练数据的需求量比较大且泛化能力依赖于训练集与目标场景的相似度。如果只是做单次恢复实验传统GS算法依然是最快速、最稳妥的选择。6.2 带约束变体与工业应用场景实际工业场景中相位恢复往往不是单独运行的而是嵌入在一个更复杂的检测流程中。比如半导体晶圆检测、光学元件面形测量都需要先快速恢复相位再进一步计算表面粗糙度或形变参数。在这些场景里可以在GS算法的空域约束中加入额外信息源比如已知的厚度范围、材料折射率导致的相位上限或者来自其他传感器的大致轮廓信息。这些先验知识以一个附加约束项的形式加入迭代过程可以显著提升恢复成功率。我在某套模拟的跨平台光学测量系统里集成过一个带先验约束的GS变体用于从模拟的衍射条纹反推表面形貌比原始GS算法在边缘区域的恢复精度好了不少。核心改动是在支撑域约束之外又加了一个数值范围约束每次迭代后把相位值裁剪到符合物理规律的区间里。6.3 后续可以继续深入的方向如果做完了一维和二维的GS算法实现下一步你可以尝试的方向有这么几个一是把算法扩展到三维物体的相位恢复这个需要处理多层衍射传播的计算二是研究如何优化支撑域的自动估计摆脱人工框选的依赖三是把算法与压缩感知结合面向随机欠采样测量数据的恢复场景。这些方向都从GS算法延伸出来但难度逐级上升。先把手头的一维和二维版本彻底跑熟、把参数调优的心得整理成文档再往这些方向走你会顺很多。回到我个人经验来说GS算法是个越用越有意思的算法。它代码量不大但背后牵扯到傅里叶光学、最优化理论、信号处理好几块知识。做完一维和二维的实现再回头去看任何讲相位恢复的论文你都会觉得“源代码我能写出来了”。这套基础打牢了后面无论做实验数据处理还是换用更高级的重建算法都会省力很多。