CN102663693A - 一种基于最小二乘法的线阵推扫式影像自适应辐射校正方法 - Google Patents

Abstract

一种基于最小二乘法的线阵推扫式影像自适应辐射校正方法，处理步骤包括：首先，定位影像中的坏扫描列并对该列左右正常列的像素灰度值进行加权平均，消除坏扫描列对后续处理的影响；其次，分析线阵推扫式影像中条带的成因和分布特点，根据用户设定的辐射校正精度值自适应判别条带列；再次，对于判别得到的条带列，提出用最小二乘法拟合计算相对辐射校正系数，并对条带列的所有像素进行相对辐射校正；最后，将相对辐射校正系数存储入库，为后续快速应用提供查询数据库。本发明的技术方案能够有效保证相对辐射校正精度和影像数据的物理意义不改变。

Description

一种基于最小二乘法的线阵推扫式影像自适应辐射校正方法

技术领域

本发明涉及一种基于最小二乘法的线阵推扫式影像自适应辐射校正方法，属于遥感影像处理领域。

背景技术

随着空间技术、信息技术和传感器技术的快速发展，遥感影像的分辨率得到很大提高。而在各类遥感成像传感器中，线阵CCD推扫式传感器得到广泛应用。例如，国外的SPOT、IKONOS、QuickBird卫星以及我国的CBERS-01/02、HJ1-A/B等卫星都搭载有线阵CCD推扫式传感器。到目前为止，线阵推扫式传感器是对地观测有效的传感器，预计在今后该类传感器仍将具有良好的应用前景。线阵CCD相机采用推扫方式成像，即CCD线阵沿垂直卫星飞行方向排列在焦面上，并随卫星运动完成飞行方向的扫描(参见《CBERS-01、02卫星CCD影像相对辐射校正研究》，中国科学E辑，2005年)。

目前，CCD相机探测器数目比较多，多则达几万个，少则也有几千个，而由于相机研制的技术和制作工艺的水平限制，可能会产生相机光路上的不均匀、各片CCD器件在光电响应上的不一致、各片CCD器件内部奇偶像素电路的不一致、各片CCD后端的处理电路特性上的不一致等现象，这些现象将导致线阵CCD传感器各个探测器之间存在一系列的偏差。要保证每个探元的响应一致是很难做到的，因此，探测器之间的响应差异导致成像系统推扫生成的影像在扫描方向存在灰度不均匀现象，出现明显的条带。条带现象是线阵CCD推扫式成像的一个固有问题，它将严重影响着影像数据的后续应用：卫星遥感数据应用越来越趋向于定量化应用，而条带现象严重影响地物信息的定量反演，因此保证数据的准确性是地面预处理算法的一个重要目的，在对数据进行定量分析应用之前，必须消除这种影响。

如果借助星上定标系统去除条带固然是一种原理上可行的相对辐射校正方法。然而，一方面，到目前为止，有些卫星已经取消了定标灯的使用；另一方面，星上定标系统基本不用，并且当卫星运行一段时间后，卫星携带的传感器性能将发生一定的变化，随着元器件的老化定标灯已经失效或不可靠。此时，基于遥感影像计算相对辐射校正系数的优势显得更加突出。

资料检索表明，基于遥感影像计算相对辐射校正系数所采用的方法主要有矩匹配、直方图匹配和均匀场景统计等方法。其中，均匀场景的均匀性直接影响着矩匹配法和均匀场景法的有效性，而现实中绝对均匀的地物场景是不存在的，所有用这些方法计算出的相对辐射校正系数存在误差，并且矩匹配方法要求影像的行数远远大于探测器的个数；直方图匹配法是根据数据的概率分布进行校正，达到对影像的相对辐射校正，直方图查找表的生成是进行直方图匹配处理的关键。

可以看出，上述方法对影像的要求比较苛刻，直方图均衡法和均匀场景法适应的前提条件有很大的局限性，需要的均匀场景、影像行数足够大等条件在实际应用中难以满足；而直方图匹配法假设整个影像区域灰度分布相同或相似的，对于包括不同地物的复杂地表不适用，所以直方图匹配方法会导致影像所反映的灰度信息分布发生畸变；并且基于统计的方法利用整幅影像进行校正，改变了原始影像的反射率分布信息，并且在去除条带噪声的同时带来影像细节信息的丢失。随着用户对遥感数据的定量化应用需求，在对影像进行相对辐射校正的同时，必须考虑校正精度，以保证原始数据的物理意义，并且避免影响后续的信息反演。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提出一种基于最小二乘法的遥感影像自适应辐射校正方法，该方法解决了利用统计模型对遥感影像进行相对辐射校正后所带来的灰度分布发生改变、影像细节信息丢失等问题，能有效保障影像像素间的平滑过渡、确保辐射校正精度，对处理影像的要求较传统统计方法宽松，不需要均匀场景，也不需要足够多的影像行数。

本发明的技术解决方案是：一种基于最小二乘法的线阵推扫式影像自适应辐射校正方法，通过以下步骤实现：

(1)计算线阵推扫式遥感影像的列灰度均值和标准差序列

推扫式线阵相机生成的待处理影像大小为M×N个像素，对于推扫式线阵相机而言，值N也相当于推扫式线阵相机的探元个数；

影像中第j列的灰度平均值m_j和灰度标准差v_j表示为：

$m_j={\mathrm{\Σ}}_{i=1}^M{\mathrm{DN}}_{\mathrm{ij}}/M---\left(1\right)$

$v_j=\sqrt{{\mathrm{\Σ}}_{i=1}^M{({\mathrm{DN}}_{\mathrm{ij}}-m_j)}^2/M}---\left(2\right)$

其中，i表示影像的行号，i＝1，2，...，M，M为影像的总行数，j表示影像的列号，j＝1，2，...，N，N为影像的总列数，且M与N都为正整数，DN_ij表示第i行第j列像素的灰度值；

计算待处理影像所有列的灰度平均值和灰度标准差，得到的列均值序列表示为{m₁，m₂，...，m_N}，列标准差序列表示为{v₁，v₂，...，v_N}；

(2)坏扫描列的定位与有效灰度值的计算

(2.1)遍历步骤(1)计算得到的列均值序列和列标准差序列，当正在遍历的列满足均值和标准差为零时，确定该列为坏扫描列，该列中的所有像素为坏像素；每个坏扫描列中共有M个坏像素，这些像素的灰度值为无效灰度值；

(2.2)坏扫描列定位之后，根据该列左右正常列的灰度值采用加权平均的方法重新计算坏像素的灰度值，得到坏扫描列的有效灰度值；

(2.3)重复步骤(2.1)和步骤(2.2)，直到列均值序列和列标准差序列遍历结束后，则得到所有坏扫描列去除后的影像；

(2.4)对坏扫描列去除后的影像，按照步骤(1)中所述影像列灰度均值的计算方法，重新计算坏扫描列去除后影像所有列的灰度平均值，得到新的列均值序列{m₁′，m₂′，...，m_N′}；

(3)影像中条带列的自适应判别

(3.1)从左到右依次遍历新的列均值序列，将正在遍历的列j称为当前列，j＝1，2，...，N，根据遥感影像中条带的分布特征，依据当前列及其临近列的均值大小关系来判别当前列是否为条带列；

(3.2)用符号Sign_j标识当前列是否为条带列，若当前列为条带列则将其标识Sign_j赋值为“T”，否则赋值为“F”，重复步骤(3.1)直到整个新的列均值序列都遍历完毕，得到长度为N的条带标识序列S＝{Sign₁，Sign₂，...，Sign_j，...，Sign_N}；

(4)对条带列中的所有像素进行相对辐射校正处理；

(4.1)遍历步骤(3.2)得到的标识序列S，若某列的标识符号为“T”且其左右相邻列的标识符号都为“F”，则该列称为单条带列，执行步骤(4.2.1)；若存在多个标识符号“T”连续分布在一起，则存在连续的多条带列，执行步骤(4.2.2)；

(4.2)依据相对辐射校正输入输出的线性关系，推扫式线阵相机中探元的输出表示为原始输入的函数：

DN_new＝DN_old×Gain+Offset    (3)

式中，Gain为探元的增益，Offset为探元的偏移量，即称之为相对辐射校正系数，DN_old为输入值，即探元的初始响应值，DN_new为输出值，即经辐射校正后的新值；

(4.2.1)若条带列为单一条带列，则利用最小二乘法拟合计算该条带列的相对辐射校正系数，继续步骤(4.3)；

(4.2.2)若存在连续的多条带列，列数为W，W为正整数，2≤W＜N，则利用最小二乘法计算多条带列中所有列的相对辐射校正系数，继续步骤(4.3)；

(4.3)利用相对辐射校正系数对条带列中的每个像素进行相对辐射校正，具体计算公式如下：

DN_{ij_new}＝DN_{ij_old}×Gain_j+Offset_j    (4)

其中，Gain_j和Offset_j是上述步骤计算得到的相对辐射校正系数，分别表示第j列的增益和偏移量，DN_{ij_old}是第i行第j列像素的初始输入值，DN_{ij_new}是辐射校正后的输出值，即经过辐射校正后的灰度值；

(4.4)影像列号继续增加，重复执行步骤(4.1)至(4.3)，直至所有列遍历一遍；至此，完成影像中所有条带列的判别与其对应像素的相对辐射校正，得到相对辐射校正后的影像；

(5)保存辐射校正后的影像，并可根据用户要求将计算的辐射校正系数入库，为后续快速应用提供查询数据库。

所述步骤(3.1)中判别条带列的步骤如下：

(3.1.1)当前列为第1列时，则利用第1、2列的灰度均值m₁′、m₂′与给定阈值λ的大小关系来判别当前列是否为条带列，若m₂′-m₁′≥λ，则第1列为条带列；若m₁′-m₂′≥λ，则第1列不是条带列；若|m₁′-m₂′|＜λ，则第1、2列都是条带列或都不是条带列，继续步骤(3.1.2)；

