基于在轨分类统计的热红外影像辐射校正方法

张炳先， 李岩， 何红艳

(北京空间机电研究所,北京 100076)



基于在轨分类统计的热红外影像辐射校正方法

张炳先， 李岩， 何红艳

(北京空间机电研究所,北京 100076)

相对辐射校正精度直接影响影像的应用效果，因此高精度的相对辐射校正是卫星影像处理的关键环节。目前热红外影像相对辐射校正的常用方法是借用星上黑体数据获取校正参数，因受到采样数据不足的限制，处理模型无法真实反映卫星载荷在轨工作状态，故无法满足高精度相对辐射校正的需要。因此，结合热红外载荷的工作特点，设计了基于在轨分类统计的热红外影像辐射校正方法，用以弥补星上校正精度低的不足。该方法利用卫星影像中地物类型丰富的特点，采用在轨获取的原始影像作为样本数据得到校正参数，弥补了星上校正样本不足的缺陷； 同时考虑到探元在不同能量下具有不同的光电响应状态，引入了分类思想进一步提高校正精度。结果表明，上述方法相对于星上校正方法能更好地抑制影像中的系统条带噪声，为提高判读精度提供了良好基础。

热红外影像； 辐射校正； 分类； 在轨统计

0 引言

热红外成像侦察能够感知目标的热红外波谱信息，大气穿透力强，可用于夜间侦察，在目标状态判断和伪装揭露等方面具有优势。由于现有制造工艺水平和材料的限制，热红外传感器中各探测单元的响应特性不完全一致，导致像元之间存在一定的非均匀性，使热红外影像中常出现深浅相间的条纹状噪声； 当影像被增强后，这种噪声会更加突出，呈百叶窗状，严重影响了热红外影像的质量和解析度。对热红外影像进行相对辐射校正，就是为了校正热红外载荷各个探元响应度差异而对卫星遥感器测量的原始数字计数值进行再量化的一种处理过程[1-3]。相对辐射校正通过对实验室校正数据进行分析，拟合出不同探元的光电响应模式，并求取相应的辐射校正系数或查找表，用于消除因探元响应不一致引起的条带效应，使这些条带的影响降低到最小程度或被彻底去除。相对辐射校正的前提和基础是精确的校正系数。目前热红外影像主要是通过星上校正数据来获取校正系数，由于热红外相机星上校正装置受到星上载荷空间和重量的限制，导致用于校正的样本数据不足，只能通过简单的线性模型来拟合探元的光电响应函数； 这与探元的实际在轨光电响应函数以及地物目标特性存在偏差，不能很好地描述复杂场景下目标温度特性，因此星上校正数据满足不了热红外影像高精度相对辐射校正的需求。在使用星上校正数据进行相对辐射校正后，某些复杂区域的影像仍存在均匀性不佳的现象，主要表现为条纹噪声。这种非均匀性降低了最终输出影像的清晰度，影响了目视判读和温度反演的精度，已成为进一步提高影像质量的瓶颈，并在一定程度上限制了热红外影像的应用。

鉴于上述情况，本文通过研究热红外载荷的工作特点，提出了一种基于在轨分类统计的热红外影像相对辐射校正方法。该方法利用卫星影像中地物类型丰富的特点，采用在轨获取的原始影像作为样本数据来获取相对辐射校正系数，以弥补星上校正数据不足的缺陷； 同时考虑到热红外成像探元在不同辐射能量下具有不同的光电响应状态[4]，本方法引入了分类思想，对不同辐射能量情况下的校正参数分别进行分类统计，以进一步提高相对辐射校正的精度； 并通过实验证明该方法相对于星上校正方法能更好地抑制影像中的系统条带噪声。

1 研究方法

通过对星上校正数据进行分析发现，星上校正数据主要存在2个方面的问题： ①不能完全覆盖探元的感光范围[5-9]，对于12 bit 量化的成像数据, 星上校正数据获取的灰度值(DN)最小值为2 692 ，最大值为3 441，高亮度区域和低亮度区域校正数据明显缺失； ②星上校正数据依赖于黑体数据，只能同时获取2个灰度值，所以只能采取简单的线性模型进行拟合，采样间隔较大等问题使得采用不能完全覆盖亮度区域的校正数据进行拟合得到的校正系数并不能正确反映出探元在不同感光区域的实际光电响应规律。因此，对采用局部数据拟合分析获取辐射校正系数的有效性存在疑问。

