摘要:相位恢复法利用光波传输中某一(或某些)截面上的光强分布来传感系统波前, 其结构简单, 不易受震动及环境干扰, 被广泛应用于光学遥感和像差检测等领域. 传统相位恢复法采用迭代计算, 很难满足实时性要求, 且在一定程度上依赖于迭代转换或迭代优化初值. 为克服上述问题, 本文提出了一种基于卷积神经网络的相位恢复方法, 该方法采用基于小波变换的图像融合技术对焦面和离焦面图像进行融合处理, 可在不损失图像信息的同时简化卷积神经网络的输入. 网络模型训练完成后可依据输入的融合图像直接输出表征波前相位的4—9阶Zernike系数, 且波前传感精度均方根(root-mean-square, RMS)可达0.015λ, λ = 632.8 nm. 研究了噪声、离焦量误差和图像采样分辨率等因素对波前传感精度的影响, 验证了该方法对噪声具有一定鲁棒性, 相对离焦量误差在7.5%内时, 波前传感精度RMS仍可达0.05λ, 且随着图像采样分辨率的提升, 波前传感精度有所改善, 但训练时间成本随之增加. 此外, 分析了实际应用中, 当系统像差阶数与网络训练阶数略有差异时, 本方法所能实现的传感精度, 并给出了解决方案.
关键词: 相位恢复/
卷积神经网络/
图像融合/
传感精度

English Abstract

全文HTML

--> --> -->

空间遥感系统在轨工作时, 由于受环境及系统自身像差的影响, 需实时对系统波前进行监测以实现高分辨成像. 在诸多波前传感方法中, 相位恢复(phase retrieval, PR)法具有传感精度高、系统结构简单、不易受震动和环境干扰等优势, 特别适宜空间应用^[1-4]. 近年来, 已成为空间遥感系统波前传感的首选方案.
1972年Gerhberg和Saxton^[5]提出G-S迭代算法, 奠定了相位恢复方法的核心理论. 随后出现了很多改进算法, 如误差减少算法(ER)^[6]、混合输入输出算法(HIO) ^[6]及杨-顾算法(Y-G) ^[7]等. 但这些迭代方法耗时长, 难以满足实时性需求, 且在一定程度上依赖于迭代转换或迭代优化过程中用到的迭代初值.
为克服传统迭代方法存在的问题, 近年来许多****将人工神经网络与波前传感结合^[8-12], 提出了基于深度学习的波前传感新方法. 2018年, Paine和Fienup^[8]考虑到PR法对迭代初值的依赖性, 提出利用卷积神经网络(convolutional neural network, CNN)对迭代初值进行预测, 提升了波前传感精度. 但该方法应用神经网络后, 仍要进行迭代计算, 并未完全克服传统迭代方法的耗时问题. 同年鞠国浩等^[11]通过对焦面和离焦面图像提取切比雪夫矩特征, 结合反向传播(back-propagation, BP)网络实现相位恢复, 但切比雪夫矩特征提取方法复杂, 编程困难. 2019年齐鑫等^[12]通过对焦面和离焦面图像提取一种与目标无关的特征, 并结合长短期记忆(long short-term memory, LSTM)网络实现了相位变更波前传感方法, 但其用到的网络层数为128层, 结构复杂.
为提高波前传感精度, 并降低系统复杂度, 本文提出了一种基于图像融合和卷积神经网络的相位恢复方法, 该方法采用基于小波变换的图像融合技术对焦面和离焦面的点扩散函数(point spread function, PSF)图像进行融合处理, 可在不损失图像信息的同时简化CNN的输入. 训练好的网络模型, 可依据输入的融合图像直接输出表征波前相位的Zernike系数, 有效地提升波前传感精度的同时简化了操作难度.

