CN102736128A - 无人机光学遥感影像数据处理方法及装置 - Google Patents

Abstract

一种无人机光学遥感影像数据处理方法，其步骤包括：A)无人机光学遥感影像数据预处理；B)对无人机光学遥感影像数据进行大气水汽反演，得到大气水汽含量产品；C)利用步骤B的大气水汽含量产品，对无人机光学遥感影像数据进行地表反射率反演，得到地表反射率产品；D)利用步骤B的大气水汽含量产品和步骤C的地表反射率产品，对无人机光学遥感影像数据进行地表温度反演，得到地表温度产品。

Description

无人机光学遥感影像数据处理方法及装置

技术领域

本发明涉及一种遥感影像数据处理方法，具体地涉及一种无人机光学遥感影像数据处理方法。

本发明还涉及一种实现上述方法的装置。

背景技术

自上世纪60年代民用空间对地观测技术发展以来，遥感载荷技术正朝着追求高分辨率和应用实用化的方向发展。无人机遥感作为一种新型的遥感数据获取手段，具有续航时间长、影像实时传输、高危地区探测、成本低、高分辨率、机动灵活等优点，是卫星遥感与有人机航空遥感的有力补充，可广泛应用于国家生态环境保护、矿产资源勘探、海洋环境监测、土地利用调查和国土资源监测等领域。因此，无人机遥感已经成为世界各国争相研究的热点课题，现已逐步从研究开发发展到实际应用阶段，将成为未来主要航空遥感技术之一。

遥感反演地表参数的精度是衡量遥感载荷应用性能的重要方面，对于地表参数反演，机载和星载的主要区别是传感器测得的辐射受到大气层的影响不一样，对于机载传感器，地物的反射辐射和发射辐射仅受到从地面到飞机飞行高度这一层大气的影响，而星载传感器受到了整个大气层的影响。目前，利用机载遥感技术反演地表反射率，难点主要在大气校正、双向反射分布函数校正、邻近效应校正三个方面，特别是邻近效应校正，由于机载数据的空间分辨率高，因而邻近效应的影响更为显著。在机载数据的大气水汽反演方面，由于仅受到从地表到飞机飞行高度这一层大气的影响，需要进行从地表到飞机飞行高度这一层的大气水汽到整个大气层大气水汽的转换。目前关于地表温度的机载遥感反演算法，可以根据传感器通道数的不同，划分为单通道算法、分裂窗算法和多通道算法。然而，利用传感器探测到的地物发射辐射推断出相应的地表温度和比辐射率却非常棘手。其中地表温度和通道比辐射率都是待求参数，但是N个通道观测的热辐射却对应N+1个未知数(N个通道比辐射率和1个地表温度)，由此导致地表温度和比辐射率的反演成为一个欠定问题。

发明内容

本发明的目的在于提供一种无人机光学遥感影像数据处理方法。

本发明的又一目的在于提供一种用于实现上述方法的装置。

为实现上述目的，本发明提供的无人机光学遥感影像数据处理方法，其步骤包括：

A)无人机光学遥感影像数据预处理；

B)对无人机光学遥感影像数据进行大气水汽反演，得到大气水汽含量产品；

C)利用步骤B的大气水汽含量产品，对无人机光学遥感影像数据进行地表反射率反演，得到地表反射率产品；

D)利用步骤B的大气水汽含量产品和步骤C的地表反射率产品，对无人机光学遥感影像数据进行地表温度反演，得到地表温度产品。

所述的无人机光学遥感影像数据处理方法，其中，对无人机光学遥感影像数据进行预处理过程是：

(1)根据影像数据本身，对无人机光学遥感影像数据进行相对辐射定标，消除由于各个探测元的响应差异在影像上产生的条纹；

(2)根据无人机传感器的辐射定标系数，对无人机光学遥感影像数据进行绝对辐射定标，将影像灰度值转换为辐射亮度。

所述的无人机光学遥感影像数据处理方法，其中，对无人机光学遥感影像数据进行大气水汽反演过程是：

(1)根据水汽通道和相邻两个大气窗口通道的辐射亮度计算比值指数；

