光子学报  2018, Vol. 47 Issue (5): 0511001  DOI: 10.3788/gzxb20184705.0511001
0

引用本文  

谭政, 相里斌, 吕群波, 方煜, 孙建颖, 赵娜. 基于像差选择性校正的光学-数字联合设计[J]. 光子学报, 2018, 47(5): 0511001. DOI: 10.3788/gzxb20184705.0511001.
TAN Zheng, XIANG-LI Bin, LV Qun-bo, FANG Yu, SUN Jian-ying, ZHAO Na. Optics/Digital Processing Co-design Method Based on Aberration Optional-correcting[J]. Acta Photonica Sinica, 2018, 47(5): 0511001. DOI: 10.3788/gzxb20184705.0511001.

基金项目

国家自然科学基金(No.61505219),中国科学院国防科技创新基金(No.CXJJ-16S045)和青岛市光电智库联合基金项目资助

第一作者

谭政(1982-), 男, 工程师, 硕士, 主要研究方向为计算光学成像技术、遥感数据增强.Email:tanzheng@aoe.ac.cn

导师

相里斌(1967-), 男, 研究员, 博士, 主要研究方向为计算光学成像技术、干涉成像光谱技术等.Email:xiangli@aoe.ac.cn

通讯作者

吕群波(1979-), 男, 研究员, 博士, 主要研究方向为计算光学成像技术.Email:lvqunbo@aoe.ac.cn

文章历史

收稿日期:2017-11-21
录用日期:2018-01-12
基于像差选择性校正的光学-数字联合设计
谭政1,2, 相里斌1, 吕群波1,2, 方煜1, 孙建颖1, 赵娜1    
(1 中国科学院计算光学成像技术重点实验室, 北京 100094)
(2 中国科学院大学, 北京 100039)
摘要:提出面向成像系统的光学和数字图像处理联合优化设计.根据泽尼克多项式模型分析了每种像差对成像清晰度的影响;在变分贝叶斯框架下,提出了基于加权双向差分先验模型的像差数字校正算法;将像差分为易于数字校正和不易于数字校正两类,建立基于像差选择性校正策略的光学-数字处理联合设计方法.采用联合设计方法设计了三镜片大像差光学系统,其奈奎斯特频率处的MTF值约为0.50,采用传统方法设计的六镜片光学系统MTF值约为0.53,两系统成像质量相当,本文方法可降低系统复杂度.
关键词成像系统    光学和数字图像处理    联合优化    像差    数字校正算法    选择性校正    
中图分类号:O439;TP751      文献标识码:A      文章编号:1004-4213(2018)05-0511001-9
Optics/Digital Processing Co-design Method Based on Aberration Optional-correcting
TAN Zheng1,2, XIANG-LI Bin1, LV Qun-bo1,2, FANG Yu1, SUN Jian-ying1, ZHAO Na1    
(1 Key Laboratory of Computation Optical Imaging Technology, Chinese Academy of Sciences, Beijing 100094, China)
(2 University of Chinese Academy of Sciences, Beijing 100039, China)
Foundation item: The National Natural Science Foundation of China (No. 61505219), and the National Defense Technology Innovation Foundation of Chinese Academy of Sciences (No. CXJJ-16S045)
Abstract: This paper is aiming at jointly optimizing the design and analysis of the optics and digital image processing for imaging systems. The influence of each aberration on the imaging sharpness is analyzed according to the Zernike polynomial model. A digital image processing algorithm is proposed for aberration correction, which is based on the variational bayesian framework and weighted bi-direction difference prior. Then the aberrations are divided into two categories:easy to digital correct and hard to digital correct, so an optics-digital processing co-design method can be established based aberration optional-correction. At last, a three-lens large aberration optical system is designed by the co-design method, its MTF value at the nyquist frequency is about 0.50; in contrast, a six-lens system is designed by traditional method, and its MTF value is about 0.53. These results show that the imaging quality of two system is similar, therefore, the proposed method would be capable of reducing the complexity of the optical system.
Key words: Optical imaging system    Optics and digital image processing    Jointly optimizing    Aberration    Correction algorithm    Optional-correction    
OCIS Codes: 110.0110;220.4830;220.1010;200.4740;100.2980
0 引言