潘志强等[4]通过对卫星在轨影像的分析发现，探元在不同的感光区域具有不同的光电响应函数，其中，在中等亮度区域探元的光电响应函数为线性模型。由于校正数据主要集中在中等亮度区, 并且在该区域探元的光电响应函数呈线性分布，因此星上校正采用线性模型计算中等亮度区域的相对辐射校正参数是较为合理的。但是，采用中等亮度区域的光电响应函数直接替代探元在低亮度区和高亮度区的光电响应函数，在理论上没有支撑，在实际中也存在问题，因为探元在不同的亮度区域的光电响应函数是不同的[5]。因此，为了正确反映探元的光电响应函数，需要对探元的不同感光范围分别进行拟合处理, 其处理流程如图1所示。

图1 基于在轨分类统计的热红外影像相对辐射校正流程图

1.1 K均值分类

考虑到探元载荷主要是对低、中和高亮度3个区间段的响应存在较大差异，本文以影像的灰度量化信息作为样本集，取k=3，用X1,X2,X3分别表示低、中及高亮度地物聚类样本集，处理流程见图2。

若d1

图2 K均值算法处理示意图

Fig.2 Sketch map of K-MEANS processing

如图2所示，假设k=3；m1,m2,m3分别为3个样本聚类的初始中心值；d1,d2,d3分别为样本mp到m1,m2,m3的距离。若d1

假设热红外遥感影像采用BN位量化，那么使用

(1)

m2=2BN-1,

(2)

(3)

计算3个样本聚类的初始中心值m1,m2和m3，式中BN为量化的位数。计算每个样本数据xj(j=1,2,3,…,N)到3个样本聚类中心m1,m2,m3的距离dti，即

dti=d(xj,mi)=‖xj-mi‖,i=1,2,3 。

(4)

取dti中的最小值所对应聚类中心，将xj分配到对应的类别中。当所有样本数据全部分配结束后，重新计算每个聚类的中心，即

(5)

式中：Ni为第i个地物聚类Xi中的样本数目；N为参与统计的样本总数目；mi为样本的均值。

多次遍历所有的样本数据，并反复计算聚类中心，直到mi(i=1,2,3)的变化小于域值即停止迭代，此时便可得到低亮度(Dl)、中亮度和高亮度(Dh)分段区间[0,Dl)，[Dl,Dh)和[Dh,Dmax]。

1.2 分段辐射校正参数获取

由于探元在低亮度区域的光电响应函数往往具有较高的非线性特征，无法用数学函数去精确地拟合，而2个DN值的误差在影像中却会被肉眼所察觉而形成噪声，因此需要对低亮度区进行逐灰度级处理。考虑到若探元载荷之间不存在光电响应差异，那么在载荷获取相似灰度地物类型的地面区域时，其对应的灰度直方图分布就会相似或者相同，因此，本文对低亮度区间[0,Dl)采用直方图匹配处理的方式，逐灰度级地得到低亮度区域的校正参数。图3为直方图匹配示意图。

(a) 综合直方图 (b) 期望直方图

图3 直方图匹配处理示意图

Fig.3 Sketch map of histogram matching

先将影像中每行数据的相同探元的直方图进行合并，获得每个探元的综合直方图； 再将每个探元的综合直方图进行合并，即可获得所有探元的综合直方图。对于直方图匹配算法来说， 所有探元的综合直方图就是期望直方图，建立查找表的原理是使匹配处理后每个探元的综合直方图的概率密度函数与期望直方图的概率密度函数相同。

因此，首先统计每个探元综合直方图的概率密度函数和期望直方图的概率密度函数。令每个探元灰度值为k时对应的概率密度函数为pk，即

(6)

(7)

式中nj为此探元的灰度值为j的像素个数。

然后计算低亮度区间的期望直方图的概率密度函数Sk，即

(8)

(9)

式中rj为所有探元中灰度值为j的像元个数。

在获取探元的概率密度函数及期望的概率密度函数后，便可计算低亮度区间的辐射校正参数，具体计算方法为： 在标准累计概率密度函数上找出灰度值L，使得SL≤pk≤SL+1； 若|pk-SL|-|pk-SL+1|≤0，则用灰度值L代替灰度值k； 若|pk-SL|-|pk-SL+1|>0，则用灰度值(L+1)代替灰度值k。

由于探元在中间亮度区域的光电响应函数往往具有较高的线性特征，因此可以对中间亮度区间[Dl,Dh)进行最小二乘拟合处理，得到中间亮度区间的辐射校正系数[10]，最终得到该区间的辐射校正查找表，具体计算步骤如下：

1)计算中间亮度区间不同辐亮度下影像的整体灰度均值Zk及每个探元的灰度均值zi,k，“k=1,2,3，…”，其中k对应的最大值为中间区间的不同辐亮度级数的个数； “i=1,2,3，…”，i的最大值为探元的个数。

