CN104992411A - 一种红外多光谱原始图像处理方法 - Google Patents

Abstract

本发明公开了一种红外多光谱原始图像处理方法，包括图像恢复、辐射校正和几何校正，图像恢复：包括缺元填补、正反扫错位校正、正反扫一致性校正、谱段内响应一致性校正、B7谱段去干扰、B8谱段缺行填补、B9波段插值、B9波段去细条纹、B9波段去条带拖尾；几何校正：包括成像投影方式引起的变形、地球曲率变化产生图像畸变、大气折射产生的影响、地球自转校正。与现有技术相比，本发明使图像处理水平有了明显的提高，本发明技术法运用使项目图像质量比以前有明显改善。

Description

一种红外多光谱原始图像处理方法

技术领域

本发明涉及图像处理技术领域，特别是一种红外多光谱原始图像处理方法。

背景技术

遥感是通过对反映地物电磁波辐射水平的灰度信息的处理分析与解译来进行地物识别和专题研究的。由于遥感成像过程中各种因素(例如，卫星速度变化、大气与地物反射与发射电磁波的相互作用、随机噪声等)的影响，实际的图像灰度值并不完全是地物辐射电磁波能量大小的反映，其中还包含着上述因素作用的结果，因此在进行遥感图像处理前，还需要进行校正处理(图像恢复处理)以消除上述因素的影响。这些处理习惯上称为预处理。

卫星遥感图像预处理主要目的是消除两大类影响：即遥感图像的辐射失真和遥感图像的几何畸变。遥感数字图像必须经过几何畸变处理(包括几何粗校正与几何精校正)、图像的辐射校正、噪声压制处理等预处理后，才能根据实际问题的需要进行其它的专门处理(例如图像的增强处理和图像的分类处理)。

预处理是图像处理中的一大类技术，通常意义上的预处理主要包含以下内容：

辐射校正：遥感的目的是利用遥感传感器有效地收集来自地物的电磁波辐射能量，然而由于电磁波在大气中的传输和传感器的测量过程中，受到遥感传感器本身灵敏特性、地物光照条件(如太阳高度角及地形变化等)以及大气作用等的影响，遥感传感器的测量值与地物实际的光谱、辐射率是不一致的，测量值存在着辐射失真。为了正确的评价地物的电磁波辐射特征，需要消除这些失真的影响，即进行辐射量校正处理。根据辐射失真的原因对遥感图像进行辐射校正处理可采取四种方法，见图1，即传感器校正、太阳高度角和地形引起的畸变校正以及大气散射校正。

几何校正：感图像在其成像过程中发生几何变形是难以避免的，引起遥感图像几何畸变的因素有三个方面，即传感器方面的原因、遥感平台(卫星)方面的原因以及地球本身的原因。

处理的红外多光谱扫描仪图像数据是原始图像数据：图像数据信号编码为8位，图像的：B6、B7、B8为1712×1536象元；B9为864×768象元。在B6、B7、B8中，每一个扫描带宽度为8个象元；在B9中，每一个扫描带宽度为4个象元。

红外多光谱扫描仪图像采集系统是通过扫描装置正向和反向扫描获取地物的反射、辐射能量。由于探测器的性能指标、扫描系统的控制精度级图像数据截取、卫星姿态等问题，红外多光谱扫描仪的原始图像数据存在一定的畸变(错位、压缩和拉伸等)和由于探测器响应的差异引起的一些条纹。针对各种不同问题产生的物理机制，进行相应的校正处理。预处理的内容主要包括数据恢复、红外多光谱扫描仪辐射校正、和几何校正四部分。

现有的技术方法往往是针对某一种缺陷设计的方法，并不够系统，不能够系统的解决图像中所存在的问题以用于实际生产中；其次，现有的技术方法对图像处理中的某些环节，如图像恢复的各个环节处理的不是很好。

发明内容

本发明的目的是要提供一种红外多光谱原始图像处理方法，针对图像中存在的缺陷，设计相应的数据恢复方法，根据传感器的设计特点，制定图像绝对定标和辐射校正以及几何校正方法，消除在成像过程中，卫星的运行、地球的自转和地球曲率以及遥感器等各种因素对图像数据造成的影响。

为达到上述目的，本发明是按照以下技术方案实施的：

一种红外多光谱原始图像处理方法，包括以下步骤：

1)图像恢复：所述图像恢复包括缺元填补、正反扫错位校正、正反扫一致性校正、谱段内响应一致性校正、B7谱段去干扰、B8谱段缺行填补、B9波段插值、B9波段去细条纹、B9波段去条带拖尾：