现代光学成像系统的传统设计方法主要由两步构成:第一步,根据像差理论设计光学系统[1],通过光线追迹来选择光学元件平衡各种像差的影响,尽可能地校正像差获得高质量的图像;第二步,根据需求设计数字图像处理算法提高第一步中成像结果的清晰度修正剩余缺陷.而这种顺序式的设计方式也会带来一定的局限性:首先,传统设计方式使光学系统的结构十分复杂,为了尽可能校正各种像差需使用较多的光学镜片或较复杂的镜片面型,但是各级次的像差数量仍远远多于光学元件结构参数的数量,导致每种像差都会有一定的残留,而且衍射对成像清晰度影响也不可避免.其次,尽管可以通过图像复原算法来提高图像的信噪比和对比度从而提高成像清晰度,削弱残余像差和衍射的影响,但提高信噪比和提高对比度相互制约,图像平滑可提高信噪比但会造成对比度下降,图像锐化能提高对比度但也会放大噪声使信噪比降低,此外,图像复原一般都采用反卷积算法,存在求解的病态性问题[2].再次,光学设计参数和数字图像处理设计参数之间关联度较弱,导致成像系统的设计成为局部最优.光学设计一般不考虑图像复原算法的设计需求,图像复原处理也不考虑光学设计参数对算法的影响.而一个最优的光学系统和一个最优的图像处理算法的简单结合,并不一定意味着从目标到图像的最优化设计.

因此,为解决以上问题,近些年也发展出了多种新型的计算成像技术[3-8],将传统的光学设计和数字图像信号处理相融合,对成像系统进行联合优化,提高系统的指标,如增大系统景深、降低光学系统复杂度等.波前编码技术[3-4]是在光学系统孔径光阑处插入一面型复杂的非球面透镜编码板,对非相干波前进行编码,使得当焦平面沿光轴在一定范围内移动时,由离焦像差引起的图像模糊一致,令成像对离焦点扩散函数的作用不敏感,再通过数字处理算法对模糊图像进行解码,增大成像的景深,但该系统编码板的面型计算十分复杂;文献[5]和文献[6]提出一种将光学传递函数和图像复原滤波器联合设计的计算成像系统,将重建图像和理想图像之间差值的均方根作为联合优化函数,这种设计方法仅从图像处理算法角度对成像结果进行优化,并没有考虑光学像差的影响,引入的物理规律不充分.文献[7]和文献[8]将光程差函数与图像复原函数联系到一起,实现光学系统和数字处理系统的联合优化,并且一定程度上简化了光学系统的结构.但是,这种联合设计方法需已知目标的功率谱密度,只能应用于黑白文档、条形码等目标,并且其仅用维纳滤波器作为数字处理系统,难以解决图像复原时的病态问题,限制了该设计方法在其它类型光学系统上的推广.

本文从每种光学像差的物理特性——调制传递函数不尽相同为出发点,研究成像系统的光学-数字联合设计方法,首先分析每种光学像差对成像的影响,再基于变分贝叶斯框架设计了加权双向变分先验的像差数字校正算法来修正像差、抑制系统噪声、提高调制度,进而提出了像差选择性校正策略,将像差数字校正算法作为光学系统设计的补偿环节,分担光学系统的部分压力,一定程度上降低光学系统的复杂度,为成像系统的简化设计理论提供一定的参考依据.

1 像差物理特性分析

在空间域上像差的物理特性反映为点扩散函数,在频率域上像差的特性反映为调制传递函数.为方便计算分析,采用波像差来分析每种像差的特性,光学系统出射光瞳中(x, y)点的出瞳函数[9]p(x, y)表示为

$ P\left( {x,y} \right) = a\left( {x,y} \right){{\rm{e}}^{{\rm{i}}kW\left( {x,y} \right)}} $ (1)

式中,a(x, y)表示光学系统孔径函数,W(x, y)表示各个光学镜片等元件引起的波像差函数,该波像差函数可用一系列正交泽尼克圆多项式的线性组合来表示为

$ W\left( {x,y} \right) = \sum\limits_n^k {\sum\limits_{m = - n}^n {W_n^mZ_n^m\left( {x,y} \right)} } $ (2)

式中,Wnm表示泽尼克多项式的系数,n为泽尼克多项式的径向频率数,表示多项式的阶数;m为泽尼克多项式的角向频率数,表示第n阶的模.由于泽尼克多项式组成一个正交的完备集,光学系统的每种像差就是这些完备的基矢所张Hilbert空间中的一个矢量,通过选择泽尼克多项式的不同参数,就能够实现将每种单色像差单独分离出来进行计算分析.

