推扫式多光谱遥感相机动态范围拓展方法

孙 武，韩诚山，吕恒毅,2，薛旭成，胡长虹

(1.中国科学院 长春光学精密机械与物理研究所，吉林 长春 130033；2.中国科学院大学，北京 100049)

1 引 言

国内外现有高分辨率航天遥感相机多为推扫式航天遥感相机[1-2]，其动态范围一般不大于80 dB，无法满足实际场景的远大于100 dB动态范围的要求。目前，提高相机动态范围的方法主要有4类：(1)针对同一场景进行多次不同量的曝光，获取多幅不同曝光量的图像，再将它们融合为一幅新的高动态范围的图像，这是较为经典的办法。这类方法于1997年由Paul E.D., Jitendra M.首次提出[3]，2003年至2005年，Robertson M.A.，Szeliski R.，Goshtasby A.A.等人针对该技术做了进一步的探索和完善[4-6]。该方法主要针对同一静态场景进行动态范围拓展，由于推扫式航天遥感相机难以在较短的时间内对同一场景进行多次曝光，该方法不适用于该场合。(2)薛旭成等人提出在推扫式航天遥感相机上应用双排TDI CCD(Time Delay Integration-Charge Coupled Device)，对同一场景成两幅或多幅积分级数不同的图像，再通过图像融合算法获取高动态范围遥感图像[7]。周继权[8]等人提出基于相机阵列的高动态范围图像合成方法也是同样的原理，可以应用于动态变化的场景。该类方法可以有效提高相机的动态范围，但是由于成本、功耗、重量、光学设计等问题，在航天领域的工程应用价值比较低。(3)吕涛、吕伟振等人提出利用数字微镜获取高动态范围图像[9-10]，但是考虑到技术和产品的性能[11-13]、成熟度、可靠性等因素，目前难以应用到航天相机领域。(4)Wu Sun等人提出采用面阵CMOS，并应用数字TDI技术，设置2个高低不同的数字积分级数，针对同一场景可以获取2幅不同曝光量的图像，再通过融合算法获取高动态范围遥感图像[14]。该技术简单易行，可在不改变现有光机设计、不增加成本的前提下有效提高动态范围，但是该技术只能在新研制的遥感相机中应用，对现有遥感相机需要需要进行硬件改造。

本文提出了一种特殊的TDI-CCD积分级数设置方法：通过设置较低的P谱段积分级数，配合较高的多光谱谱段积分级数，提高动态范围。具体做法是：通过多谱段图像融合获取一幅低积分级数全色图像；该图像经过上采样后，与高积分级数、高分辨率全色图像融合，以获得一幅高动态范围、高分辨率的全色遥感图像。该方法可以作为一种特殊的高动态范围成像模式，广泛应用于推扫式多光谱航天遥感相机中。

2 积分级数设置及动态范围拓展方法

为实现推扫式多光谱遥感相机的动态范围拓展，可结合TDI-CCD的像元尺寸、谱段分布等特点，将P谱段积分级数调高、B谱段积分级数调低，以分别获取长、短曝光时间的图像，再通过图像融合获取高动态范围的图像。

2.1 推扫式多光谱遥感相机

推扫式多光谱遥感相机常采用TDI-CCD作为探测器，将光信号转换为数字图像信号。TDI-CCD中分为全色谱段(记作P谱段)和多光谱谱段(统称为B谱段，分别记作B1, B2, B3, B4,…,Bn谱段)。一般而言，P谱段地面像元分辨率为δp，B谱段地面像元分辨率为δb。表1 列出了几颗常见商业高分辨率遥感卫星的探测器像元尺寸。一般而言，δb是δp的k倍，k为整数。B谱段图像中第i行第j列的像元Bi，j对应k×k个全色图像的像元。与同一光学系统配合使用时，像元尺寸与地面像元分辨率成正比，Bx(x≤n)谱段图像像元与P谱段图像像元的对应关系如图1 所示。

表1 商业高分辨率遥感卫星地面像元分辨率Tab.1 Ground pixel resolution of commercial high-resolution remote sensing satellites(nadir)

图1 Bx谱段图像与P谱段图像的像元对应关系 Fig.1 Corresponding relationship of pixels betweene Bx band image and P band image