缺元填补：由于数据传输系统存在的问题，虽然探测元已将地表辐射信息接收下来，但是在向地面传送数据的时候有时会出现缺元的现象，即红外多光谱扫描仪的原始图像数据的某些像元的值为0；而在B8波段的图像则存在缺少第四和第七行探测元的数据的情况，使得这两行的数据值为0，将图像按8个条带分段，对每一段图像条带中的每个象素及其周围8行2列的象素与0灰度值进行比较，象素值为0的表示缺元，确定缺元的位置，然后利用左右插值的方法解决缺元；由于B9谱段为4个探测元，所以在B9的原始图像上表现出来的缺元的形式为4行2列共8个像元的灰度值为0，将图像按4个条带分段，对每一段图像条带中的每个象素及其周围4行2列的象素与0灰度值进行比较，象素值为0的表示缺元，确定缺元的位置，然后利用左右插值的方法进行缺元填补；

正反扫错位校正：由于红外多光谱扫描仪的图像数据在地面截取时正向和反向扫描数据的起始点选取不一致，因此B6、B7、B8中，正反扫描带中有n个像元的错位；在B9中，正反扫描带中有m个像元错位，以正向扫描带为基准，对反向的扫描带进行相应的平移以消除正向扫描带与反向扫描带之间的错位现象；在红外多光谱扫描仪的B6、B7、B8中8个探测元、B9为4个探测元之间的排列不是在一条直线上的，而是在奇偶像元之间有两个像元的错位，采用平移的方法进行消除；

正反扫一致性校正：由于扫描装置的正向和反向扫描的速度不完全一致，经过以上处理之后的图像在正反扫描带之间还存在一定的错位，采用以正向扫描带为基准，在相邻的反向扫描带上选取一定的控制点，把反向扫描带与正向配准，这组控制点可以通过大量图像的阅读得到的经验值，但是，由于扫描成像的特征，正向扫描带的扫描速度也不是匀速的，为了使图像不产生新的变形，采用CCD相机这种特定时刻一次成像的图像对正向和反向扫描带同时进行配准；

谱段内响应一致性校正：在B6、B7、B8谱段采用傅里叶变换，将各谱段的图像由时域转换为频域，然后将8个探测元一次扫描图像的数据在零频域上取平均，解决在像素灰度变化很小时8个探测元响应不一致的问题，最后再利用傅里叶逆变换恢复到时域，使图像的条带问题得到改善；

B7谱段去干扰：由于B7谱段中存在较为明显的周期性噪声，对图像的质量有所影响，所以应去除噪声，采用二维傅里叶变换在傅立叶频谱中提取尖峰状噪声分量进行去除，再进行傅立叶反变换；

B8谱段缺行填补：B8波段的图像缺少第四和第七行探测元的数据，采用缺行插值的方法弥补，用相邻的上下两行的数据值进行平均差值处理；

B9波段插值：将B9图像扩展为原始宽度*2×原始高度*2大小；

B9波段去细条纹：由于温度的变化和探测器材料的不同，探测器的传递函数也互不相同，而且探测器具有非线性传递特性，这些原因使得图像具有规则条纹，在B9波段上表现得尤为严重，采用图像中每一条数据带与其相邻的上下两条数据条带求其平均值后替代原有的数据，以减少条带；

B9波段去条带拖尾：条带产生的原因是探测器与前置放大器是交流耦合，由RC组成的高通交流耦合电路，在对目标进行扫描时会产生图像缺陷，采用在空间域上利用差分的方法探测条带及海陆边缘，并控制海面的灰度阈值，利用正扫条带对反扫条带进行平滑和利用无拖尾条带对拖尾条带去拖尾；

2)辐射校正：

利用实验室绝对定标数据进行辐射校正：实验中谱段6、7、8的定标源采用的是500W的太阳模拟照射漫反射版，漫反射板置于准直仪的焦面，经准直仪产生平行光照射扫描仪。改变太阳模拟器输出的辐亮度，测量扫描仪输出的数据量，从而对扫描仪进行定标；

谱段9的定标方法是利用置于准直仪的焦平面的黑体，经准直仪产生平行辐射照射扫描仪主体，改变黑体温度，测量扫描仪输出的数字量，对扫描仪谱段9进行定标；

卫星遥感器接收到的地－气系统反射辐射信号或辐射信号是以电信号方式并以计数器输出的，这些计数器与卫星遥感器所接收的辐亮度成正比：

DN＝G·L_e+DN。

式中DN―遥感器的计数值；

G―遥感器的辐射响应度，也称定标系数；

DN₀―遥感器的零位计数值，也称定标系数；

L_e―遥感器入瞳处的等效辐亮度：

$L_e=\frac{\underset{{\mathrm{\λ}}_1}{\overset{{\mathrm{\λ}}_2}{\∫}}{L}_{\mathrm{\λ}}{R}_{\mathrm{\λ}}{d}_{\mathrm{\λ}}}{\underset{{\mathrm{\λ}}_1}{\overset{{\mathrm{\λ}}_2}{\∫}}{R}_{\mathrm{\λ}}{d}_{\mathrm{\λ}}}$

