CN103268618B - 一种多光谱遥感数据真彩色定标方法 - Google Patents

Abstract

本发明公开了一种多光谱遥感数据真彩色定标方法，属于遥感图像处理技术领域，可用于航空、航天传感器平台获取的多光谱图像真彩色定标。根据实际地物光谱辐射信息，利用颜色匹配函数计算出地物真实颜色的RGB值，将真实RGB颜色与多光谱遥感图像中相同地物的红、绿、蓝通道数值相比较，建立图像真彩色转换矩阵，利用转换矩阵对整幅图像进行真彩色定标，得到的定标结果更加接近地物真实颜色，符合人眼视觉特性。

Description

一种多光谱遥感数据真彩色定标方法

技术领域

本发明涉及一种多光谱遥感数据真彩色定标方法，属于遥感图像处理技术领域，可用于航空、航天探测器平台获取的多光谱图像真彩色定标。

背景技术

多光谱遥感图像的真彩色合成是指彩色合成后的图像上地物色彩与实际地物色彩接近一致，且合成图像适合人眼视觉习惯。在虚拟现实、环境仿真、资源调查、影响解译尤其是非专业人士在使用遥感数据方面，真彩色遥感影像有广泛的应用。目前一般的多光谱彩色合成方法多是直接将多光谱图像的红、绿、蓝波段对应R、G、B通道进行合成。由于多光谱成像时，探测器的红、绿、蓝波段通道的光谱响应及波段范围的影响，直接用红绿蓝三通道合成彩色图像与地物的真实色彩有一定差异，影响遥感图像的解译与处理。

目前遥感图像的色彩定标方法大多以解决多景遥感图像镶嵌时色彩不均问题为目的，通过在景间重叠部分的特征点中抽取样本并计算统计参数，利用重叠区域统计线性模型定标方法使各景图像数据在统计分布上达到基本一致。此类方法归根结底是图像的非均匀性定标，由于脱离了地物色彩形成机理，仅仅利用统计手段使不同图像上地物色彩尽可能接近，并不能反映地物的真实颜色。有研究根据典型地物特征，将多光谱遥感图像各个波段进行加权处理再合成，如针对植被覆盖度高的地区，利用植被指数在植被区域对绿波段进行加权运算，再和红色、蓝色波段进行真彩色合成得到比较理想的效果，但这种方法容易造成图像中其他地物颜色的失真。

发明内容

本发明的目的是针对上述背景技术描述的不足，提出了一种多光谱遥感数据真彩色定标方法。本发明方法具体包含如下四个步骤：

步骤一、模拟地物在探测器入瞳出的辐亮度信息，计算出对应的RGB真彩色值；

步骤二、提取步骤一中地物在多光谱遥感图像中的RGB颜色值，建立同一个地物多光谱图像RGB值与真彩色RGB值之间的图像真彩色转换矩阵；

步骤三、利用步骤二中得到的图像真彩色转换矩阵对真彩色图像进行定标，得到更为接近地物真实色彩的影像。

对比现有技术，本发明的有益效果在于：本发明方法利用地物真实光谱信息合成真彩色RGB值，并与多光谱遥感图像红、绿、蓝波段灰度值相比较，建立真彩色转换矩阵，对整幅图像进行真彩色定标，主要有以下两点优势：

1.本方法采用典型地物真实光谱信息作为真彩色定标依据，得到的定标结果更加接近地物真实颜色。

2.本方法利用颜色匹配函数进行真彩色合成，符合人眼视觉特性。

附图说明

图1为多光谱遥感数据真彩色定标方法流程图。

具体实施方式

下面结合图1对本发明进行解释。

本发明方法包含如下四个步骤：

步骤一、模拟地物在探测器入瞳出的辐亮度信息，计算出对应的RGB真彩色值，具体过程包括如下步骤：

(1)模拟地物在探测器入瞳出的辐亮度信息

基于朗伯均一地表假定，探测器所接收的辐射亮度L_sensor表示为：