(3.1.2)当前列j不等于第1列或最后一列N时，从左到右遍历新的列均值序列，利用当前列及其左右两列的灰度均值，m_j-1′、m_j′、m_j+1′，之间的差值与λ的关系确定条带列：①若m_j-1′-m_j′≥λ且m_j+1′-m_j′≥λ，则当前列为条带列，左右两列不是条带列；②若m_j′-m_j-1′≥λ且m_j′-m_j+1′≥λ，则当前列不是条带列，而左右两列为条带列；③若m_j′-m_j-1′≥λ且m_j+1′-m_j-1′≥λ，则当前列不是条带列，而左边列为条带列；④若m_j′-m_j+1′≥λ且m_j+1′-m_j-1′≥λ，则当前列不是条带列，而右边列为条带列；⑤若m_j-1′-m_j′≥λ且m_j-1′-m_j+1′≥λ，则当前列j及其右侧列j+1为条带列，此时为多条带列分布在一起；⑥若m_j+1′-m_j′≥λ且m_j+1′-m_j-1′≥λ，则当前列j与左侧列j-1都为条带列，此时为多列条带分布在一起；⑦若|m_j+1′-m_j′|＜λ且|m_j-1′-m_j′|＜λ，则该列和左右两列全是条带或全不是条带；

(3.1.3)当前列为第N列时，利用第N-1、N列的灰度均值m_N-1′、m_N′判别第N列是否为条带列，若m_N-1′-m_N′≥λ，则第N列为条带列；若m_N′-m_N-1′≥λ，则第N列不是条带列；若|m_N-1′-m_N′|＜λ，则第N列与第N-1列一致，即当第N-1列为条带列时，则第N列也是条带列，否则当第N-1列不是条带列时，则第N列也不是条带列。

所述阈值λ也称作相对辐射校正精度，可根据用户对精度的需求进行配置，依据经验λ的取值为0.5≤λ≤2之间的实数，当要求相对辐射校正精度优于0.5个灰度值，则可设定λ为0.5，...，当要求相对辐射校正精度优于2个灰度值，则可设定λ为2。

所述步骤(4.2.1)的单一条带列的相对辐射校正系数计算步骤如下：

(4.2.1.1)当第1列是单一条带列时，对该列的每个像素建立一个线性方程，方程的输入为条带列像素的初始值，方程的输出为第2列相应像素的加权平均值；

该列第1个像素的灰度值DN₁₁作为方程一的输入值，第2列中第1、2像素的加权均值作为方程一的输出值，加权系数为[2，1]；该列第2个像素的灰度值DN₂₁作为方程二的输入值，第2列中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[1，2，1]；依次类推，...；该列第M个像素的灰度值DN_M1作为方程M的输入值，第2列中第M-1、M像素的加权均值作为方程M的输出值，加权系数为[1，2]；因此，该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{11}\×{\mathrm{Gain}}_1+{\mathrm{Offset}}_1=[2{\mathrm{DN}}_{12}+{\mathrm{DN}}_{22}]/3\\ {\mathrm{DN}}_{21}\×{\mathrm{Gain}}_1+{\mathrm{Offset}}_1=[{\mathrm{DN}}_{12}+2{\mathrm{DN}}_{22}+{\mathrm{DN}}_{32}]/4\\ ...\\ {\mathrm{DN}}_{M1}\×{\mathrm{Gain}}_1+{\mathrm{Offset}}_1=[{\mathrm{DN}}_{(M-1)2}+2{\mathrm{DN}}_{M2}]/3\end{array}\right.---\left(5\right)$

式中，Gain₁和Offset₁分别为第1列的增益和偏移量，即相对辐射校正系数；

(4.2.1.2)当第k列，k＝2，3，...，N-1，为单一条带列时，则条带列的相邻两列第k-1和第k+1列为正常列，针对条带列中的每个像素建立一个线性方程，利用条带列像素的初始值作为输入值，以条带列像素的周围正常像素的加权平均作为输出值；

该列第1个像素的灰度值DN_1k作为方程一的输入值，正常列第k-1和第k+1列中第1、2像素的加权均值作为方程一的输出值，第k-1和第k+1列中对应像素的加权系数均为[2，1]；该列第2个像素的灰度值DN_2k作为方程二的输入值，第k-1和第k+1列中第1、2、3像素的加权均值作为方程二的输出值，第k-1和第k+1列中对应像素的加权系数均为[1，2，1]；依次类推，...；该列第M个像素的灰度值DN_Mk作为方程M的输入值，第k-1和第k+1列中第M-1、M像素的加权均值作为方程M的输出值，第k-1和第k+1列中对应像素的加权系数均为[1，2]；因此，该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1k}\×{\mathrm{Gain}}_k+{\mathrm{Offset}}_k=[2{\mathrm{DN}}_{1(k-1)}+{2\mathrm{DN}}_{1(k+1)}+{\mathrm{DN}}_{2(k-1)}+{\mathrm{DN}}_{2(k+1)}]/6\\ {\mathrm{DN}}_{2k}\×{\mathrm{Gain}}_k+{\mathrm{Offset}}_k=[{\mathrm{DN}}_{1(k-1)}+{\mathrm{DN}}_{1(k+1)}+2{\mathrm{DN}}_{2(k-1)}+\\ 2{\mathrm{DN}}_{2(k+1)}+{\mathrm{DN}}_{3(k-1)}+{\mathrm{DN}}_{3(k+1)}]/8\\ ...\\ {\mathrm{DN}}_{\mathrm{Mk}}\×{\mathrm{Gain}}_k+{\mathrm{Offset}}_k=[{\mathrm{DN}}_{(M-1)(k-1)}{+\mathrm{DN}}_{(M-1)(k+1)}+{2\mathrm{DN}}_{M(k-1)}+{2\mathrm{DN}}_{M(k+1)}]/6\end{array}\right.---\left(6\right)$

式中，Gain_k和Offset_k分别为第k列的增益和偏移量，即相对辐射校正系数；

(4.2.1.3)当最后一列N是条带列时，针对该列中的每个像素建立一个线性方程，第N-1列为正常列，利用条带列像素的初始值作为方程输入值，利用第N-1列的正常像素的加权平均作为方程输出值；

该列第1个像素的灰度值DN_1N作为方程一的输入值，第N-1中第1、2像素的加权均值作为方程一的输出值，加权系数为[2，1]；该列第2个像素的灰度值DN_2N作为方程二的输入值，第N-1中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[1，2，1]；依次类推，...；该列第M个像素的灰度值DN_MN作为方程M的输入值，第N-1列中第M-1、M像素的加权均值作为方程M的输出值，加权系数均为[1，2]；因此，该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1N}\×{\mathrm{Gain}}_N+{\mathrm{Offset}}_N=[2{\mathrm{DN}}_{1(N-1)}+{\mathrm{DN}}_{2(N-1)}]/3\\ {\mathrm{DN}}_{2N}\×{\mathrm{Gain}}_N+{\mathrm{Offset}}_N=[{\mathrm{DN}}_{1(N-1)}+2{\mathrm{DN}}_{2(N-1)}+{\mathrm{DN}}_{3(N-1)}]/4\\ ...\\ {\mathrm{DN}}_{\mathrm{MN}}\×{\mathrm{Gain}}_N+{\mathrm{Offset}}_N=[{\mathrm{DN}}_{(M-1)(N-1)}+2{\mathrm{DN}}_{M(N-1)}]/3\end{array}\right.---\left(7\right)$

式中，Gain_N和Offset_N分别为第N列的增益和偏移量，即相对辐射校正系数；

(4.2.1.4)最小二乘法拟合计算辐射校正系数；

用最小二乘法拟合上述M个方程的系数，即可得到相对辐射校正系数。

所述步骤(4.2.2)的连续多条带列的相对辐射校正系数具体计算步骤如下：

(4.2.2.1)当从第1列开始连续W列为条带列，即第1列～第W列组成连续的多条带，右侧第W+1列为正常列，从左到右逐列计算相对辐射校正系数；

对于第α列，1≤α≤W，针对每个像素建立一个线性方程，利用条带列像素的初始值作为方程的输入值，以正常列即第W+1列中像素的加权均值作为方程的输出值，加权系数与距离成反比；

该列第1个像素的灰度值DN_1α作为方程一的输入值，第W+1列中第1、2像素的加权均值作为方程一的输出值，加权系数为[1+α，≈]；该列第2个像素的灰度值DN_2α作为方程二的输入值，第W+1列中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[α，1+α，α]；依次类推，...；该列第M个像素的灰度值DN_Mα作为方程M的输入值，第W+1列中第M-1、M像素的加权值作为方程M的输出值，加权系数为[α，1+α]；因此，该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1\mathrm{\α}}\×{\mathrm{Gain}}_{\mathrm{\α}}+{\mathrm{Offset}}_{\mathrm{\α}}=[(1+\mathrm{\α}){\mathrm{DN}}_{1(W+1)}+{\mathrm{\αDN}}_{2(W+1)}]/(2\mathrm{\α}+1)\\ {\mathrm{DN}}_{2\mathrm{\α}}\×{\mathrm{Gain}}_{\mathrm{\α}}+{\mathrm{Offset}}_{\mathrm{\α}}=[{\mathrm{\αDN}}_{1(W+1)}+(1+\mathrm{\α}){\mathrm{DN}}_{2(W+1)}+{\mathrm{\αDN}}_{3(W+1)}]/(3\mathrm{\α}+1)\\ ...\\ {\mathrm{DN}}_{\mathrm{M\α}}\×{\mathrm{Gain}}_{\mathrm{\α}}+{\mathrm{Offset}}_{\mathrm{\α}}=[{\mathrm{\αDN}}_{(M-1)(W+1)}+(1+\mathrm{\α}){\mathrm{DN}}_{M(W+1)}]/(2\mathrm{\α}+1)\end{array}\right.---\left(8\right)$

式中，Gain_α和Offset_α分别为第α列的增益和偏移量，即相对辐射校正系数；

(4.2.2.2)当第Q+1～Q+W共W列为连续多条带列，其中，2≤Q+1＜Q+W≤N-1，则条带列左右正常两列为第Q与Q+W+1列，每列有M个像素；从左到右逐列计算相对辐射校正系数：

对于第Q+β列条带，1≤β≤W，针对该列中的每个像素建立一个线性方程，利用条带列像素的初始值作为方程的输入值，以第Q与Q+W+1正常列像素的距离加权平均值作为方程的输出值；

