CN106940219B - 一种在轨运行的宽波段卫星遥感器的光谱响应获取方法 - Google Patents

Abstract

遥感器光谱响应是评价遥感器光谱性能的必要参数，是遥感器把接收到的地物辐射量转换到遥感信号的关键参数，在遥感数据定量处理和应用中必不可少。本发明从理论上提出、并形成一种不依赖于实验室测试数据和星上光谱定标数据的在轨运行宽波段卫星遥感器的光谱响应获取方法。本发明通过分析遥感器成像机理和辐射传输过程，结合优化计算方法建立宽波段卫星遥感器光谱响应解算模型；利用地物光谱库和模拟图像来分析模型解算对数据要求与解算规则；利用同步高光谱遥感图像获取满足模型解算要求的解算样本，并最终实现宽波段卫星遥感器光谱响应的数学解算。

Description

一种在轨运行的宽波段卫星遥感器的光谱响应获取方法

技术领域

本发明是一种对在轨运行的宽波段卫星遥感器的光谱响应获取方法，属于″摄影测量与遥感″学科中″定量遥感″技术领域。

背景技术

随着遥感及其相关技术发展日趋成熟，遥感应用也从定性走向定量，定量遥感是遥感发展的需求和必然结果。定量遥感的前提和基础是遥感数据的定量化，即建立遥感器输出的遥感信号和地表物质反射、发射的辐射量之间的定量关系，在这个过程中，只有通过遥感器光谱响应作用后的能量才最终被量化成遥感信号。因此，遥感器光谱响应是遥感器把接收到的地物辐射量转换到遥感信号的关键参数，是定量遥感的″定量″的前提和基础。另外，由于光谱响应是遥感器对入射能量进行有效转换的一个转换关系参数，光谱响应在遥感数据后续的定量反演和对反演结果的定量化分析中也是不可或缺的必要参数。

遥感器光谱响应一般只能在实验室条件下进行精确测量，在卫星上天后不能直接进行精确测量，目前主要通过星上定标法和光谱吸收带定标法这两种方法对光谱响应进行监测和修正。星上定标法是在卫星上安置稳定的、对光谱吸收或反射敏感的特殊材料和设备，对遥感器进行定期监测和修正。这种方法对星上定标设备的性能要求很高，在实际应用中很难保证精度，且星上定标设备昂贵，因而没有得到广泛应用。光谱吸收特征定标法的前提是遥感器光谱分辨率足够高，小于某些特殊物质在特定的波长位置上的光谱吸收特征，如氧气在760nm的吸收特征，通过比较图像光谱特征和固有光谱特征来评价和修正遥感器光谱响应。

对于波段响应范围在50nm-100nm左右的宽波段卫星遥感器来说，比如Landsat的TM、SPOT卫星的HRV等，由于光谱响应范围大于一般物质的光谱吸收特征，不能用光谱吸收定标法对光谱响应进行标定，现在通用的做法是在卫星发射前在实验室对遥感器进行光谱定标，精确测量每一个波段的光谱响应，假设卫星在发射前后遥感器光谱响应不变化或者其变化可以忽略，直接使用实验室测定的光谱响应。

对于宽波段卫星遥感器而言，假设其光谱响应在其整个生命周期不变化的做法在实际应用中存在一些不足和问题：首先，在严格意义上这种假设不成立。事实上，卫星发射后由于仪器元件的老化和外界因素的干扰，遥感器的性能和灵敏度不断衰减和下降。其次，在特殊情况下，一些遥感器没有在卫星发射前全面测定其光谱响应，极大限制了遥感数据定量化处理和应用。因此，无论是对宽波段卫星遥感器光谱性能的在轨监测和修正，还是为了没有光谱响应的在轨卫星提供光谱响应，都有必要寻求一种新的、不依赖于实验室测试的光谱响应获取方法。

发明内容

本发明针对以上现有技术中存在的不足，提出了一种不依赖实验室测试数据和星上光谱定标数据的宽波段卫星遥感器光谱响应获取方法，实现对在轨运行宽波段遥感器光谱响应的实际测量和更新，解决目前仅依靠卫星发射前在实验室测量的光谱响应作为卫星遥感器整个生命周期内唯一可用值的缺陷和不足。本发明具体步骤包括：

A.根据遥感器成像原理，确定光谱响应在遥感器入瞳辐亮度和遥感图像DN值之间转换过程中的能量传递关系；