$L_{\mathrm{sensor}}=L_p+\frac{{\mathrm{\μ}}_s{E}_0{\mathrm{\ρ}}_{s1}[e^{-\mathrm{\τ}{/\mathrm{\μ}}_s}+t_d\left({\mathrm{\μ}}_s\right)][e^{-\mathrm{\τ}/{\mathrm{\μ}}_v}+t_d\left({\mathrm{\μ}}_v\right)]}{\mathrm{\π}(1-{\mathrm{\ρ}}_{s1}S)}---\left(1\right)$

式中：E₀为大气层顶的太阳辐照度；S为大气半球反射率；为从太阳到地表路径上的直射透过率；μ_s＝cosθ_s为太阳天顶角θ_s的余弦值；为从地表到探测器路径上的直射透过率；μ_v＝cosθ_v为观测天顶角θ_v的余弦值；τ为大气光学厚度；t_d(μ_s)为从太阳到地表路径上的散射透过率；t_d(μ_v)为从地表到探测器路径上的散射透过率；L_p为大气路径辐射；ρ_s1为朗伯均一地表反射率。

通过公式(1)可以基于朗伯均一地表反射率ρ_s1推演得到探测器所接收的辐射亮度L_sensor需要知道6个大气参数t_d(μ_s)、t_d(μ_v)、L_p和S。这6个大气参数与地表无关，只与大气状况有关。大气层顶的太阳辐照度E₀可以通过大气层顶的太阳辐照度曲线与探测器的通道响应函数卷积计算得到。

给定地表反射率ρ_s1＝0，ρ_s1＝0.5，ρ_s1＝1，代入大气辐射传输模型软件MODTRAN运算3次可以得到以下参数：

·从地表到探测器路径上的直射透过率

从地表到探测器路径上的直射透过率可以通过MODTRAN输出的结果与探测器的通道响应函数卷积得到。

·从太阳到地表路径上的直射透过率

根据公式从太阳到地表路径上的直射透过率表示为：

$e^{-{\mathrm{\τ}/u}_s}=\frac{{\mathrm{\πL}}_{\mathrm{drct}}\left({\mathrm{\ρ}}_3\right)}{{\mathrm{\μ}}_s{E}_0{\mathrm{\ρ}}_3{e}^{-\mathrm{\τ}/u_v}}---\left(2\right)$

式中，L_drct为探测器接收的来自地表的直射辐射，可以通过MODTRAN输出的结果与探测器的通道响应函数卷积得到。

·大气路径辐射L_p

大气路径辐射L_p可直接由地表反射率ρ₁＝0的路径辐射得到：

L_p＝L_path(ρ₁) (3)

式中，L_path为探测器接收的路径辐射，可以通过MODTRAN输出的结果与探测器的通道响应函数卷积得到。

·从太阳到地表路径上的散射透过率t_d(μ_s)

给定地表反射率ρ₂＝0.5或者ρ₃＝1.0，根据公式和从太阳到地表路径上的散射透过率t_d(μ_s)表示为：

$t_d\left({\mathrm{\μ}}_s\right)=(\frac{L_{\mathrm{grnd}}\left({\mathrm{\ρ}}_3\right)(1-{\mathrm{\ρ}}_{3}S)}{L_{\mathrm{drct}}\left({\mathrm{\ρ}}_3\right)}-1)e^{-{\mathrm{\τ}/u}_s}---\left(4\right)$

式中，L_grnd为探测器接收的来自像元本身的总辐射，可以通过MODTRAN输出的结果与探测器的通道响应函数卷积得到。

·从地表到探测器路径上的散射透过率t_d(μ_v)

给定地表反射率ρ₂＝0.5或者ρ₃＝1.0，根据公式 $L_{\mathrm{path}}\left(\mathrm{\ρ}\right)=\frac{{\mathrm{\μ}}_s{E}_0}{\mathrm{\π}}\frac{(e^{{-\mathrm{\τu}}_s}+t_d\left({\mathrm{\μ}}_s\right))t_d\left({\mathrm{\μ}}_v\right)\mathrm{\ρ}}{1-\mathrm{\ρS}}+L_p$ 和 $L_{\mathrm{grnd}}\left(\mathrm{\ρ}\right)=\frac{{\mathrm{\μ}}_s{E}_0}{\mathrm{\π}}\frac{(e^{{-\mathrm{\τu}}_s}+t_d\left({\mathrm{\μ}}_s\right))e^{-\mathrm{\τ}/u_v}\mathrm{\ρ}}{1-\mathrm{\ρS}},$ 从地表到探测器路径上的散射透过率t_d(μ_v)表示为：

