1 引言

时间延迟积分(TDI)技术通过对同一目标多次曝光,在与推扫速度严格匹配时可大幅提高成像系统的灵敏度和信噪比^[1],因此被广泛应用于高分辨率航天遥感、机器视觉、微光成像等领域。无论是电荷域、模拟域还是数字域TDI,TDI成像模式都存在像移失配,出现这种现象的原因是卫星在轨的姿态、轨道参数在实时变化,星下点的经纬度也不断改变,导致像平面的像移速度矢不固定;且在敏捷成像模式下像移变化现象更为显著^[2-4]。因此,在航天相机推摆扫成像过程中必须有精密的像移补偿机制,尤其在敏捷成像时还要求有很高的卫星姿态稳定度,否则像移失配现象将引起系统传递函数下降,甚至使图像发生畸变和模糊^[5]。

目前国内在相移补偿方面普遍采用机械式调偏流和电子式调像移速度大小相结合的方法,即根据像移速度矢的计算结果调整调偏流机构的运动来补偿像移速度矢的角度偏差,通过调整TDI行转移频率来补偿像移速度矢的大小偏差^[5-6]。

对地遥感卫星为了获取高分辨率图像,通常运行于太阳同步轨道。受稀薄大气阻力的影响,卫星轨道会逐渐衰变。为了使卫星维持在标称轨道,每隔一段时间就需要燃料推进。然而,当携带燃料耗尽时,卫星轨道只能逐渐衰退并最终坠入大气层烧毁,这个过程称为卫星轨道末期。由于太空中的空气甚为稀薄,卫星轨道末期比较漫长,短则几个月,长则几年^[7-8]。

在卫星轨道末期,成像系统通常会受自身参数限制而不能进行正常成像。尤其对于目前普遍采用的TDI航天遥感相机而言,轨道末期像移速度要求的行转移时间已超出CCD或CMOS(complementary metal oxide semiconductor)器件允许的工作范围,因此只能得到区域性混叠的欠采样图像。然而,轨道衰变是一个漫长的时期,此时若能进行一些非常规模式的成像,则对于提高卫星利用率意义重大。目前,研究较多的方法大多是图像复原算法。如文献[ 8]阐述了一种通过建立轨道末期成像模型进行图像复原的方法,但该算法的计算量大。文献[ 9]提出一种TDI行点扩展函数计算和相邻行合并复原的方法,并采用该方法对TDI模糊图像进行复原,虽然提升了计算速度,但却牺牲了图像的复原质量。上述图像复原算法均是基于欠采样混叠图像的后处理,图像信息已经丢失。

本课题组在前期工作中开展了大量逐行累加的数字域TDI通用算法的研究^[10-12]。为了解决上述问题,本文利用数字域TDI操作灵活的优势,在前期数字域TDI通用算法研究的基础上进行进一步改进,提出一种适用于沿扫描方向任意像移大小的数字TDI算法。该算法结合传感器的工作特点,利用图像插值和配准像素对位累加思想,构建电子补偿和图像补偿相结合的数字积分模型,在标称轨道和轨道末期均可直接获取无区域性混叠图像。相对于现有TDI相机在轨道末期有限的像移可补偿范围,本文算法可适用于标称轨道及轨道末期沿扫描方向广域像移变化的情况。

2 数字域TDI通用算法工作原理

数字域TDI通用算法工作原理如图1所示。在推扫成像过程中,目标像素将自上而下依次扫过各级TDI像元,而每级像元对目标的凝视时间为一个行周期T。在第一个行周期,利用第一级像元对目标进行曝光。在第二个行周期,目标恰好移动到第二级像元的成像区域,利用第二级像元对同一目标进行曝光,并在前级转入信号的基础上继续收集该次曝光信号,当积分结束时,累加信号再转入第三级像元。依次类推,直到第M个行周期,第M级像元信号与前面共M-1级像元对同一目标的信号累加后,读出积分信号。显然,M级积分信号是单级信号的M倍,若响应率为R_e,则其输出信号为

S=EAReTM,(1)

