一种MSG2-SEVIRI数据估算地表温度日较差的方法

摘要

本发明公开了一种MSG2-SEVIRI数据估算地表温度日较差的方法，包括以下步骤：获取研究区域的MSG2-SEVIRI数据，对MSG2-SEVIRI数据进行辐射定标，将DN值转换为辐亮度；根据普朗克方程，将辐亮度转换为星上亮温；利用地表温度分裂窗算法计算任意两个时刻的地表温差；耦合地表温度日变化模型和分裂窗算法，通过最小二乘拟合法，计算地表温度日变化模型参数；根据步骤3获取的地表温度日变化模型参数计算地表温度日较差。本发明的目的是充分利用静止卫星数据高时间分辨率的特点，通过耦合分裂窗算法和地表温度日变化模型，在不需要反演地表温度和比辐射率的情况下，从静止卫星数据直接估算地表温度日较差。

说明书

技术领域

本发明属于遥感技术领域，具体涉及一种MSG2-SEVIRI数据估算地表温度日较差的方法。

背景技术

地表温度日较差，也称为地表温度日振幅，是一天中地表温度最高值与最低值之差。其大小和纬度、季节、地表性质及天气情况有关。地表温度日较差在气候变化和农情监测等研究方面具有重要意义。传统的地表温度日较差测量是通过地面观测获取。常规的地面测量只能代表点尺度的数据，难以获取区域尺度上的数据，而卫星遥感是获取区域尺度上地表温度日较差的唯一手段。要获取地表温度日较差，首先需要从卫星数据中反演得到地表温度。但是，受大气和地表比辐射率的影响，从卫星数据中准确地反演地表温度目前仍然比较困难，因此从卫星数据获取区域尺度上的地表温度日较差也比较困难。静止卫星具有较高的时间分辨率，例如MSG-SEVIRI每15分钟获取一次数据，GOES每30分钟获取一次数据，因此静止卫星数据可以描述地表温度日变化。

发明内容

本发明的目的是提供一种MSG2-SEVIRI数据估算地表温度日较差的方法，解决了现有技术中存在的需要已知地表温度和比辐射率估算地表温度日较差的问题。

本发明所采用的技术方案是，一种MSG2-SEVIRI数据估算地表温度日较差的方法，包括以下步骤：

步骤1、获取研究区域的MSG2-SEVIRI数据，对MSG2-SEVIRI数据进行辐射定标，将DN值转换为辐亮度；根据普朗克方程，将辐亮度转换为星上亮温；

步骤2、利用地表温度分裂窗算法计算任意两个时刻的地表温差；

步骤3、耦合地表温度日变化模型和分裂窗算法，通过最小二乘拟合法，计算地表温度日变化模型参数；

步骤4、根据步骤3获取的地表温度日变化模型参数计算地表温度日较差。

进一步地，步骤1中获取研究区域的MSG2-SEVIRI数据，对MSG2-SEVIRI数据进行辐射定标，将DN值转换为辐亮度；根据普朗克方程，将辐亮度转换为星上亮温具体为：

对MSG2-SEVIRI数据进行辐射定标，将MSG2-SEVIRI获取的DN值转换为辐亮度；

L＝CAL_gain*DN+CAL_offset(1)

式中，L为MSG2-SEVIRI测量的星上辐亮度，DN为计数值，CAL_gain和CAL_offset为辐射定标系数，对于MSG2-SEVIRI第9通道：CAL_gain＝0.205034，CAL_offset＝-10.4568；对于第10通道：CAL_gain＝0.222311，CAL_offset＝-11.3379；

根据普朗克方程，进一步将辐亮度转换为星上亮温：

$>T_i=[C_2{v}_c/ln(\frac{C_1{v}_c^3}{L}+1)-D]/E---\left(2\right);$>

式中，T_i为星上亮温，C₁和C₂为辐射常量，C₁＝1.19104e5mWm^-2sr^-1(cm^-1)^-4,C₂＝1.43877K(cm^-1)^-1；v_c为中心波数，D和E为转换系数；对于MSG2-SEVIRI第9通道：v_c＝931.700cm^-1,D＝0.640K,E＝0.9983；对于第10通道：v_c＝836.445cm^-1，D＝0.408K,E＝0.9988。

进一步地，步骤2中地表温度分裂窗算法计算任意两个时刻的地表温差具体为：

地表温度分裂窗算法表示为：