$t_d\left({\mathrm{\μ}}_v\right)=\frac{L_{\mathrm{path}}\left({\mathrm{\ρ}}_3\right)-L_p}{L_{\mathrm{grnd}}\left({\mathrm{\ρ}}_3\right)}{e}^{-{\mathrm{\τ}/u}_v}---\left(5\right)$

·大气半球反照率S

给定地表反射率ρ₂＝0.5和ρ₃＝1.0，根据公式

$L_{\mathrm{grnd}}\left(\mathrm{\ρ}\right)=\frac{{\mathrm{\μ}}_s{E}_0(e^{{-\mathrm{\τ}/u}_s}+t_d\left({\mathrm{\μ}}_s\right))e^{{-\mathrm{\τ}/u}_v}\mathrm{\ρ}}{1-\mathrm{\ρS}}$

大气半球反照率S表示为：

$S=\frac{{\mathrm{\ρ}}_2{L}_{\mathrm{grnd}}\left({\mathrm{\ρ}}_3\right)-{\mathrm{\ρ}}_3{L}_{\mathrm{grnd}}\left({\mathrm{\ρ}}_2\right)}{{\mathrm{\ρ}}_2{\mathrm{\ρ}}_3[L_{\mathrm{grnd}}\left({\mathrm{\ρ}}_3\right)-L_{\mathrm{grnd}}\left({\mathrm{\ρ}}_2\right)]}---\left(6\right)$

由公式1～6可计算得到地物在探测器入瞳处的辐亮度值。

(2)将步骤(1)得到的地物在探测器入瞳处的辐亮度值合成XYZ颜色值；

在380-780nm可见光波长范围按式(7)对高光谱遥感图像做加权积分。

$\left\{\begin{array}{c}X=k{\∫}_{380}^{780}{L}_{\mathrm{\λ}}{\stackrel{\‾}{x}}_{\mathrm{\λ}}\mathrm{d\λ}\\ Y=k{\∫}_{380}^{780}{L}_{\mathrm{\λ}}{\stackrel{\‾}{y}}_{\mathrm{\λ}}\mathrm{d\λ}\\ Z=k{\∫}_{380}^{780}{L}_{\mathrm{\λ}}{\stackrel{\‾}{z}}_{\mathrm{\λ}}\mathrm{d\λ}\end{array}\right.---\left(7\right)$

$k=1/{\∫}_{380}^{780}{S}_{\mathrm{\λ}}{\stackrel{\‾}{y}}_{\mathrm{\λ}}\mathrm{d\λ}---\left(8\right)$

其中，为对应波长的颜色匹配函数值，S_λ为入射光源能量，L_λ为地物在λ波长处的辐亮度值，它与地物反射率的关系如下：

L_λ＝E_λρ_λ (9)其中，ρ_λ为地面靶标反射率，E_λ是地面太阳入射光辐照度。

由于应探测器通道的中心波长与CIE1931颜色匹配函数表中心波长有可能不会一一对应，因此需要利用公式(10)将探测器各个通道的辐亮度值映射到与CIE1931颜色匹配函数相对应的中心波长上，若λ₁＜λ＜λ₂则插补公式如下：

${\mathrm{\φ}}_{\mathrm{\λ}}=\frac{[L_{\mathrm{\λ}1}{({\mathrm{\λ}}_2-\mathrm{\λ})}^2+L_{\mathrm{\λ}2}{\left({\mathrm{\λ}-\mathrm{\λ}}_1\right)}^2]}{{({\mathrm{\λ}}_2-\mathrm{\λ})}^2+{\left(\mathrm{\λ}{-\mathrm{\λ}}_1\right)}^2}---\left(10\right)$

λ1、λ2为探测器通道中心波长，λ为任意波长，求得的Φ_λ在λ处的辐照度，可作为一个新的通道辐射能量参与下一步运算。

(3)将步骤(2)中得到的XYZ值转为真彩色RGB值

XYZ三色刺激值并不是人眼对短、中和长波的真实反应，因此需要将合成后的XYZ彩色图像按式(11)转换为RGB真彩色图像。采用CIE1931标准在接近太阳光的D65光源下的XYZ-RGB色彩空间转换矩阵：

$\left[\begin{array}{c}R\\ G\\ B\end{array}\right]{=M}_{\mathrm{XYZ}->\mathrm{RGB}}\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]---\left(11\right)$