(2)根据比值指数和拟合系数计算大气斜程水汽；

(3)根据大气斜程水汽和大气水汽比计算大气垂直水汽。

所述的无人机光学遥感影像数据处理方法，其中，对无人机光学遥感影像数据进行地表反射率反演过程是：

(1)根据气溶胶光学厚度、大气水汽含量、飞行高度、地表高程、太阳天顶角、观测天顶角和相对方位角参数及大气参数查找表，计算得到朗伯均一地表反射率；

(2)根据朗伯均一地表反射率计算得到周围环境平均地表反射率；

(3)根据朗伯均一地表反射率和周围环境平均反射率逐像元计算得到朗伯非均一地表反射率；

(4)根据步骤3得到的朗伯非均一地表反射率再计算周围环境平均地表反射率，重复步骤2和步骤3进行迭代，直到两次计算的朗伯非均一地表反射率之差小于设定的阈值，迭代结束。

所述的无人机光学遥感影像数据处理方法，其中，对无人机光学遥感影像数据进行地表温度反演过程是：

(1)根据普朗克函数将辐射亮度转换为亮度温度；

(2)根据地表反射率计算归一化植被指数，并计算植被覆盖度；

(3)根据归一化植被指数和植被覆盖度计算地表比辐射率；

(4)根据地表比辐射率、大气水汽含量和亮度温度计算地表温度。

本发明提供的实现上述方法的装置，其包括：

数据预处理模块进行相对辐射定标和绝对辐射定标；以及

大气水汽反演模块计算大气水汽含量；以及

地表反射率反演模块计算地表反射率；以及

地表温度反演模块计算地表温度。

本发明的无人机光学遥感影像数据处理方法是通过大气水汽比来实现从地表到飞机飞行高度这一层的大气水汽到整个大气层的大气水汽的转换。

本发明的无人机光学遥感影像数据处理方法是通过对朗伯均一地表反射率进行迭代来实现邻近效应校正，得到精度更高，更接近真实值的朗伯非均一地表反射率。

本发明的无人机光学遥感影像数据处理方法是通过建立普适性分裂窗算法来实现不需要辅助数据只从数据本身对热红外数据进行大气校正。

本发明的无人机光学遥感影像数据处理装置是通过数据预处理模块、大气水汽反演模块、地表反射率反演模块和地表温度反演模块来实现从可见光到热红外光谱范围的大气和地表参数的一体化反演，具有操作简单、实用性强、精度高及可扩展性的特点。

附图说明

图1是本发明无人机光学遥感影像数据处理流程示意图。

图2是本发明无人机大气水汽反演算法系数拟合流程示意图。

图3是本发明无人机大气水汽含量反演流程示意图。

图4是本发明无人机可见光及近红外地表反射率反演流程示意图。

图5是本发明无人机热红外地表比辐射率计算流程示意图。

图6是本发明采用的普适性分裂窗算法系数拟合流程示意图。

图7是本发明无人机热红外地表温度反演流程示意图。

具体实施方式

本发明提供的无人机光学遥感影像数据处理方法包括：数据预处理模块、大气水汽反演模块、地表反射率反演模块、地表温度反演模块。

1、数据预处理模块

(1)相对辐射定标

相对辐射定标是通过调整每个探测元的响应值，消去偏移值带来的影响，将各个探测元的输出值调整到同一个基准上，使各个探测元对相同的输入产生相同的输出值：

Y＝(X-μ_i)σ_r/σ_i+μ_r    (1)

式中，X和Y为变换前后的影像灰度值；μ_i为第i列影像灰度值的平均值；σ_i为第i列影像灰度值的标准差；μ_r为整幅图像影像灰度值的平均值；σ_r为整幅图像影像灰度值的标准差。

(2)绝对辐射定标

无人机传感器获取的数据为影像灰度值，需要经过辐射定标将影像灰度值转换为辐射亮度：

L_sensor＝DN*gain+offset    (2)

式中，DN为影像灰度值；L_sensor为辐射亮度，单位为Wm^-2sr^-1μm^-1；gain和offset分别为增益值和偏移值，单位为Wm^-2sr^-1μm^-1。