式中L_λ―遥感器入瞳处的光谱辐亮度；

R_λ―遥感器的光谱响应函数；

λ₁、λ₂―某谱段的起始波长和终止波长，

辐射定标就是由不同的Le和DN用最小二乘法计算得出定标系数G和DN；

3)几何校正：所述几何校正包括成像投影方式引起的变形、地球曲率变化产生图像畸变、大气折射产生的影响、地球自转校正：

成像投影方式引起的变形：全景投影是一种特殊形式的线中心投影，全景投影的影像面不是一个平面，而是一个圆柱面，地物点P在全景面上的像点p在扫描方向上的坐标为y_p

y_p＝f·θ

其中：f为焦距，θ为成像角(以度为单位)

设L是一个等效的中心投影成像面，P点在L上的像点p’的坐标为y’_p为：

y'_p＝f·tgθ

从而根据以上公式可以导出全景变形公式为：

dy＝y_p-y'_p＝f·(θ-tgθ)；

地球曲率变化产生的畸变：地球曲率将引起像点位移，使得图像边缘的比例尺比中心的比例尺小，经过推算，一幅红外遥感图像大致将要损失10个像元，其表达式如下：

dx＝0

$\mathrm{dy}=-\frac{f\·\mathrm{tg\θ}\·{\sin }^2\mathrm{\θ}}{2R\·f^2}\·h$

其中，R为地球半径，H为传感器航高，f为焦距；

大气折射产生的影响：对于光波和电磁波的传播而言，大气层并非是一个均匀的介质，因为它的密度是随离地面的高度增加而递减的，所以光波、电磁波在大气层中的折射率也是随高度而变的，从而使电磁波传播的路径不是一条直线而变成了曲线，进而引起了像点位移，象点的位移公式如下，K是一个与传感器航高和地物点高程有关大气的条件常数，其表达式为：

Δ_r＝K(r+r³/f²)

$K=\frac{2410H}{H^2-6H+250}$

其中的H为传感器航高，在x,y两个方向上考虑大气折射时，将r分别用图像坐标x和y来代替处理，

d_x＝k*x*(1+r²/f)

d_y＝k*y*(1+r²/f)；

4)地球自转的影响：地球自转主要是对动态传感器的图像产生变形影响，特别是对卫星遥感影像，当卫星由北向南运行的同时，地球表面也在由西向东自转，由于卫星图像每条扫描线的成像时间不同，因而造成扫描线在地面上的投影依次向西平移，最终使得图像发生扭曲，计算公式如下：

$\mathrm{\Δx}=\left[\frac{m_y}{m_x}\right]\·\left[\frac{{\mathrm{\ω}}_e}{{\mathrm{\ω}}_s}\right]\·\sin \mathrm{\ϵ}\·x$

$\mathrm{\Δy}=\frac{{\mathrm{\ω}}_e}{{\mathrm{\ω}}_s}]\·\sqrt{{\cos }^2\mathrm{\φ}-{\sin }^2\mathrm{\ϵ}}\·y$

$\mathrm{\θ}=\left[\frac{{\mathrm{\ω}}_e}{{\mathrm{\ω}}_s}\right]\·\sqrt{{\cos }^2\mathrm{\φ}-{\sin }^2\mathrm{\ϵ}}$

其中，△x、△y、θ分别为图像底边中点的坐标位移，以及平均偏航角，1/m_x和1/m_y分别为影像x和y方向的比例尺，ω_e和ω_s分别为地球自转角速度和卫星轨道角速度。

与现有技术相比，本发明使图像处理水平有了明显的提高，本发明技术法运用使项目图像质量比以前有明显改善；而对于图像数据结构的分析研究，在生产中解决了遥感数据用常规软件方法无法读取的问题，可以明显节约数据重新购置费用，保证项目的按时进展。

附图说明

图1是本发明的流程图；

图2是本发明的条带拖尾的耦合信号波形图，a为耦合前的信号波形图，b为耦合后的信号波形图；

图3是本发明的谱段6的定标曲线；

图4是本发明的谱段7的定标曲线；

图5是本发明的谱段8的定标曲线；

图6是本发明的谱段9的主份的定标曲线；

图7是本发明的谱段6的备份的定标曲线；

图8是本发明全景投影时的物像关系图；

图9是本发明地球曲率的影响图；

图10是本发明大气折射差关系图；

图11：a是缺元填补前的图像，b是缺元填补后的图像；

图12：a是正反扫错位校正前的图像，b是正反扫错位校正后的图像；

图13：a是正反扫一致性校正前的图像，b是正反扫一致性校正后的图像；

图14：a是谱段内响应一致性校正前的图像，b是谱段内响应一致性校正后的图像；

图15：a是傅里叶变换前的图像，b是傅里叶变换后的图像；