B.依据光谱响应在遥感信号转换过程中的作用和光谱响应具有的曲线特性，将其分解为高斯函数形式的快变函数和待求的慢变函数，建立光谱响应的物理模型和解算方程；

C.在Philips光滑约束条件和在上述步骤B中光谱响应解算方程解的2阶导数平方和最小的条件下，综合得到光谱响应解的表达式；根据终得到光谱响应解的表达式，采用迭代方法计算其最优解，该最优解即为遥感器的光谱响应解算值；

D.依据上述步骤A-C的原理和方法，用已知光谱响应遥感器为实验对象，从光谱库中选取在待解算遥感器波段范围内差异较大的不同地物光谱；依据辐射传输方程，模拟得到上述C中光谱响应解算方程的观测样本；通过改变方程观测样本条件，分析建立观测样本个数、观测样本之间的相关性、以及待求解光谱响应被离散的个数与最终解算获取的光谱响应函数之间的内在联系，并形成利用同步卫星影像进行光谱响应解算的先验知识；

E.获取待解算宽波段卫星遥感器卫星影像、同步高光谱遥感影像、同步测量大气参数、待解算卫星和参考高光谱卫星的观测条件参数，依据上述步骤D中获取的先验知识，采用交叉辐射定标的方法，利用同步的高光谱影像计算得到待解算遥感器卫星影像的入瞳辐亮度；以计算得到的入瞳辐亮度和待解算遥感器图像DN值为解算样本，按照上述步骤A的方法计算得到待解算遥感器的光谱响应。

作为一种解算方法，上述步骤A-C中，基于遥感器成像原理，采用数学优化的方法推导建立基于同步高光谱图像的光谱响应获取方法，具体步骤包括：

1)根据遥感器成像原理，建立遥感器入瞳辐亮度g(λ)和遥感图像DN值之间的定量函数关系，该函数关系中包含了遥感器入瞳辐亮度g(λ)与遥感器归一化光谱响应s′(λ)之间的积分运算；

2)依据光谱响应的特性，将s′(λ)分解为已知的高斯函数形式的快变函数h(λ)和待求的慢变函数f(λ)，即s′(λ)＝h(λ)f(λ)；

3)综合上述步骤1)、2)中的函数关系，得到光谱响应解算模型，D_i＝A_ijf+ε，其中ε为误差向量，D_i＝DN_i是DN值向量，f是慢变函数f(λ)的离散化向量；

4)按照Philips光滑约束条件，使上述步骤3)得到的光谱响应解算模型D_i＝A_ijf+ε的解的2阶导数平方和最小，该模型解Q通过最小化性能函数得到：其中γ是Lagrange光滑因子，为非负值；

5)求Q对f的导数，可得，-A^TC^-1ε+γHf＝0，其中C_ij是观测值的协方差矩阵C的元素，如果观测值之间相互独立，C是单位矩阵；H是光滑矩阵；

6)综合上述步骤4)和5)所得公式，最终得到光谱响应解的表达式为f(λ)＝(A^TC^-1A+γH)^-1A^TC^-1D；

7)采用迭代方法计算上述步骤6)所述的光谱响应解的最优解。

进一步地，上述步骤7)中，采用迭代法计算最优解的具体步骤包括：

a)根据遥感器光谱响应一般形状满足高斯函数分布的假定，确定快变函数其中μ₀＝(λ₂+λ₁)/2，σ₀＝(λ₂-λ₁)/4，上标(0)表示初始值，λ₁和λ₂分别是待解算遥感器光谱响应的最小、最大波长；

b)把h⁽⁰⁾(λ)带入公式中计算出A的元素；

c)利用公式f(λ)＝(A^TC^-1A+γH)^-1A^TC^-1D计算出f⁽¹⁾(λ)，上标(1)表示第一次计算结果；

d)h⁽¹⁾(λ)＝f⁽¹⁾(λ)h⁽⁰⁾(λ)，重复上述步骤b)和c)，直到得到稳定解f⁽ⁿ⁾(λ)，上标(n)表示第n次迭代计算结果，下同；

e)最终求得归一化光谱响应s′(λ)＝f⁽ⁿ⁾(λ)h⁽ⁿ⁾(λ)。

作为一种解算方法，其特征在于所述步骤D中，采用和待解算遥感器卫星同步或准同步的高光谱卫星图像获得真实的方程解算样本，解决了在卫星同步的地面光谱采集不足问题，具体步骤包括：