PR是一种间接测量波前方法, 即利用可直接测量的强度分布来恢复波前相位^[13]. 由于传统G-S算法仅从一幅焦面图像恢复波前, 无法确保解的唯一性, 故本文采用引入了离焦相位变更的相位变更相位恢复(phase-diversity phase retrieval, PDPR)法^[1,14]. 通过引入已知的离焦相位变更, 消除了光瞳面相位分布的不确定性. 图1为典型的离焦型PDPR法工作原理.
图 1 离焦型PDPR法工作原理图
Figure1. Schematic diagram of PDPR (defocus diversity).

依据标量衍射理论, 焦面光场分布$U(x, y)$可表示为

$\begin{split} &U(x,y) = \frac{{\exp ({\rm{j}}kf)}}{{{\rm{j}}\lambda f}}\exp \Big[{\rm{j}}\frac{k}{{2f}}({x^2} + {y^2})\Big] \iint\limits_{( - \infty, + \infty )} {{U_0}({x_0},{y_0})\exp \left[ { - {\rm{j}} \cdot 2{\rm{\pi }} \cdot \frac{{({x_0}x + {y_0}y)}}{{\lambda f}}} \right]}{\rm{d}}{x_0}{\rm{d}}{y_0}. \end{split} $

同理, 离焦面光场分布$U'(x, y)$可表示为

$\begin{split} U'(x,y) =\; & \frac{{\exp [{\rm{j}}k(f + d)]}}{{{\rm{j}}\lambda (f + d)}}\exp \Big[{\rm{j}}\frac{k}{{2(f + d)}}({x^2} + {y^2})\Big] \iint\limits_{( - \infty, + \infty )} {{U_0}({x_0},{y_0}) \cdot \exp \left[ { - {\rm{j}} \cdot k \cdot d\frac{{(x_0^2 + y_0^2)}}{{2(f + d) \cdot f}}} \right]}\\& \times \exp \left[ { - {\rm{j}} \cdot 2{\rm{\pi }} \cdot \frac{{({x_0}x + {y_0}y)}}{{\lambda (f + d)}}} \right]{\rm{d}}{x_0}{\rm{d}}{y_0} , \\[-15pt] \end{split}$

其中$\lambda $为光波波长; $f$为焦距; $d$为已知离焦距离; ${\rm j^2} = - 1$; $k$为波矢数. 显然, 焦面及离焦面光场均可看作广义光瞳函数的傅里叶变换, 所不同的是, 对于离焦面光场, 光瞳函数的相位因子中引入了一离焦相位. 本文选用Zernike多项式表征相位$\varphi (x, y)$, 如(3)式所示:

$\varphi (x,y) = \sum\limits_{i = 1}^N {{\alpha _i} \cdot {Z_i}(x,y)}, $

其中$i$为阶数; ${Z_i}(x, y)$为第$i$阶Zernike多项式的表达式; ${\alpha _i}$为第$i$阶多项式的系数. 可以看出, 焦面、离焦面上的光场分布与系统波前存在非线性关系, 故可利用CNN建立两者联系, 由已知的光强信息反演出波前相位.

如前所述, 应用于波前传感的神经网络模型主要有BP, CNN, LSTM等. 由于相位恢复法是基于图像信息来实现波前传感, 故本文采用更擅长图片处理的CNN进行相位恢复训练^[15-17].
本文提出了利用CNN拟合焦面、离焦面PSF融合图像与系统波前的非线性关系, 实现由已知的光强信息预测波前的传感方法. 其中, CNN的输入为焦面、离焦面PSF图像经小波变换图像融合方法所得的融合图像, 输出为采用Zernike多项式表征的系统波前信息(即各阶Zernike系数值). 由于Zernike多项式前三项分别表示平移, x, y方向的倾斜对图像像质没有影响, 故在相位恢复过程中不予考虑. 同时, 针对空间遥感系统以低空间频率像差为主的特点, 本文选择了4—9阶Zernike多项式表征系统波前, 分别对应离焦、与轴成0°或90°的像散、与轴成45°的像散、X轴的三级彗差、Y轴的三级彗差和三级球差^[18].
2

-->

3.1.CNN结构