图16：a是B7谱段去干扰前的图像，b是B7谱段去干扰后的图像；

图17：a是B8谱段的缺行填补前的图像，b是谱段的缺行填补后的图像；

图18：a是B9谱段的去细条纹前的图像，b是B9谱段的去细条纹后的图像：

图19：a是B9谱段的去条带拖尾前的图像，b是B9谱段的去条带拖尾后的图像；

图20：a是辐射校正前的图像，b是辐射校正后的图像；

图21：a是几何校正前的原始图像，b是几何校正后的图像。

具体实施方式

下面结合具体实施例对本发明作进一步描述，在此发明的示意性实施例以及说明用来解释本发明，但并不作为对本发明的限定。

如图1所示的红外多光谱原始图像处理方法，包括以下步骤：

1)图像恢复：所述图像恢复包括缺元填补、正反扫错位校正、正反扫一致性校正、谱段内响应一致性校正、B7谱段去干扰、B8谱段缺行填补、B9波段插值、B9波段去细条纹、B9波段去条带拖尾：

缺元填补：由于数据传输系统存在的问题，虽然探测元已将地表辐射信息接收下来，但是在向地面传送数据的时候有时会出现缺元的现象，即红外多光谱扫描仪的原始图像数据的某些像元的值为0；而在B8波段的图像则存在缺少第四和第七行探测元的数据的情况，使得这两行的数据值为0，缺元填补前的图像如图11a，将图像按8个条带分段，对每一段图像条带中的每个象素及其周围8行2列的象素与0灰度值进行比较，象素值为0的表示缺元，确定缺元的位置，然后利用左右插值的方法解决缺元；由于B9谱段为4个探测元，所以在B9的原始图像上表现出来的缺元的形式为4行2列共8个像元的灰度值为0，将图像按4个条带分段，对每一段图像条带中的每个象素及其周围4行2列的象素与0灰度值进行比较，象素值为0的表示缺元，确定缺元的位置，然后利用左右插值的方法进行缺元填补，缺元填补后的图像如图11b所示，缺元填补的运行代码：

正反扫错位校正：由于红外多光谱扫描仪的图像数据在地面截取时正向和反向扫描数据的起始点选取不一致，因此B6、B7、B8中，正反扫描带中有n个像元的错位，正反扫错位校正前的图像如图12a所示；在B9中，正反扫描带中有m个像元错位，以正向扫描带为基准，对反向的扫描带进行相应的平移以消除正向扫描带与反向扫描带之间的错位现象；在红外多光谱扫描仪的B6、B7、B8中8个探测元、B9为4个探测元之间的排列不是在一条直线上的，而是在奇偶像元之间有两个像元的错位，采用平移的方法进行消除，正反扫错位校正后的图像如图12b所示，正反扫错位校正的代码：

正反扫一致性校正：如图13a所示，由于扫描装置的正向和反向扫描的速度不完全一致，经过以上处理之后的图像在正反扫描带之间还存在一定的错位，采用以正向扫描带为基准，在相邻的反向扫描带上选取一定的控制点，把反向扫描带与正向配准，这组控制点可以通过大量图像的阅读得到的经验值，但是，由于扫描成像的特征，正向扫描带的扫描速度也不是匀速的，为了使图像不产生新的变形，采用CCD相机这种特定时刻一次成像的图像对正向和反向扫描带同时进行配准，正反扫一致性校正后的图像如图13b所示，正反扫一致性校正的代码：

谱段内响应一致性校正：在B6、B7、B8谱段采用傅里叶变换，将各谱段的图像由时域转换为频域，然后将8个探测元一次扫描图像的数据在零频域上取平均，解决在像素灰度变化很小时8个探测元响应不一致的问题，如图14中a的谱段内响应一致性校正前的图像中条带现象明显，如图14的b的谱段内响应一致性校正的图像中条带现象消失；最后再利用傅里叶逆变换恢复到时域，使图像的条带问题得到改善，傅里叶变换前的图像如图15a，傅里叶变换后的图像如图15b；谱段内响应一致性校正的代码：

B7谱段去干扰：由于B7谱段中存在较为明显的周期性噪声，对图像的质量有所影响，去干扰前谱域图像如图16a所示，所以应去除噪声，采用二维傅里叶变换在傅立叶频谱中提取尖峰状噪声分量进行去除，再进行傅立叶反变换，去干扰后谱域图像如图16b；B7谱段去干扰的代码：

B8谱段缺行填补：B8波段的图像缺少第四和第七行探测元的数据如图17a所示，采用缺行插值的方法弥补，用相邻的上下两行的数据值进行平均差值处理，得到缺行处理之后的图像如图17b；B8谱段的缺行填补的代码：

B9波段插值：将B9图像扩展为原始宽度*2×原始高度*2大小；B9波段插值的