T_s＝T₉+a(T₉-T₁₀)+b(T₉-T₁₀)²+c+B(ε,W)(3)

式中，T_s为地表温度；T₉和T₁₀分别为MSG2-SEVIRI第9和10通道的星上亮温；a、b和c为拟合系数，只与大气状况有关，而与地表条件无关；系数B为地表比辐射率和大气水汽含量的函数，表示为：

B(ε,W)＝(α₀+α₁W)(1-ε)-(β₀+β₁W+β₂W²)Δε(4)

其中，ε＝(ε₉+ε₁₀)/2为MSG2-SEVIRI第9和10通道比辐射率的平均值；Δε＝ε₉-ε₁₀为MSG2-SEVIRI第9和10通道比辐射率的差值；W为大气水汽含量；α₀、α₁、β₀、β₁和β₂为拟合系数。

一天中任意两个时刻的大气水汽含量差值导致系数B的变化表示为：

$>\Delta B=[{\alpha }_1(1-ϵ)-{\beta }_1\Delta ϵ](W_i-W_j)-{\beta }_2\Delta ϵ(W_i^2-W_j^2)---\left(5\right)$>

式中，W_i和W_j分别为i和j时刻的大气水汽含量。

为了通过公式(5)计算大气水汽含量的变化导致系数B的变化，取一天内任意两个时刻的大气水汽含量之和为10gcm^-2，之差为1gcm^-2，同时，给定2个像元的地表比辐射率：沙地中，ε＝0.96和Δε＝-0.02，植被，ε＝0.99和Δε＝0)。

根据上述分析，系数B在一天内可以假定为常数。由此，任意时刻t的地表温度可以表示为：

T_s(t)＝T₉(t)+a[T₉(t)-T₁₀(t)]+b[T₉(t)-T₁₀(t)]²+c+B(6)

根据公式(6)，任意时刻t和参考时刻t_r(t_r设定为13:00UTC)之间的地表温度差值表示为：

ΔT_s(t)＝T_s(t)-T_s(t_r)

＝T₉(t)+a[T₉(t)-T₁₀(t)]+b[T₉(t)-T₁₀(t)]²-(7)

{T₉(t_r)+a[T₉(t_r)-T₁₀(t_r)]+b[T₉(t_r)-T₁₀(t_r)]²}

在公式(7)中，a和b是与地表条件无关的大气系数。由于MSG2-SEVIRI观测天顶角的范围较大，为了得到任意角度下的系数a和b，系数a和b要表示成观测天顶角的函数，即a(θ)＝a₀+a₁[sec(θ)-1]和b(θ)＝b₀+b₁[sec(θ)-1]。a₀、a₁、b₀和b₁的值分别为1.07、-0.04、0.259和0.165。当MSG2-SEVIRI第9和10通道星上亮温获得后，通过公式(7)可以计算得到任意时刻的地表温差。

进一步地，步骤3中耦合地表温度日变化模型和分裂窗算法，通过最小二乘拟合法，计算地表温度日变化模型参数具体为：

地表温度日变化模型的公式表达如下：

$>\begin{array}{c}{T}_s^d\left(t\right)=T_0+T_a\cos \left(\frac{\pi }{\omega }\right(t-t_m\left)\right),t<t_s\\ T_s^n\left(t\right)=T_0+\delta T+[T_a\cos \left(\frac{\pi }{\omega }\right(t_s-t_m\left)\right)-\delta T]\frac{k}{(k+t-t_s)},t\ge t_s\end{array}---\left(8\right)$>

其中，

$>\omega =\frac{4}{3}(t_m-t_{sr})---\left(9\right)$>

$>k=\frac{\omega }{\pi }[{\tan }^{-1}\left(\frac{\pi }{\omega }\right(t_s-t_m\left)\right)-\frac{\delta T}{T_a}{\sin }^{-1}\left(\frac{\pi }{\omega }\right(t_s-t_m\left)\right)]---\left(10\right)$>

式中，T_s^d和T_sⁿ分别为白天和晚上部分的地表温度；T₀为日出时刻附近的初始温度；T_a为温度的幅度；ω为余弦函数半周期宽；t_m为最大温度出现的时刻；t_s为温度开始衰减的时刻；δT为T₀和T(t→∞)之间的温差；k为衰减系数；该地表温度日变化模型具有5个自由参数，即T₀、T_a、δT、t_m、t_s；参数ω和k分别通过公式(9)和(10)计算得到；t_sr通过太阳几何计算得到；