该列第1个像素的灰度值DN_1(Q+β)作为方程一的输入值，第Q与Q+W+1列中第1、2像素的加权均值作为方程一的输出值，列Q中两个像素的加权系数为[W+2-β，W+1-β]，列Q+W+1中两个像素的加权系数为[1+β，β]；该列第2个像素的灰度值DN_2(Q+β)作为方程二的输入值，第Q与Q+W+1列中第1、2、3像素的加权均值作为方程二的输出值，列Q中三个像素的加权系数为[W+1-β，W+2-β，W+1-β]，列Q+W+1中三个像素的加权系数为[β，1+β，β]；依次类推，...；该列第M个像素的灰度值DNM_(Q+β)作为方程M的输入值，第Q与Q+W+1列中第M-1、M像素的加权均值作为方程M的输出值，列Q中两个像素的加权系数为[W+1-β，W+2-β]，列Q+W+1中两个像素的加权系数为[1+β，β]；因此，该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1(Q+\mathrm{\β})}\×{\mathrm{Gain}}_{(Q+\mathrm{\β})}+{\mathrm{Offset}}_{(Q+\mathrm{\β})}=[(W+2-\mathrm{\β}){\mathrm{DN}}_{1Q}+(W+1-\mathrm{\β}){\mathrm{DN}}_{2Q}+\\ (1+\mathrm{\β}){\mathrm{DN}}_{1(Q+W+1)}+\mathrm{\β}{\mathrm{DN}}_{2(Q+W+1)}]/(2W+4)\\ {\mathrm{DN}}_{2(Q+\mathrm{\β})}\×{\mathrm{Gain}}_{(Q+\mathrm{\β})}+{\mathrm{Offset}}_{(Q+\mathrm{\β})}=[(W+1-\mathrm{\β}){\mathrm{DN}}_{1Q}+(W+2-\mathrm{\β}){\mathrm{DN}}_{2Q}+\\ (W+1-\mathrm{\β}){\mathrm{DN}}_{3Q}+\mathrm{\β}{\mathrm{DN}}_{1(Q+W+1)}+(1+\mathrm{\β}){\mathrm{DN}}_{2(Q+W+1)}+\mathrm{\β}{\mathrm{DN}}_{3(Q+W+1)}]/(3W+5)\\ ...\\ {\mathrm{DN}}_{M(Q+\mathrm{\β})}\×{\mathrm{Gain}}_{(Q+\mathrm{\β})}+{\mathrm{Offset}}_{(Q+\mathrm{\β})}=[(W+1-\mathrm{\β}){\mathrm{DN}}_{(M-1)Q}+(W+2-\mathrm{\β}){\mathrm{DN}}_{\mathrm{MQ}}+\\ (1+\mathrm{\β}){\mathrm{DN}}_{(M-1)(Q+W+1)}+\mathrm{\β}{\mathrm{DN}}_{M(Q+W+1)}]/(2W+4)\end{array}\right.---\left(9\right)$

式中，Gain_(Q+β)和Offset_(Q+β)分别为第(Q+β)列的增益和偏移量，即相对辐射校正系数；

(4.2.2.3)当从第N-W+1列开始到最后一列N，共连续W列为条带列，即第N-W+1、N-W+2、...、N共W列组成连续的多条带，左侧第N-W列为正常列，从左到右逐列计算相对辐射校正系数；

对于第N-γ列，N-W+1≤N-γ≤N，针对该列中的每个像素建立一个线性方程，利用条带列像素的初始值作为方程的输入值，以正常列像素的加权均值作为方程的输出值，加权系数与距离成反比；

该列第1个像素的灰度值DN_1(N-γ)作为方程一的输入值，正常列第N-W列中第1、2像素的加权均值作为方程一的输出值，加权系数为[W-γ，W-γ-1]；该列第2个像素的灰度值DN_2(N-γ)作为方程二的输入值，第N-W列中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[W-γ-1，1 W-γ，W-γ-1]；依次类推，...；该列第M个像素的灰度值DN_M(N-γ)作为方程M的输入值，第N-W列中第M-1、M像素的加权值作为方程M的输出值，加权系数为[W-γ-1，W-γ]；因此，该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1(N-\mathrm{\γ})}\×{\mathrm{Gain}}_{(N-\mathrm{\γ})}+{\mathrm{Offset}}_{(N-\mathrm{\γ})}=[(W-\mathrm{\γ}){\mathrm{DN}}_{1(N-W)}+\\ (W-\mathrm{\γ}-1){\mathrm{DN}}_{2(N-W)}]/(2W-2\mathrm{\γ}-1)\\ {\mathrm{DN}}_{2(N-\mathrm{\γ})}\×{\mathrm{Gain}}_{(N-\mathrm{\γ})}+{\mathrm{Offset}}_{(N-\mathrm{\γ})}=[(W-\mathrm{\γ}-1){\mathrm{DN}}_{1(N-W)}+(W+\mathrm{\γ}){\mathrm{DN}}_{2(N-W)}+\\ (W-\mathrm{\γ}-1){\mathrm{DN}}_{3(N-W)}]/(3W-2\mathrm{\γ}-2)\\ ...\\ {\mathrm{DN}}_{M(N-\mathrm{\γ})}\×{\mathrm{Gain}}_{(N-\mathrm{\γ})}+{\mathrm{Offset}}_{(N-\mathrm{\γ})}=[(W-\mathrm{\γ}-1){\mathrm{DN}}_{(M-1)(N-W)}+\\ (W-\mathrm{\γ}){\mathrm{DN}}_{M(N-W)}]/(2W-2\mathrm{\γ}-1)\end{array}\right.---\left(10\right)$

式中，Gain_(N-γ)和OFfset_(N-γ)分别为第(N-γ)列的增益和偏移量，即相对辐射校正系数；

所述步骤(5)的影像数据与辐射校正系数的存储方式如下：将影像中所有条带列的相对辐射校正结果以实数形态单独存储成一个附加文件，在后续应用过程中可以直接读取附加文件中的数据代替对应列的影像数据，为后续处理提供更精确的数据；同时，将计算得到的相对辐射校正系数存储到数据库中，在后续应用过程中可利用事先存储的辐射校正系数对输入的影像数据进行处理。

本发明与现有技术相比的有益效果是：

(1)本发明针对目前地面处理系统中遥感影像相对辐射校正处理采用固定校正模型的缺点。根据卫星成像过程中线阵相机的物理特性，分析并确定条带分布情况，并根据分布情况提出自动检测列条带方法；较之现有地面应用系统中传统的相对辐射校正方法，能够根据用户对精度的需求对影像进行自适应辐射校正，只改变条带所在列像素的灰度值，对正常列的灰度值不做变化。

(2)本发明面向用户对遥感数据定量化应用的需求，为保证遥感数据辐射校正处理的精度以及信息反演的有效性，需要确保遥感数据的物理意义不发生畸变。本发明针对条带列中所有像素，利用最小二乘法拟合输入与输出的线性关系，计算增益和偏移量校正系数并对条带列像素进行校正。一方面，充分利用条带列的探测器成像信息作为输入，在此基础上进行校正；另一方面，从提高辐射校正精度角度考虑，克服了传统地面系统的相对辐射处理方法更多关注影像处理，而忽视辐射校正所带来的对影像灰度信息的改变，进一步导致遥感数据物理意义的扭曲。

(3)本发明的辐射校正以提高辐射精度为基本目的，一方面辐射校正精度可配置，另一方面利用实数数据类型保存影像数据，避免或减少了相对辐射校正处理结果的取整运算所带来的误差，降低了信息量丢失对影像质量的影响。

附图说明

图1是本发明基于最小二乘法的线阵推扫式影像自适应辐射校正流程；

图2(a)是本发明坏扫描列分布情况的示意图，图2(b)是对坏像素处理的加权平均模板；

图3是本发明线阵推扫式遥感影像中条带的可能分布情况，其中，图3(a)多列影像中可能存在的条带，图3(b1)～(b8)为影像中条带可能的分布情况；

图4(a)是本发明中单像素宽度条带中各个像素的示意图，图4(b)是本发明中连续多像素宽度条带中各个像素的示意图，图4(c)是对条带像素处理的加权平均模板；

图5为本发明算法、矩匹配法处理结果。其中，图5(a1)～(a3)分别为原始影像、本发明算法的处理结果、整幅影像的差值影像，图5(b1)～(b3)与图5(c1)～(c3)分别为区域A与区域B对应的原始影像、本发明算法的处理结果、差值影像，图5(d1)与图5(d2)表示矩匹配算法处理结果及其差值影像。

图6是原始影像列均值曲线图，其中，图6(a)表示列灰度均值曲线，图6(b)表示探测器序号从2000～3000所对应的影像列灰度均值曲线，而6(c)表示相邻两列均值差值的绝对值曲线；

图7是基于本发明算法处理后的列均值曲线图，其中，图7(a)表示相对辐射校正后的列灰度均值曲线，图7(b)表示探测器序号从2000～3000所对应的校正影像列灰度均值曲线，而7(c)表示校正后相邻两列均值差值的绝对值曲线。

具体实施方式

本发明基于最小二乘法的线阵推扫式影像自适应辐射校正方法，具体处理流程如图1所示。

1、计算线阵推扫式遥感影像的列灰度均值和标准差序列

推扫式线阵相机生成的待处理影像大小为M×N个像素，M为影像的总行数，N为影像的总列数，且M与N都为正整数，对于推扫式线阵相机而言，值N也相当于推扫式线阵相机的探元个数。影像中第j列的灰度平均值m_j和灰度标准差v_j可表示为：

$m_j={\mathrm{\Σ}}_{i=1}^M{\mathrm{DN}}_{\mathrm{ij}}/M---\left(1\right)$

$v_j=\sqrt{{\mathrm{\Σ}}_{i=1}^M{({\mathrm{DN}}_{\mathrm{ij}}-m_j)}^2/M}---\left(2\right)$

其中，i表示影像的行号，i＝1，2，...，M，j表示影像的列号，j＝1，2，...，N，DN_ij表示第i行第j列像素的灰度值。

计算待处理影像所有列的灰度平均值和灰度标准差，得到的列均值序列表示为{m₁，m₂，...，m_N}，列标准差序列表示为{v₁，v₂，...，v_N}。