B9波段去细条纹：由于温度的变化和探测器材料的不同，探测器的传递函数也互不相同，而且探测器具有非线性传递特性，这些原因使得图像具有规则条纹，在B9波段上表现得尤为严重，参见图18a，采用图像中每一条数据带与其相邻的上下两条数据条带求其平均值后替代原有的数据，以减少条带，经去细条纹处理之后的B9图像如图18b；B9去细条纹的代码：

B9波段去条带拖尾：条带产生的原因是探测器与前置放大器是交流耦合，由RC组成的高通交流耦合电路，在对目标进行扫描时会产生图像缺陷。对中等温度差异的大目标来说，耦合前的信号为矩形波，耦合后的信号失真产生直流下跌和负尖峰，因此图像目标产生畸变，耦合前的信号波形图如图2a)，对高温小目标来说，由于电路输出的平均值很低，所以输出信号在正信号响应之后将伴随一个振幅较低但持续时间较长的负信号响应，从图像上看会产生较严重的黑色拖尾现象如图19a所示。处理的方法是采用在空间域上利用差分的方法探测条带及海陆边缘，并控制海面的灰度阈值，利用正扫条带对反扫条带进行平滑和利用无拖尾条带对拖尾条带去拖尾，经去条带拖尾处理后的结果图像如图19b所示。这种方法可以很好的去除海面上的条带和拖尾，但是阈值需要经验值给定，经过处理后的图像随着阈值的不同，可能会对海陆边缘部分或者海面较弱信息有所削弱，耦合后的信号波形图如图2b)；B9去条带拖尾的代码：

4)辐射校正：

利用传感器观测目标的反射和折射能量时，传感器由于在测量时得到的测量值包含了太阳位置和角度条件、薄雾等大气条件所引起的失真，使得测量值与目标的光谱反射率或光谱辐射度等物理量是不一致的，从而影响了图像的质量和应用，必须进行消除。即辐射校正。通过反演地物的光谱反射特性和光谱辐射特性，从而可正确地判读和有效地利用卫星遥感图像。利用实验室绝对定标数据进行辐射校正：卫星发射前，星载遥感器进行实验室辐射定标。红外多光谱扫描仪真空辐射定标的目的是在模拟空间环境的条件下确定扫描仪各通道的定标系数。测量的方法是在已知辐射输入条件下测量扫描仪的输出信号，并在不同的温度条件下进行测试，获得不同温度时的定标数据，实验中谱段6、7、8的定标源采用的是500W的太阳模拟照射漫反射版，漫反射板置于准直仪的焦面，经准直仪产生平行光照射扫描仪。改变太阳模拟器输出的辐亮度，测量扫描仪输出的数据量，从而对扫描仪进行定标；

谱段9的定标方法是利用置于准直仪的焦平面的黑体，经准直仪产生平行辐射照射扫描仪主体，改变黑体温度，测量扫描仪输出的数字量，对扫描仪谱段9进行定标；

卫星遥感器接收到的地－气系统反射辐射信号或辐射信号是以电信号方式并以计数器输出的，这些计数器与卫星遥感器所接收的辐亮度成正比：

DN＝G·L_e+DN。

式中DN―遥感器的计数值；

G―遥感器的辐射响应度，也称定标系数；

DN₀―遥感器的零位计数值，也称定标系数；

L_e―遥感器入瞳处的等效辐亮度：

$L_e=\frac{\underset{{\mathrm{\λ}}_1}{\overset{{\mathrm{\λ}}_2}{\∫}}{L}_{\mathrm{\λ}}{R}_{\mathrm{\λ}}{d}_{\mathrm{\λ}}}{\underset{{\mathrm{\λ}}_1}{\overset{{\mathrm{\λ}}_2}{\∫}}{R}_{\mathrm{\λ}}{d}_{\mathrm{\λ}}}$

式中L_λ―遥感器入瞳处的光谱辐亮度；

R_λ―遥感器的光谱响应函数；

λ₁、λ₂―某谱段的起始波长和终止波长，

辐射定标就是由不同的Le和DN用最小二乘法计算得出定标系数G和DN，并绘制出定标曲线，如图3所示是谱段6的定标曲线，谱段6的定标方程为：

y＝1.1710x-0.44  (3.2-9)

y＝1.1700x-0.96  (3.2-10)

y＝1.1755x-0.62  (3.2-11)

y＝1.1718x-1.45  (3.2-12)

y＝1.1799x-0.63  (3.2-13)

y＝1.1872x-1.08  (3.2-14)

y＝1.1572x+1.07  (3.2-15)

y＝1.1773x-1.34  (3.2-16)；

图4是谱段7定标曲线，谱段7的定标方程为：

y＝4.4981x-0.67  (3.2-17)

y＝4.4256x-1.03  (3.2-18)

y＝4.5261x-1.07  (3.2-19)

y＝4.5015x-1.01  (3.2-20)