2、大气水汽反演模块

(1)大气水汽反演算法

通过940nm处的水汽通道及其两侧的大气窗口通道来计算比值指数R_index：

$R_{\mathrm{index}}=\frac{L_m}{w_{r1}{L}_{r1}+w_{r2}{L}_{r2}}---\left(3\right)$

其中， $w_{r1}=\frac{{\mathrm{\λ}}_{w_{r2}}-{\mathrm{\λ}}_m}{{\mathrm{\λ}}_{w_{r2}}-{\mathrm{\λ}}_{w_{r1}}}$ 和 $w_{r2}=\frac{{\mathrm{\λ}}_m-{\mathrm{\λ}}_{w_{r1}}}{{\mathrm{\λ}}_{w_{r2}}-{\mathrm{\λ}}_{w_{r1}}}.$

式中，L_m为水汽通道的辐射亮度；L_r1和L_r2为大气窗口通道的辐射亮度；w_r1和w_r2为权重系数；λ_m为水汽通道的中心波长；λ_r1和λ_r2为大气窗口通道的中心波长。

通过比值指数R_index计算大气水汽含量W_p：

R_index＝exp(-a+b(W_p)^c)    (4)

式中，a、b和c为拟合系数。

(2)拟合系数的确定

从ASTER波谱库(URL：http://speclib.jpl.nasa.gov/)中选取6条地物波谱，包括：裸土、植被、雪和沙，给定不同的地表和大气条件以及太阳和观测几何，利用大气辐射传输模型MODTRAN，模拟得到不同大气和地表条件下水汽通道和大气窗口通道的辐射亮度，通过公式(3)计算得到比值指数R_index。根据比值指数R_index和模拟时输入的大气水汽含量，通过公式(4)拟合得到系数a，b和c(如图2所示)。

(3)斜程水汽到垂直水汽的转换

从公式(4)得到的大气水汽含量为斜程水汽，要得到垂直水汽，需要经过斜程水汽到垂直水汽的转换。与星载传感器不同，无人机传感器接收到的来自地物的反射辐射仅受到从地表到飞机飞行高度这一层大气的影响，而星载传感器受到了整个大气层的影响。因此，无人机传感器的斜程水汽W_p表示为：

$W_p=\frac{W_v}{\cos {\mathrm{\θ}}_s}+\frac{W_h}{\cos {\mathrm{\θ}}_v}---\left(5\right)$

式中，W_v为从太阳到地表路径上的垂直水汽；W_h为从地表到无人机传感器路径上的垂直水汽；θ_s为太阳天顶角；θ_v为观测天顶角。

大气水汽比R_w定义为从大气层顶到某一高度上的垂直水汽与从大气层顶到地表的垂直水汽之比：

$R_w=\frac{W_s}{W_v}=\frac{W_v-W_h}{W_v}---\left(6\right)$

式中，W_s为从大气层顶到无人机传感器高度上的垂直水汽。

根据斜程水汽W_p和大气水汽比R_w，从太阳到地表路径上的垂直水汽W_v表示为：

$W_v=\frac{W_p}{\frac{1}{\cos {\mathrm{\θ}}_s}+\frac{1-R_w}{\cos {\mathrm{\θ}}_v}}---\left(7\right)$

大气水汽含量反演流程如图3所示，其步骤为：

1)根据公式(3)，通过无人机传感器940nm处的水汽通道和相邻两个大气窗口通道的辐射亮度计算比值指数R_index。

2)根据公式(4)，通过比值指数R_index和拟合系数a，b和c计算斜程水汽W_p。

3)根据公式(7)，通过斜程水汽W_p和水汽比R_w计算大气垂直水汽W_v。

3、地表反射率反演模块

(1)大气参数查找表

从无人机传感器获得的辐射亮度转换为地表反射率需要经过大气校正。利用大气辐射传输模型逐像元进行大气校正是非常复杂的计算，需要占用大量的计算机时间和资源，所以要在业务化运行系统中实时地调用大气辐射传输模型几乎是不可能的，查找表可用于解决这个矛盾。根据无人机机载传感器的特性，基于大气辐射传输模型MODTRAN，建立以气溶胶光学厚度、大气水汽含量、飞行高度、地表高程、太阳天顶角、观测天顶角和相对方位角为索引的多维大气参数查找表。表1为大气参数查找表的模型输入参数的设置。