2、坏扫描列的定位与有效灰度值的计算

坏扫描列的产生是由于线阵的探测器失效，对于线阵相机来说主要表现在推扫式影像上会出现整列方向没有光谱信息，其灰度值为零。这是个比较严重的问题，因为没有获取的数据是不能恢复的。但是，为提高影像辐射质量必须对坏扫描列进行灰度插值处理。

2.1、待处理影像中的坏扫描列的定位

遍历步骤1计算得到的列均值序列和列标准差序列，当正在遍历的列满足均值和标准差为零时，可确定该列为坏扫描列，对应的探测器为无效探测器，从而该列中的所有像素为坏像素，每个坏扫描列中共有M个坏像素，这些像素的灰度值为无效灰度值。

需要说明的是，地面场景可能存在例如水之类的均匀物体，理论上整条影像列只输出一个灰度值；但是实际情况是，一方面地面均匀场景几乎不存在；另一方面，这种情况即使存在，经过成像系统后影像中会叠加加性和乘性噪声、以及成像过程大气的影响，因此导致不可能存在列灰度值全是零的实际影像。

2.2、坏扫描列中有效灰度值的计算

坏扫描列一旦定位之后，就可以根据该列左右正常列的像素灰度值确定坏扫描列的灰度值。为更好地利用景物在空间分布的相关性，本发明采用对坏像素周边的正常像素进行加权平均的方法重新计算坏扫描列中所有像素的灰度值，加权系数的大小与坏像素与正常像素的距离成反比。

(1)当坏扫描列不是影像的第1列和最后一列即第N列时，如图2(a)所示，假设通过上述方法探测到第L列～第L+D列共D+1列连续的坏扫描列，其中，2≤L＜L+D≤N-1，则坏扫描列左右正常两列为第L-1与L+D+1列。对于L+δ列坏扫描列中，0≤δ≤D，对应的坏像素计算方法如下：

①当该列的坏像素既不在第一行又不在最后一行时，即1＜K＜M，则加权模板如图2(b)所示，正常列L-1与L+D+1中的系数分别为[D+1-δ，D+2-δ，D+1-δ]和[1+δ，2+δ，1+δ]。如图2(a)所示，第L-1和第L+D+1列中对应的第K-1，K，K+1行像素的灰度值分别为DN₁、DN₂、DN₃和DN₄、DN₅、DN₆，因此，坏像素DN_K(L+δ)的值等于模板中的系数与对应像素灰度值DN₁～DN₆加权求和，再除以加权系数总和(3D+8)后得到的结果。

②当K＝1或M时，即在处理第一行或最后一行的坏像素时，则模板中的第一行或最后一行将不在影像内，因此，此时仅利用模板的后两行或前两行系数与对应像素灰度值加权求和。

(2)当坏扫描列是影像的第一列或最后一列时，则模板中的第一列或最后一列将不在图像内，因此，此时仅利用模板的第三列(对应扫描列是第一列)或第一列(对应扫描列是最后一列)系数与对应像素值加权求和。

2.3、重复步骤2.1和步骤2.2，直到列均值序列和列标准差序列遍历结束后，则得到所有坏扫描列去除后的影像。

2.4、针对坏扫描列去除后的影像，重新计算坏扫描列去除后影像所有列的灰度平均值，得到新的列均值序列{m₁′，m₂′，...，m_N′}。

3、影像中条带列的自适应判别

推扫式相机在成像过程中由于相机光路的不均匀、以及各片CCD器件的光电响应、内部奇偶像素电路与后端的处理电路特性上会存在不一致现象，这些现象都是导致推扫式影像产生条带的原因。依据条带产生的原因以及在轨推扫式相机影像数据的经验分析，具有如下结论：

a.当内部采用奇偶像素电路的情况下，推扫式影像中的条带往往呈现奇偶分布，因此，该种情况下条带往往表现为单列宽度；

b.各个探测器及放大增益存在不均匀也会导致条带的产生，因此可能导致影像中存在多列连续条带的情况；

c.遥感影像中条带所在列的灰度值一般偏低。

推扫式影像中可能存在的条带情况如示意图3(a)所示，可以看出，图3(b1)中条带分布在中间列，图3(b2)中条带分布在左右两列，图3(b3)中条带分布在左列，图3(b4)中条带分布在右列，图3(b5)中右边两列为条带，图3(b6)中左边两列为条带，图3(b7)中无条带，图3(b8)中三列全是条带列。

3.1、针对待处理的影像，从左到右依次遍历影像新的列均值序列，将正在遍历的列j称为当前列，j＝1，2，...，N，根据遥感影像中条带的分布特征，依据当前列及其临近列的均值大小来判别当前列是否为条带列，并用符号Sign_j标识，若当前列为条带列则其标识Sign_j赋值为“T”，否则赋值为“F”。判别条带列的具体步骤如下：

(1)当前列为第1列时，利用第1、2列的灰度均值m₁、m₂判别当前列是否为条带列，若m₂′-m₁′≥λ，则第1列为条带列；若m₁′-m₂′≥λ，则第1列不是条带列；若|m₁′-m₂′|＜λ，则第1、2列都是条带列或都不是条带列，继续利用步骤(2)判别。

其中，阈值λ表示相对辐射校正精度，可根据用户需求进行配置，依据经验λ的取值一般为0.5≤λ≤2之间的实数，当要求相对辐射校正精度优于0.5个灰度值，则可设定λ为0.5，...，要求相对辐射校正精度优于2个灰度值，则可设定λ为2依次类推，以下λ的取值与此相同。

(2)当前列j不等于第1列或最后一列N时，从左到右遍历新的列灰度均值序列，利用当前列及其左右两列的灰度均值，m_j-1′、m_j′、m_j+1′，之间的差值与给定λ的关系确定条带列：①若m_j-1′-m_j′≥λ且m_j+1′-m_j′≥λ，则当前列为条带列，左右两列不是条带列，对应图3(b1)；②若m_j′-m_j-1′≥λ且m_j′-m_j+1′≥λ，则当前列不是条带列，而左右两列为条带列，对应图3(b2)；③若m_j′-m_j-1′≥λ且m_j+1′-m_j-1′≥λ，则当前列不是条带列，而左边列为条带列，对应图3(b3)；④若m_j′-m_j+1′≥λ且m_j+1′-m_j-1′≥λ，则当前列不是条带列，而右边列为条带列，对应图3(b4)；⑤若m_j-1′-m_j′≥λ且m_j-1′-m_j+1′≥λ，则当前列j及其右侧列j+1为条带列，此时为多条带列分布在一起，对应图3(b5)；⑥若m_j+1′-m_j′≥λ且m_j+1′-m_j-1′≥λ，则当前列j与左侧列j-1都为条带列，此时为多列条带分布在一起，对应图3(b6)；⑦若|m_j+1′-m_j′|＜λ且|m_j-1′-m_j′|＜λ，则该列和左右两列全是条带或全不是条带，对应图3(b7)和(b8)。

(3)当前列为第N列时，利用第N-1、N列的灰度均值m_N-1′、m_N′判别当前列是否为条带列，若m_N-1′-m_N′≥λ，则第N列为条带列；若m_N′-m_N-1′≥λ，则第N列不是条带列；若|m_N-1′-m_N′|＜λ，则第N列与第N-1列一致，即当第N-1列为条带列时，则第N列也是条带列，否则当第N-1列不是条带列时，则第N列也不是条带列。

3.2、重复步骤3.1直到整个新的列均值序列都遍历完毕，则得到长度为N的条带标识序列S＝{Sign₁，Sign₂，...，Sign_j，...，Sign_N}。

4、对条带列中的像素进行相对辐射校正处理

当通过上述判别方式识别得到条带列之后，则需要对此列中的所有像素进行处理，以消除与两侧正常列像素的差异。

根据线阵相机CCD制造厂商提供的特性指标，以及实际实验数据的分析可知，线阵CCD探测器的光电响应在很大的输入范围为线性关系，因此，确定线阵CCD探测器的模型为一次线性模型，对该探测器对应的列像素需要建立辐射校正模型进行校正。依据线性模型的结论，则推扫式线阵相机中探元的输出DN_new可写成原始输入DN_old的函数：

DN_new＝DN_old×Gain+Offset    (3)

式中，Gain为探元的增益，Offset为探元的偏移量，即称之为相对辐射校正系数，DN_old为输入值，即探元的初始响应值，DN_new为输出值，即经辐射校正后的新值。

如果条带列中的所有像素都按照上述公式进行计算处理，则完成条带列像素的相对辐射校正。基于这种方法，只要使用一定的方法确定出相对辐射校正系数即可。

4.1、遍历步骤3.2得到的标识序列S，若某列的标识符号为“T”且其左右列的标识符号都为“F”，则该列称为单一条带列，执行步骤4.2.1；若存在多个标识符号“T”连续分布在一起，则存在连续的多条带列，执行步骤4.2.2。

4.2、计算条带列的相对辐射校正系数

本发明充分利用条带列周边的正常列像素的灰度信息计算相对辐射校正系数，从而根据校正系数对条带列中的像素作变换以消除条带现象。根据一次线性辐射校正模型，可对每个需要校正的像素建立线性方程。

4.2.1、若条带列为单一条带列，则利用最小二乘法拟合计算该条带列的相对辐射校正系数，具体处理步骤如下：

(1)当第1列是条带列时，对该列的每个像素建立一个线性方程，方程的输入为条带列像素的初始值，方程的输出为第2列对应像素的加权平均值：

该列第1个像素的灰度值DN₁₁作为方程一的输入值，第2列中第1、2像素的加权均值作为方程一的输出值，加权系数为[2，1]；该列第2个像素的灰度值DN₂₁作为方程二的输入值，第2列中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[1，2，1]；依次类推，...；该列第M个像素的灰度值DNM1作为方程M的输入值，第2列中第M-1、M像素的加权均值作为方程M的输出值，加权系数为[1，2]；因此，该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{11}\×{\mathrm{Gain}}_1+{\mathrm{Offset}}_1=[2{\mathrm{DN}}_{12}+{\mathrm{DN}}_{22}]/3\\ {\mathrm{DN}}_{21}\×{\mathrm{Gain}}_1+{\mathrm{Offset}}_1=[{\mathrm{DN}}_{12}+2{\mathrm{DN}}_{22}+{\mathrm{DN}}_{32}]/4\\ ...\\ {\mathrm{DN}}_{M1}\×{\mathrm{Gain}}_1+{\mathrm{Offset}}_1=[{\mathrm{DN}}_{(M-1)2}+2{\mathrm{DN}}_{M2}]/3\end{array}\right.---\left(4\right)$