y＝4.5423x-1.33  (3.2-21)

y＝4.5220x-1.65  (3.2-22)

y＝4.5205x-0.98  (3.2-23)

y＝4.5131x-0.90  (3.2-24)；

图5是谱段8的定标曲线，谱段8的定标方程：

y＝12.117x-1.39  (3.2-25)

y＝12.102x-1.56  (3.2-26)

y＝12.334x-1.68  (3.2-27)

y＝12.295x-2.05  (3.2-28)

y＝12.273x-1.54  (3.2-29)

y＝11.958x-1.83  (3.2-30)

y＝12.298x-1.70  (3.2-31)

y＝12.101x-1.51  (3.2-32)；

图6是谱段9的定标曲线，谱段9的定标方程：

y＝23.246x-84.80  R²＝1  (3.2-33)

y＝23.536x-86.91  R²＝0.9999  (3.2-34)

y＝22.622x-79.73  R²＝0.9999  (3.2-35)

y＝23.723x-88.40  R²＝1  (3.2-36)；图7是谱段9定标曲线，谱段9定标方程：

y＝22.43x-86.34  R²＝0.9999  (3.2-37)

y＝24.119x-91.71  R²＝0.9999  (3.2-38)

y＝22.886x-81.88  R²＝0.9999  (3.2-39)

y＝23.81x-89.16  R²＝0.9999  (3.2-40)；

以Band6谱段为例具体处理过程为：

具体效果见图20，其中a是辐射校正前的图像，b是辐射校正后的图像。

5)几何校正：所述几何校正包括成像投影方式引起的变形、地球曲率变化产生图像畸变、大气折射产生的影响、地球自转校正：

成像投影方式引起的变形：全景投影是一种特殊形式的线中心投影，全景投影的影像面不是一个平面，而是一个圆柱面，如图8所示，地物点P在全景面上的像点p在扫描方向上的坐标为y_p

y_p＝f·θ

其中：f为焦距，θ为成像角(以度为单位)

设L是一个等效的中心投影成像面，P点在L上的像点p’的坐标为y’_p为：

y'_p＝f·tgθ

从而根据以上公式可以导出全景变形公式为：

dy＝y_p-y'_p＝f·(θ-tgθ)；

地球曲率变化产生的畸变：地球曲率将引起像点位移，使得图像边缘的比例尺比中心的比例尺小，经过推算，一幅红外遥感图像大致将要损失10个像元，其表达式如下：

dx＝0

$\mathrm{dy}=-\frac{f\·\mathrm{tg\θ}\·{\sin }^2\mathrm{\θ}}{2R\·f^2}\·h$

其中，R为地球半径，H为传感器航高，f为焦距；

如图9所示：s为卫星成像时的位置，H为卫星到地面的高度，p为地物点坐标，p0为地物点p在全景面上的坐标。

大气折射产生的影响：对于光波和电磁波的传播而言，大气层并非是一个均匀的介质，因为它的密度是随离地面的高度增加而递减的，所以光波、电磁波在大气层中的折射率也是随高度而变的，从而使电磁波传播的路径不是一条直线而变成了曲线，进而引起了像点位移。中心投影属于方向投影，其成像点的位置取决于地物点的入射光线的方向。如图10所示，在无大气折射影响的时候，地物点A通过直线光线AS成像于a0点，而有大气干扰时，A点通过曲线AS成像于a1点，因此而引起的像点位移△r＝a1 a0。象点的位移公式如下，K是一个与传感器航高和地物点高程有关大气的条件常数，其表达式为：

Δ_r＝K(r+r³/f²)

$K=\frac{2410H}{H^2-6H+250}$

其中的H为传感器航高，在x,y两个方向上考虑大气折射时，将r分别用图像坐标x和y来代替处理，

d_x＝k*x*(1+r²/f)

d_y＝k*y*(1+r²/f)；

4)地球自转的影响：地球自转主要是对动态传感器的图像产生变形影响，特别是对卫星遥感影像，当卫星由北向南运行的同时，地球表面也在由西向东自转，由于卫星图像每条扫描线的成像时间不同，因而造成扫描线在地面上的投影依次向西平移，最终使得图像发生扭曲，计算公式如下：

$\mathrm{\Δx}=\left[\frac{m_y}{m_x}\right]\·\left[\frac{{\mathrm{\ω}}_e}{{\mathrm{\ω}}_s}\right]\·\sin \mathrm{\ϵ}\·x$

$\mathrm{\Δy}=\frac{{\mathrm{\ω}}_e}{{\mathrm{\ω}}_s}]\·\sqrt{{\cos }^2\mathrm{\φ}-{\sin }^2\mathrm{\ϵ}}\·y$

$\mathrm{\θ}=\left[\frac{{\mathrm{\ω}}_e}{{\mathrm{\ω}}_s}\right]\·\sqrt{{\cos }^2\mathrm{\φ}-{\sin }^2\mathrm{\ϵ}}$