对于圆形孔径,式(2)表示为极坐标形式为

$ W\left( {\rho ,\theta } \right){\rm{ = }}\sum\limits_j {{W_j}{Z_j}\left( {\rho ,\theta } \right)} \left| {_{x = \rho \cos \theta ,y = \rho \sin \theta }} \right. $ (3)

式中,圆内泽尼克多项式的极坐标形式为

$ Z_n^m\left( {\rho ,\theta } \right){\rm{ = }}\left\{ {\begin{array}{*{20}{c}} {\sqrt {2\left( {n + 1} \right)} R_n^m\left( \rho \right)\cos \left( {m\theta } \right)}&{m \ne 0\;且\;j\;为奇数}\\ {\sqrt {2\left( {n + 1} \right)} R_n^m\left( \rho \right)\sin \left( {m\theta } \right)}&{m \ne 0\;且\;j\;为偶数}\\ {\sqrt {\left( {n + 1} \right)} R_n^0\left( \rho \right)}&{m = 0} \end{array}} \right. $ (4)
$ R_n^{\left| m \right|}\left( \rho \right) = \sum\limits_{s = 0}^{\left( {n - \left| m \right|} \right)/2} {\frac{{{{\left( { - 1} \right)}^s}\left( {n - s} \right)!}}{{s!\left[ {0.5 \times \left( {n + \left| m \right|} \right) - s} \right]!\left[ {0.5 \times \left( {n - \left| m \right|} \right) - s} \right]!}}} {\rho ^{n - 2s}} $ (5)

由于光学系统受初级像差影响较大,因此,在本文分析中先不考虑二级以上像差的影响,另外,像散可表示为子午场曲和弧矢场曲的差值[10],故可认为场曲的与像散的特性近似一致.则需分析畸变、像散、离焦、彗差、球差的特性,它们对应的泽尼克多项式参数[11]为:畸变,n=1,m=±1;像散,n=2,m=±2;离焦,n=2,m=0;彗差,n=3,m=±1;球差,n=4,m=0.

像面上像点的点扩散函数表示为

$ {\rm{PSF}} = {\rm{F}}{{\rm{T}}^{ - 1}}\left[ {P\left( {{s_x},{s_y}} \right) \cdot {P^ * }\left( {{s_x},{s_y}} \right)} \right] $ (6)

式中,FT-1表示逆傅里叶变换,P(sx, sy)表示出瞳函数p(x, y)的傅里叶变换,符号“*”表示共轭.像差的光学调制传递函数(Modulation Transfer Function,MTF)与点扩散函数的关系为

$ {\rm{MTF}}\left( {{s_x},{s_y}} \right) = \left\| {\frac{{{\rm{FT}}\left\{ {{\rm{PSF}}} \right\}}}{{{\rm{FT}}\left\{ {{\rm{PSF}}} \right\}\left| {_{{s_x} = 0,{s_y} = 0}} \right.}}} \right\| $ (7)

式中,FT表示逆傅里叶变换.

结合式(3)、(6)、(7)即可计算出每种像差的调制传递函数,计算时取波长为可见光波段的中心波长为570 nm,以瑞利判据作为每种像差的计算容限,即波像差为1/4波长.

图 1中可以看出,尽管波像差大小相同,但不同的像差的调制传递函数曲线并不相同,在调制传递函数曲线中,高频部分反映了该像差对成像目标细节、边缘的影响,中频部分反映了对成像目标层次的影响,低频部分反映了对成像目标轮廓的影响.同一频率下调制传函数的值越低说明此种像差对成像的影响越大,即导致的清晰度损失越大.

图 1 每种像差的调制传递函数曲线 Fig.1 Modulation transfer function curves of each aberration
2 像差数字校正算法

成像系统对目标的观测就是在随机性噪声的影响下,目标场景与像差点扩散函数进行卷积运算的结果,观测方程为

$ \mathit{\boldsymbol{y}} = \sum\limits_k {\mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{h}}_k}} \right)} \otimes \mathit{\boldsymbol{x}} + \mathit{\boldsymbol{n}} $ (8)

式中,x为观测目标,hk为像差点扩散函数,下标k代表畸变、像散、离焦、彗差、球差、场曲等像差,B为Toeplitz矩阵转化算子,n为由探测器和外界环境导致的加性噪声,y为在噪声和各种像差综合作用下目标的成像结果.像差数字校正就是利用图像复原理论从该卷积运算中尽量消除点扩散函数对成像的影响,削弱反卷积的病态性对求解的稳定性和收敛性的作用,权衡对比度增强和信噪比提升的关系,提高像差校正的精度,恢复出清晰的目标图像.