式中，Gain₁和Offset₁分别为第1列的增益和偏移量，即相对辐射校正系数。

(2)当第k列为条带列时，且k＝2，3，...，N-1，如图4(a)所示，条带列的相邻两列第k-1和第k+1列为正常列，针对条带列中的每个像素建立一个线性方程，利用条带列像素的初始值作为输入值，以条带列像素的周围正常像素的加权平均结果作为输出值：

该列第1个像素的灰度值DN_1k作为方程一的输入值，正常列第k-1和第k+1列中第1、2像素的加权均值作为方程一的输出值，第k-1和第k+1列中对应像素的加权系数均为[2，1]；该列第2个像素的灰度值DN_2k作为方程二的输入值，第k-1和第k+1列中第1、2、3像素的加权均值作为方程二的输出值，第k-1和第k+1列中对应像素的加权系数均为[1，2，1]；依次类推，...；该列第M个像素的灰度值DN_MK作为方程M的输入值，第k-1和第k+1列中第M-1、M像素的加权均值作为方程M的输出值，第k-1和第k+1列中对应像素的加权系数均为[1，2]；因此，M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1k}\×{\mathrm{Gain}}_k+{\mathrm{Offset}}_k=[2{\mathrm{DN}}_{1(k-1)}+{2\mathrm{DN}}_{1(k+1)}+{\mathrm{DN}}_{2(k-1)}+{\mathrm{DN}}_{2(k+1)}]/6\\ {\mathrm{DN}}_{2k}\×{\mathrm{Gain}}_k+{\mathrm{Offset}}_k=[{\mathrm{DN}}_{1(k-1)}+{\mathrm{DN}}_{1(k+1)}+2{\mathrm{DN}}_{2(k-1)}+\\ 2{\mathrm{DN}}_{2(k+1)}+{\mathrm{DN}}_{3(k-1)}+{\mathrm{DN}}_{3(k+1)}]/8\\ ...\\ {\mathrm{DN}}_{\mathrm{Mk}}\×{\mathrm{Gain}}_k+{\mathrm{Offset}}_k=[{\mathrm{DN}}_{(M-1)(k-1)}{+\mathrm{DN}}_{(M-1)(k+1)}+{2\mathrm{DN}}_{M(k-1)}+{2\mathrm{DN}}_{M(k+1)}]/6\end{array}\right.---\left(5\right)$

式中，Gain_k和Offset_k分别为第k列的增益和偏移量，即相对辐射校正系数。

(3)当最后一列N是条带列时，针对该列中的每个像素建立一个线性方程，第N-1列为正常列，利用条带列像素的初始值作为方程输入值，利用第N-1列的正常像素的加权平均作为方程输出值：

该列第1个像素的灰度值DN_1N作为方程一的输入值，第N-1中第1、2像素的加权均值作为方程一的输出值，加权系数均为[2，1]；该列第2个像素的灰度值DN_2N作为方程二的输入值，第N-1中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[1，2，1]；依次类推，...；该列第M个像素的灰度值DN_MN作为方程M的输入值，第N-1列中第M-1、M像素的加权均值作为方程M的输出值，加权系数均为[1，2]；因此，该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1N}\×{\mathrm{Gain}}_N+{\mathrm{Offset}}_N=[2{\mathrm{DN}}_{1(N-1)}+{\mathrm{DN}}_{2(N-1)}]/3\\ {\mathrm{DN}}_{2N}\×{\mathrm{Gain}}_N+{\mathrm{Offset}}_N=[{\mathrm{DN}}_{1(N-1)}+2{\mathrm{DN}}_{2(N-1)}+{\mathrm{DN}}_{3(N-1)}]/4\\ ...\\ {\mathrm{DN}}_{\mathrm{MN}}\×{\mathrm{Gain}}_N+{\mathrm{Offset}}_N=[{\mathrm{DN}}_{(M-1)(N-1)}+2{\mathrm{DN}}_{M(N-1)}]/3\end{array}\right.---\left(6\right)$

式中，Gain_N和Offset_N分别为第N列的增益和偏移量，即相对辐射校正系数。

(4)最小二乘法拟合计算辐射校正系数

用最小二乘法拟合上述M个方程的系数，即可得到相对辐射校正系数。

4.2.2、若存在连续的多条带列，列数为W，W为正整数，2≤W＜N，则利用最小二乘法计算多条带列中所有列的相对辐射校正系数

(1)当从第1列开始连续W列为条带列，即第1列～第W列组成连续的多条带，右侧第W+1列为正常列。

从左到右逐列计算相对辐射校正系数的方法如下：

对于第α列，1≤α≤W，针对该列中的每个像素建立一个线性方程，利用条带列像素的初始值作为方程的输入值，以正常列像素的加权均值作为方程的输出值，加权系数与距离成反比：

该列第1个像素的灰度值DN_1α作为方程一的输入值，第W+1列中第1、2像素的加权均值作为方程一的输出值，加权系数为[1+α，α]；该列第2个像素的灰度值DN_2α作为方程二的输入值，第W+1列中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[α，1+α，α]；依次类推，...；该列第M个像素的灰度值DN_Mα作为方程M的输入值，第W+1列中第M-1、M像素的加权值作为方程M的输出值，加权系数为[α，1+α]；因此，该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1\mathrm{\α}}\×{\mathrm{Gain}}_{\mathrm{\α}}+{\mathrm{Offset}}_{\mathrm{\α}}=[(1+\mathrm{\α}){\mathrm{DN}}_{1(W+1)}+{\mathrm{\αDN}}_{2(W+1)}]/(2\mathrm{\α}+1)\\ {\mathrm{DN}}_{2\mathrm{\α}}\×{\mathrm{Gain}}_{\mathrm{\α}}+{\mathrm{Offset}}_{\mathrm{\α}}=[{\mathrm{\αDN}}_{1(W+1)}+(1+\mathrm{\α}){\mathrm{DN}}_{2(W+1)}+{\mathrm{\αDN}}_{3(W+1)}]/(3\mathrm{\α}+1)\\ ...\\ {\mathrm{DN}}_{\mathrm{M\α}}\×{\mathrm{Gain}}_{\mathrm{\α}}+{\mathrm{Offset}}_{\mathrm{\α}}=[{\mathrm{\αDN}}_{(M-1)(W+1)}+(1+\mathrm{\α}){\mathrm{DN}}_{M(W+1)}]/(2\mathrm{\α}+1)\end{array}\right.---\left(7\right)$

式中，Gain_α和Offset_α分别为第α列的增益和偏移量，即相对辐射校正系数。

(2)当第Q+1～Q+W共W列为连续多条带列，其中，2≤Q+1＜Q+W≤N-1，如图4(b)所示，条带列左右正常两列为第Q与Q+W+1列，每列有M个像素。

从左到右逐列计算相对辐射校正系数的方法如下：

对于第Q+β列条带，1≤β≤W，针对该列中的每个像素建立一个线性方程，利用条带列像素的初始值作为方程的输入值，以第Q与Q+W+1正常列像素的距离加权平均值作为方程的输出值，加权模板与系数如图4(c)所示：

该列第1个像素的灰度值DN_1(Q+β)作为方程一的输入值，第Q与Q+W+1列中第1、2像素的加权均值作为方程一的输出值，列Q中两个像素的加权系数为[W+2-β，W+1-β]，列Q+W+1中两个像素的加权系数为[1+β，β]；该列第2个像素的灰度值DN_2(Q+β)作为方程二的输入值，第Q与Q+W+1列中第1～3像素的加权均值作为方程二的输出值，列Q中三个像素的加权系数为[W+1-β，W+2-β，W+1-β]，列Q+W+1中三个像素的加权系数为[β，1+β，β]；依次类推，...；该列第M个像素的灰度值DN_M(Q+β)作为方程M的输入值，第Q与Q+W+1列中第M-1、M像素的加权均值作为方程M的输出值，列Q中两个像素的加权系数为[W+1-β，W+2-β]，列Q+W+1中两个像素的加权系数为[1+β，β]；因此，该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1(Q+\mathrm{\β})}\×{\mathrm{Gain}}_{(Q+\mathrm{\β})}+{\mathrm{Offset}}_{(Q+\mathrm{\β})}=[(W+2-\mathrm{\β}){\mathrm{DN}}_{1Q}+(W+1-\mathrm{\β}){\mathrm{DN}}_{2Q}+\\ (1+\mathrm{\β}){\mathrm{DN}}_{1(Q+W+1)}+\mathrm{\β}{\mathrm{DN}}_{2(Q+W+1)}]/(2W+4)\\ {\mathrm{DN}}_{2(Q+\mathrm{\β})}\×{\mathrm{Gain}}_{(Q+\mathrm{\β})}+{\mathrm{Offset}}_{(Q+\mathrm{\β})}=[(W+1-\mathrm{\β}){\mathrm{DN}}_{1Q}+(W+2-\mathrm{\β}){\mathrm{DN}}_{2Q}+\\ (W+1-\mathrm{\β}){\mathrm{DN}}_{3Q}+\mathrm{\β}{\mathrm{DN}}_{1(Q+W+1)}+(1+\mathrm{\β}){\mathrm{DN}}_{2(Q+W+1)}+\mathrm{\β}{\mathrm{DN}}_{3(Q+W+1)}]/(3W+5)\\ ...\\ {\mathrm{DN}}_{M(Q+\mathrm{\β})}\×{\mathrm{Gain}}_{(Q+\mathrm{\β})}+{\mathrm{Offset}}_{(Q+\mathrm{\β})}=[(W+1-\mathrm{\β}){\mathrm{DN}}_{(M-1)Q}+(W+2-\mathrm{\β}){\mathrm{DN}}_{\mathrm{MQ}}+\\ (1+\mathrm{\β}){\mathrm{DN}}_{(M-1)(Q+W+1)}+\mathrm{\β}{\mathrm{DN}}_{M(Q+W+1)}]/(2W+4)\end{array}\right.---\left(8\right)$