CNN结构有很多变体, 但其基本结构大多包括输入层、卷积层、池化层、全连接层及输出层^[19]. 笔者基于这些基本结构进行多次搭建尝试, 依据训练结果的不同表现, 对网络结构及参数进行调整, 最终确定了如图2所示的网络结构. 该结构共包含22层, 包括1层输入层, 13层卷积层, 6层池化层, 1层平坦层和1层全连接层(即输出层). 处理后的焦面和离焦面融合图像作为网络的输入, 经过卷积核大小为3 × 3, 步长为1的第一层卷积层后得到32个大小为32 × 32的特征图; 接着进入池化层, 采用最大池化, 池化核尺寸为2 × 2, 输出特征图尺寸减半, 即此时特征图大小为16 × 16, 特征图数量仍为32个; 之后经过5次相同的卷积层和池化层的交替结构后, 得到64个大小为1 × 1的特征图; 再经过7层卷积核大小为3 × 3, 步长为1的卷积层, 1层平坦层及1层包含6个神经元的全连接层, 最终输出长度为6的一维数组, 对应4—9阶Zernike多项式系数. 训练过程中应用的学习率为0.0001, 损失函数为均方误差(mean square error, MSE), 表达式为
图 2 CNN结构示意图
Figure2. The schematic of CNN structure.

${\rm{MSE}} = \frac{1}{N}\sum\limits_{i = 1}^N {{{\left( {{\alpha _{{\rm{true}}}} - {\alpha _{{\rm{est}}}}} \right)}^2}},$

其中${\alpha _{{\rm{true}}}}$为Zernike多项式系数真值; ${\alpha _{{\rm{est}}}}$为多项式系数网络预测值; MSE的值越小, 说明预测模型精确度越好. 应用的优化函数为Adam, 该优化算法是一种对随机梯度下降法的扩展, 结合了自适应梯度算法(AdaGrad)和均方根传播算法(RMSProp)的优点^[20]. 应用的激活函数为ReLU (rectified linear unit), 表达式为

