相位恢复（Phase Retrieval）方法详解

相位恢复（Phase Retrieval）方法详解

相位恢复是信号处理中的一个基础问题，它涉及从信号或波的振幅信息中重建其相位信息。在许多实际应用中，我们只能测量信号的强度或振幅，而相位信息却无法直接获取。然而，相位信息往往包含了关于信号的重要特征，因此如何从有限的振幅信息中恢复相位成为了一个关键挑战。

问题背景与数学表述

在物理学和工程学中，信号通常可以用复数形式表示：

$$f(x)=|f(x)|e^{i\varphi (x)}$$

其中$|f(x)|$是振幅，$\varphi (x)$是相位。在很多情况下，我们只能测量到强度$|f(x){|}^2$，而失去了相位信息$\varphi (x)$。相位恢复问题就是要从已知的强度测量值重建相位信息。在傅里叶域中，这个问题可以表述为：给定信号的傅里叶振幅$|\hat{f}(\omega )|$，求解原始信号$f(x)$。数学上，这涉及到求解以下方程：

$$|\mathcal{F}[f(x)]{|}^2=|\hat{f}(\omega ){|}^2$$

其中$\mathcal{F}[\cdot ]$表示傅里叶变换操作。

在离散域中，我们可以用线性算子$\mathbf{A}\in {\mathbb{C}}^{m\times n}$表示测量过程，则相位恢复问题可以表述为：

$$\mathbf{y}=|\mathbf{Ax}{|}^2+\mathbf{n}$$

其中$\mathbf{y}\in {\mathbb{R}}^m$是测量值，$\mathbf{x}\in {\mathbb{C}}^n$是待恢复的信号，$\mathbf{n}$是测量噪声。

数学基础：唯一性与不适定性

相位恢复问题的一个核心挑战是解的唯一性。对于一维信号，如果没有额外约束，相位恢复通常存在多解。这可以通过零点分析来理解：若$\hat{f}(\omega )$是$f(x)$的傅里叶变换，则$\hat{f}(\omega )$的零点与${\hat{f}}^{\ast }(-\omega )$的零点的镜像结合可以产生无穷多个有效解。

定理1（一维信号的双胞胎歧义性）：设$f(x)$和$g(x)$是具有有限支撑的一维函数，则$|f(x)|=|g(x)|$且$|\hat{f}(\omega )|=|\hat{g}(\omega )|$当且仅当以下之一成立：

1. $g(x)=e^{i\alpha }f(x)$，其中$\alpha$是常数；
2. $g(x)=e^{i\alpha }{f}^{\ast }(-x)$，其中$\alpha$是常数；
3. $g(x)=e^{i\alpha }f(x-a)$，其中$\alpha$和$a$是常数；
4. $g(x)=e^{i\alpha }{f}^{\ast }(-x+a)$，其中$\alpha$和$a$是常数。

然而，在二维或更高维空间中，相位恢复问题的解在几乎所有情况下都是唯一的（忽略平凡的模糊性如全局相位和平移）。这是因为高维空间中的零点集合是余维数为2的，这使得零点集合的随机组合产生有效解的概率趋近于零。

定理2（二维信号的唯一性）：对于几乎所有的具有有限支撑的二维函数$f(x,y)$，如果$g(x,y)$满足$|f(x,y)|=|g(x,y)|$且$|\hat{f}(u,v)|=|\hat{g}(u,v)|$，则$g(x,y)=e^{i\alpha }f(x\pm a,y\pm b)$，其中$\alpha$，$a$和$b$是常数。

高级变分方法：全变分正则化

在处理含噪声的相位恢复问题时，全变分（Total Variation，TV）正则化是一种有效的方法，特别是当原始信号具有分段平滑特性时。TV正则化的相位恢复问题可以表述为：

$$\underset{\mathbf{x}}{\min }\Vert \mathbf{y}-|\mathbf{Ax}{|}^2{\Vert }_2^2+\lambda \text{TV}(\mathbf{x})$$

其中，对于二维信号$\mathbf{X}\in {\mathbb{C}}^{n_1\times n_2}$（将$\mathbf{x}$重塑为矩阵），全变分定义为：

$$\text{TV}(\mathbf{X})=\sum _{i=1}^{n_1-1}\sum _{j=1}^{n_2-1}\sqrt{|{\mathbf{X}}_{i+1,j}-{\mathbf{X}}_{i,j}{|}^2+|{\mathbf{X}}_{i,j+1}-{\mathbf{X}}_{i,j}{|}^2}$$