表1  大气参数查找表的模型输入参数设置

(2)大气校正

基于朗伯均一地表假定，传感器所接收的辐射亮度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(9\right)$

式中：E₀为大气层顶的太阳辐照度；S为大气半球反射率；

为从太阳到地表路径上的直射透过率；μ_s＝cosθ_s为太阳天顶角θ_s的余弦值；

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

通过公式(9)可以看出，从传感器所接收的辐射亮度L_sensor反演得到朗伯均一地表反射率ρ_s1需要知道6个大气参数

t_d(μ_s)、t_d(μ_v)、L_p和S。这6个大气参数与地表无关，只与大气状况有关。大气层顶的太阳辐照度E₀可以通过大气层顶的太阳辐照度曲线与传感器的通道响应函数卷积计算得到。

根据6个大气参数

t_d(μ_s)、t_d(μ_v)、L_p和S，朗伯均一地表反射率ρ_s1表示为：

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

(3)大气参数计算

给定地表反射率ρ₁＝0、ρ₂＝0.5和ρ₃＝1.0，运行3次MODTRAN可以得到以下参数：

·从地表到传感器路径上的直射透过率

从地表到传感器路径上的直射透过率

可以通过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(11\right)$

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

·大气路径辐射L_p

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

L_p＝L_path(ρ₁)    (12)

式中，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(13\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{\τ}/{\mathrm{\μ}}_v}---\left(14\right)$

·大气半球反照率S

给定地表反射率ρ₂＝0.5和ρ₃＝1.0，根据公式 $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}},$ 大气半球反照率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(15\right)$

(4)邻近效应校正

基于朗伯非均一地表假定，传感器所接收的辐射亮度L_sebsor表示为：

$L_{\mathrm{sensor}}=L_p+\frac{{\mathrm{\&mu;}}_s{E}_0[e^{-\mathrm{\&tau;}/{\mathrm{\&mu;}}_s}+t_d\left({\mathrm{\&mu;}}_s\right)][e^{-\mathrm{\&tau;}{/\mathrm{\&mu;}}_v}{\mathrm{\&rho;}}_{s2}+t_d\left({\mathrm{\&mu;}}_v\right)<\mathrm{\&rho;}>]}{\mathrm{\&pi;}(1-<\mathrm{\&rho;}>S)}---\left(16\right)$

根据6个大气参数

t_d(μ_s)、t_d(μ_s)、L_p和S，朗伯非均一地表反射率ρ_s2表示为：

${\mathrm{\&rho;}}_{s2}=\frac{\mathrm{\&pi;}(L_{\mathrm{sensor}}-L_p)(1-<\mathrm{\&rho;}>S)}{{\mathrm{\&mu;}}_s{E}_0{e}^{-\mathrm{\&tau;}/{\mathrm{\&mu;}}_v}[e^{-\mathrm{\&tau;}/{\mathrm{\&mu;}}_s}+t_d\left({\mathrm{\&mu;}}_s\right)]}-\frac{t_d\left({\mathrm{\&mu;}}_v\right)<\mathrm{\&rho;}>}{e^{-\mathrm{\&tau;}/{\mathrm{\&mu;}}_v}}---\left(17\right)$

式中，<ρ>为周围环境平均地表反射率。根据朗伯均一地表反射率ρ_s1，由下式计算得到：

$<\mathrm{\&rho;}>=\frac{1}{{(2N+1)}^2}\underset{x=-N}{\overset{N}{\mathrm{\&Sigma;}}}\underset{y=-N}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_{s1}(x,y)---\left(18\right)$

地表反射率反演流程如图4所示，其步骤为：