在贝叶斯框架下,式(8)的解表示为

$ p\left( {\mathit{\boldsymbol{x}}\left| \mathit{\boldsymbol{y}} \right.} \right) = p\left( {\mathit{\boldsymbol{y}}\left| \mathit{\boldsymbol{x}} \right.} \right)p\left( \mathit{\boldsymbol{x}} \right)/p\left( \mathit{\boldsymbol{y}} \right) $ (9)

令噪声服从零均值的高斯分布,则

$ p\left( {\mathit{\boldsymbol{y}}\left| \mathit{\boldsymbol{x}} \right.,\beta } \right) \propto {\beta ^{P/2}}\exp \left[ { - \frac{\beta }{2}{{\left\| {\mathit{\boldsymbol{y}} - \sum\limits_k {\mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{h}}_k}} \right)\mathit{\boldsymbol{x}}} } \right\|}^2}} \right] $ (10)

式中,P表示图像的像素总数,未知参数β为超参数,表示高斯分布方差的倒数.

图像中每个像素点可以用马尔科夫模型描述,则p(x)也可表示为吉布斯分布的形式

$ p\left( {\mathit{\boldsymbol{x}}\left| \alpha \right.} \right) = Z\left( \alpha \right)\exp \left[ { - \alpha U\left( \mathit{\boldsymbol{x}} \right)} \right] $ (11)

式中,Z(α)为关于常数α的配分函数,U(x)为能量函数.那么求解方程(9)的问题就转化为寻找一个恰当的能量函数,作为方程指标函数收敛的引导,通过迭代使方程的解尽量逼近真解.

采用全变分先验(Total Variation Prior, TV)模型作为能量函数,因其较强的空间约束性,使算法在边缘保持和噪声滤除中都取得了较好的图像复原结果[12-15],但是从反卷积算法逆问题的角度,TV先验仍旧会存在求解不稳定和局部最优解的问题,导致方程解还存在一定的病态性.为解决这个问题,本文采用式(12)形式的加权双向差分模型作为能量函数.

$ U\left( \mathit{\boldsymbol{x}} \right) = \sum\limits_i {\sqrt {{{\left[ {\Delta _{hi}^ \leftarrow \left( \mathit{\boldsymbol{x}} \right)} \right]}^2} + {{\left[ {\Delta _{vi}^ \uparrow \left( \mathit{\boldsymbol{x}} \right)} \right]}^2} + \lambda \left\{ {{{\left[ {\Delta _{hi}^ \to \left( \mathit{\boldsymbol{x}} \right)} \right]}^2} + {{\left[ {\Delta _{vi}^ \downarrow \left( \mathit{\boldsymbol{x}} \right)} \right]}^2}} \right\}} } @\sum\limits_i {\sqrt {{\mathit{\Lambda} _i}\left( \mathit{\boldsymbol{x}} \right)} } $ (12)

式中,“箭头”符号表示沿着该方向向点xi相邻的像素点取一阶差分,Δhi和Δvi分别表示水平和垂直方向的差分算子,为防止过平滑,引入常数权值,因此,该能量函数相当于沿着像素xi水平和垂直的正负两个方向做差分.与TV先验相比,该先验模型通过建立每个像素点和周围像素点的不一致关系对高频噪声加以约束,进一步增强了先验信息的空间约束性,更能削弱反卷积本身造成的噪声放大和传播等缺点,使算法在提高图像调制度的同时能够保证较高的信噪比,从而进一步的减轻成像方程求解的病态性,获得更优的解.

定义未知量的集合为Θ={x, α, β},则式(9)变为对p(Θ|y)的求解,而该分布不易计算,引入变分贝叶斯推论[16]解决该问题,即用易于计算的分布q(Θ)来近似,使q(Θ)与后验分布p(Θ|y)之间的Kullback-Leibler(KL)距离最小为