式中:E为输入光照度;A为传感器的光敏元面积;T为积分时间。

根据上述描述,数字域TDI通用算法除了将信号累加操作置于数字域外,其基本原理与TDI CCD、模拟域TDI并无区别,均采用同列像素逐行累加策略^[13-16]。即在推扫成像过程中,要求每个行周期目标恰好移动一个像元尺寸的距离,否则在固定的逐行累加操作下,将会引起图像畸变和模糊。

3 广域像移变化的数字TDI算法

根据光学系统设计,遥感相机分辨率随着轨道高度的降低而提高。地面像元分辨率D_GS求解公式为

DGS=Ha/f,(2)

式中:a为CMOS像元尺寸;f为光学系统焦距;H为轨道高度。在保证逐行叠加的数字TDI算法能清晰成像的前提下,随着地面分辨率的提高,要求的行转移时间T_need会显著缩短,T_need满足

vgTneed=DGS,(3)

图 1. 数字域TDI通用算法工作原理

Fig. 1. Working principle of digital TDI general algorithm

下载图片 查看所有图片

式中:v_g为卫星地速,与卫星轨道速度具有相同的角速率。根据万有引力提供向心力理论,可得v_g与H的关系为

vg=[R/(R+H)]μ/(R+H),(4)

式中:R为地球半径;μ为地球引力常数。

将(4)式代入(3)式,得到T_need的计算式为

Tneed=(Ha/f)/{[R/(R+H)]μ/(R+H)}。(5)

由(5)式可知,T_need∝H^5/2,因此要求的T_need会随着H降低而急剧下降。然而,任何传感器的读出频率都有一个上限值,当达到其上限读出频率时,行周期将不再继续降低,因此在轨道末期受传感器读出频率的限制,像移速度失配将会越来越严重。

为解决上述问题,本文提出一种适用于广域像移变化条件下的数字TDI算法。算法流程如图2所示。

首先,对拍照时刻地物点在像面上的像移速度进行时间积分,得到第i+1帧图像的像移量计算式为

SP(i)=∫0iTfrmvpdt,(6)

式中:T_frm为帧周期。为使数字域像素可以配准叠加,像移量对应的像素个数应为

Δn(i)=SP(i)/a。(7)

图 2. 广域像移变化的数字TDI算法流程图

Fig. 2. Flowchart of digital TDI algorithm for wide-range image motion variation

下载图片 查看所有图片

将第i帧像素矩阵用P_i表示,其第m行第n列像素表示为P_i(m,n)。为了实现对应地物的亚像素级累加,首先需要对像移个数为非整数的图像进行插值:

P'i+1(m,n)=α(i)P'i[m+Δβ(i),n]+[1-α(i)]Pi+1(m,n),m∈R[1,nlines-Δβ(i)]α(i)Pi+1(m-1,n)+[1-α(i)]Pi+1(m,n),m∈R[nlines-Δβ(i)+1,nlines],(8)

其中

α(i)=rem[Δn(i)],(9)β(i)=int[Δn(i)],(10)Δβ(i)=β(i)-β(i-1),(11)

式中:rem(·)为取余函数;int(·)为取整函数;n_lines为图像行数。将插值后的结果与原有像素矩阵叠加,得到TDI后第i帧像素矩阵P″_i+1的第m行第n列像素为

P″i+1(m,n)=P″i[m+Δβ(i),n]+P'i+1(m,n),m∈R[1,nlines-Δβ(i)]P'i+1(m,n),m∈R[nlines-Δβ(i)+1,nlines]。(12)

结合(8)式与(12)式得到插值和叠加操作的简化公式为

P″i+1(m,n)=[1+α(i)]P″i[m+Δβ(i),n]+[1-α(i)]Pi+1(m,n),m∈R[1,nlines-Δβ(i)]α(i)Pi+1(m-1,n)+[1-α(i)]Pi+1(m,n),m∈R[nlines-Δβ(i)+1,nlines]。(13)