对于复值信号，TV正则化可以分别应用于实部和虚部，或者直接应用于幅度：

$${\text{TV}}_{\text{complex}}(\mathbf{X})=\text{TV}(\text{Re}(\mathbf{X}))+\text{TV}(\text{Im}(\mathbf{X}))$$

复杂度与信息理论界限

从信息理论角度分析，成功恢复$n$维复值信号需要至少$4n-o(n)$个实值测量，这是由信号的自由度决定的。实际上，对于随机高斯测量矩阵$\mathbf{A}$，理论证明当测量数$m\ge 4n-o(n)$时，相位恢复问题在高维空间中几乎确定有唯一解。更精确地说，当测量算子满足受限等距性质（Restricted Isometry Property, RIP）时，我们有以下定理：

定理3（相位恢复的采样复杂度）：设$\mathbf{A}\in {\mathbb{C}}^{m\times n}$满足$(2k,\delta )$-RIP，即对任意$2k$-稀疏向量$\mathbf{z}$，有：

$$(1-\delta )\Vert \mathbf{z}{\Vert }_2^2\le \Vert \mathbf{Az}{\Vert }_2^2\le (1+\delta )\Vert \mathbf{z}{\Vert }_2^2$$

则对于任意$k$-稀疏信号${\mathbf{x}}_0$，如果测量值$\mathbf{y}=|\mathbf{A}{\mathbf{x}}_0{|}^2$，则${\mathbf{x}}_0$可以通过求解以下优化问题精确恢复（忽略全局相位）：

$$\underset{\mathbf{x}}{\min }\Vert \mathbf{x}{\Vert }_0\text{ subject to }|\mathbf{Ax}{|}^2=\mathbf{y}$$

当$m\ge Ck\log (n/k)$时，随机高斯矩阵$\mathbf{A}$满足上述RIP条件的概率趋近于1，其中$C$是一个常数。

高级算法：交替方向乘子法（ADMM）

交替方向乘子法（ADMM）是解决复杂约束优化问题的强大工具，适用于相位恢复。考虑以下带约束的形式：

$$\underset{\mathbf{x},\mathbf{z}}{\min }f(\mathbf{x})+g(\mathbf{z})\text{ subject to }\mathbf{x}=\mathbf{z}$$

其中$f(\mathbf{x})=\Vert \mathbf{y}-|\mathbf{Ax}{|}^2{\Vert }_2^2$，$g(\mathbf{z})$是正则化项（如${\ell }_1$范数或TV范数）。

ADMM迭代步骤为：

$${\mathbf{x}}_{k+1}=\arg \underset{\mathbf{x}}{\min }\left\{f(\mathbf{x})+\frac{\rho }{2}\Vert \mathbf{x}-{\mathbf{z}}_k+{\mathbf{u}}_k{\Vert }_2^2\right\}$$

$${\mathbf{z}}_{k+1}=\arg \underset{\mathbf{z}}{\min }\left\{g(\mathbf{z})+\frac{\rho }{2}\Vert {\mathbf{x}}_{k+1}-\mathbf{z}+{\mathbf{u}}_k{\Vert }_2^2\right\}$$

$${\mathbf{u}}_{k+1}={\mathbf{u}}_k+{\mathbf{x}}_{k+1}-{\mathbf{z}}_{k+1}$$

其中$\mathbf{u}$是标度拉格朗日乘子，$\rho >0$是惩罚参数。尽管$f(\mathbf{x})$是非凸的，但ADMM在实践中仍表现良好。对于第一个子问题，可以使用Wirtinger梯度下降；对于第二个子问题，若$g(\mathbf{z})=\lambda \Vert \mathbf{z}{\Vert }_1$，则有闭式解：

$${\mathbf{z}}_{k+1}={\text{Shrink}}_{\lambda /\rho }({\mathbf{x}}_{k+1}+{\mathbf{u}}_k)$$

其中${\text{Shrink}}_{\tau }(\mathbf{v})=\text{sign}(\mathbf{v})\max (|\mathbf{v}|-\tau ,0)$是软阈值算子。

高维非均匀傅里叶变换与多测量相位恢复

在多测量相位恢复问题中，我们有一组相关的测量值：

$${\mathbf{y}}_j=|{\mathbf{A}}_j\mathbf{x}{|}^2,j=1,2,\dots ,J$$

这种设置在多角度X射线成像等应用中很常见。若${\mathbf{A}}_j$表示不同角度的非均匀傅里叶变换（NUFT）算子，则可以写成：