$ {C_{{\rm{KL}}}}\left( {q\left( \Theta \right)\left\| {p\left( {\Theta \left| \mathit{\boldsymbol{y}} \right.} \right)} \right.} \right) = \int {q\left( \Theta \right)\log \left( {\frac{{q\left( \Theta \right)}}{{p\left( {\Theta ,\mathit{\boldsymbol{y}}} \right.}}} \right){\rm{d}}\Theta } + {\rm{const}} $ (13)

式中,q(Θ)=q(x)q(α)q(β),联合分布p(Θ, y)=p(y|x, βp(x|αp(αp(β),对式(13)应用最大最小方法,即要使KL距离最小就是使联合分布p(Θ, y)的下边界最大,具体推导过程可参见文献[16],以下对该过程做简要描述.

对式(11)应用均值不等式

$ p\left( {\mathit{\boldsymbol{x}}\left| \alpha \right.} \right) \propto \exp \left[ { - \alpha \sum\limits_i {\sqrt {{\mathit{\Lambda }_i}\left( \mathit{\boldsymbol{x}} \right)} } } \right] \ge \exp \left[ { - \frac{\alpha }{2}\sum\limits_i {\frac{{{\mathit{\Lambda } _i}\left( \mathit{\boldsymbol{x}} \right) + {\omega _i}}}{{\sqrt {{\omega _i}} }}} } \right] $ (14)

式中,令Υ$(\alpha, \mathit{\boldsymbol{x}}, \mathit{\boldsymbol{\omega}})@{\rm{exp}}\left[\left(-\alpha /2 \right)\sum\limits_{i}{\left( {{\mathit{\Lambda }}_{i}}\left( \mathit{\boldsymbol{x}} \right)+{{\omega }_{i}} \right)}~/\sqrt{{{\omega }_{i}}} \right]~$,则联合分布p(Θ, y)的下边界的最大值P(Θ, y, ω)=p(y|x, β)Υ(α, x, ω)p(α)p(β),此时,KL距离

$ {C_{{\rm{KL}}}}\left( {q\left( \Theta \right)\left\| {p\left( {\Theta \left| \mathit{\boldsymbol{y}} \right.} \right)} \right.} \right) \le \int {q\left( \Theta \right)\log \left( {\frac{{q\left( \Theta \right)}}{{p\left( {\Theta ,\mathit{\boldsymbol{y}},\mathit{\boldsymbol{\omega }}} \right)}}} \right){\rm{d}}\Theta } $ (15)

令式(15)右边对q(x)的导数为零,可得KL距离最小时的q(x)估计值$\mathit{\hat{q}}\left( \mathit{\boldsymbol{x}} \right)$

$ \mathop q\limits^ \wedge \left( \mathit{\boldsymbol{x}} \right) = c \cdot \exp \left[ {{E_{q\left( \mathit{\boldsymbol{x}} \right)}}\left( {\log P\left( {\Theta ,\mathit{\boldsymbol{y}},\mathit{\boldsymbol{\omega }}} \right)} \right)} \right] $ (16)

Eq(x)表示q(x)=$\mathit{\hat{q}}\left( \mathit{\boldsymbol{x}} \right)$时的x均值.$\mathit{\hat{q}}\left( \mathit{\boldsymbol{x}} \right)$服从多变量高斯分布,其协方差为

$ {\rm{Co}}{{\rm{v}}_{q\left( \mathit{\boldsymbol{x}} \right)}} = \beta {\left[ {\sum\limits_k {B\left( {{\mathit{\boldsymbol{h}}_k}} \right)} } \right]^{\rm{T}}}\left[ {\sum\limits_k {B\left( {{\mathit{\boldsymbol{h}}_k}} \right)} } \right] + \alpha \sum\limits_{d = 1}^4 {{\lambda _{\rm{d}}}{{\left( {{\Delta _{\rm{d}}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{W}}\left( {{\Delta _{\rm{d}}}} \right)} $ (17)

式中,Δd表示式(12)中四个方向的差分矩阵,λd表示权值,$\mathit{\boldsymbol{W=}}{\rm{diag}}{{\left( {{w}_{i}} \right)}^{-\frac{1}{2}}}$.均值Eq(x)表达式为

$ {E_{q\left( \mathit{\boldsymbol{x}} \right)}} = \beta {\left[ {{\rm{Co}}{{\rm{v}}_{q\left( \mathit{\boldsymbol{x}} \right)}}} \right]^{ - 1}}{\left[ {\sum\limits_k {\mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{h}}_k}} \right)} } \right]^{\rm{T}}}y $ (18)

在求解超参数αβ时,同样可令式(15)右边分别对q(α)和q(β)的导数为零,分别得到,$\alpha =\frac{P}{2}\sum\limits_{i}{{{\left( {{\omega }_{i}} \right)}^{-\frac{1}{2}}};\beta =P/{{\left\| \mathit{\boldsymbol{y}}-\sum{\mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{h}}}_{k}} \right)\mathit{\boldsymbol{x}}} \right\|}^{2}}}$,对式(18)采用通过共轭梯度法即可解得x.为了说明本算法的有效性,对图像处理常用的Lena图像同时施以半径为9、方差为3的高斯模糊,和均值为0、方差为3的随机噪声;比较本文的加权(λ=1.3)双向差分先验和TV先验对该模糊图像的恢复效果,采用与标准无降质图像的峰值信噪比(Peak Signal to Noise Ratio, PSNR)作为评价指标,PSNR值越大表明复原结果和无降质的原图越接近,其调制度和信噪比就越高.