定义像移量变化率ε=β(i+j)/β(i)。以图3为例,得到各帧周期输出的行数据如下。其中当ε>1时,为避免丢行,需在每帧叠加完的第一行输出前,从前一叠加帧输出缺失行:

Nsignal(i)=P″i(1),i=1,2,3,4P″i-1(2)P″i(1),i=5,6。(14)

以此类推,得到第i帧输出信号的表达式为

Nsignal(i)=P″i-1[β(i-1)-β(i-2)]P″i-1[β(i-1)-β(i-2)-1]︙P″i-1(2)P″i(1)。 (15)

以上插值方法完成了像素位移个数向下取整的上溯插值。该方法操作简便,但是会存在某些行数据叠加次数小于积分级数M的问题。因此,为了避免级数差异引起的失真,也为了使插值数据更接近真实值,当ε>1时,将插值改为向上追溯4个邻近帧插值缺失的1行数据,即

P'i+j+1(m,n)=∑k=03Pi+j-k[1+∑x=1k+1β(i+j-x),n]/4,M≥4Pi+j-k[1+∑x=2Mβ(i+j-x),n]/M,M<4。(16)

因为相机在数字域TDI模式下,每个行周期将第一行累加数据输出作为M级积分值,所以第i个行周期第n个输出像素为

OM,i(j,n)=∑j=0β(i-M+1){Pi-M+1[Mβ(i-M+1)-j,n]+∑x=1M-1P'i-M+1+x[(M-x)β(i-M+1)-j,n]}。(17)

图 3. 上溯插值示意图

Fig. 3. Retrospective interpolation diagram

下载图片 查看所有图片

算法设置一个最小行周期阈值(对于CMOS传感器,该阈值为最大读出频率时的帧周期),当通过姿轨参数和(5)式计算得到的所需行周期大于该阈值时,通过电子补偿像移;当所需行周期小于该阈值时,将行转移周期锁定为该阈值,然后利用图像法补偿像移。

在所提数字TDI算法中,电子像移补偿模型为

nlines=(Tneed-Tnicheng)/(trb+npixels/fs),(18)Tnicheng=Tneed-{(trb+npixels/fs)×rem[Tneed/(trb+npixels/fs)]},(19)

式中:t_rb为行空白时间,CMOS传感器需利用此时间进行相关双采样;f_s为像素读出频率;n_pixels为每行像元数;T_nicheng为逆程时间。当像移速度不匹配且行周期大于阈值时,可以通过调整开窗行数(粗调)或逆程时间(精调)来改变行转移速度,从而改善像移失配引起的成像模糊。

4 实验结果与分析

首先选取合适的样本图像,然后用计算机模拟数字TDI成像,继而利用像移传递函数方法和互相关相似性测度方法对仿真图像进行分析比较。

1) 像移传递函数方法

假设将成像系统看作一个线性不变系统,那么景物经成像系统转化为积分后的数字图像的过程,可视为景物经过成像系统的传递过程。传递的频率不变,变化量是对比度和相位。随着对比度下降和相位推移,某一频率处的对比度会变为零。此时该频率的光强分布已经不存在亮度变化,即该频率被截止。一般而言,高频部分反映物体的细节传递情况,中频部分反映物体的层次传递情况,低频部分反映物体的轮廓传递情况,而表明各种频率传递情况的则是调制传递函数(MTF)。因此可以将像移对应的MTF用于评价像移匹配误差对成像质量的影响。

2) 互相关相似性测度方法

若原始图像表示为T(m,n),仿真图像表示为S(m,n),数字TDI CMOS成像系统扫描图像对应像素表示为S^i,j(m,n),其中(i,j)代表对应像素的位置,M×N表示选取图像的大小。则其归一化互相关相似性测度为

σ=∑m=1M∑n=1N[Si,j(m,n)×T(m,n)]/{∑m=1M∑n=1N[Si,j(m,n)]2×∑m=1M∑n=1N[T(m,n)]2}。(20)

3) 仿真分析