式中，Gain_(Q+β)和Offset_(Q+β)分别为第(Q+β)列的增益和偏移量，即相对辐射校正系数。

(3)当从第N-W-1列开始到最后一列N共连续W列为条带列，即第N-W+1、N-W+2、...、N共W列组成连续的多条带，左侧第N-W列为正常列，从左到右逐列计算相对辐射校正系数；

对于第N-γ列，N-W+1≤N-γ≤N，针对该列中的每个像素建立一个线性方程，利用条带列像素的初始值作为方程的输入值，以正常列像素的加权均值作为方程的输出值，加权系数与距离成反比：

该列第1个像素的灰度值DN_1(N-γ)作为方程一的输入值，正常列第N-W列中第1、2像素的加权均值作为方程一的输出值，加权系数为[W-γ，W-γ-1]；该列第2个像素的灰度值DN_2(N-γ)作为方程二的输入值，第N-W列中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[W-γ-1，1 W-γ，W-γ-1]；依次类推，...；该列第M个像素的灰度值DN_M(N-γ)作为方程M的输入值，第N-W列中第M-1、M像素的加权值作为方程M的输出值，加权系数为[W-γ-1，W-γ]；因此，该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1(N-\mathrm{\γ})}\×{\mathrm{Gain}}_{(N-\mathrm{\γ})}+{\mathrm{Offset}}_{(N-\mathrm{\γ})}=[(W-\mathrm{\γ}){\mathrm{DN}}_{1(N-W)}+\\ (W-\mathrm{\γ}-1){\mathrm{DN}}_{2(N-W)}]/(2W-2\mathrm{\γ}-1)\\ {\mathrm{DN}}_{2(N-\mathrm{\γ})}\×{\mathrm{Gain}}_{(N-\mathrm{\γ})}+{\mathrm{Offset}}_{(N-\mathrm{\γ})}=[(W-\mathrm{\γ}-1){\mathrm{DN}}_{1(N-W)}+(W+\mathrm{\γ}){\mathrm{DN}}_{2(N-W)}+\\ (W-\mathrm{\γ}-1){\mathrm{DN}}_{3(N-W)}]/(3W-2\mathrm{\γ}-2)\\ ...\\ {\mathrm{DN}}_{M(N-\mathrm{\γ})}\×{\mathrm{Gain}}_{(N-\mathrm{\γ})}+{\mathrm{Offset}}_{(N-\mathrm{\γ})}=[(W-\mathrm{\γ}-1){\mathrm{DN}}_{(M-1)(N-W)}+\\ (W-\mathrm{\γ}){\mathrm{DN}}_{M(N-W)}]/(2W-2\mathrm{\γ}-1)\end{array}\right.---\left(9\right)$

式中，Gain_(N-γ)和Offset_(N-γ)分别为第(N-γ)列的增益和偏移量，即相对辐射校正系数。

4.3、利用相对辐射校正系数对条带列中的所有像素进行处理

利用最小二乘法拟合的校正系数，对条带列中的每个像素进行相对辐射校正，达到恢复影像目的，具体计算公式如下：

DN_{ij_new}＝DN_{ij_old}×Gain_j+Offset_j     (10)

其中，Gain_j和Offset_j是上述步骤计算得到的相对辐射校正系数，分别表示第j列的增益和偏移量，DN_{ij_old}是第i行第j列像素的初始输入值，DN_{ij_new}是辐射校正后的输出值，即经过辐射校正后的灰度值。

本发明的辐射校正以高精度为基本目的，将计算得到的灰度值以实数数据类型保存，避免或减少了相对辐射校正处理过程中的四舍五入取整运算所带来的误差，降低了信息量丢失对影像质量的影响。

4.4、影像列号继续增加，重复执行步骤4.1至4.3，直至所有列遍历一遍；至此，完成影像中所有条带列的判别与其对应像素的相对辐射校正，得到相对辐射校正后的影像。

5、保存辐射校正后的影像，并可根据用户要求将计算的辐射校正系数入库

影像数据的每个像素一般以正整数形式存储，因此，原则上需要将各个像素的相对辐射校正结果取整。为了为后续处理步骤提供更精确的数据，本发明中将所有条带列的相对辐射校正结果以实数形态单独存储成一个附加文件，在后续应用过程中可以用附加文件中的数据代替对应列的影像数据。

另一方面，将计算的相对辐射校正系数存储到数据库中，在后续应用过程中可利用事先生成的辐射校正模型对输入数据进行处理，可以实现对由于相机CCD及其电路特性的不均匀性造成的影像列方向条带现象进行校正。

下面结合具体实例对本发明算法的处理效果做说明。

以某推扫式成像卫星为例，线阵CCD阵列共有12000个探测器，如图5(a)所示为大小为12000×12000像素的推扫影像(显示比例3％)。下面以如图5(a)所示的遥感影像为例，说明本发明的相对辐射校正过程。

1、列灰度均值与方差计算

为了突出显示效果，在图5(a1)中选择了两个区域A和B，大小100×100像素，其中，A区域为薄云覆盖的区域(左上角坐标为(4450，6150))，而B区域为包含丰富地物信息的区域(左上角坐标为(10800，6500))，图5(b1)和图5(c1)分别为A和B区域缩放比例达到300％时的放大图。首先，通过目视效果可以看出，影像中条带现象明显；计算影像列灰度均值，图6(a)所示为图5(a)影像的列均值曲线图，图6(b)所示为局部放大曲线图，而图6(c)所示为相邻两列均值之差绝对值对应的曲线图，图中显示，大部分差值主要集中在1.0与2.0之间。可以看出，该影像中存在明显的条带噪声，影像的列方向均值变化剧烈，列之间的灰度均值存在很大的偏差，影响了影像的辐射质量。

2、坏扫描列的定位与去除

根据本发明提出的坏扫描列定位方法以及计算出的标准差序列，可得知图5(a)影像中不存在坏扫描列。

3、影像中条带列的自适应判别

根据本发明中提出的条带列判别原则，针对影像5(a1)，从左到右以“每三列”为判别单位，依次遍历列均值序列。对本实施例，设定λ为0.5。

4、影像的相对辐射校正

对判别出的条带列利用本发明提出的最小二乘法拟合相对辐射校正系数，并完成每个像素的校正。以图5(a1)中区域A为例具体说明校正效果。表1和表2所示分布为区域A中的局部区域在校正前后灰度信息：

表1局部区域校正前灰度信息

表2局部区域校正后的灰度信息

从表1与表2可以看出，非条带列信息没发生变化，条带列通过最小二乘法拟合出的辐射校正系数校正后其值发生了变化，如表3所示，可以看出校正后数据输出基本一致，有效地消除了探测器的响应不一致现象，说明校正数据及其校正算法有效。

图5(a2)是经过本发明算法处理后的效果图，图5(b2)和图5(c2)分别对应处理后的区域A和B；而图5(a3)为原始影像与处理结果的差值影像(注：为突出显示效果，本发明提供的所有差值结果图都是在实际差值基础上加上一个大小为20的附加值，以下同)，图5(b3)和图(c3)分别对应差值影像中的区域A和区域B。图5(d1)为矩匹配算法处理后的效果图，而图5(d2)为原始影像与矩匹配方法结果的差值图。通过图5中的目视效果可以看出：本发明方法在处理前后信息量丢失少，仅仅去掉了影像中的条带噪声，而其他部分灰度信息并没发生大的变化，从而保证影像细节未丢失；而矩匹配法虽然去掉了影像条带噪声，但是大量的原始影像中的信息被丢失。

图7为利用本发明方法处理后的影像列方向灰度均值和相邻列灰度差值两类曲线，其中，通过图7(a)是处理后的列灰度均值，图7(b)对应局部范围的放大图，图7(c)为相邻两列灰度差值曲线。可以看出，本发明方法处理后，其对应的列方向灰度均值比较平滑，其中，相邻两列之间的灰度均值之差小于0.5个像素，而本实施例设定的λ为0.5，符合处理精度的要求。

本发明未详细说明部分属本领域技术人员公知常识。

Claims (6)

1.一种基于最小二乘法的线阵推扫式影像自适应辐射校正方法，其特征在于实现步骤如下：

(1)计算线阵推扫式遥感影像的列灰度均值和标准差序列

推扫式线阵相机生成的待处理影像大小为M×N个像素，对于推扫式线阵相机而言，值N也相当于推扫式线阵相机的探元个数；

影像中第j列的灰度平均值m_j和灰度标准差v_j表示为：

$m_j={\mathrm{\Σ}}_{i=1}^M{\mathrm{DN}}_{\mathrm{ij}}/M---\left(1\right)$

$v_j=\sqrt{{\mathrm{\Σ}}_{i=1}^M{({\mathrm{DN}}_{\mathrm{ij}}-m_j)}^2/M}---\left(2\right)$

其中，i表示影像的行号，i＝1，2，...，M，M为影像的总行数，j表示影像的列号，j＝1，2，...，N，N为影像的总列数，且M与N都为正整数，DN_ij表示第i行第j列像素的灰度值；

计算待处理影像所有列的灰度平均值和灰度标准差，得到的列均值序列表示为{m₁，m₂，...，m_N}，列标准差序列表示为{v₁，v₂，...，v_N}；

(2)坏扫描列的定位与有效灰度值的计算

(2.1)遍历步骤(1)计算得到的列均值序列和列标准差序列，当正在遍历的列满足均值和标准差为零时，确定该列为坏扫描列，该列中的所有像素为坏像素；每个坏扫描列中共有M个坏像素，这些像素的灰度值为无效灰度值；

(2.2)坏扫描列定位之后，根据该列左右正常列的灰度值采用加权平均的方法重新计算坏像素的灰度值，得到坏扫描列的有效灰度值；

(2.3)重复步骤(2.1)和步骤(2.2)，直到列均值序列和列标准差序列遍历结束后，则得到所有坏扫描列去除后的影像；

(2.4)对坏扫描列去除后的影像，按照步骤(1)中所述影像列灰度均值的计算方法，重新计算坏扫描列去除后影像所有列的灰度平均值，得到新的列均值序列{m₁′，m₂′，...，m_N′}；