(1)依据先验知识，确定可用的高光谱遥感器，获取和待解算遥感器同步或准同步的影像以及待解算遥感器图像；并获取同步测量的大气参数；

(2)依据所述的先验知识，确定选取的观测样本个数，从参考的高光谱影像获取同名地物的辐亮度值，从待解算遥感影像获取同名地物的DN值；

(3)使用同步获取的大气参数和卫星观测条件，对参考的光谱影像进行大气校正；

(4)使用同步获取的大气参数和卫星观测条件，对经过大气校正的辐亮度进行辐射传输计算，计算为待解算遥感器的入瞳辐亮度；

(5)依据上述步骤1)-7)或a)-e)所述方法对光谱响应进行解算，得到待解算遥感器的光谱响应。

具体实施方式

在轨运行宽波段卫星遥感器的光谱响应获取过程分为三步，(1)基于遥感器成像原理建立解算模型，并根据数学方法推导出模型的解算方法和步骤；(2)利用光谱库数据，基于辐射传输过程模拟得到满足解算需求的数值，通过对模拟数据的分析建立实际解算的先验知识；(3)利用同步的高光谱卫星图像，基于交叉定标的方法，获得真实的解算样本，并最终求得待求遥感器的光谱响应。三部分工作具体说明如下：

一、模型建立与推导

基于线性遥感***的假设，遥感器接收到的辐射能量表达为入射辐亮度和遥感器光谱响应的函数，即：

其中，λ₁和λ₂光谱响应范围，λ₁是最小值，λ₂是最大值；g(λ)是遥感器入瞳的辐亮度；s(λ)遥感器光谱响应。

由于原始图像提供的是DN₀值，辐亮度和DN₀值，以及扣除暗电流的DN值的关系为：

L＝A·(DN₀-B)＝Gain·DN (2)

公式1代入2中得：

由于定标系数Gain是一个固定常量，由于遥感器的光谱响应通常要做归一化处理，常数Gain对归一化的光谱响应没有影响，因此，

令为归一化光谱响应，公式3可以改写为：

如果s′(λ)符合某种函数分布则有利于计算s′(λ)。众所周知，对于高光谱分辨率波段来说，其光谱响应可以近似用高斯函数代替，对于宽波段卫星遥感器来说，其光谱响应s′(λ)有高斯函数分布特征，但函数形式随波长变化较大，如果按照高斯函数进行反演容易导致解的不稳定。为此，把宽波段卫星遥感器光谱响应s′(λ)分解为已知的高斯函数形式的快变函数h(λ)和待求的慢变函数f(λ)，即s′(λ)＝h(λ)f(λ)，则方程4变为：

则可以把方程5写成线性方程式，即光谱响应解算模型：

D＝Af+ε (6)

其中ε为误差向量；

D_i＝DN_i，i＝1，2，3，...，p，p是第i个DN值，p是参与运算的观测样本(图像DN值个数)。

j＝1，2，3，...，q，q是宽波段卫星遥感器光谱响应被离散的个数；f是慢变函数f(λ)离散化的向量；

对于方程6，如果直接通过最小化误差矢量求解f，由于f只能被离散为有限几个未知数，另外观测样本g之间也存在一定相关性，所以导致f的解不稳定，不能得到满意解。Philips引入了光滑约束条件，在方程6解的2阶导数平方和最小的条件下得到解Q表达式：

其中γ是Lagrange光滑因子，为非负值。

求Q对f的导数，可得：

-A^TC^-1ε+γHf＝0 (9)

C_ij是观测值的协方差矩阵C的元素，如果观测值之间相互独立，C是单位矩阵；H是光滑矩阵，其形式为：

联合方程8和9可得f的解为

f(λ)＝(A^TC^-1A+γH)^-1A^TC^-1D (11)

对方程11采用迭代方法解算，具体步骤如下：

a)根据遥感器光谱响应一般形状满足高斯函数分布的假定，确定快变函数h⁽⁰⁾(λ)

式中：

μ₀＝(λ₂+λ₁)/2，

σ₀＝(λ₂-λ₁)/4

b)把h⁽⁰⁾(λ)带入公式7中计算出A的元素，上标(0)表示初始值；

c)利用公式11计算出f⁽¹⁾(λ)，上标(1)表示第一次计算结果；