$$[{\mathbf{A}}_j\mathbf{x}{]}_m=\sum _{n=0}^{N-1}{x}_n{e}^{-i{\omega }_{j,m}n},m=0,1,\dots ,M-1$$

其中${\omega }_{j,m}$是第$j$个测量中第$m$个频率点。在此情况下，可以利用测量之间的相关性设计更高效的算法。例如，通过结合各个测量的梯度信息：

$$\nabla f(\mathbf{x})=\sum _{j=1}^J\nabla f_j(\mathbf{x})$$

其中$f_j(\mathbf{x})=\Vert {\mathbf{y}}_j-|{\mathbf{A}}_j\mathbf{x}{|}^2{\Vert }_2^2$，每个梯度的Wirtinger形式为：

$$\nabla f_j(\mathbf{x})=2{\mathbf{A}}_j^H\left((|{\mathbf{A}}_j\mathbf{x}{|}^2-{\mathbf{y}}_j)\odot {\mathbf{A}}_j\mathbf{x}\right)$$

其中$\odot$表示逐元素乘法，${\mathbf{A}}_j^H$是${\mathbf{A}}_j$的共轭转置。

微分相位恢复与传输矩阵方法

在光学中，传输矩阵方法将输入场与输出场之间的关系表示为线性映射：

$${\mathbf{E}}_{\text{out}}=\mathbf{T}{\mathbf{E}}_{\text{in}}$$

其中$\mathbf{T}\in {\mathbb{C}}^{m\times n}$是传输矩阵，${\mathbf{E}}_{\text{in}}$和${\mathbf{E}}_{\text{out}}$分别是输入和输出光场。在某些情况下，传输矩阵可以表示为微分算子。例如，在近场衍射中，传输矩阵可以用菲涅尔衍射公式表示：

$$E_{\text{out}}(x^{\prime },y^{\prime })=\frac{e^{ikz}}{i\lambda z}\iint E_{\text{in}}(x,y)e^{\frac{ik}{2z}[(x^{\prime }-x{)}^2+(y^{\prime }-y{)}^2]}dxdy$$

或者在傅里叶光学中，传输矩阵可以表示为傅里叶变换和相位调制的组合：

E out ( x ′ , y ′ ) = F − 1 { H ( f x , f y ) ⋅ F { E in ( x , y ) } } E_{\text{out}}(x', y') = \mathcal{F}^{-1}\{H(f_x, f_y) \cdot \mathcal{F}\{E_{\text{in}}(x, y)\}\} Eout​(x′,y′)=F−1{H(fx​,fy​)⋅F{Ein​(x,y)}}

其中$H(f_x,f_y)$是光学系统的传递函数。在这种情况下，相位恢复问题可以表述为：给定输出场的强度$|E_{\text{out}}(x^{\prime },y^{\prime }){|}^2$和光学系统的传递函数$H(f_x,f_y)$，恢复输入场$E_{\text{in}}(x,y)$。

高维投影梯度法与模拟退火

对于高维相位恢复问题，投影梯度法是一种有效的方法。其迭代步骤为：

$${\mathbf{x}}_{k+1}={\mathcal{P}}_{\mathcal{C}}\left({\mathbf{x}}_k-{\eta }_k\nabla f({\mathbf{x}}_k)\right)$$

其中${\mathcal{P}}_{\mathcal{C}}$是投影到约束集$\mathcal{C}$上的算子，${\eta }_k$是步长。对于相位恢复，目标函数通常为：

$$f(\mathbf{x})=\frac{1}{2}\sum _{j=1}^m{\left(|⟨{\mathbf{a}}_j,\mathbf{x}⟩{|}^2-y_j\right)}^2$$

其梯度为：

$$\nabla f(\mathbf{x})=\sum _{j=1}^m\left(|⟨{\mathbf{a}}_j,\mathbf{x}⟩{|}^2-y_j\right)\cdot ⟨{\mathbf{a}}_j,\mathbf{x}⟩{\mathbf{a}}_j$$

为了避免局部极小值，可以结合模拟退火技术，引入温度参数$T$：

$${\mathbf{x}}_{k+1}={\mathcal{P}}_{\mathcal{C}}\left({\mathbf{x}}_k-{\eta }_k\nabla f({\mathbf{x}}_k)+{\sigma }_k{\mathbf{\xi }}_k\right)$$

其中${\mathbf{\xi }}_k$是随机噪声，${\sigma }_k=\sqrt{T_k}$是与温度相关的标准差。温度按指数衰减：$T_k=T_0{\alpha }^k$，其中$0<\alpha <1$是衰减因子。