其中，M_XYZ-＞RGB为CIE1931标准在接近太阳光的D65光源下的XYZ-RGB色彩空间转换矩阵：

$M_{\mathrm{XYZ}->\mathrm{RGB}}=\left[\begin{array}{ccc}3.241& -1.5374& -0.4986\\ -0.9692& 1.876& 0.0416\\ 0.0556& -0.204& 1.057\end{array}\right]---\left(12\right)$

步骤二、提取步骤一中地物在多光谱遥感图像中的RGB颜色值，建立同一个地物多光谱图像RGB值与真彩色RGB值之间的图像真彩色转换矩阵；

提取步骤一中地物在多光谱遥感图像中的红、绿、蓝通道灰度值作为其真实RGB颜色值并与真彩色RGB建立如式(13)的对应关系。

$\left(\begin{array}{c}{R}_{\mathrm{real}}\\ G_{\mathrm{real}}\\ B_{\mathrm{real}}\end{array}\right)=M_{\mathrm{multi}\→\mathrm{real}}\left(\begin{array}{c}{R}_{\mathrm{multi}}\\ G_{\mathrm{multi}}\\ B_{\mathrm{multi}}\end{array}\right)---\left(13\right)$

其中， $\left(\begin{array}{c}{R}_{\mathrm{real}}\\ G_{\mathrm{real}}\\ B_{\mathrm{real}}\end{array}\right)$ 和 $\left(\begin{array}{c}{R}_{\mathrm{multi}}\\ G_{\mathrm{multi}}\\ B_{\mathrm{multi}}\end{array}\right)$ 分别为地物的真实RGB值及在多光谱图像上的RGB值。M_multi→real即表征了二者转换的定标矩阵，可利用求解超定方程组的数学方法可以求解该矩阵。

步骤三、利用步骤二中得到的图像真彩色转换矩阵对真彩色图像进行定标，得到更为接近地物真是色彩的影像。

实施例：

下面结合一个实例对本发明做进一步说明。

本例选取2011年9月获取的无人机多光谱影像作为实验数据，影像具有4个通道，其中前三个通道为蓝、绿、红颜色通道。影像获取时在试验场地铺设了多种反射率和颜色的靶标11处。

按照本发明阐述的实施方式，首先提取上述11处靶标真实的RGB值：

步骤一：根据现场测量的样本靶标1～11光谱反射率信息及大气参数，带入MORTRAN4.0大气辐射传输模型中，求得样本靶标在探测器入瞳处的光谱辐亮度信息。

步骤二：根据CIE1931颜色匹配函数公式(7)，对上一步骤求得的各靶标入瞳辐射亮度值进行加权积分，得到样本靶标的XYZ彩色值。公式如下：

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}k\underset{\mathrm{\λ}=380}{\overset{780}{\mathrm{\Σ}}}{{\mathrm{\φ}}_{\mathrm{\λ}}}^{*}{\stackrel{\‾}{x}}_{\mathrm{\λ}}\\ k\underset{\mathrm{\λ}=380}{\overset{780}{\mathrm{\Σ}}}{{\mathrm{\φ}}_{\mathrm{\λ}}}^{*}{\stackrel{\‾}{y}}_{\mathrm{\λ}}\\ k\underset{\mathrm{\λ}=380}{\overset{780}{\mathrm{\Σ}}}{{\mathrm{\φ}}_{\mathrm{\λ}}}^{*}{\stackrel{\‾}{z}}_{\mathrm{\λ}}\end{array}\right]---\left(14\right)$