d)h⁽¹⁾(λ)＝f⁽¹⁾(λ)h⁽⁰⁾(λ)，重复上述步骤b，c，直到得到稳定解f⁽ⁿ⁾(λ)，上标(n)表示第n次计算结果；

e)最终求得光谱响应s′(λ)＝f⁽ⁿ⁾(λ)h⁽ⁿ⁾(λ)。

二、模型解算先验知识获取

由于解算模型中观测样本之间存在一定相关性，而且受限于高光谱图像的光谱分辨率，连续的宽波段光谱响应只能离散为有限的几个值，因此，解算模型是一个病态的，虽然解算方法采用了优化算法，但是仍然不能完全消除病态效果，因此，模型的解算样本的先验知识就尤为重要。本发明采用光谱库数据，基于辐射传输方程，模拟得到待解算遥感器的入瞳辐亮度和图像DN值，通过调整观测样本个数、观测样本之间的光谱差异、光谱响应被离散的个数来获取模型解算的先验知识。具体步骤如下：

1)用已知光谱响应遥感器为实验对象；

2)从光谱库中选取在待解算遥感器波段范围内差异较大的不同地物光谱；

3)依据辐射传输方程，模拟得到在一定成像条件下待解算遥感器入瞳辐亮度数值和输出图像DN值，即方程解算的观测样本；

4)通过改变方程观测样本个数、观测样本之间的相关性、待求解光谱响应被离散的个数，通过比较分析数学解算的光谱响应和已知光谱响应之间的差别，建立光谱响应在观测样本个数、观测样本之间的相关性、待求解光谱响应被离散的个数三方面的先验知识；

三、真实光谱响应解算数据获取

由于光谱响应解算的观测样本个数多于被离散的光谱响应的未知数个数，观测样本的获取应该是和待解算遥感器图像同步的地面光谱，考虑到在实际应用中很难在卫星过顶时获取足够的同步观测光谱。为此，本发明提出利用和待解算卫星同步的高光谱遥感图像，采用交叉定标的方法，从高光谱图像上获取待解算遥感器图像入瞳辐亮度。具体步骤如下：

1)依据模型解算的先验知识，确定可用的高光谱遥感器，获取和待解算遥感器同步或准同步的影像以及待解算遥感器图像；并获取同步测量的大气参数；

2)依据模型解算的的先验知识，确定选取的观测样本个数，从参考的高光谱影像获取同名地物的辐亮度值，从待解算遥感影像获取同名地物的DN值；

3)使用同步获取的大气参数和卫星观测条件，对参考的光谱影像进行大气校正；

4)使用同步获取的大气参数和卫星观测条件，对经过大气校正的辐亮度进行辐射传输计算，计算为待解算遥感器的入瞳辐亮度；

5)依据光谱响应解算模型的步骤行解算，得到待解算遥感器的光谱响应。

Claims (4)

1.一种在轨运行宽波段卫星遥感器的光谱响应获取方法，其特征在于，不依赖于卫星发射前的实验室测试数据和卫星发射后的星上光谱定标数据，利用卫星自身和同步影像及观测条件参数来计算获取卫星遥感器光谱响应，实现对在轨运行宽波段遥感器光谱响应的实际测量和更新，具体包括以下步骤：

A.根据遥感器成像原理，确定光谱响应在遥感器入瞳辐亮度和遥感图像DN值之间转换过程中的能量传递关系；

B.依据光谱响应在遥感信号转换过程中的作用和光谱响应具有的曲线特性，将其分解为高斯函数形式的快变函数和待求的慢变函数，建立光谱响应的物理模型和解算方程；

C.在Philips光滑约束条件和在上述步骤B中光谱响应解算方程解的2阶导数平方和最小的条件下，综合得到光谱响应解的表达式；根据终得到光谱响应解的表达式，采用迭代方法计算其最优解，该最优解即为遥感器的光谱响应解算值；

D.依据上述A－C的原理和方法，用已知光谱响应遥感器为实验对象，从光谱库中选取在待解算遥感器波段范围内差异较大的不同地物光谱；依据辐射传输方程，模拟得到上述步骤C中光谱响应解算方程的观测样本；通过改变方程观测样本条件，分析建立观测样本个数、观测样本之间的相关性、以及待求解光谱响应被离散的个数与最终解算获取的光谱响应函数之间的内在联系，并形成利用同步卫星影像进行光谱响应解算的先验知识；