1)根据气溶胶光学厚度、大气水汽含量、飞行高度、地表高程、太阳天顶角、观测天顶角和相对方位角7个输入参数及大气参数查找表，逐像元插值计算得到6个大气参数t_d(μ_s)、t_d(μ_v)、L_p和S，通过公式(10)计算得到朗伯均一地表反射率ρ_s1。

2)根据朗伯均一地表反射率ρ_s1，通过公式(18)计算得到周围环境平均地表反射率<ρ>。

3)根据朗伯均一地表反射率ρ_s1和周围环境平均反射率<ρ>，结合6个大气参数

t_d(μ_s)、t_d(μ_v)、L_p和S，通过公式(17)逐像元计算得到朗伯非均一地表反射率ρ_s2。

4)将第3步得到的朗伯非均一地表反射率ρ_s2代入公式(18)再计算周围环境平均地表反射率<ρ>，重复步骤2和3，进行迭代，直到两次计算的朗伯非均一地表反射率ρ_s2之差小于设定的阈值，迭代结束。

4、地表温度反演模块

(1)亮度温度计算

由传感器所接收的辐射亮度L_sensor和普朗克函数求出传感器的亮度温度T_sensor：

$T_{\mathrm{sensor}}=\frac{C_2}{{\mathrm{\&lambda;}}_c\ln (\frac{C_1}{{\mathrm{\&lambda;}}_c^5{L}_{\mathrm{sensor}}}+1)}---\left(19\right)$

式中，C₁为第一辐射常数，等于1.191×10⁸W(μm)⁴sr^-1m^-2；C₂为第二辐射常数，等于1.439×10⁴μmK；λ_c为中心波长，单位为μm。

(2)归一化植被指数和植被覆盖度计算

归一化植被指数(NDVI)根据地表反射率计算得到：

$\mathrm{NDVI}=\frac{{\mathrm{\&rho;}}_{\mathrm{nir}}-{\mathrm{\&rho;}}_{\mathrm{red}}}{{\mathrm{\&rho;}}_{\mathrm{nir}}+{\mathrm{\&rho;}}_{\mathrm{red}}}---\left(20\right)$

式中，ρ_nir为近红外波段的地表反射率；ρ_red为可见光红波段的地表反射率。

植被覆盖度P_v根据NDVI由以下公式计算得到：

$P_v={\left[\frac{{\mathrm{NDVI}-\mathrm{NDVI}}_s}{{\mathrm{NDVI}}_v-{\mathrm{NDVI}}_s}\right]}^2---\left(21\right)$

式中，NDVI_s为裸土的NDVI值，取值0.2；NDVI_v为植被的NDVI值，取值0.5。

(3)地表比辐射率的确定

要从无人机热红外遥感数据中反演得到地表温度，必须要先知道地表比辐射率。无人机热红外遥感数据的地表比辐射率可以根据NDVI和P_v计算得到(如图5所示)：

$\mathrm{\&epsiv;}=\left\{\begin{array}{cc}{\mathrm{\&epsiv;}}_s,& \mathrm{NDVI}<{\mathrm{NDVI}}_s\\ {\mathrm{\&epsiv;}}_s+({\mathrm{\&epsiv;}}_v-{\mathrm{\&epsiv;}}_s)P_v,& {\mathrm{NDVI}}_s\&le;\mathrm{NDVI}\&le;{\mathrm{NDVI}}_v\\ {\mathrm{\&epsiv;}}_v,& \mathrm{NDVI}>{\mathrm{NDVI}}_v\end{array}\right.---\left(22\right)$

式中，ε_v为植被的比辐射率，取值0.99；ε_s为裸土的平均比辐射率。从ASTER波谱库(URL：http://speclib.jpl.nasa.gov/)中选取49条裸土波谱，通过与无人机热红外每个通道的通道响应函数进行卷积，计算每个通道裸土的平均比辐射率。

(4)普适性分裂窗地表温度反演算法

根据大气辐射传输方程，利用10-13um大气窗口内，无人机传感器两个相邻热红外通道(10.5-11.5um、11.5-12.5um)对大气吸收作用的不同，通过这两个通道的测量值来剔除大气的影响，进行大气和地表比辐射率的修正：

$T_s=C+(A_1+A_2\frac{1-\mathrm{\&epsiv;}}{\mathrm{\&epsiv;}}+A_3\frac{\mathrm{\&Delta;\&epsiv;}}{{\mathrm{\&epsiv;}}^2})\frac{T_{11}+T_{12}}{2}+(B_1+B_2\frac{1-\mathrm{\&epsiv;}}{\mathrm{\&epsiv;}}+B_3\frac{\mathrm{\&Delta;\&epsiv;}}{{\mathrm{\&epsiv;}}^2})\frac{T_{11}-T_{12}}{2}---\left(23\right)$

式中，T_s为地表温度；ε＝(ε₁₁+ε₁₂)/2为11μm和12μm通道的地表比辐射率平均值；Δε＝ε₁₁-ε₁₂为11μm和12μm通道的地表比辐射率差值；ε₁₁和ε₁₂分别为11μm和12μm通道的地表比辐射率；T₁₁和T₁₂分别为11μm和12μm通道的亮度温度；A1、A2、A3、B1、B2、B3和C为拟合系数。

(5)拟合系数的确定

从公式(23)可以看出，要想反演地表温度，必须要先确定拟合系数A1、A2、A3、B1、B2、B3和C。拟合系数的确定分为两步：

1)利用MODTRAN模拟得到模拟数据(如图6所示)；

2)利用模拟数据确定拟合系数。