2)将单一探元在不同辐亮度下的灰度均值与不同辐亮度下整景影像的灰度均值通过最小二乘进行线性拟合，得到每个探元的增益Ki和偏置Ci，即

Zk=zi,kKi+Ci。

(10)

通过对实际在轨影像进行分析发现，在高亮度感光区域，由于能量充足，探元在这段区域内的光电响应函数趋于1∶1 的关系。因此，对高亮度感光区域，本文采用斜率为1、截距为0 的线性模型，以保持输入与输出1∶1 的关系不变。

1.3 过渡区处理

上述处理模型是基于分段进行处理的，在处理过程中难免会遇到辐射校正参数在边缘处断裂的现象。然而考虑到地物的灰度分布是一个连续的函数，因此需要在获取分段辐射校正参数的分段处进行过渡区处理，以保证辐射校正参数灰度分布的连续性。具体处理步骤如下。

(11)

(12)

式中：DNDl+5为低端区间在灰度值(Dl+5)校正后所对应的灰度值；DNDl-5为低端区间在灰度值(Dl-5)校正后所对应的灰度值。然后建立低亮度过渡区间的辐射校正查找表。

(13)

(14)

式中：DNDh+5为高端区间在灰度值(Dh+5)校正后所对应的灰度值；DNDh-5为高端区间在灰度值(Dh-5)校正后所对应的灰度值。然后建立高端过渡区间的辐射校正查找表。

2 实验与分析

为了验证本文算法的有效性，从中国资源卫星应用中心提供的热红外成像数据中选取若干区域的数据进行实验。各实验区域的影像大小都为512像元×512像元。采用本文算法与已有的基于星上校正数据的校正算法分别对该数据进行处理，并对处理结果进行了目视比较分析和定量分析。具体实验结果如图4和图5所示。

(a) 原始影像(b) 使用星上校正数据校正的影像 (c) 使用本文方法校正的影像

图4 陆地区域热红外影像不同方法辐射校正结果对比

Fig.4 Comparison of radiometric correction results of thermal-infrared image covered land area by using different methods

(a) 原始影像(b) 使用星上校正数据校正的影像 (c) 使用本文方法校正的影像

图5 水域区域热红外影像不同方法辐射校正结果对比

Fig.5 Comparison of radiometric correction results of thermal-infrared image covered water area by using different methods

通过目视可以看出，图4和图5中2种校正方法获得的校正影像的质量相比于原始影像均有了很大的提高，但在图4中，采用星上校正数据校正的影像仍存在大量的残余条带噪声，且分布较为分散，说明利用星上校正数据无法很好地消除能量较低区域中存在的系统条带噪声。图5中的水域区域同样也说明了上述情况。然而，对于图5中的能量充足区域，无论是采用星上校正的方法，还是采用本文方法都能够获得较好的校正效果。这正好验证了本文开头的假设，即探元在不同亮度区域的光电响应函数存在明显差异。为了进一步验证本文算法的有效性，本文选取广义噪声对校正前、后影像进行定量比较，广义噪声的计算公式为[11]

(15)

(16)

Re=E/Ave,

(17)

式中：m′和n′分别为影像的列数和行数；Ave为整景影像的均值；DNi′j′为影像第j′行第i′列的DN值；E为各列DN 值的平均值与整景影像均值差绝对值的平均；Re为通过该景影像计算得到的广义噪声。若相对辐射校正的精度较差, 则计算出来的广义噪声值越大。对图4和图5的广义噪声计算结果对比如表1所示。

表1 本文算法与已有算法处理结果的广义噪声值

从表1可以看出，图4中的地物类型较为均一，因此其广义噪声值整体偏小； 而图5中的地物类型对比差异较大，因此其广义噪声值整体偏大。通过表1的对比可以看出，本文方法无论是对低对比度区域的校正，还是对高对比度区域的校正，都具有良好的效果。

3 结论

本文根据实际生产实验中的经验和原则，分析了已有基于星上校正的热红外影像辐射校正方法的不足和局限性，提出了一种基于在轨分类统计的热红外影像辐射校正方法。该方法针对探元在不同能量段具有不同的光电响应函数的情况，引入了K均值分类的思想，对原始影像的地物进行粗分类； 然后对不同分类段的结果进行了辐射校正参数求解。实验结果表明，本文方法能够很好地去除影像中的系统条带噪声，为热红外影像的后续应用(尤其是目标判读精度的提高)提供良好的基础。