(3)影像中条带列的自适应判别

(3.1)从左到右依次遍历新的列均值序列，将正在遍历的列j称为当前列，j＝1，2，...，N，根据遥感影像中条带的分布特征，依据当前列及其临近列的均值大小关系来判别当前列是否为条带列；

(3.2)用符号Sign_j标识当前列是否为条带列，若当前列为条带列则将其标识Sign_j赋值为“T”，否则赋值为“F”，重复步骤(3.1)直到整个新的列均值序列都遍历完毕，得到长度为N的条带标识序列S＝{Sign₁，Sign₂，...，Sign_j，...，Sign_N}；

(4)对条带列中的所有像素进行相对辐射校正处理；

(4.1)遍历步骤(3.2)得到的标识序列S，若某列的标识符号为“T”且其左右相邻列的标识符号都为“F”，则该列称为单条带列，执行步骤(4.2.1)；若存在多个标识符号“T”连续分布在一起，则存在连续的多条带列，执行步骤(4.2.2)；

(4.2)依据相对辐射校正输入输出的线性关系，推扫式线阵相机中探元的输出表示为原始输入的函数：

DN_new＝DN_old×Gain+Offset    (3)

式中，Gain为探元的增益，Offset为探元的偏移量，即称之为相对辐射校正系数，DN_old为输入值，即探元的初始响应值，DN_new为输出值，即经辐射校正后的新值；

(4.2.1)若条带列为单一条带列，则利用最小二乘法拟合计算该条带列的相对辐射校正系数，继续步骤(4.3)；

(4.2.2)若存在连续的多条带列，列数为W，W为正整数，2≤W＜N，则利用最小二乘法计算多条带列中所有列的相对辐射校正系数，继续步骤(4.3)；

(4.3)利用相对辐射校正系数对条带列中的每个像素进行相对辐射校正，具体计算公式如下：

DN_{ij_new}＝DN_{ij_old}×Gain_j+Offset_j    (4)

其中，Gain_j和Offset_j是上述步骤计算得到的相对辐射校正系数，分别表示第j列的增益和偏移量，DN_{ij_old}是第i行第j列像素的初始输入值，DN_{ij_new}是辐射校正后的输出值，即经过辐射校正后的灰度值；

(4.4)影像列号继续增加，重复执行步骤(4.1)至(4.3)，直至所有列遍历一遍；至此，完成影像中所有条带列的判别与其对应像素的相对辐射校正，得到相对辐射校正后的影像；

(5)保存辐射校正后的影像，并可根据用户要求将计算的辐射校正系数入库，为后续快速应用提供查询数据库。

2.根据权利要求1所述基于最小二乘法的线阵推扫式影像自适应辐射校正方法，其特征在于：所述步骤(3.1)中判别条带列的步骤如下：

(3.1.1)当前列为第1列时，则利用第1、2列的灰度均值m₁′、m₂′与给定阈值λ的大小关系来判别当前列是否为条带列，若m₂′-m₁′≥λ，则第1列为条带列；若m₁′-m₂′≥λ，则第1列不是条带列；若|m₁′-m₂′|＜λ，则第1、2列都是条带列或都不是条带列，继续步骤(3.1.2)；

(3.1.2)当前列j不等于第1列或最后一列N时，从左到右遍历新的列均值序列，利用当前列及其左右两列的灰度均值，m_j-1′、m_j′、m_j+1′，之间的差值与λ的关系确定条带列：①若m_j-1′-m_j′≥λ且m_j+1′-m_j′≥λ，则当前列为条带列，左右两列不是条带列；②若m_j′-m_j-1′≥λ且m_j′-m_j+1′≥λ，则当前列不是条带列，而左右两列为条带列；③若m_j′-m_j-1′≥λ且m_j+1′-m_j-1′≥λ，则当前列不是条带列，而左边列为条带列；④若m_j′-m_j+1′≥λ且m_j+1′-m_j-1′≥λ，则当前列不是条带列，而右边列为条带列；⑤若m_j-1′-m_j′≥λ且m_j-1′-m_j+1′≥λ，则当前列j及其右侧列j+1为条带列，此时为多条带列分布在一起；⑥若m_j+1′-m_j′≥λ且m_j+1′-m_j-1′≥λ，则当前列j与左侧列j-1都为条带列，此时为多列条带分布在一起；⑦若|m_j+1′-m_j′|＜λ且|m_j-1′-m_j′|＜λ，则该列和左右两列全是条带或全不是条带；

(3.1.3)当前列为第N列时，利用第N-1、N列的灰度均值m_N-1′、m_N′判别第N列是否为条带列，若m_N-1′-m_N′≥λ，则第N列为条带列；若m_N′-m_N-1′≥λ，则第N列不是条带列；若|m_N-1′-m_N′|＜λ，则第N列与第N-1列一致，即当第N-1列为条带列时，则第N列也是条带列，否则当第N-1列不是条带列时，则第N列也不是条带列。

3.根据权利要求2所述基于最小二乘法的线阵推扫式影像自适应辐射校正方法，其特征在于：所述阈值λ也称作相对辐射校正精度，根据用户对精度的需求进行配置，依据经验λ的取值为0.5≤λ≤2之间的实数，当要求相对辐射校正精度优于0.5个灰度值，则可设定λ为0.5，...，当要求相对辐射校正精度优于2个灰度值，则可设定λ为2。

4.根据权利要求1所述基于最小二乘法的线阵推扫式影像自适应辐射校正方法，其特征在于：所述步骤(4.2.1)的单一条带列的相对辐射校正系数计算步骤如下：

(4.2.1.1)当第1列是单一条带列时，对该列的每个像素建立一个线性方程，方程的输入为条带列像素的初始值，方程的输出为第2列相应像素的加权平均值；

该列第1个像素的灰度值DN₁₁作为方程一的输入值，第2列中第1、2像素的加权均值作为方程一的输出值，加权系数为[2，1]；该列第2个像素的灰度值DN₂₁作为方程二的输入值，第2列中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[1，2，1]；依次类推，...；该列第M个像素的灰度值DN_M1作为方程M的输入值，第2列中第M-1、M像素的加权均值作为方程M的输出值，加权系数为[1，2]；该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{11}\×{\mathrm{Gain}}_1+{\mathrm{Offset}}_1=[2{\mathrm{DN}}_{12}+{\mathrm{DN}}_{22}]/3\\ {\mathrm{DN}}_{21}\×{\mathrm{Gain}}_1+{\mathrm{Offset}}_1=[{\mathrm{DN}}_{12}+2{\mathrm{DN}}_{22}+{\mathrm{DN}}_{32}]/4\\ ...\\ {\mathrm{DN}}_{M1}\×{\mathrm{Gain}}_1+{\mathrm{Offset}}_1=[{\mathrm{DN}}_{(M-1)2}+2{\mathrm{DN}}_{M2}]/3\end{array}\right.---\left(5\right)$

式中，Gain₁和Offset₁分别为第1列的增益和偏移量，即相对辐射校正系数；

(4.2.1.2)当第k列，k＝2，3，...，N-1，为单一条带列时，则条带列的相邻两列第k-1和第k+1列为正常列，针对条带列中的每个像素建立一个线性方程，利用条带列像素的初始值作为输入值，以条带列像素的周围正常像素的加权平均作为输出值；

该列第1个像素的灰度值DN_1k作为方程一的输入值，正常列第k-1和第k+1列中第1、2像素的加权均值作为方程一的输出值，第k-1和第k+1列中对应像素的加权系数均为[2，1]；该列第2个像素的灰度值DN_2k作为方程二的输入值，第k-1和第k+1列中第1、2、3像素的加权均值作为方程二的输出值，第k-1和第k+1列中对应像素的加权系数均为[1，2，1]；依次类推，...；该列第M个像素的灰度值DN_Mk作为方程M的输入值，第k-1和第k+1列中第M-1、M像素的加权均值作为方程M的输出值，第k-1和第k+1列中对应像素的加权系数均为[1，2]；该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1k}\×{\mathrm{Gain}}_k+{\mathrm{Offset}}_k=[2{\mathrm{DN}}_{1(k-1)}+{2\mathrm{DN}}_{1(k+1)}+{\mathrm{DN}}_{2(k-1)}+{\mathrm{DN}}_{2(k+1)}]/6\\ {\mathrm{DN}}_{2k}\×{\mathrm{Gain}}_k+{\mathrm{Offset}}_k=[{\mathrm{DN}}_{1(k-1)}+{\mathrm{DN}}_{1(k+1)}+2{\mathrm{DN}}_{2(k-1)}+\\ 2{\mathrm{DN}}_{2(k+1)}+{\mathrm{DN}}_{3(k-1)}+{\mathrm{DN}}_{3(k+1)}]/8\\ ...\\ {\mathrm{DN}}_{\mathrm{Mk}}\×{\mathrm{Gain}}_k+{\mathrm{Offset}}_k=[{\mathrm{DN}}_{(M-1)(k-1)}{+\mathrm{DN}}_{(M-1)(k+1)}+{2\mathrm{DN}}_{M(k-1)}+{2\mathrm{DN}}_{M(k+1)}]/6\end{array}\right.---\left(6\right)$

式中，Gain_k和Offset_k分别为第k列的增益和偏移量，即相对辐射校正系数；

(4.2.1.3)当最后一列N是条带列时，针对该列中的每个像素建立一个线性方程，第N-1列为正常列，利用条带列像素的初始值作为方程输入值，利用第N-1列的正常像素的加权平均作为方程输出值；