其中，△x、△y、θ分别为图像底边中点的坐标位移，以及平均偏航角，1/m_x和1/m_y分别为影像x和y方向的比例尺，ω_e和ω_s分别为地球自转角速度和卫星轨道角速度。

几何校正采用的图像为系统接收的tif图像即几何校正前图像如图21a(在几何校正处理中，认为成像时卫星运行轨道与设计相同)，几何校正后的图像具体效果见图21b。

本发明的技术方案不限于上述具体实施例的限制，凡是根据本发明的技术方案做出的技术变形，均落入本发明的保护范围之内。

Claims (1)

1.一种红外多光谱原始图像处理方法，其特征在于，包括以下步骤：

1)图像恢复：所述图像恢复包括缺元填补、正反扫错位校正、正反扫一致性校正、谱段内响应一致性校正、B7谱段去干扰、B8谱段缺行填补、B9波段插值、B9波段去细条纹、B9波段去条带拖尾：

缺元填补：由于数据传输系统存在的问题，虽然探测元已将地表辐射信息接收下来，但是在向地面传送数据的时候有时会出现缺元的现象，即红外多光谱扫描仪的原始图像数据的某些像元的值为0；而在B8波段的图像则存在缺少第四和第七行探测元的数据的情况，使得这两行的数据值为0，将图像按8个条带分段，对每一段图像条带中的每个象素及其周围8行2列的象素与0灰度值进行比较，象素值为0的表示缺元，确定缺元的位置，然后利用左右插值的方法解决缺元；由于B9谱段为4个探测元，所以在B9的原始图像上表现出来的缺元的形式为4行2列共8个像元的灰度值为0，将图像按4个条带分段，对每一段图像条带中的每个象素及其周围4行2列的象素与0灰度值进行比较，象素值为0的表示缺元，确定缺元的位置，然后利用左右插值的方法进行缺元填补；

正反扫错位校正：由于红外多光谱扫描仪的图像数据在地面截取时正向和反向扫描数据的起始点选取不一致，因此B6、B7、B8中，正反扫描带中有n个像元的错位；在B9中，正反扫描带中有m个像元错位，以正向扫描带为基准，对反向的扫描带进行相应的平移以消除正向扫描带与反向扫描带之间的错位现象；在红外多光谱扫描仪的B6、B7、B8中8个探测元、B9为4个探测元之间的排列不是在一条直线上的，而是在奇偶像元之间有两个像元的错位，采用平移的方法进行消除；

正反扫一致性校正：由于扫描装置的正向和反向扫描的速度不完全一致，经过以上处理之后的图像在正反扫描带之间还存在一定的错位，采用以正向扫描带为基准，在相邻的反向扫描带上选取一定的控制点，把反向扫描带与正向配准，这组控制点可以通过大量图像的阅读得到的经验值，但是，由于扫描成像的特征，正向扫描带的扫描速度也不是匀速的，为了使图像不产生新的变形，采用CCD相机这种特定时刻一次成像的图像对正向和反向扫描带同时进行配准；

谱段内响应一致性校正：在B6、B7、B8谱段采用傅里叶变换，将各谱段的图像由时域转换为频域，然后将8个探测元一次扫描图像的数据在零频域上取平均，解决在像素灰度变化很小时8个探测元响应不一致的问题，最后再利用傅里叶逆变换恢复到时域，使图像的条带问题得到改善；

B7谱段去干扰：由于B7谱段中存在较为明显的周期性噪声，对图像的质量有所影响，所以应去除噪声，采用二维傅里叶变换在傅立叶频谱中提取尖峰状噪声分量进行去除，再进行傅立叶反变换；

B8谱段缺行填补：B8波段的图像缺少第四和第七行探测元的数据，采用缺行插值的方法弥补，用相邻的上下两行的数据值进行平均差值处理；

B9波段插值：将B9图像扩展为原始宽度*2×原始高度*2大小；

B9波段去细条纹：由于温度的变化和探测器材料的不同，探测器的传递函数也互不相同，而且探测器具有非线性传递特性，这些原因使得图像具有规则条纹，在B9波段上表现得尤为严重，采用图像中每一条数据带与其相邻的上下两条数据条带求其平均值后替代原有的数据，以减少条带；

B9波段去条带拖尾：条带产生的原因是探测器与前置放大器是交流耦合，由RC组成的高通交流耦合电路，在对目标进行扫描时会产生图像缺陷，采用在空间域上利用差分的方法探测条带及海陆边缘，并控制海面的灰度阈值，利用正扫条带对反扫条带进行平滑和利用无拖尾条带对拖尾条带去拖尾；

2)辐射校正：