选取某高分辨率原始遥感图像(如图4所示)作为假想景物样本。由像移速度匹配残余误差对MTF的影响可知,积分级数越大,MTF下降得越严重,图像越模糊,因此为了更好地验证算法性能,实验中应尽量选取高积分级数。参考IT-EN-06144-00-R型航天TDI CCD可设置积分级数为8,16,32,48,64,96,最终选取96级积分来验证所提算法,若96级积分仍可清晰成像,则其他低积分级数下也可满足要求。利用建立的广域像移大小数字TDI算法模型,仿真生成不同像移失配率情况下的96级数字域TDI图像,如图5(a)、图6(a)所示。为了对比分析,根据TDI通用算法工作原理,模拟同样条件下其扫描所得图像,结果如图5(b)、图6(b)所示。

图 4. 原始遥感图像

Fig. 4. Original remote sensing image

下载图片 查看所有图片

图 5. 像移失配率为0.5%时两种算法获取的图像。(a) 本文算法;(b) TDI通用算法

Fig. 5. Images obtained by two algorithms when image mismatch rate is 0.5%. (a) Proposed algorithm; (b) TDI general algorithm

下载图片 查看所有图片

在本实验中,影响最终图像传递函数的因素为目标自身对比度特性以及像移失配,图像传递函数f_MT-IMG、目标传递函数f_MT-OBJ和像移传递函数f_MT-MAT三者之间满足关系f_MT-IMG=f_MT-OBJf_MT-MAT。由于实验采用同一个样本图像作为基准,因此在各幅图像中由目标自身对比度特性引起的传递函数属于固定传递函数,此部分传递函数可通过由无像移失配(f_MT-MAT=1)的基准样本图像计算的图像传递函数得到。采用刃边法对图像传递函数进行计算。分别在样本图像、两种算法获取的图像中选取同一区域的刃边图像,首先采用刃边法计算上述各图像中选取区域的图像传递函数,然后以样本图像的图像传递函数作为其他图像的目标传递函数,并利用3种传递函数关系式求取像移传递函数,算法越优,像移传递函数值越大。归一化互相关相似性测度中样本图像为S(m,n),将其与两种算法在不同像移失配率下获取的图像进行比较,得到计算结果,测度越高,说明图像失真越小。采用像移传递函数和互相关相似性测度方法对图像进行分析,结果见表1。

图 6. 像移失配率为2%时两种算法获取的图像。(a)本文算法;(b) TDI通用算法

Fig. 6. Images obtained by two algorithms when image mismatch rate is 2%. (a) Proposed algorithm; (b) TDI general algorithm

下载图片 查看所有图片

观察图5、图6以及表1数据,可以发现:当96级积分的总像移量远小于像元尺寸时,本文算法与逐行累加的TDI通用算法的成像结果差不多,但因为本文算法可以进行相邻像素的加权修正,因此其成像质量略优于TDI通用算法;当96级积分的总像移量大于像元尺寸时,总像移量超出电子式像移可补偿的范围,因此TDI通用算法扫描图像发生严重混叠,成像质量急剧下降,像移传递函数和互相关性测度分别只有0.8201和0.8009;而本文算法利用匹配像素叠加,可以有效地控制成像质量不发生急剧恶化,图像像移传递函数和互相关性测度分别达到0.9386和0.9109,相较于TDI通用算法这两项指标至少提高0.11,但因为本文算法上溯插值的非整数像移量的像素由多像素插值得到,或者某些行像素由多帧图像获得,因此像移失配率为2.0%的成像效果相比于0.5%时会略微差一点。

5 结论

为了弥补传统电子式像移补偿机制存在的不足,以轨道末期遥感相机成像为例展开研究,提出了适用于广域像移大小变化的数字TDI算法。该算法以数字TDI算法为基础,同时将像移补偿与数字TDI算法结合在一起,扩大了传统电子式像移补偿的可补偿范围和TDI遥感相机可清晰成像的工作范围。设计了不同像移失配率的数字TDI算法的验证实验,实验结果表明:所提设计算法在不同像移失配率下均可获得清晰的图像,正确可行;尤其在2%的大像移失配率条件下仍可保持较高的传递函数值,与通用TDI算法的成像结果对比鲜明。