该列第1个像素的灰度值DN_1N作为方程一的输入值，第N-1中第1、2像素的加权均值作为方程一的输出值，加权系数为[2，1]；该列第2个像素的灰度值DN_2N作为方程二的输入值，第N-1中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[1，2，1]；依次类推，...；该列第M个像素的灰度值DN_MN作为方程M的输入值，第N-1列中第M-1、M像素的加权均值作为方程M的输出值，加权系数均为[1，2]；该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1N}\×{\mathrm{Gain}}_N+{\mathrm{Offset}}_N=[2{\mathrm{DN}}_{1(N-1)}+{\mathrm{DN}}_{2(N-1)}]/3\\ {\mathrm{DN}}_{2N}\×{\mathrm{Gain}}_N+{\mathrm{Offset}}_N=[{\mathrm{DN}}_{1(N-1)}+2{\mathrm{DN}}_{2(N-1)}+{\mathrm{DN}}_{3(N-1)}]/4\\ ...\\ {\mathrm{DN}}_{\mathrm{MN}}\×{\mathrm{Gain}}_N+{\mathrm{Offset}}_N=[{\mathrm{DN}}_{(M-1)(N-1)}+2{\mathrm{DN}}_{M(N-1)}]/3\end{array}\right.---\left(7\right)$

式中，Gain_N和Offset_N分别为第N列的增益和偏移量，即相对辐射校正系数；

(4.2.1.4)最小二乘法拟合计算辐射校正系数；

用最小二乘法拟合上述M个方程的系数，即可得到相对辐射校正系数。

5.根据权利要求1所述基于最小二乘法的线阵推扫式影像自适应辐射校正方法，其特征在于，所述步骤(4.2.2)的连续多条带列的相对辐射校正系数具体计算步骤如下：

(4.2.2.1)当从第1列开始连续W列为条带列，即第1列～第W列组成连续的多条带，右侧第W+1列为正常列，从左到右逐列计算相对辐射校正系数；

对于第α列，1≤α≤W，针对每个像素建立一个线性方程，利用条带列像素的初始值作为方程的输入值，以正常列即第W+1列中像素的加权均值作为方程的输出值，加权系数与距离成反比；

该列第1个像素的灰度值DN_1α作为方程一的输入值，第W+1列中第1、2像素的加权均值作为方程一的输出值，加权系数为[1+α，α]；该列第2个像素的灰度值DN_2α作为方程二的输入值，第W+1列中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[α，1+α，α]；依次类推，...；该列第M个像素的灰度值DN_Mα作为方程M的输入值，第W+1列中第M-1、M像素的加权值作为方程M的输出值，加权系数为[α，1+α]；该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1\mathrm{\α}}\×{\mathrm{Gain}}_{\mathrm{\α}}+{\mathrm{Offset}}_{\mathrm{\α}}=[(1+\mathrm{\α}){\mathrm{DN}}_{1(W+1)}+{\mathrm{\αDN}}_{2(W+1)}]/(2\mathrm{\α}+1)\\ {\mathrm{DN}}_{2\mathrm{\α}}\×{\mathrm{Gain}}_{\mathrm{\α}}+{\mathrm{Offset}}_{\mathrm{\α}}=[{\mathrm{\αDN}}_{1(W+1)}+(1+\mathrm{\α}){\mathrm{DN}}_{2(W+1)}+{\mathrm{\αDN}}_{3(W+1)}]/(3\mathrm{\α}+1)\\ ...\\ {\mathrm{DN}}_{\mathrm{M\α}}\×{\mathrm{Gain}}_{\mathrm{\α}}+{\mathrm{Offset}}_{\mathrm{\α}}=[{\mathrm{\αDN}}_{(M-1)(W+1)}+(1+\mathrm{\α}){\mathrm{DN}}_{M(W+1)}]/(2\mathrm{\α}+1)\end{array}\right.---\left(8\right)$

式中，Gain_α和Offset_α分别为第α列的增益和偏移量，即相对辐射校正系数；

(4.2.2.2)当第Q+1～Q+W共W列为连续多条带列，其中，2≤Q+1＜Q+W≤N-1，则条带列左右正常两列为第Q与Q+W+1列，每列有M个像素；从左到右逐列计算相对辐射校正系数：

对于第Q+β列条带，1≤β≤W，针对该列中的每个像素建立一个线性方程，利用条带列像素的初始值作为方程的输入值，以第Q与Q+W+1正常列像素的距离加权平均值作为方程的输出值；

该列第1个像素的灰度值DN_1(Q+β)作为方程一的输入值，第Q与Q+W+1列中第1、2像素的加权均值作为方程一的输出值，列Q中两个像素的加权系数为[W+2-β，W+1-β]，列Q+W+1中两个像素的加权系数为[1+β，β]；该列第2个像素的灰度值DN_2(Q+β)作为方程二的输入值，第Q与Q+W+1列中第1、2、3像素的加权均值作为方程二的输出值，列Q中三个像素的加权系数为[W+1-β，W+2-β，W+1-β]，列Q+W+1中三个像素的加权系数为[β，1+β，β]；依次类推，...；该列第M个像素的灰度值DN_M(Q+β)作为方程M的输入值，第Q与Q+W+1列中第M-1、M像素的加权均值作为方程M的输出值，列Q中两个像素的加权系数为[W+1-β，W+2-β]，列Q+W+1中两个像素的加权系数为[1+β，β]；该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1(Q+\mathrm{\β})}\×{\mathrm{Gain}}_{(Q+\mathrm{\β})}+{\mathrm{Offset}}_{(Q+\mathrm{\β})}=[(W+2-\mathrm{\β}){\mathrm{DN}}_{1Q}+(W+1-\mathrm{\β}){\mathrm{DN}}_{2Q}+\\ (1+\mathrm{\β}){\mathrm{DN}}_{1(Q+W+1)}+\mathrm{\β}{\mathrm{DN}}_{2(Q+W+1)}]/(2W+4)\\ {\mathrm{DN}}_{2(Q+\mathrm{\β})}\×{\mathrm{Gain}}_{(Q+\mathrm{\β})}+{\mathrm{Offset}}_{(Q+\mathrm{\β})}=[(W+1-\mathrm{\β}){\mathrm{DN}}_{1Q}+(W+2-\mathrm{\β}){\mathrm{DN}}_{2Q}+\\ (W+1-\mathrm{\β}){\mathrm{DN}}_{3Q}+\mathrm{\β}{\mathrm{DN}}_{1(Q+W+1)}+(1+\mathrm{\β}){\mathrm{DN}}_{2(Q+W+1)}+\mathrm{\β}{\mathrm{DN}}_{3(Q+W+1)}]/(3W+5)\\ ...\\ {\mathrm{DN}}_{M(Q+\mathrm{\β})}\×{\mathrm{Gain}}_{(Q+\mathrm{\β})}+{\mathrm{Offset}}_{(Q+\mathrm{\β})}=[(W+1-\mathrm{\β}){\mathrm{DN}}_{(M-1)Q}+(W+2-\mathrm{\β}){\mathrm{DN}}_{\mathrm{MQ}}+\\ (1+\mathrm{\β}){\mathrm{DN}}_{(M-1)(Q+W+1)}+\mathrm{\β}{\mathrm{DN}}_{M(Q+W+1)}]/(2W+4)\end{array}\right.---\left(9\right)$

式中，Gain_(Q+β)和Offset_(Q+β)分别为第(Q+β)列的增益和偏移量，即相对辐射校正系数；

(4.2.2.3)当从第N-W+1列开始到最后一列N，共连续W列为条带列，即第N-W+1、N-W+2、...、N共W列组成连续的多条带，左侧第N-W列为正常列，从左到右逐列计算相对辐射校正系数；

对于第N-γ列，N-W+1≤N-γ≤N，针对该列中的每个像素建立一个线性方程，利用条带列像素的初始值作为方程的输入值，以正常列像素的加权均值作为方程的输出值，加权系数与距离成反比；

该列第1个像素的灰度值DN_1(N-γ)作为方程一的输入值，正常列第N-W列中第1、2像素的加权均值作为方程一的输出值，加权系数为[W-γ，W-γ-1]；该列第2个像素的灰度值DN_2(N-γ)作为方程二的输入值，第N-W列中第1、2、3像素的加权均值作为方程二的输出值，加权系数为[W-γ-1，1 W-γ，W-γ-1]；依次类推，...；该列第M个像素的灰度值DN_M(N-γ)作为方程M的输入值，第N-W列中第M-1、M像素的加权值作为方程M的输出值，加权系数为[W-γ-1，W-γ]；该列的M个像素可建立M个方程：

$\left\{\begin{array}{c}{\mathrm{DN}}_{1(N-\mathrm{\γ})}\×{\mathrm{Gain}}_{(N-\mathrm{\γ})}+{\mathrm{Offset}}_{(N-\mathrm{\γ})}=[(W-\mathrm{\γ}){\mathrm{DN}}_{1(N-W)}+\\ (W-\mathrm{\γ}-1){\mathrm{DN}}_{2(N-W)}]/(2W-2\mathrm{\γ}-1)\\ {\mathrm{DN}}_{2(N-\mathrm{\γ})}\×{\mathrm{Gain}}_{(N-\mathrm{\γ})}+{\mathrm{Offset}}_{(N-\mathrm{\γ})}=[(W-\mathrm{\γ}-1){\mathrm{DN}}_{1(N-W)}+(W+\mathrm{\γ}){\mathrm{DN}}_{2(N-W)}+\\ (W-\mathrm{\γ}-1){\mathrm{DN}}_{3(N-W)}]/(3W-2\mathrm{\γ}-2)\\ ...\\ {\mathrm{DN}}_{M(N-\mathrm{\γ})}\×{\mathrm{Gain}}_{(N-\mathrm{\γ})}+{\mathrm{Offset}}_{(N-\mathrm{\γ})}=[(W-\mathrm{\γ}-1){\mathrm{DN}}_{(M-1)(N-W)}+\\ (W-\mathrm{\γ}){\mathrm{DN}}_{M(N-W)}]/(2W-2\mathrm{\γ}-1)\end{array}\right.---\left(10\right)$

式中，Gain_(N-γ)和Offset_(N-γ)分别为第(N-γ)列的增益和偏移量，即相对辐射校正系数。

6.根据权利要求1所述基于最小二乘法的线阵推扫式影像自适应辐射校正方法，其特征在于：所述步骤(5)的影像数据与辐射校正系数的存储方式如下：将影像中所有条带列的相对辐射校正结果以实数形态单独存储成一个附加文件，在后续应用过程中能够直接读取附加文件中的数据代替对应列的影像数据，为后续处理提供更精确的数据；同时，将计算得到的相对辐射校正系数存储到数据库中，在后续应用过程中利用事先存储的辐射校正系数对输入的影像数据进行处理。

Images