$y = \left\{ {\begin{array}{*{20}{c}} x, \\ 0, \end{array}} \right.\,\,\,\,\,\,\,\begin{array}{*{20}{c}} {(x > 0){\rm{,}}} \\ {(x \leqslant 0){\rm{ }}{\rm{.}}} \end{array}$

将(5)式中y对x求导, 则得到

$\frac{{\partial y}}{{\partial x}} = \left\{ {\begin{array}{*{20}{c}} 1, \\ 0, \end{array}} \right.\,\,\,\,\,\begin{array}{*{20}{c}} {(x > 0){\rm{,}}} \\ {(x \leqslant 0){\rm{ }}{\rm{.}}} \end{array}$

可见数据在网络内传播过程中, 正向输入x > 0时, 反向传播会将上游的值原封不动传给下游; 正向输入$ x \leqslant 0 $时, 反向传播中传给下游的信号将停在此处. 因此ReLU函数会使一部分神经元输出为0, 致网络稀疏, 从而减少参数之间的依存性, 防止过拟合, 且与其他激活函数相比不会出现梯度消失问题.
2 -->

3.2.基于小波变换图像融合方法生成CNN输入数据

为简化CNN输入, 降低系统复杂度, 本文采用基于小波变换的极值法对焦面、离焦面PSF图像进行融合. 其核心思想为挑选待融合图像中相应各像素点中像素值的极大或极小值, 作为融合图像在该点的像素值$F(i, j)$^[21], 如(7)式、(8)式所示, 分别称为极大值法和极小值法:

$F(i,j) = \max \{ A(i,j),B(i,j)\}, $

$F(i,j) = \min \{ A(i,j),B(i,j)\}, $

其中$A(i, j)$与$B(i, j)$分别表示待融合图像A和B在$(i, j)$点的像素值.
基于小波变换机理的图像融合原理如图3所示. 首先, 采用长度为4的多贝西小波对焦面、离焦面PSF图像进行1层小波分解^[22], 即采用低通(L)、高通(H)滤波器分别对两幅源图像的水平、垂直方向进行滤波, 得到一个低频带(水平低频垂直低频LL)和三个高频带(水平低频垂直高频LH、水平高频垂直低频HL及水平高频垂直高频HH)子图像, 且能量主要集中在低频部分^[23]; 接着, 对高频带和低频带图像分别采用极大值法和极小值法, 将不同频率信息分别进行融合, 既可不损失图像信息, 又可将焦面、离焦面PSF图像细节信息融合, 减少冗余数据; 最后, 采用第一步小波变换的重构算法对处理后的小波系数进行反变换, 实现图像重构^[24,25].
图 3 基于小波变换机理的图像融合原理
Figure3. Principle of image fusion based on wavelet transform.

该融合方法计算简单、适宜实时处理. 本质上, 融合图像的过程可看作是数据压缩过程, 如图4所示, 在本文中, 两幅64 × 64的焦面、离焦面图像(图4(a)和图4(b)), 经小波变换图像融合处理后合成了一幅64 × 64的融合图像(图4(c)), 同时, 为了最大化排除冗余信息, 仅将融合图像中心32 × 32的部分(图4(d))作为CNN的输入. 相比于将全部图像信息作为输入的传统方法, 采用图像融合技术可有效地减少冗余数据, 加快网络训练速度. 仿真结果表明, 若图像采样率为64 × 64, 以融合图像作为网络的输入数据, 训练时间会缩短1/2, 且随着图像采样率的增加时间优势愈发显著.
图 4 融合图像示意图　(a) 焦面PSF图像; (b) 离焦面PSF图像; (c) 融合图像; (d) CNN的输入图像(即图(c)中红框内部分)
Figure4. Schematic diagram of fusion image: (a) PSF image of the focal plane; (b) PSF image of the defocal plane; (c) the fusion image; (d) the input image of CNN.

建立了网络结构及输入数据生成、处理方式后, 本文开展了大量的仿真实验, 对所提方法的波前传感精度及其影响因素进行了分析.
2

-->

4.1.CNN训练

首先确定系统参数, 生成训练数据. 本文仿真过程中使用单透镜成像系统, 具体系统参数如表1所列.

透镜焦距/mm	入瞳直径/mm	波长/nm	离焦距离/mm
150	10	632.8	4

表1仿真系统参数
Table1.Simulation system parameters.

在系统中引入10000组由4—9阶Zernike多项式表征的随机波前, 其均方根(root-mean-square, RMS)值分布如图5所示. 可以看出, 波前RMS值在$[0.1\lambda, 1.1\lambda ]$范围内, 且近似符合高斯分布.
图 5 10000组随机波前RMS分布图
Figure5. RMS distributions of 10000 random wavefront groups.

依据傅里叶光学原理生成相应的10000对焦面、离焦面PSF图, 并在图中引入信噪比为35 dB的噪声来模拟实际噪声环境. 然后利用小波变换图像融合方法将10000对PSF图像融合, 得到10000张融合图像, 作为网络训练的输入数据. 而输出数据则为表征波前畸变的10000组对应的4—9阶Zernike系数. 训练过程中将10000组训练数据集分成三部分: 80%用作训练集数据, 用以拟合网络权值; 15%用作验证集数据, 用以调整网络最终架构; 5%用作测试集数据, 用以评估网络性能. 网络训练环境为python3.7, tensorflow-gpu2.1.0及keras-gpu2.3.1, 训练所用计算机处理器为Intel(R)Core(TM)i7-8750H CPU @2.20 GHz, 显卡型号为NVIDIA GeForce GTX 1060.
2 -->

4.2.训练结果及分析

完成网络训练后, 可得到一训练完成的网络模型. 采用4.1节所述5%的测试集数据对该网络模型进行测试, 预测4—9阶Zernike系数值, 实现波前传感.
3 -->

4.2.1.网络预测系数与真实系数比较

500组测试数据的真实系数与网络预测系数间的标准差统计直方图如图6所示, 从图中数据可见, 系数间的标准差不大于0.025的占90%. 故利用本文方法得到的Zernike系数与真实系数极为接近.
图 6 真实系数与网络预测系数间的标准差统计直方图
Figure6. Statistical histogram of standard deviation between the real coefficient and the network prediction coefficient.

3 -->

4.2.2.波前传感精度

根据网络输出的各阶Zernike系数值可拟合出系统波前畸变${\varphi _{{\rm{est}}}}$, 将其与引入系统的真实波前畸变${\varphi _{{\rm{true}}}}$做差得到残余波面, 并以该波面RMS值作为评价波前传感精度的指标. 图7所示为500组测试数据波前传感精度RMS统计直方图. 可见, 本文提出的方法传感精度最高可达RMS 0.0016λ, 最低为RMS 0.113λ, 均值为RMS 0.0104λ, 且${\rm{RMS}} \leqslant 0.015\lambda $占比90%, 相比于文献[11]给出的传感精度0.025λ有所提高.
图 7 测试数据波前传感精度RMS统计直方图
Figure7. RMS statistical histogram of wavefront sensing accuracy of test data.

2 -->

4.3.影响本方法波前传感精度的因素分析

3 -->

4.3.1.像差阶数的影响

如前文所述, 遥感系统在轨工作时, 系统像差以低空间频率误差为主, 故本文采用了4—9阶Zernike多项式表征波前畸变, 由于系统实际像差阶数未必局限于4—9阶, 考虑到CNN网络的普适性, 研究了系统存在略低或略高阶数像差时, 本文方法的波前传感精度.
假设系统实际像差阶数分别为4—7阶、4—8阶、4—10阶及4—11阶. 其中第4—9阶像差系数范围仍为$[ - 0.5\lambda, 0.5\lambda ]$, 第10, 11阶的像差系数范围为$[ - 0.1\lambda, 0.1\lambda ]$, 基于表1的系统参数, 产生各500组验证数据, 然后运用训练好的网络, 预测4—9阶Zernike系数值, 实现波前传感. 得到的传感精度如表2所列. 由表2可以看出, 当实际像差阶数小于9阶时, 本方法的传感精度RMS呈上升趋势; 当实际像差阶数在10阶以内时, 传感精度RMS仍可控制在0.015λ以内; 当实际像差阶数达到11阶时, 传感精度RMS降至0.1λ, 通常情况下, 遥感系统的传感精度需达到0.05λ, 故已很难满足遥感系统波前传感要求. 鉴于CNN网络的普适性, 此时, 可依据遥感系统像差空间频率特性, 基于本文所提方法增加网络输出层单元数, 并调整网络结构, 训练适用于更高阶像差的网络.

像差阶数	4—7阶	4—8阶	4—10阶	4—11阶
传感精度RMS/$\lambda $	0.010	0.012	0.015	0.100

表2系统存在不同阶像差时本文方法的波前传感精度
Table2.The wavefront sensing accuracy of the proposed method when the system has different order of aberration.

3 -->

4.3.2.离焦量误差的影响

受实际图像采集过程中探测器定位精度的限制, 应用相位恢复法进行波前传感时, 很难保证预设的离焦量被精确引入系统. 为此, 笔者在传感系统中引入了一定量的相对离焦量误差, 研究了该误差对传感精度的影响.
相对离焦量误差$D$的定义为

$D({\text{%}}) = \frac{{\delta ({\rm{mm}})}}{{d({\rm{mm}})}} \times 100{\text{%}}, $

其中$\delta $为离焦量误差; $d$为离焦量. 仿真中, 将相对离焦量误差设置为2.5%, 5.0%, 7.5%和10.0%, 应用本文方法进行波前传感, 所得传感精度如表3所列.

相对离焦量误差/%	2.5	5.0	7.5	10.0
传感精度RMS/$\lambda $	0.022	0.035	0.050	0.065

表3不同离焦量误差下, 本文方法的传感精度
Table3.The sensing accuracy of the proposed method under different defocusing errors.

可以看出, 当相对离焦量误差在2.5%—7.5%时, 传感精度RMS值介于0.022λ—0.050λ之间, 离焦量误差对传感精度影响不大. 当相对离焦量误差达到10%时, 传感精度RMS降至0.065λ, 已无法满足通常情况下遥感系统的传感要求. 可见, 本文所提方法对系统相对离焦量误差的容限为7.5%, 随着离焦量误差的增大, 传感精度随之降低.
3 -->

4.3.3.噪声的影响

采用探测器获取焦面、离焦面图像时, 图像噪声会影响波前传感精度, 为此在PSF图像上叠加了不同程度的高斯噪声, 研究了不同信噪比情况下, 本文所提方法的传感精度. 其中图像信噪比SNR定义为

${\rm{SNR}} = 10 \times \log \left(\frac{{\sigma _{\rm{s}}^2}}{{\sigma _{\rm{n}}^2}}\right), $

其中$\sigma _{\rm{s}}^{\rm{2}}$和$\sigma _{\rm{n}}^{\rm{2}}$分别表示图像和噪声的方差.
研究结果如表4所列. 由表4可以看出, 当SNR为30—50时, 该方法的传感精度介于$0.015\lambda$—$0.020\lambda $之间. 当SNR降至25时, 传感精度RMS降为$0.06\lambda $, 已不能满足通常情况下遥感系统波前传感要求.

信噪比/dB	50	40	35	30	25
传感精度RMS/$\lambda $	0.015	0.015	0.015	0.020	0.060

表4噪声对传感精度的影响
Table4.The influence of noise on the sensing accuracy.

3 -->

4.3.4.图像采样分辨率的影响

当采样分辨率提升时, 获取的PSF图像信息将更为丰富, 为此, 采用128 × 128的图像训练网络并进行测试, 同时与上文采用的64 × 64图像训练结果进行对比. 不同采样率下, 500组验证数据所得的波前传感精度RMS对比图如图8所示. 此时, 波前传感精度最高达RMS 0.0004λ, 最小为 RMS0.1500λ, 均值为RMS 0.0071λ且${\rm{RMS}} \leqslant 0.015\lambda$占比92%. 可以看出, 随着图像采样分辨率的提升, 波前传感精度有所改善, 但网络输入数据亦随着采样率的升高而增加, 训练时间成本随之增加.
图 8 采用不同分辨率图像训练网络时, 波前传感精度RMS统计直方图对比　(a) 图像采样率 64 × 64; (b) 图像采样率128 × 128
Figure8. When image training networks with different resolutions are adopted, the RMS statistical histogram of wavefront sensing precision are compared: (a) Image sampling rate 64 × 64; (b) image sampling rate 128 × 128

本文提出了一种基于小波变换图像融合和卷积神经网络的相位恢复方法, 采用小波变换图像融合技术对焦面、离焦面图像进行处理, 在不损失图像信息的同时简化网络输入. 将融合所得图像作为CNN的输入数据, 以相应的4—9阶Zernike系数作为CNN的输出数据进行训练. 训练结果表明当波前RMS值在1.1λ内时, 该方法波前传感精度RMS值可达0.015λ.
研究了噪声、离焦量误差、图像采样分辨率等因素对波前传感精度的影响, 验证了该方法对噪声具有一定的鲁棒性, 且相对离焦量误差在7.5%内时, 波前传感精度RMS仍可达0.05λ. 同时, 随着图像采样分辨率的提升, 波前传感精度有所改善, 但网络输入数据亦随着采样率的升高而增加, 训练时间成本随之增加.
此外, 在实际应用中, 为实现较高精度波前传感, 可依据系统像差特性, 基于本文方法, 改变网络输出层单元个数, 并调整网络结构, 训练适用于更高阶像差的新网络.
该方法与目前现有的相位恢复方法相比, 无需进行迭代处理, 更能满足实时性需求. 与现有的结合深度学习的相位恢复方法相比, 数据处理方法简单易行, 且精度更高, 除用于空间遥感系统外, 亦可用于天文望远镜波前传感, 大型光学镜面面形检测等领域.