光学成像系统有限维本征频域分析方法
【技术领域】
[0001]本发明属于光学系统信息光学与成像分析技术领域，涉及一种光学成像系统有限 维本征频域分析方法。
【背景技术】
[0002] 根据傅里叶光学的分析方法，光学系统成像空域成像的卷积过程可由傅里叶频域 中的频谱乘积关系描述。傅里叶光学频域成像分析方法的基本原理是将物、像及光学系统 的点扩散函数(PSF)作傅里叶变换分别得到物频谱、像频谱及光学系统的光学传递函数 (0TF)，像频谱等于物频谱与0TF的乘积。但傅里叶光学的频域成像分析方法在理论与数值 计算中均存在不完善的方面。
[0003] 在理论方面，根据傅里叶变换的性质可知，有界函数的频谱是无界函数，而光学系 统的像面在空域是有界的，所以，从傅里叶变换的角度看，像频谱应该是无界函数。但由于 光学系统的口径不可能是无限大的，有限口径光学系统的0TF必然存在截止频率，而像频谱 等于物频谱与0TF的乘积，所以，从光学系统成像的角度看，像频谱应该是有界函数。傅里叶 变换的性质得出的有界像频谱结果与0TF性质得出的无界像频谱果是互相矛盾的，这是傅 里叶频域分析方法存在的问题。
[0004]在数值计算方面，由于傅里叶变换的积分变换核是复指数函数(余弦函数和正弦 函数的组合形式），而复指数函数是具有无穷数据量的无界函数，无法用计算机进行数值计 算。目前普遍采用快速傅里叶变换方法(FFT)来实现傅里叶变换的理论公式的数值计算，而 FFT的计算结果与傅里叶变换公式的计算结果不一致，因而无法获得精确的数值计算结果。
[0005] 由于傅里叶变换在理论与数值计算方面存在上述问题，研究一种精确的光学成像 系统有限维频域分析方法具有重要的理论意义与实用价值。
【发明内容】
[0006] 本发明的目的是提供一种光学成像系统有限维本征频域分析方法，该方法克服了 傅里叶光学的频域成像分析方法在理论与数值计算中存在的缺点，给出了物、像本征频谱 向量的求解方法，实现了物本征频谱向量乘以光学系统特征值向量等于像本征频谱向量的 光学成像系统有限维本征频域分析方法。
[0007] 本发明的目的是通过以下技术方案实现的：
[0008] -种光学成像系统有限维本征频域分析方法，包括如下步骤：
[0009] 第一步、二维光强传输矩阵的构建：
[0010]根据光学系统的PSF矩阵和一维物向量求解光学系统的二维光强传输矩阵；
[0011] 第二步、二维光强传输矩阵的特征值向量和本征函数向量组的求解：
[0012] 利用QR分解法求解二维光强传输矩阵的特征值向量，利用幂法求解二维光强传输 矩阵的本征函数向量组；
[0013] 第三步、物、像的有限维本征频谱向量计算：
[0014] 通过物向量、像向量与本征函数向量组的矩阵运算求解物、像的有限维本征频谱 向量。
[0015] 本发明通过建立光强传输矩阵实现了用矩阵运算方法取代了卷积的数值计算，利 用光强传输矩阵的本征函数向量组与特征值向量取代了傅里叶分析方法中的复指数函数 积分变换核与0TF，进而利用本征函数向量组求解出物、像空域向量的本征频谱向量，实现 了物本征频谱向量乘以特征值向量等于像本征频谱向量的光学成像系统有限维本征频域 分析方法。
【附图说明】
[0016] 图1为光学系统空域成像原理图；
[0017] 图2为将像面绕光轴旋转180度的成像原理图；
[0018] 图3为二维物矩阵转化为一维物向量的一维化过程原理图；
[0019] 图4为二维像矩阵转化为一维像向量的一维化过程原理图；
[0020] 图5为一维物向量、一维像向量与二维光强传输矩阵的关系图；
[00211图6为光学系统的光路图；
[0022]图7为光学系统的PSF曲线；
[0023]图8为光学系统的MTF曲线；
[0024] 图9为光学系统的第0~8阶本征函数向量；
[0025] 图10为光学系统的特征值向量；
[0026] 图11为有限维余弦函数物向量；
[0027] 图12为该物的有限维本征频谱向量；
[0028]图13为该物的FFT频谱向量的模；
[0029] 图14为有限维余弦函数的像向量
[0030] 图15为像的有限维本征频谱向量；
[0031 ]图16为像的FFT频谱向量的模。
【具体实施方式】
[0032] 下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本 发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖 在本发明的保护范围中。
[0033] 本发明提供了一种光学成像系统有限维本征频域分析方法，具体上述步骤如下：
[0034] 第一步:二维光强传输矩阵的构建。
[0035] 1)二维空域成像情况下的二维光强传输矩阵构建
[0036] 此部分内容与CN104574315A及CN104360481A的【具体实施方式】中的第一、二步所述 原理相似，但物理量的意义及处理方法是有细微差别的。
[0037] 如图1所示，二维物面1上的某一物点2的光强经光学系统3在二维像面4上成像为 一个弥散斑5,弥散斑5的光强分布可由以相应像点为中心的点扩散函数(PSF)6的数据精确 地描述。在未加旋转光路的光学系统中，像面相对于物面是倒立的，即像面相对于物面旋转 了 180度角。为了便于分析和处理，可将像面以光轴为中心旋转180度，如图2所示，这种处理 方法不改变光学系统的成像规律，在实际光学系统中，可通过加旋转光路实现。如图2所示， 二维物面1上各点的光强数据可由图3中的二维物矩阵7表示，同样，二维像面4上各点的光 强数据可由图4中的二维像矩阵9表示。二维物矩阵的某个元素 am,n到二维像矩阵的某个元 素 bt,^光强传输系数ckn,t,A、须是四维参数，则由dm, n,t,w构成的矩阵必然是四维矩阵，而 四维矩阵是不便于分析和处理的，也不便于绘图表示。而对于一维物向量中的元素^到一 维像向量中的元素比的光强传输系数PW是二维参数，由PW构成的光强传输矩阵是二维矩 阵，且二维矩阵便于分析、处理，也便于绘图表示。
[0038]如图3所示，物的一维化处理就是将二维物矩阵7各行元素中的每一行依次首尾相 接依次排列为一个行向量。如图4所示，像的一维化处理就是将二维像矩阵9各行元素中的 每一行依次首尾相接依次排列为一个行向量。对物向量与像向量中的元素重新按前后顺序 编号后，便形成了物向量8与一维像向量10。
[0039]如图3、图4所示，在对物矩阵、像矩阵进行一维化后，二维物矩阵的某个元素 am,n转 变为一维物向量中的相应元素 ai，二维像矩阵的某个元素 bt,w转变为一维像向量中的相应 元素氏，应满足如下关系：
[0040]
[0041]
[0042] 其中，m、n为二维物矩阵7中元素的序号，i为一维物向量8中元素的序号，t、w为二 维像矩阵9中元素的序号，j为一维像向量10中元素的序号，M、N为二维物矩阵7的总行数、总 列数，T、W为二维像矩阵9的总行数、总列数，且1、」、!11、11、^1、11'、1均为正整数。通常情况 下，物的数据量与像的数据量相等，即Μ与T是相等的、N与W是相等的。
[0043] 二维物面上的每一个物点发出的光被光学系统会聚在像面上而形成一个覆盖多 个像点的弥散斑，弥散斑中心位于与物点相对应的像点（即图2中与物点二维位置相同的像 点），该弥散斑的光强分布可由以该像点为中心的PSF精确描述。需要说明的是，在线性空变 光学系统中以每个像点为中心的PSF均是不同的，而线性空不变光学系统中以每个像点为 中心的PSF均是相同的。
[0044] 如式(3)所示，PSF数据可由矩阵Η表示，其中，k、l为二维PSF矩阵中某个元素的序 号，K、L为二维PSF矩阵的总行数、总列数，k、l、K、L均为正整数。
[0045]
(3)fl
[0046] 由二维物矩阵的某个元素 am,n在像面形成的像矩阵Γ'η*: (4)。
[0047]
[0048] 则二维物矩阵的某个元素 am,n传输到二维像矩阵的某个元素心^的光强传输系数 dm, n, t, w^J ：
[0049]
(5)0
[0050] 由式（1)、（2)、（5)可知，一维物向量中的元素€4到像向量中的元素氏的光强传输系 数PU(即二维光强传输矩阵11中的元素)是相等的，BP:
[0051 ]
(6)。
[0052] 且参数^1\1、1(、1^满足如下关系：
[0053]
(7)。
[0054] 因此，只要记录好二维物矩阵、像矩阵各元素间的四维光强传输系数dm,n,t, w，便可 获得一维物向量各元素到一维像向量各元素的二维光强传输系数PU，该系数即图5所示的 二维光强传输矩阵11中的相应元素。按照这种方法，依次求出一维物向量各元素到一维像 向量各元素的光强传输系数，便可得到二维光强传输矩阵11。由于二维光强传输矩阵中的 传输系数PU均为PSF中的数据。
[0055] 求解出二维光强传输矩阵11后，一维物向量8乘以二维光强传输矩阵11等于一维 像向量10的关系便确立了，如图5所示。这种乘积运算关系可写成式(8)的形式：
[0056] Α·Ρ = Β (8)；
[0057] 其中，Α为一维物向量，意义等同于图5中的一维物向量8;Ρ为二维光强传输矩阵， 意义等同于图5中的二维光强传输矩阵11 ;Β为一维像向量，意义等同于图5中的像向量10。
[0058] 2)-维空域成像情况下的二维光强传输矩阵构建
[0059] -维空域成像情况下，物与像均是一维的，所以，物、像均可直接由物向量8、像向 量10表示，无需进行一维化处理。PSF也可由一维向量形式表示。此时，直接根据式(3)~(7) 即可求出一维空域成像情况下二维光强传输矩阵11。
[0060]第二步、二维光强传输矩阵的特征值向量和本征函数向量组的求解。
[0