E.获取待解算宽波段卫星遥感器卫星影像、同步高光谱遥感影像、同步测量大气参数、待解算卫星和参考高光谱卫星的观测条件参数，依据上述D中获取的先验知识，采用交叉辐射定标的方法，利用同步的高光谱影像计算得到待解算遥感器卫星影像的入瞳辐亮度；以计算得到的入瞳辐亮度和待解算遥感器图像DN值为解算样本，按照上述步骤A的方法计算得到待解算遥感器的光谱响应。

2.根据权利要求1所述的在轨运行宽波段卫星遥感器的光谱响应获取方法，其特征在于所述步骤A－C中，基于遥感器成像原理，采用数学优化的方法推导建立了一种基于同步高光谱图像的光谱响应获取方法，具体步骤包括：

1)根据遥感器成像原理，建立遥感器入瞳辐亮度g(λ)和遥感图像DN值之间的定量函数关系，该函数关系中包含了遥感器入瞳辐亮度g(λ)与遥感器归一化光谱响应s'(λ)之间的积分运算；

2)依据光谱响应的特性，将s'(λ)分解为已知的高斯函数形式的快变函数h(λ)和待求的慢变函数f(λ)，即s'(λ)＝h(λ)f(λ)；

3)综合上述步骤1)、2)中的函数关系，得到光谱响应解算模型，D_i＝A_ijf+ε，其中ε为误差向量，D_i＝DN_i是DN值向量，f是慢变函数f(λ)的离散化向量；

4)按照Philips光滑约束条件，使上述3)得到的光谱响应解算模型D_i＝A_ijf+ε解的2阶导数平方和最小，该模型解Q通过最小化性能函数得到其中γ是Lagrange光滑因子，为非负值；

5)求Q对f的导数，可得，-A^TC^-1ε+γHf＝0，其中C_ij是观测值的协方差矩阵C的元素，如果观测值之间相互独立，C是单位矩阵；H是光滑矩阵；

6)综合上述步骤4)和5)所得公式，最终得到光谱响应解的表达式为f＝(A^TC^-1A+γH)^{- 1}A^TC^-1D；

7)采用迭代方法计算上述步骤6)所述的光谱响应解的最优解。

3.根据权利要求2所述的在轨运行宽波段卫星遥感器的光谱响应获取方法，其特征在于，所述步骤7)中，采用迭代法计算最优解，具体步骤包括：

a)根据遥感器光谱响应一般形状满足高斯函数分布的假定，确定快变函数其中μ₀＝(λ₂+λ₁)/2，σ₀＝(λ₂-λ₁)/4，上标(0)表示初始值，λ₁和λ₂分别是待解算遥感器光谱响应的最小、最大波长；

b)把h⁽⁰⁾(λ)带入公式中计算出A的元素；

c)利用公式f(λ)＝(A^TC^-1A+γH)^-1A^TC^-1D计算出f⁽¹⁾(λ)，上标(1)表示第一次计算结果；

d)h⁽¹⁾(λ)＝f⁽¹⁾(λ)h⁽⁰⁾(λ)，重复上述步骤b)和c)，直到得到稳定解f⁽ⁿ⁾(λ)，上标(n)表示第n次迭代计算结果，下同；

e)最终求得归一化光谱响应s'(λ)＝f⁽ⁿ⁾(λ)h⁽ⁿ⁾(λ)。

4.根据权利要求3所述的在轨运行宽波段卫星遥感器的光谱响应获取方法，其特征在于，所述步骤E中，采用和待解算遥感器卫星同步的高光谱卫星图像获得真实的方程解算样本，解决了在卫星同步的地面光谱采集不足问题，具体步骤包括：

1)依据先验知识，确定可用的高光谱遥感器，获取和待解算遥感器同步的影像以及待解算遥感器图像；并获取同步测量的大气参数；

2)依据所述的先验知识，确定选取的观测样本个数，从参考的高光谱影像获取同名地物的辐亮度值，从待解算遥感影像获取同名地物的DN值；

3)使用同步获取大气参数和卫星观测条件，对参考光谱影像进行大气校正；

4)使用同步获取的大气参数和卫星观测条件，对经过大气校正的辐亮度进行辐射传输计算，计算为待解算遥感器的入瞳辐亮度；

5)依据权利要求2或3所述方法步骤对光谱响应进行解算，得到待解算遥感器的光谱响应。