从公式(7)看出，任意时刻t和参考时刻t_r之间的地表温度差值表示为：

ΔT_s^d(t)＝T_s^d(t)-T_s^d(t_r),t<t_s

(11)

ΔT_sⁿ(t)＝T_sⁿ(t)-T_s^d(t_r),t≥t_s

结合公式(7)和(8)，任意时刻t和参考时刻t_r之间的地表温差进一步表示为：

$>\begin{array}{c}{\mathrm{\Delta T}}_s^d\left(t\right)=T_a\cos \left(\frac{\pi }{\omega }\right(t-t_m\left)\right)-T_a\cos \left(\frac{\pi }{\omega }\right(t_r-t_m\left)\right),t<t_s\\ {\mathrm{\Delta T}}_s^n\left(t\right)=\delta T+[T_a\cos \left(\frac{\pi }{\omega }\right(t_s-t_m\left)\right)-\delta T]\frac{k}{(k+t-t_s)}\\ -T_a\cos \left(\frac{\pi }{\omega }\right(t_r-t_m\left)\right),t\ge t_s\end{array}---\left(12\right)$>

在公式(12)中，公式左边由公式(7)或者公式(11)从MSG-SEVIRI第9和10通道的星上亮温计算得到；在不需要反演地表温度和比辐射率情况下，通过最小二乘拟合，从公式(12)直接得到参数T_a、δT、t_m和t_s的值；在最小二乘拟合时，参数T_a、δT、t_m和t_s设定初值为：T_a＝ΔT_max-ΔT_min，其中，ΔT_max和ΔT_min分别为公式(7)计算的地表温差的最大值和最小值，δT＝0.5K，t_m＝12.5h，t_s＝17h。

进一步地，步骤4中根据步骤3获取的地表温度日变化模型参数计算地表温度日较差具体为：

从公式(8)看出，当t＝t_m时，地表温度达到最大值，即T_max＝T₀+T_a；当t＝∞时，地表温度达到最小值，即T_min≈T₀+δT，根据地表温度日较差DTR的定义，地表温度日较差表示为：

DTR＝T_max-T_min≈T_a-δT(13)

式中，T_a和δT从步骤3中可以获取得到。

本发明的有益效果是：本发明的目的是充分利用静止卫星数据高时间分辨率的特点，通过耦合分裂窗算法和地表温度日变化模型，在不需要反演地表温度和比辐射率的情况下，从静止卫星数据直接估算地表温度日较差。

本发明不需要对MSG2-SEVIRI数据进行大气校正。本发明利用热红外两个相邻通道对大气水汽吸收的差异实现大气校正，消除了MSG2-SEVIRI数据中的大气影响。

本发明不需要计算MSG2-SEVIRI地表比辐射率。本发明假定一天中地表比辐射率变化不大，通过计算任意两个时刻的地表温差，而不是计算每个时刻的地表温度，消除了地表比辐射率对地表温差的影响。

本发明不需要计算MSG2-SEVIRI地表温度。本发明利用分裂窗算法计算任意两个时刻的地表温差，利用地表温度日变化模型拟合地表温差日变化数据，计算得到地表温度日变化模型参数，从而进一步估算地表温度日较差，解决了现有方法中需要计算地表温度和比辐射率来估算地表温度日较差的问题。

附图说明

图1是MSG2-SEVIRI地表温度日较差计算流程图；

图2是MSG2-SEVIRI第9通道的星上亮温图(13:00UTC)；

图3是MSG2-SEVIRI第10通道的星上亮温图(13:00UTC)；

图4是地表温差的拟合及拟合误差图；

图5是估算的MSG2-SEVIRI地表温度日较差图。

具体实施方式

下面结合具体实施方式对本发明进行详细说明。

本发明提供一种MSG2-SEVIRI数据估算地表温度日较差的方法，如图1所示，包括以下步骤：

步骤1：获取研究区域的MSG2-SEVIRI数据，对MSG2-SEVIRI数据进行辐射定标，将DN值转换为辐亮度；根据普朗克方程，将辐亮度转换为星上亮温；

对MSG2-SEVIRI数据进行辐射定标，将MSG2-SEVIRI获取的DN值转换为辐亮度；

L＝CAL_gain*DN+CAL_offset(1)