本文方法与TV先验方法[13-15]复原结果PSNR值随迭代次数的收敛情况如图 2所示,可以看出两种算法在相同的迭代停止条件(‖xn+1-xn‖/‖xn‖≤1×10-4)下,基于TV先验的复原结果的最优值并不是最终收敛值,而本文算法的迭代收敛性相对较好,复原结果的PSNR值也更高,因此,本文算法更稳定有效,进一步减轻了反卷积算法的病态性.

图 2 加权双向差分先验和TV先验复原峰值信噪比随迭代次数的收敛情况对比 Fig.2 Comparison of convergence of PSNR with iterative times under weighted bi-directional difference priors and TV prior restoration
3 光学/数字联合设计——像差选择性校正策略

根据前述的各像差点扩散函数,以及像差数字校正算法,就可以分析各种像差对数字校正算法的敏感程度,若不考虑噪声的影响,用四分之一波长波像差下畸变、像散、离焦、彗差、球差的点扩散函数分别作用于标准无降质的Lena图像,再用上节提出的数字校正算法进行像差的数字校正,分析每种像差校正前后的PSNR值,如表 1所示.

表 1 每种像差校正前后图像的PSNR值对比 Tab.1 Comparison of PSNR between pre and post aberration digital correction

表 1中第一行代表像差导致的模糊图像相对于原图的PSNR值,第二行代表经数字校正后复原图像相对于原图的PSNR值;对于畸变而言,由于其不会使成像产生模糊,因此在数字校正前后图像的PSNR值相差不大;对于彗差、像散、离焦、球差而言,由于它们的波像差大小一致,均为四分之一波长,因此这四种像差造成图像模糊的PSNR值大小差别不大;而它们数字校正结果的PSNR值相差较大,这就说明这几种像差对数字校正算法的敏感程度并不一样,从PSNR值看来,PSNR值越大,说明这种像差对数字校正算法更敏感,就更容易数字校正.

由于彗差引起的模糊经数字校正后改善最大,因此彗差对数字处理校正的敏感度最高,其次为像散和场曲,而对离焦和球差导致模糊的改善都相对较小,因此,离焦和球差对数字处理校正的敏感度相对较低.图 3中,从目视效果来看,彗差和像散的数字校正后的图像质量也要优于离焦和球差的数字校正结果.