各个谱段可以分别获取一幅图像，由于P谱段与B谱段地面像元尺寸大小不同，一般来说，P谱段可以获取高分辨率的图像，B1,B2,B3,B4,…,Bn谱段可以分别获取各自波长范围的较低分辨率的图像。通过全色图像与多光谱图像的融合，可以获取高分辨率彩色或假彩色图像。

2.2 特殊积分级数设置方法及相应动态范围

常规的遥感相机在设置积分级数时，仅考虑信噪比这一参数，一般选择满足TDI CCD不饱和的最大积分级数[15]。

为了获取长短曝光的遥感图像，进一步达到动态范围拓展的目的，本文提出一种新的积分级数设置方法，将各谱段积分级数设置偏离原有最佳数值：P谱段积分级数从原值MP调高μ1倍至NMP；B1,B2,B3,B4,…,Bn各谱段积分级数从原值MB1，MB2，MB3，MB4，…，MBn，调低μ2倍至NMB1，NMB2，NMB3，NMB4，…,NMBn。

P谱段直接成像，得到一幅全色图像，记作ΓP。B谱段多幅图像融合成一幅全色灰度图像ΓB，插值后记作ΓBIN。在完成积分级数调整，提高ΓP的曝光量，降低ΓB的曝光量之后，将ΓP和ΓBIN融合在一起，可以获取一幅高动态范围的全色图像Γ。

新旧积分级数设置的关系如公式(1)～(2)所示。

NMP=μ1·MP,

(1)

NMBx=μ2·MBx,

(2)

其中x=1,2,3,…,n，μ1>1,0<μ2<1。

在原积分级数下，P谱段和B谱段单个像元的曝光量为E0。改变积分级数后，全色谱段的曝光量变为

EP=E0·μ1,

(3)

多光谱谱段曝光量变为

EB=E0·μ2,

(4)

探测器原有动态范围是：

(5)

其中，E指探测器可探测到的曝光量大小。

积分级数调整之后，考虑到曝光量与积分级数成正比，拓展得到的新的动态范围为：

(6)

进而可知，动态范围拓展量为：

(7)

由式(7)可知，当μ1=2,μ2=0.25的时候，动态范围提高量为18.06 dB。

3 融合与插值算法

3.1 多谱段图像融合算法

观察各多光谱遥感卫星图像可知，各谱段图像的相关性很强，在多光谱图像数据集合中，存在相当多的信息冗余[16]。主成分分析算法(Principal Component Analysis,PCA)是一种统计分析方法，主要用于减少数据的维数[17-21]。在遥感领域，其常被用于多光谱与全色谱段图像融合：一般是将配准后的多光谱图像ΓB1，ΓB2，ΓB3，…，ΓBn，作为输入，经过PCA变换后，得到PC1，PC2，PC3，…,PCn共n个主成分，这些主成分对应n个从大到小的特征值γ1，γ2，γ3，…，γn。选取前3个主成分，用经过直方图匹配的全色谱段、高分辨率图像替换PCA变换结果中的第一主成分，然后再进行PCA反变换，得到高分辨率彩色图像的RGB成分。

主成分分析算法过程中，第一主成分图像包含了各输入谱段中的相同信息，而独特的光谱信息则体现在其它输出成分中[16]。各主成分分量对方差的贡献大小可以用各自对应的特征值来表示，从PC1到PCn逐渐减小。以特征值为权重，由各个主成分分量的加权平均值得到一幅灰度图像，该图像既包含了各个谱段的信息，又避免了信息冗余，与全色谱段的图像最为接近。

Γ1，Γ2，Γ3，…，Γn分别对应PC1，PC2，PC3，…,PCn共n个主成分，其像元灰度值可用Γ1(x,y)，Γ2(x,y)，Γ3(x,y)，…，Γn(x,y)表示。

经过PCA变换后的融合过程可以用公式(8)表示：

(8)

3.2 图像插值方法

对ΓB进行插值可以获取全色、高分辨率、长曝光量的遥感图像ΓBIN。

现有的各种插值方法，包括最近邻域法、双线性内插法和立方卷积内插法等方法，都是基于周边像元的灰度值进行推断的。本文提出充分利用对应区域内高积分级数、高分辨率全色谱段图像ΓP的像素分布特点，对ΓB进行插值得到ΓBIN。

根据图1所示的像元对应关系，与一个B谱段图像像元Bi，j所对应的是k×k个P谱段图像像元，利用这k×k个P谱段图像像元的灰度值比例关系估算插值后的B谱段k×k个像元的像素值，具体可以表示为：