利用MODTRAN模拟数据时，大气廓线、大气水汽含量、观测角度、地表温度和比辐射率的设置如下：

大气廓线：选择6条MODTRAN标准大气廓线和12条从TIGR大气廓线数据库中挑选的大气廓线。这12条TIGR大气廓线是参照MODTRAN的6条标准大气廓线的大气水汽含量和大气底层温度挑选的，每2条TIGR大气廓线对应一条MODTRAN标准大气廓线。

大气水汽含量：为了使大气廓线更具有代表性，对每条大气廓线的大气水汽含量进行倍数变化，倍数从0.5倍到1.5倍之间变化，以0.1倍为间隔。

观测角度：给定6个观测角度0.0°、33.6°、44.4°、51.3°、56.3°、60.0°。

地表温度：地表温度的设置与大气廓线底层温度Ta有关，低温大气状况下(Ta≤280K)，地表温度从Ta-10K到Ta+15K之间变化，以5K为间隔；高温大气状况下(Ta＞280K)，地表温度从Ta-5K到Ta+20K之间变化，以5K为间隔。

地表比辐射率：两个通道的地表比辐射率通过地表比辐射率平均值ε和地表比辐射率差值Δε计算得到，ε₁₁＝(2ε+Δε)/2，ε₁₂＝(2ε-Δε)/2。其中，地表比辐射率平均值从0.9到1.0之间变化，间隔为0.02；地表比辐射率差值从-0.025到0.015之间变化，间隔为0.005。

根据无人机两个热红外通道的通道响应函数，利用MODTRAN模拟6个给定的观测角度下两个通道的大气上下行辐射和透过率，结合地表温度和两个通道的地表比辐射率，根据大气辐射传输方程计算得到两个通道的辐射亮度，再利用普朗克函数将两个通道的辐射亮度转换为亮度温度。

利用模拟数据确定拟合系数时，为了提高地表温度的反演精度，在每个给定的观测角度下，对地表温度、大气水汽含量和地表比辐射率平均值3个参数进行子区间的划分，再利用每个子区间内和每个观测角度下的模拟数据确定拟合系数。使用每个子区间内和每个观测角度下的拟合系数使得地表温度的反演更加准确。考虑到子区间的连续，在进行子区间划分时，对每个参数的子区间设置一定的重叠区。地表温度划分5个子区间，每个子区间有5K的重叠区：Ts≤280K，275K≤Ts≤295K，290K≤Ts≤310K，305K≤Ts≤325K，Ts≥320K。大气水汽含量划分6个子区间，每个子区间有0.5g cm^-2的重叠区：0g cm^-2≤CWV≤1.5g cm^-2，1.0g cm^-2≤CWV≤2.5g cm^-2，2.0g cm^-2≤CWV≤3.5g cm^-2，3.0g cm^-2≤CWV≤4.5gcm^-2，4.0g cm^-2≤CWV≤5.5g cm^-2，5.0g cm^-2≤CWV≤6.5g cm^-2。地表比辐射率平均值划分2个区间，每个子区间有0.02的重叠区：0.9≤ε≤0.96，0.94≤ε≤1.0。

