一种针对星载相机的调制传递函数测算方法

摘要

一种针对星载相机的调制传递函数测算方法，涉及空间光学技术领域，解决了现有的点光源测算调制传递函数方法存在的需要用参数模型近似和拟合的问题。首先在卫星过境前铺设多个点光源靶标，然后在捕获的图像中截取对应点光源的图像，去除背景亮度数值后，根据成像数学模型，采用正则化求解方法直接求解二维点扩散函数，最后由点扩散函数计算调制传递函数。本发明的测算方法实施简单，数据处理过程简便、精确，对于星载相机的成像质量评估意义重大，有着广泛的应用前景。

说明书

技术领域

本发明涉及空间光学技术领域，具体涉及星载相机的点扩散函数及调制传递函数测算方法。

背景技术

星载相机的调制传递函数是评估成像系统质量的重要指标，可以很大程度上反应图像的清晰程度，故而星载相机都对调制传递函数有指标要求，尤其关注调制传递函数在奈奎斯特频率处的数值。将相机成像系统作为一个线性移不变系统，点扩散函数描述了系统对点光源激励的响应，可以表征系统特性。作为点扩散函数傅里叶变换的归一化模值，调制传递函数则描述了该系统对于不同空间频率的响应，它受到空间相机光学系统、探测器、相对运动、电子学处理系统等众多因素的影响。虽然相机入轨之前会在实验室对相机的调制传递函数进行测试，但是由于入轨之后光机结构的变化、成像环境等因素的影响，相机的调制传递函数会发生一定的变化，故而需要对相机在轨时的调制传递函数进行重新评测。

目前针对星载相机的调制传递函数测量方法，常用的有点光源法和倾斜刃边法。倾斜刃边法通过分别测算沿轨方向和垂轨方向的线扩散函数来计算对应的调制传递函数，虽然对校直问题不敏感，但实施过程也受刃边的倾斜角度、对比度等因素影响。而点光源法是一种简单直接的测量方法，但也存在以下技术问题：

现有方法在对点光源的成像结果进行处理时，需要对点光源图像进行拟合，分析峰值位置以及点扩散函数形状，其数据处理流程如图1所示。数据拟合时的采用的参数模型为高斯函数，即假设点扩散函数呈二维高斯形状。而实际的成像系统点扩散函数通常都不能够用高斯函数模型描述，故而近似为高斯函数模型存在着较大的误差，影响调制传递函数的测量精度。同时，由于随机因素的存在，实际应用中的星载相机的点扩散函数很难用一个参数模型来表达，而且当未知参数较多、参数模型复杂时，还存在着拟合过程不稳定、拟合结果误差大的问题。

发明内容

本发明为解决现有点光源测算调制传递函数方法存在的需要用参数模型近似和拟合的问题，提供了一种针对星载相机的调制传递函数测量方法。

一种针对星载相机的调制传递函数测量方法，该方法由以下步骤实现：

步骤一、在卫星过境之前，在具有相同背景的场景中铺设规格相同的凸镜，由所述凸镜反射的太阳光形成多个点光源靶标；

步骤二、根据星载相机捕获的多个点光源靶标的图像，截取对对应各点光源的图像数据；

步骤三、估算步骤一中场景的背景亮度，并采用步骤二中获得的图像数据减去背景亮度的数值，获得去除背景亮度后的图像数据；

步骤四、根据成像数学模型，对步骤三中获得的去除背景亮度后的图像数据采用正则化求解方法计算点扩散函数；并根据调制传递函数与点扩散函数的数学关系，计算调制传递函数。

本发明的有益效果：本发明较好地避免了现有基于点光源的测算方法需要估计点扩散函数参数模型的问题，避免了使用二维高斯函数等模型来拟合点扩散函数的操作，而且通过合理的约束条件提高了点扩散函数的求解准确度，从而提高了调制传递函数的测算精度。在正则化求解过程中，由于需要求解的目标尺寸不大，故而求解过程快速方便，具有很强的工程应用价值。

附图说明

图1为现有的点光源方法测试调制传递函数的数据处理流程。

图2为本发明针对星载相机的调制传递函数测算方法实施流程。

图3为本发明所述的针对星载相机的调制传递函数测量方法中捕获的点光源靶标的图像。

图4为本发明所述的针对星载相机的调制传递函数测量方法中采用正则化求解方法计算得到的二维点扩散函数。

图5为本发明所述的针对星载相机的调制传递函数测量方法中不同空间频率下的调制传递函数。

具体实施方式

具体实施方式一、结合图2至图5说明本实施方式，一种针对星载相机的调制传递函数测量方法，该方法由以下步骤实现：

步骤一、卫星过境之前，在具有一致背景的场景中铺设一组规格相同的凸镜，两两间隔不低于六个地面像元间隔，由凸镜反射太阳光形成多个点光源靶标。所述的地面像元间隔是已知的，但是由于卫星轨道形状的缘故，地面像元间隔有可能在一个确定的范围内变化的。