式中，L为MSG2-SEVIRI测量的星上辐亮度，DN为计数值，CAL_gain和CAL_offset为辐射定标系数，对于MSG2-SEVIRI第9通道：CAL_gain＝0.205034，CAL_offset＝-10.4568；对于第10通道：CAL_gain＝0.222311，CAL_offset＝-11.3379。

根据普朗克方程，进一步将辐亮度转换为星上亮温：

$>T=[C_2{v}_c/ln(\frac{C_1{v}_c^3}{L}+1)-D]/E---\left(2\right)$>

式中，T_i为星上亮温，C₁和C₂为辐射常量，C₁＝1.19104e5mWm^-2sr^-1(cm^-1)^-4,C₂＝1.43877K(cm^-1)^-1；v_c为中心波数，D和E为转换系数；对于MSG2-SEVIRI第9通道：v_c＝931.700cm^-1,D＝0.640K,E＝0.9983；对于第10通道：v_c＝836.445cm^-1,D＝0.408K,E＝0.9988。

步骤2：利用地表温度分裂窗算法计算任意两个时刻的地表温差，

地表温度分裂窗算法表示为：

T_s＝T₉+a(T₉-T₁₀)+b(T₉-T₁₀)²+c+B(ε,W)(3)

式中，T_s为地表温度；T₉和T₁₀分别为MSG2-SEVIRI第9和10通道的星上亮温；a、b和c为拟合系数，只与大气状况有关，而与地表条件无关；系数B为地表比辐射率和大气水汽含量的函数，表示为：

B(ε_,W)＝(α₀+α₁W)(1-ε)-(β₀+β₁W+β₂W²)Δε(4)

其中，ε＝(ε₉+ε₁₀)/2为MSG2-SEVIRI第9和10通道比辐射率的平均值；Δε＝ε₉-ε₁₀为MSG2-SEVIRI第9和10通道比辐射率的差值；W为大气水汽含量；α₀、α₁、β₀、β₁和β₂为拟合系数。图2和图3分别显示了13:00UTC时刻的MSG2-SEVIRI第9和10通道的星上亮温图。

假定地表比辐射率在一天内变化不大(LiandBecker,1993)，一天中任意两个时刻的大气水汽含量差值导致系数B的变化可以表示为：

$>\Delta B=[{\alpha }_1(1-ϵ)-{\beta }_1\Delta ϵ](W_i-W_j)-{\beta }_2\Delta ϵ(W_i^2-W_j^2)---\left(5\right)$>

式中，W_i和W_j分别为i和j时刻的大气水汽含量。

为了通过公式(5)计算大气水汽含量的变化导致系数B的变化，假定一天内任意两个时刻的大气水汽含量之和为10gcm^-2，之差为1gcm^-2(SobrinoandElKharrza,1999)。同时，给定2个像元的地表比辐射率：沙地(ε＝0.96和Δε＝-0.02)和植被(ε＝0.99和Δε＝0)(SobrinoandElKharrza,1999)。

通过上述假定，利用公式(5)计算结果表明，对于沙地和植被像元，大气水汽含量的变化导致系数B的变化分别为-0.37K和-0.05K。因此，对于沙地和植被像元，系数B的变化导致地表温度的变化分别为-0.37K和-0.05K，这些误差几乎可以忽略。

根据上述分析，系数B在一天内可以假定为常数。由此，任意时刻t的地表温度可以表示为：

T_s(t)＝T₉(t)+a[T₉(t)-T₁₀(t)]+b[T₉(t)-T₁₀(t)]²+c+B(6)

根据公式(6)，任意时刻t和参考时刻t_r(t_r设定为13:00UTC)之间的地表温度差值表示为：

ΔT_s(t)＝T_s(t)-T_s(t_r)

＝T₉(t)+a[T₉(t)-T₁₀(t)]+b[T₉(t)-T₁₀(t)]²-(7)

{T₉(t_r)+a[T₉(t_r)-T₁₀(t_r)]+b[T₉(t_r)-T₁₀(t_r)]²}

在公式(7)中，a和b是与地表条件无关的大气系数。由于MSG2-SEVIRI观测天顶角的范围较大，为了得到任意角度下的系数a和b，系数a和b要表示成观测天顶角的函数，即a(θ)＝a₀+a₁[sec(θ)-1]和b(θ)＝b₀+b₁[sec(θ)-1]。a₀、a₁、b₀和b₁的值分别为1.07、-0.04、0.259和0.165。当MSG2-SEVIRI第9和10通道星上亮温获得后，通过公式(7)可以计算得到任意时刻的地表温差。