在每个特定的参数子区间内和每个给定的观测角度下，利用地表温度、比辐射率和相应的亮度温度，通过多元回归拟合确定拟合系数。对于每个子区间内，其它观测角度下的拟合系数可以通过给定的观测角度的拟合系数按照1/cos(θ_v)线性插值得到。

地表温度反演流程如图7所示，其步骤为：

1)根据普朗克函数将辐射亮度转换为亮度温度。

2)根据地表反射率，通过公式(20)计算归一化植被指数NDVI，通过公式(21)计算植被覆盖度P_v。

3)根据NDVI和P_v，通过公式(22)计算地表比辐射率。

4)根据地表比辐射率、大气水汽含量和亮度温度，利用拟合系数A₁、A₂、A₃、B₁、B₂、B₃和C，通过公式(23)计算地表温度。

Claims (6)

1.一种无人机光学遥感影像数据处理方法，其步骤包括：

A)无人机光学遥感影像数据预处理；

B)对无人机光学遥感影像数据进行大气水汽反演，得到大气水汽含量产品；

C)利用步骤B的大气水汽含量产品，对无人机光学遥感影像数据进行地表反射率反演，得到地表反射率产品；

D)利用步骤B的大气水汽含量产品和步骤C的地表反射率产品，对无人机光学遥感影像数据进行地表温度反演，得到地表温度产品。

2.根据权利要求1所述的无人机光学遥感影像数据处理方法，其中，对无人机光学遥感影像数据进行预处理过程是：

(1)根据影像数据本身，对无人机光学遥感影像数据进行相对辐射定标，消除由于各个探测元的响应差异在影像上产生的条纹；

(2)根据无人机传感器的辐射定标系数，对无人机光学遥感影像数据进行绝对辐射定标，将影像灰度值转换为辐射亮度。

3.根据权利要求1所述的无人机光学遥感影像数据处理方法，其中，对无人机光学遥感影像数据进行大气水汽反演过程是：

(1)根据水汽通道和相邻两个大气窗口通道的辐射亮度计算比值指数；

(2)根据比值指数和拟合系数计算大气斜程水汽；

(3)根据大气斜程水汽和大气水汽比计算大气垂直水汽。

4.根据权利要求1所述的无人机光学遥感影像数据处理方法，其中，对无人机光学遥感影像数据进行地表反射率反演过程是：

(1)根据气溶胶光学厚度、大气水汽含量、飞行高度、地表高程、太阳天顶角、观测天顶角和相对方位角参数及大气参数查找表，计算得到朗伯均一地表反射率；

(2)根据朗伯均一地表反射率计算得到周围环境平均地表反射率；

(3)根据朗伯均一地表反射率和周围环境平均反射率逐像元计算得到朗伯非均一地表反射率；

(4)根据步骤3得到的朗伯非均一地表反射率再计算周围环境平均地表反射率，重复步骤2和步骤3进行迭代，直到两次计算的朗伯非均一地表反射率之差小于设定的阈值，迭代结束。

5.根据权利要求1所述的无人机光学遥感影像数据处理方法，其中，对无人机光学遥感影像数据进行地表温度反演过程是：

(1)根据普朗克函数将辐射亮度转换为亮度温度；

(2)根据地表反射率计算归一化植被指数，并计算植被覆盖度；

(3)根据归一化植被指数和植被覆盖度计算地表比辐射率；

(4)根据地表比辐射率、大气水汽含量和亮度温度计算地表温度。

6.一种实现权利要求1所述方法的装置，其包括：

数据预处理模块进行相对辐射定标和绝对辐射定标；以及

大气水汽反演模块计算大气水汽含量；以及

地表反射率反演模块计算地表反射率；以及

地表温度反演模块计算地表温度。

Images