(9)

3.3 高、低积分级数图像融合方法

根据公式(1)～(4)对现有各谱段积分级数关系的描述，可以将ΓBIN的像元灰度值γ1(x,y)和ΓP的像元灰度值γ2(x,y)按照以下方式进行融合：

由于

μ1<μ2,

(10)

一般有

γ1(x,y)≤γ2(x,y) ,

(11)

满量程灰度值表示为:

ΦFS=2β.

(12)

在高积分级数图像中，较暗区域显示比较清晰，低积分级数图像较亮区域显示比较清晰，为充分利用上述特征，可以作以下融合，得到高动态范围图像ΓHDR各像元灰度值γHDR(x,y)的表达式：

(13)

4 实验与结果

4.1 图像采集

选取一台TDI CCD相机进行成像实验。首先将TDI CCD积分级数分别按照表2第一组、第二组数据设置，获取相应图像并配准。

表2 两组积分级数设置Tab.2 Two groups of TDI stages

图2 两个积分级数各自对应的P谱段图像 Fig.2 Two images of band P with different values of MP

图3 按照第二组积分级数设置所得到的8个B谱段图像 Fig.3 Eight images of band B when selecting the second group of TDI stages

按照第一组积分级数设置，成像并配准后的P谱段图像为ΓP0，如图2(a)所示；按照第二组积分级数设置，成像并配准后的P谱段图像为ΓP，如图2(b)所示，B谱段图像ΓB1，ΓB2，ΓB3，…，ΓB8如图3所示。

4.2 多光谱图像的融合与上采样

将ΓB1，ΓB2，ΓB3，…，ΓB8作为输入，进行PCA变换，得到8个主成分分量，并按照3.1节中的方法得到融合后的低积分级数、低分辨率的灰度图像ΓB，如图4所示。

利用公式(9)对ΓB进行上采样，得到新的低积分级数、高分辨率的灰度图像ΓBIN，如图5所示。

图4 ΓB：由8个B谱段图像融合所得到的灰度图像 Fig.4 ΓB：gray image fused from the 8 images of band B

图5 ΓBIN：上采样后的ΓB Fig.5 ΓBIN:ΓB image after upsampled

4.3 高动态范围图像获取

依据公式(13)将低积分级数、高分辨率的图像ΓBIN与高积分级数、高分辨率图像ΓP进行融合，得到新的高动态范围、高分辨率灰度图像ΓHDR，如图6所示。

为了获取更好的显示效果，可以对灰度进行调整，使之更符合人眼的视觉特性，达到图像增强的目的[22-23]。对ΓHDR进行灰度调整，得到Γadjust，如图7所示。

图6 ΓHDR：由ΓBIN和ΓP融合得到的高动态范围图像 Fig.6 ΓHDR: the HDR image fused by ΓBIN and ΓP

图7 灰度调整后的高动态范围图像 Fig.7 HDR image after a grayscale adjustment

比较ΓBIN、ΓP和Γadjust，尤其是椭圆圈出的部分可以看出，ΓBIN中的较亮区域细节较为清晰，ΓP较暗区域细节较为清晰，而在高动态范围图像ΓHDR中，这两个区域的细节均可以较好地观察到。

显然，由于μ1=2,μ2=0.25，由公式(7)可得动态范围拓展量为18.06 dB。

5 结 论

根据推扫式航天遥感相机一般采用单排多谱段TDI CCD探测器，相机与场景之间存在高速相对运动的特点，本文提出了通过设置特殊积分级数提高相机动态范围的方法，并介绍了应用PCA变换实现多光谱图像到全色图像的融合方法，以及根据高分辨率全色谱段图像实现低分辨率图像上采样的方法。

实验与计算结果表明：该方法能有效提高推扫式航天遥感相机的动态范围，所获取的HDR图像能够涵盖比原图像更大的动态范围。当P谱段积分级数提高为原来的2倍， B谱段积分级数降低到原有积分级数的1/4时，动态范围可以提高18.06 dB。与常规遥感相机成像方法相比，利用该方法进行动态范围拓展时，遥感图像处理的时间复杂度有所增加，但图像处理的工作可完全在地面实现，即可以在不增加相机硬件成本、不改变原有光机结构的前提下，实现推扫式TDI CCD航天遥感相机动态范围拓展的目的，具有较高的工程应用价值。