图 3 彗差、像散、离焦、球差校正前后图像对比 Fig.3 Comparison of four kinds of aberration`s digital correction results

综合以上分析,提出一种成像系统的光学/数字联合设计思路,即对用数字校正算法补偿敏感度相对较低的像差留给光学设计校正,对数字处理校正敏感度相对较高的像差用数字算法校正,这就是像差选择性校正策略.

4 仿真实验对比

由于在实际应用中光学系统通常都是对可见光等宽谱段入射光成像,此时光学像差点扩散函数是所有入射光波长对应的各种像差的点扩散函数综合作用的结果,在传统设计中为获得较好的成像质量,都通过增加光学镜片数量或提高镜片面型复杂度等方式,将像差导致的点扩散效应优化至衍射艾利斑范围内,也就是使光线追迹得到的各个波长综合的点列图RMS(Root Mean Square, 均方根)半径小于艾利斑半径,此时, 衍射艾利斑导致的点扩散效应对光学系统成像质量的影响要大于像差的点扩散效应.本文所述的像差选择性校正方法就是先放宽光学系统的约束,利用镜片数量相对较少或面型复杂度相对较低的光学系统,使像差对像质的影响大于衍射对像质的影响,并且以彗差、像散等较易于数字校正的像差为主,先获得一个成像质量较模糊的结果,然后使用像差数字校正算法对该成像结果进行处理,最终获得一幅较为清晰的目标图像.

仿真实验采用Zemax软件设计了一个六镜片透射式成像系统和一个三镜片透射式成像系统,两系统除镜片数量外参数相同,分别为:焦距100 mm、F数4、视场12°、成像谱段范围480~680 nm.六镜片成像系统采用传统的方法设计,光学系统使用6个镜片对像差进行控制,使全谱段内成像的光线追迹点列图RMS半径为1.46 μm,艾利斑半径2.86 μm,因此,该系统的像质主要受衍射影响,点扩散效应以衍射艾利斑为主,其归一化的点扩散函数如图 4(a)所示,右侧柱状图表示能量相对强度.三镜片系统采用本文光学/数字联合的方法设计,光学系统采用3个镜片,放宽对光学像差的约束,在全谱段内成像的点列图RMS半径为6.75 μm,艾利斑半径与六镜成像系统相同为2.86 μm,因此,该成像系统的像质主要受像差影响,其归一化的像差点扩散函数如图 4(b)所示,主要以像差导致的点扩散为主.

图 4 两系统点扩散函数 Fig.4 Point spread function of 6 lens system and 3 lens system

由于点扩散函数反映了物面一点的光能在像面的分布,为进一步说明三镜片成像系统的像差组成,在成像谱段范围480~680 nm内,分别对480 nm、530 nm、580 nm、630 nm、680 nm五种波长光线追迹点列图进行采样,结果如图 5(a)~(f)所示,各图中纵向长度单位为微米.

图 5 三镜片成像系统光线追迹点列图 Fig.5 Spot diagram through ray trace of Three-lens imaging system

根据图 5(b)~(f)光线追迹点列图的形状,以及像散、离焦、彗差、球差的定义[17]可以看出,在表 2所示的各波长光线追迹点列图的RMS半径范围内,光线的分布主要以像散、彗差、场曲为主,离焦和球差主要表现在点列图RMS半径范围以外的光线分布中,而像差点扩散函数由RMS半径范围内的光线决定,范围外的光线对能量分布影响较小.

表 2 图 5(a)~(f)的点列图RMS半径和艾利斑半径 Tab.2 The spot diagram RMS radius and ariy disk radius of Fig. 5(a)~(f)

仿真实验中,将两个成像系统都赋以均值为0、方差为0.5的随机加性噪声,采用Lena图像和成像质量测试专用的ISO12233标准靶标图像[18]中刃边部分作为成像目标进行成像仿真,并测量[9]两个系统调制传递函数曲线.对三镜片系统成像结果使用数字校正算法修正残余像差,对六镜片系统成像结果使用本文的数字处理算法进行增强.两系统的MTF曲线和结果对比如图 6.

图 6 三镜片系统MTF和六镜片系统MTF结果对比 Fig.6 Comparison of MTF results between Three-lens imaging system and Six-lens system

通过图 6图 7可以看出,三镜片成像系统在像差数字校正前后的MTF值提高明显,并且三镜片成像系统的光学/数字联合设计成像结果与六镜片成像系统的成像质量相比,仅在部分细节纹理方面略有不足,三镜片系统奈奎斯特频率处的MTF值约为0.50,六镜片系统奈奎斯特频率处的MTF值约为0.53,两系统成像质量相当,而三镜片成像系统由于使用像差数字校正分担了光学设计校正像差的部分压力,使系统在镜片数量上减小了一倍,这些结果除证明了本文提出的像差选择性校正方法外,还说明通过本文的光学/数字联合设计方法,能够在一定程度上简化光学系统的结构.

图 7 三镜片系统联合设计成像结果和六镜片系统传统设计成像结果对比 Fig.7 Comparison of imaging results between Three-lens system by co-designed method and Six-lens system by traditional-designed method
5 结论

本文面向光学/数字处理联合优化设计,首先,根据波像差理论,分析了像散、离焦、彗差、球差、畸变、场曲等单色像差对光学信息传递的影响;其次,基于贝叶斯框架和加权双向差分先验模型提出了像差数字校正算法,获得了与全变分先验模型相比更稳定的图像复原结果,进一步削弱了反卷积求解的病态性,能够同时提高系统成像的调制度和信噪比;最终,结合像差点扩散函数的物理特性和像差数字校正算法,提出了像差选择性校正的成像系统光学/数字处理联合设计方法,并设计像散和彗差较大的三镜片成像系统,和像差校正较好的六镜片成像系统,通过比较证明三镜片成像系统的联合设计结果与六镜片系统的传统成像结果成像质量相当,从而也说明了本文提出的基于像差选择性校正的联合设计方法能够为光学系统的简化设计提供一定的理论依据.由于泽尼克圆多项式中没有光学色差的解析表达,因此本文在进行像差分析时并未根据色差的定义来分析其影响,这在今后进一步的研究中作中将予以考虑,以求获得更优的联合设计结果.

参考文献
[1] KINGSLAKE R. Optical system design[M]. ACADEMIC Press, 2012: 7-22.
[2] ZOU Mou-yan. Deconvolution and signal recovery[M]. Beijing: National Defend Industry Press, 2001: 3-5.
邹谋炎. 反卷积和信号复原[M]. 北京: 国防工业出版社, 2001: 3-5.
[3] ZHANG Ji-yan, HUANG Yuan-qing, XIONG Fei-bing. Iris acquiring optical system design with wavefront coding[J]. Acta Photonica Sinica, 2016, 45(10): 1022001.
张继艳, 黄元庆, 熊飞兵. 带有波前编码板的虹膜采集光学系统的设计[J]. 光子学报, 2016, 45(10): 1022001.
[4] SOMAYAJI M, BHAKTA R, CHRISTENSEN M. Experimental evidence of the theoretical spatial frequency response of cubic phase mask wavefront coding imaging systems[J]. Optics Express, 2012, 20(2): 1878-1895. DOI:10.1364/OE.20.001878
[5] MIRANI T, CHRISTENSEN M, DOUGLAS S, et al. Optimal co-design of computational imaging systems[C]. IEEE International Conference on Acoustics, Speech, and Signal Processing, 2005, 2: 597-600.
[6] MIRANI T, RAJAN D, CHRISTENSEN M, et al. Computational imaging systems:joint design and end-to-end optimality[J]. Applied Optics, 2008, 47: 86-103. DOI:10.1364/AO.47.000B86
[7] ROBINSON D, STORK G. Joint design of lens system and digital image processing[C]. Proceedings of the International Optical Design Conference., 2006, 6342: 63421G1-63421G10.
[8] ROBINSON D, STORK G. Joint digital-optical design of imaging systems for grayscale objects[C]. SPIE, 2008, 7100: 710011.
[9] GEARY M. Introduction to lens design:with practical Zemax examples[M]. Willman-Bell Press, 2002.
[10] AN Lian-sheng. Applied optics[M]. 3rd ed, Beijing: Beijing Institute of Technology Press, 2009: 150-155.
安连生. 应用光学[M]. 第三版, 北京: 北京理工大学出版社, 2009: 150-155.
[11] MALACARA D. Optical shop testing[M]. John Willey & Sons Press, 2006.
[12] CAI J, DONG B, OSHER S, et al. Image restoration:total variation, wavelet frames, and beyond[J]. Journal of the American Mathematical Society, 2012, 25(4): 1033-1089. DOI:10.1090/jams/2012-25-04
[13] BABACAN D, MOLINA R, KATSAGGELOS K. Parameter estimation in TV image restoration using variational distribution approximation[J]. IEEE Transactions on Image Processing, 2008, 17(23): 326-339.
[14] AMIZIC B, BABACAN D, MOLINA R, et al. Fast total variation image restoration with parameter estimation using bayesian inference[C] International Conference on Acoustics, Speech, and Signal Processing, 2010, Dallas, TX.
[15] RUIZ P, ZHOU X, MATEOS J, et al. Variational bayesian blind image deconvolution:a review[J]. Digital Signal Processing, 2015, 47: 116-127. DOI:10.1016/j.dsp.2015.04.012
[16] FOX C, ROBERTS S. A tutorial on variational Bayes[R]. Artificial Intelligence Review. Springer, 2011.
[17] LI Xiao-tong, CEN Zhao-feng. Geometrical optics, aberrations and optical design 3rd ed[M]. ZheJiang: Zhejiang University Press, 2014: 93-125.
李晓彤, 岑兆丰. 几何光学像差光学设计[M]. 浙江: 浙江大学出本社, 2014: 93-125.
[18] ISO 12233: 2014, Resolution and spatial frequency responses[DB/OL]. [2014-02]. https://www.iso.org/standard/59419.html.
[19] Sfrmat3: SFR evaluation for digital cameras and scanners[CP/OL]. [2015-05-12]. http://losburns.com/imaging/software/sfrmat3_post/index.html.