利用实验室绝对定标数据进行辐射校正：实验中谱段6、7、8的定标源采用的是500W的太阳模拟照射漫反射版，漫反射板置于准直仪的焦面，经准直仪产生平行光照射扫描仪。改变太阳模拟器输出的辐亮度，测量扫描仪输出的数据量，从而对扫描仪进行定标；

谱段9的定标方法是利用置于准直仪的焦平面的黑体，经准直仪产生平行辐射照射扫描仪主体，改变黑体温度，测量扫描仪输出的数字量，对扫描仪谱段9进行定标；

卫星遥感器接收到的地－气系统反射辐射信号或辐射信号是以电信号方式并以计数器输出的，这些计数器与卫星遥感器所接收的辐亮度成正比：

DN＝G·L_e+DN_o

式中DN―遥感器的计数值；

G―遥感器的辐射响应度，也称定标系数；

DN₀―遥感器的零位计数值，也称定标系数；

L_e―遥感器入瞳处的等效辐亮度：

$L_e=\frac{\underset{{\mathrm{\λ}}_1}{\overset{{\mathrm{\λ}}_2}{\∫}}{L}_{\mathrm{\λ}}{R}_{\mathrm{\λ}}{d}_{\mathrm{\λ}}}{\underset{{\mathrm{\λ}}_1}{\overset{{\mathrm{\λ}}_2}{\∫}}{R}_{\mathrm{\λ}}{d}_{\mathrm{\λ}}}$

式中L_λ―遥感器入瞳处的光谱辐亮度；

R_λ―遥感器的光谱响应函数；

λ₁、λ₂―某谱段的起始波长和终止波长，

辐射定标就是由不同的Le和DN用最小二乘法计算得出定标系数G和DN；

3)几何校正：所述几何校正包括成像投影方式引起的变形、地球曲率变化产生图像畸变、大气折射产生的影响、地球自转校正：

成像投影方式引起的变形：全景投影是一种特殊形式的线中心投影，全景投影的影像面不是一个平面，而是一个圆柱面，地物点P在全景面上的像点p在扫描方向上的坐标为y_p

y_p＝f·θ

其中：f为焦距，θ为成像角(以度为单位)

设L是一个等效的中心投影成像面，P点在L上的像点p’的坐标为y’_p为：

y'_p＝f·tgθ

从而根据以上公式可以导出全景变形公式为：

dy＝y_p-y'_p＝f·(θ-tgθ)；

地球曲率变化产生的畸变：地球曲率将引起像点位移，使得图像边缘的比例尺比中心的比例尺小，经过推算，一幅红外遥感图像大致将要损失10个像元，其表达式如下：

dx＝0

$\mathrm{dy}=-\frac{f\·\mathrm{tg\θ}\·{\sin }^2\mathrm{\θ}}{2R\·f^2}\·h$

其中，R为地球半径，H为传感器航高，f为焦距；

大气折射产生的影响：对于光波和电磁波的传播而言，大气层并非是一个均匀的介质，因为它的密度是随离地面的高度增加而递减的，所以光波、电磁波在大气层中的折射率也是随高度而变的，从而使电磁波传播的路径不是一条直线而变成了曲线，进而引起了像点位移，象点的位移公式如下，K是一个与传感器航高和地物点高程有关大气的条件常数，其表达式为：

Δ_r＝K(r+r³/f²)

$K=\frac{2410H}{H^2-6H+250}$

其中的H为传感器航高，在x,y两个方向上考虑大气折射时，将r分别用图像坐标x和y来代替处理，

d_x＝k*x*(1+r²/f)

d_y＝k*y*(1+r²/f)；

4)地球自转的影响：地球自转主要是对动态传感器的图像产生变形影响，特别是对卫星遥感影像，当卫星由北向南运行的同时，地球表面也在由西向东自转，由于卫星图像每条扫描线的成像时间不同，因而造成扫描线在地面上的投影依次向西平移，最终使得图像发生扭曲，计算公式如下：

$\mathrm{\Δx}=\left[\frac{m_y}{m_x}\right]\·\left[\frac{{\mathrm{\ω}}_e}{{\mathrm{\ω}}_s}\right]\·\sin \mathrm{\ϵ}\·x$

$\mathrm{\Δy}=\frac{{\mathrm{\ω}}_e}{{\mathrm{\ω}}_s}]\·\sqrt{{\cos }^2\mathrm{\φ}{-\sin }^2\mathrm{\ϵ}}\·y$

$\mathrm{\θ}=\left[\frac{{\mathrm{\ω}}_e}{{\mathrm{\ω}}_s}\right]\·\sqrt{{\cos }^2\mathrm{\φ}-{\sin }^2\mathrm{\ϵ}}$

其中，△x、△y、θ分别为图像底边中点的坐标位移，以及平均偏航角，1/m_x和1/m_y分别为影像x和y方向的比例尺，ω_e和ω_s分别为地球自转角速度和卫星轨道角速度。