步骤3：耦合地表温度日变化模型和分裂窗算法，通过最小二乘拟合法，计算地表温度日变化模型参数；

Duanetal.(2013)建立了一个地表温度日变化模型，其可以用下述公式表达：

$>\begin{array}{c}{T}_s^d\left(t\right)=T_0+T_a\cos \left(\frac{\pi }{\omega }\right(t-t_m\left)\right),t<t_s\\ T_s^n\left(t\right)=T_0+\delta T+[T_a\cos \left(\frac{\pi }{\omega }\right(t_s-t_m\left)\right)-\delta T]\frac{k}{(k+t-t_s)},t\ge t_s\end{array}---\left(8\right)$>

其中，

$>\omega =\frac{4}{3}(t_m-t_{sr})---\left(9\right)$>

$>k=\frac{\omega }{\pi }[{\tan }^{-1}\left(\frac{\pi }{\omega }\right(t_s-t_m\left)\right)-\frac{\delta T}{T_a}{\sin }^{-1}\left(\frac{\pi }{\omega }\right(t_s-t_m\left)\right)]---\left(10\right)$>

式中，T_s^d和T_sⁿ分别为白天和晚上部分的地表温度；T₀为日出时刻附近的初始温度；T_a为温度的幅度；ω为余弦函数半周期宽；t_m为最大温度出现的时刻；t_s为温度开始衰减的时刻；δT为T₀和T(t→∞)之间的温差；k为衰减系数。该地表温度日变化模型具有5个自由参数，即T₀、T_a、δT、t_m、t_s。参数ω和k分别通过公式(9)和(10)计算得到。t_sr通过太阳几何计算得到。

从公式(7)可以看出，任意时刻t和参考时刻t_r之间的地表温度差值也可以表示为：

ΔT_s^d(t)＝T_s^d(t)-T_s^d(t_r),t<t_s

(11)

ΔT_sⁿ(t)＝T_sⁿ(t)-T_s^d(t_r),t≥t_s

结合公式(7)和(8)，任意时刻t和参考时刻t_r之间的地表温差可以进一步表示为：

$>\begin{array}{c}{\mathrm{\Delta T}}_s^d\left(t\right)=T_a\cos \left(\frac{\pi }{\omega }\right(t-t_m\left)\right)-T_a\cos \left(\frac{\pi }{\omega }\right(t_r-t_m\left)\right),t<t_s\\ {\mathrm{\Delta T}}_s^n\left(t\right)=\delta T+[T_a\cos \left(\frac{\pi }{\omega }\right(t_s-t_m\left)\right)-\delta T]\frac{k}{(k+t-t_s)}\\ -T_a\cos \left(\frac{\pi }{\omega }\right(t_r-t_m\left)\right),t\ge t_s\end{array}---\left(12\right)$>

在公式(12)中，公式左边可由公式(7)或者公式(11)从MSG-SEVIRI第9和10通道的星上亮温计算得到。在不需要反演地表温度和比辐射率情况下，通过最小二乘拟合，从公式(12)可以直接得到参数T_a、δT、t_m和t_s的值。在最小二乘拟合时，参数T_a、δT、t_m和t_s设定初值为：T_a＝ΔT_max-ΔT_min(ΔT_max和ΔT_min分别为公式(7)计算的地表温差的最大值和最小值)，δT＝0.5K，t_m＝12.5h，t_s＝17h。图4显示了MSG2-SEVIRI地表温差的拟合误差图。从图4可以看出，公式(12)很好地拟合了MSG2-SEVIRI地表温差数据，其均方根误差为0.32K。

步骤4：根据步骤3获取的地表温度日变化模型参数计算地表温度日较差。

从公式(8)可以看出，当t＝t_m时，地表温度达到最大值，即T_max＝T₀+T_a；当t＝∞时，地表温度达到最小值，即T_min≈T₀+δT。根据地表温度日较差DTR的定义，地表温度日较差表示为：

DTR＝T_max-T_min≈T_a-δT(13)

式中，T_a和δT从步骤3中可以获取得到。图5显示了利用本发明估算的MSG2-SEVIRI地表温度日较差图。