$k=1/\underset{\mathrm{\λ}=380}{\overset{780}{\mathrm{\Σ}}}{{\mathrm{\φ}}_{\mathrm{\λ}}}^{*}{\stackrel{\‾}{y}}_{\mathrm{\λ}}$

为了与人眼视觉颜色感知相同，利用公式(11)，将XYZ彩色值转换为RGB真彩色值矩阵。

$\left[\begin{array}{c}{R}_{\mathrm{real}}\\ G_{\mathrm{real}}\\ B_{\mathrm{real}}\end{array}\right]=\left[\begin{array}{ccccccccccc}56& 93& 109& 132& 152& 255& 0& 18& 171& 16& 15\\ 54& 83& 103& 120& 15& 167& 70& 97& 162& 17& 16\\ 55& 82& 105& 119& 25& 15& 43& 183& 165& 22& 20\end{array}\right]---\left(14\right)$

步骤三：选取多光谱数据的第3，2，1波段作为红、绿、蓝通道。

在多光谱图像中提取11个样本靶标DN值，转换为字节型的RGB彩色值：

$\left[\begin{array}{c}{R}_{\mathrm{multi}}\\ G_{\mathrm{multi}}\\ B_{\mathrm{multi}}\end{array}\right]=\left[\begin{array}{ccccccccccc}124& 168& 193& 213& 247& 252& 88& 115& 252& 74& 72\\ 152& 200& 235& 252& 127& 252& 139& 186& 255& 89& 86\\ 92& 110& 124& 129& 71& 92& 107& 50& 151& 67& 63\end{array}\right]---\left(15\right)$

根据公式(13)，最为优选使用最小二乘法进行拟合，获得多光谱图像RGB值与真实RGB值之间的转换矩阵。

$M_{\mathrm{multi}-\mathrm{real}}=\left[\begin{array}{ccc}{M}_R& M_G& M_B\end{array}\right]=\left[\begin{array}{ccc}0.8037& -0.3087& -0.1917\\ 0.9244& 1.2272& -0.1442\\ -1.7802& -0.6611& 1.4644\end{array}\right]---\left(16\right)$

步骤四：根据真彩色转换矩阵(16)，对同一时间的地区获取的农田地物多光谱图像进行真彩色校正运算。

在未经过真彩色校正的多光谱图像中，虽然可以根据形状和颜色大致判断出农田区域，但是图像颜色整体偏绿，区别度不大，无法明显的辨别农田、沙地以及不同的农作物之间的区别。

而根据真彩色转换矩阵(16)对原始多光谱图像进行真彩色校正运算后，图像整体颜色会更加自然和丰富，沙地和道路的颜色能够得到较好的校正效果。可见本专利阐述的方法对于遥感图像的真彩色校正具有一定的效果。本发明未说明部分属于本领域技术人员公知常识。

CIE1931XYZ颜色匹配函数列表

Claims (1)

1.一种多光谱遥感数据真彩色定标方法，其特征在于，包括以下步骤：

步骤一、模拟地物在探测器入瞳出的辐亮度信息，计算出对应的RGB真彩色值，具体过程包括如下步骤：

(1)模拟地物在探测器入瞳处的辐亮度信息

基于朗伯均一地表假定，探测器所接收的辐射亮度L_sensor表示为：

式中：E₀为大气层顶的太阳辐照度；S为大气半球反射率；为从太阳到地表路径上的直射透过率；μ_s＝cosθ_s为太阳天顶角θ_s的余弦值；为从地表到探测器路径上的直射透过率；μ_v＝cosθ_v为观测天顶角θ_v的余弦值；τ为大气光学厚度；t_d(μ_s)为从太阳到地表路径上的散射透过率；t_d(μ_v)为从地表到探测器路径上的散射透过率；L_p为大气路径辐射；ρ_s1为朗伯均一地表反射率；

通过公式(1)可以基于朗伯均一地表反射率ρ_s1推演得到探测器所接收的辐射亮度L_sensor需要知道6个大气参数t_d(μ_s)、t_d(μ_v)、L_p和S ；这6个大气参数与地表无关，只与大气状况有关 ；大气层顶的太阳辐照度E₀可以通过大气层顶的太阳辐照度曲线与探测器的通道响应函数卷积计算得到；