步骤二、从星载相机捕获的图像中，对应各点光源，截取相同尺寸的图像数据。

步骤三、估算场景的背景亮度，并将背景亮度数值从步骤二的图像数据中减去，得到扣除背景亮度后的图像。

步骤四、对步骤三中得到的图像，根据成像数学模型，采用正则化求解方法直接计算点扩散函数。

步骤五、根据调制传递函数与点扩散函数的数学关系，计算调制传递函数。

其中，步骤三中的背景亮度数值可以通过计算背景区域的像素平均值得到。执行步骤一至步骤三后，得到了不包含背景亮度的P幅图像z_k(k＝1,2,...,P)，每幅图像的尺寸均相同。

步骤四中，对于星载相机，在不考虑畸变时，其成像过程可以用如下公式描述：

z(i,j)＝D[h(x,y)*o(x,y)]+n(i,j)(1)

其中，z(i,j)表示捕获的数字图像，D表示抽样算子，h(x,y)为点扩散函数，o(x,y)为地面场景，n(i,j)为加性噪声，*表示卷积运算。

当地面场景为点光源时，

o(x,y)＝δ(x-ξ,y-ζ)(2)

其中，δ为狄拉克冲击函数，(ξ,ζ)为点光源所处的位置坐标。

在本实施方式中，场景中包含P个点光源目标，记为o_k，它们的位置坐标依次为(ξ_k,ζ_k)(k＝1,2,...,P)，即满足：

o_k(x,y)＝δ(x-ξ_k,y-ζ_k)(3)

将公式3代入到公式1中，即有：

$>\left(\begin{array}{c}{z}_k(i,j)=D[h(x,y)*o_k(x,y)]+n_k(i,j)\\ =D[h(x,y)*\delta (x-{\xi }_k,y-{\zeta }_k)]+n_k(i,j)\\ =D\left[h(x-{\xi }_k,y-{\zeta }_k)\right]+n_k(i,j)\end{array}\right)---\left(4\right)$>

在公式4中，要求解h，采用一种正则化求解方式，即最小化如下的代价函数E：

$>E(h,\{{\xi }_k,{\zeta }_k\left\}\right)=\stackrel{P}{\Sigma }{\left|\right|D\left[h(x-{\xi }_k,y-{\zeta }_k)\right]-z_k\left|\right|}^2+\lambda \int |▿h|---\left(5\right)$>

公式5中，第一项为保真项，第二项为平滑项，||||表示l₂范数，▽是哈密顿算子，λ为正则化系数，用于平衡保真项与平滑项之间的折衷。

当h固定时，函数E关于{ξ_k,ζ_k}(k＝1,2,...,P)是凸函数；当{ξ_k,ζ_k}固定时，函数E关于是凸函数，故而可以使用交替最小化的方法求解公式5对应的最小化问题，即第t次迭代过程表示为：

$>\left(\begin{array}{c}{\{{\xi }_k,{\zeta }_k\}}^{\left(t\right)}=\underset{\{{\xi }_k,{\zeta }_k\}}{\arg \min }E(h^{(t-1)},\{{\xi }_k,{\zeta }_k\left\}\right)\\ h^{\left(t\right)}=\underset{h}{\arg \min E}(h,\{{{\xi }_k,{\zeta }_k\}}^{\left(t\right)})\end{array}\right)---\left(6\right)$>

在求解过程中，可以借助Matlab软件最优化工具箱(OT)中的fminsearch函数进行搜索最小值点。

步骤五中，点扩散函数与调制传递函数满足的数学关系为：

MTF＝|FT[h]| (7)

其中，MTF为调制传递函数，FT为傅里叶变换算子，||表示取模。也可以借助Matlab软件的图像处理工具箱(IPT)中的psf2oft函数，将其结果取模即可。

具体实施方式二、结合图3至图5说明本实施方式，本实施方式为将具体实施方式一所述的一种针对星载相机的调制传递函数测量方法的应用到QuickBirdII铺设的凸镜靶标成像实验中，QuickBird的地面像元间隔，可以为0.6m，也可以为0.61-0.72m。所述地面像元间隔取决于像元尺寸、光学系统焦距、卫星轨道高度。

首先从捕获的图像中截取对应各点光源的图像，结合图3，使用图3中虚线区域的背景数据计算均值作为整个区域的背景亮度，并将其从图3的各个点光源图像中减去。

令公式5中λ＝0.5，利用公式6进行重复迭代，当相邻两次迭代之间h的变化小于10^-6时，停止计算，得到点扩散函数如图4所示。

使用图4所示的点扩散函数，依据公式7计算其调制传递函数，得到的对于不同空间频率的响应如图5所示，其中空间频率按照图像采样频率进行了归一化。在奈奎斯特频率处，即归一化频率为0.5时，垂轨方向的调制传递函数数值为0.183，沿轨方向的调制传递函数数值为0.171。

通过实验，得到的奈奎斯特频率处调制传递函数的数值与使用倾斜刃边法计算结果偏差小于3％，表明该方法是准确有效的。