[1] Horn B K P,Woodham R J.Destriping Landsat MSS images by histogram modification[J].Computer Graphics and Image Processing,1979,10(1):69-83.

[2] Singh A.Thematic Mapper radiometric correction research and development results and performance[J].Photogrammetric Engineering and Remote Sensing,1985,51:1379-1383.

[3] 郭建宁,于晋,曾湧,等.CBERS-01/02卫星CCD图像相对辐射校正研究[J].中国科学E辑:信息科学,2005,35(s1):11-25. Guo J N,Yu J,Zeng Y,et al.Study on the relative radiometric correction of CBERS satellite CCD image[J].Science in China Series E Engineering & Materials Science,2005,48(s2):12-28.

[4] 潘志强,顾行发,刘国栋,等.基于探元直方图匹配的CBERS-01星CCD数据相对辐射校正方法[J].武汉大学学报:信息科学版,2005,30(10):925-927. Pan Z Q,Gu X F,Liu G D,et al.Relative radiometric correction of CBERS-01 CCD data based on detector histogram matching[J].Geomatics and Information Science of Wuhan University,2005,30(10):925-927.

[5] Weinreb M P,Xie R,Lienesch J H,et al.Destriping GOES images by matching empirical distribution functions[J].Remote Sensing of Environment,1989,29(2):185-195.

[6] Wegener M.Destriping multiple sensor imagery by improved histogram matching[J].International Journal of Remote Sensing,1990,11(5):859-875.

[7] Weinreb M,Hamilton G,Brown S,et al.Nonlinearity corrections in calibration of advanced very high resolution radiometer infrared channels[J].Journal of Geophysical Research,1990,95(c5):7381-7388.

[8] Wack E C.SWIR nonlinear radiometric response due to errors in relative spectral response[C]//Proceedings of SPIE 3117,Earth Observing Systems II.San Diego,CA,USA:SPIE,1997,3117:217-224.

[9] Raytheon S T X.CBESR Radiometric Algorithm Theoretical Basis Document[M].Beijing:National Association of Local Remote Sensing Application,2002:87-93.

[10]王小燕,龙小祥.资源一号02B星相机相对辐射校正方法分析[J].航天返回与遥感,2008,29(2):29-34. Wang X Y,Long X X.Analysis of relative radiometric correction method of CBERS-02B camera[J].Spacecraft Recovery & Remote Sensing,2008,29(2):29-34.

[11]王密,张炳先,潘俊.基于分段辐射校正的星载TDI-CCD成像数据辐射处理方法[J].中国科学:信息科学,2011,41(s):32-41. Wang M,Zhang B X,Pan J.Radiometric correction method of TDI-CCD imaging data based on segmentation[J].Science China Series E Engineering & Materials Science,2011,41(s):32-41.

(责任编辑： 邢宇)

Radiometric calibration method of thermal-infrared images based on on-orbit classification and statistics

ZHANG Bingxian， LI Yan， HE Hongyan

(BeijingInstituteofSpaceMechanicsandElectricity，Beijing100076，China)

Relative radiometric calibration precision will affect the application of satellite images, and hence high precision relative radiometric is very important. Nowadays, the ordinary radiometric calibration method of the thermal-infrared images is the on-orbit radiometric calibration which has fewer sample data and lower accuracy, therefore the result of on-orbit radiometric calibration can't satisfy the application requirement. In view of such a situation, a new radiometric calibration method of the thermal-infrared images based on classification and statistics is proposed in this paper. The new method adopts the original satellite images as sample data to calculate the parameters based on rich types of surface feature in satellite images to solve the problem of insufficient sample data of on-orbit radiometric calibration, meanwhile it introduces classification into the construction of the new calibration model by considering the different characteristics of photoelectric response function that the satellite payload will have in different radiometric energy so as to improve the precision of radiometric calibration. The experimental results show that the proposed method performs better than the ordinary on-orbit radiometric calibration method.

thermal-infrared images； radiometric calibration； classification； on-orbit statistics

10.6046/gtzyyg.2016.04.04

张炳先,李岩,何红艳.基于在轨分类统计的热红外影像辐射校正方法[J].国土资源遥感,2016,28(4)：24-29.(Zhang B X,Li Y,He H Y.Radiometric calibration method of thermal-infrared images based on on-orbit classification and statistics[J].Remote Sensing for Land and Resources,2016,28(4)：24-29.)

2015-05-19；

2015-06-17

TP 751.1; P 236

A

1001-070X(2016)04-0024-06

张炳先(1986-),男,博士,工程师,主要从事卫星影像处理方面的研究。Email： zbx@whu.edu.cn。