给定地表反射率ρ_s1＝0，ρ_s1＝0.5，ρ_s1＝1，代入大气辐射传输模型软件MODTRAN运算3次可以得到以下参数：

·从地表到探测器路径上的直射透过率

从地表到探测器路径上的直射透过率可以通过MODTRAN输出的结果与探测器的通道响应函数卷积得到；

·从太阳到地表路径上的直射透过率

根据公式从太阳到地表路径上的直射透过率表示为：

式中，L_drct为探测器接收的来自地表的直射辐射，可以通过MODTRAN输出的结果与探测器的通道响应函数卷积得到；

·大气路径辐射L_p

大气路径辐射L_p可直接由地表反射率ρ₁＝0的路径辐射得到：

L_p＝L_path(ρ₁) (3)

式中，L_path为探测器接收的路径辐射，可以通过MODTRAN输出的结果与探测器的通道响应函数卷积得到；

·从太阳到地表路径上的散射透过率t_d(μ_s)

给定地表反射率ρ₂＝0.5或者ρ₃＝1.0，根据公式和从太阳到地表路径上的散射透过率t_d(μ_s)表示为：

式中，L_grnd为探测器接收的来自像元本身的总辐射，可以通过MODTRAN输出的结果与探测器的通道响应函数卷积得到；

·从地表到探测器路径上的散射透过率t_d(μ_v)

给定地表反射率ρ₂＝0.5或者ρ₃＝1.0，根据公式和从地表到探测器路径上的散射透过率t_d(μ_v)表示为：

·大气半球反照率S

给定地表反射率ρ₂＝0.5和ρ₃＝1.0，根据公式

大气半球反照率S表示为：

由公式1～6可计算得到地物在探测器入瞳处的辐亮度值；

(2)将步骤(1)得到的地物在探测器入瞳处的辐亮度值合成XYZ颜色值；

在380-780nm可见光波长范围按式(7)对高光谱遥感图像做加权积分；

其中为对应波长的颜色匹配函数值，S_λ为入射光源能量，地物在λ波长处的辐亮度值，它与地物反射率的关系如下：

L_λ＝E_λρ_λ (9)

其中，ρ_λ为地面靶标反射率，E_λ是地面太阳入射光辐照度；

由于应探测器通道的中心波长与CIE1931颜色匹配函数表中心波长有可能不会一一对应，因此需要利用公式(10)将探测器各个通道的辐亮度值映射到与CIE1931颜色匹配函数相对应的中心波长上，若λ₁＜λ＜λ₂则插补公式如下：

λ1、λ2为探测器通道中心波长，λ为任意波长，求得的Φ_λ在λ处的辐照度，可作为一个新的通道辐射能量参与下一步运算；

(3)将步骤(2)中得到的XYZ值转为真彩色RGB值

XYZ三色刺激值并不是人眼对短、中和长波的真实反映，因此需要将合成后的XYZ彩色图像按式(11)转换为RGB真彩色图像 ；采用CIE1931标准在接近太阳光的D65光源下的XYZ-RGB色彩空间转换矩阵：

其中，M_XYZ-＞RGB为CIE1931标准在接近太阳光的D65光源下的XYZ-RGB色彩空间转换矩阵：

步骤二、提取步骤一中地物在多光谱遥感图像中的RGB颜色值，建立同一个地物多光谱图像RGB值与真彩色RGB值之间的图像真彩色转换矩阵；

提取步骤一中地物在多光谱遥感图像中的红、绿、蓝通道灰度值作为其真实RGB颜色值并与真彩色RGB建立如式(13)的对应关系 ；

其中，和分别为地物的真实RGB值及在多光谱图像上的RGB值；M_multi→real即表征了二者转换的定标矩阵，可利用求解超定方程组的数学方法可以求解该矩阵；

步骤三、利用步骤二中得到的图像真彩色转换矩阵对真彩色图像进行定标，得到更为接近地物真是色彩的影像。