技术领域
[0001] 本
发明涉及定量遥感领域，尤其涉及一种基于大气尺度因子的地表温度估算方法。
背景技术
[0002] 在热红外遥感地表温度反演领域，大气校正和
辐射传输方程的高度耦合特性一直是困扰研究学者的核心难题。近几十年来，研究学者们基于不同的出发点，陆续提出了众多沿用至今的经典
算法，如分裂窗算法，TES算法，
光谱平滑算法，日夜算法，NDVI算法等。然而，现有算法中，或多或少的存在一定的局限性。如被广泛用于多光谱热红外地表温度的通用分裂窗算法(SW)，其算法由于不需要借助大气廓线信息，因此反演过程简单，且在满足算法假设与一定观测条件下能够得到较高的反演
精度。但正如所述的一样，SW算法最大的局限性往往也在于其算法自身的诸多假设，如需提供精确的通道地表发射率，需满足观测辐射亮温、大气等效亮温与地表亮温近似相等，以便更好的满足普朗克函数在不同温度下的泰勒级数展开；另一方面，由于SW算法是基于观测通道亮温的线性或非线性函数，因此其拟合数据的选择也直接影响最后的反演精度，这也是导致其在大
水汽条件/高观测天顶
角下反演精度下降的部分原因。综合而言，SW算法在缺乏精准大气廓线的前提下，是一种有效的地表温度反演方法，目前也是MODIS地表产品的官方反演算法。虽然，针对SW中普朗克函数在温度展开下的合理线性问题，Wang等人提出了一种基于普朗克函数的泰勒展开，其研究结果也表明反演精度可以实现一定程度的提升，但在仍对发射率的使用较为敏感。
[0003] 不同于SW利用相邻热红外通道消除大气辐射的原理方法，Gillespie等人通过将3种不同的经典算法进行融合的方式，提出了基于ASTER卫星热红外影像的地表温度与发射率分离(TES)算法。然而，TES算法的局限性在于，其不仅需要首先进行精准的大气校正，且算法自身仍依赖经验关系，如其算法核心算法MMD模
块，就是建立在MMD与最小发射率之间的经验拟合关系等。另一方面，Pivovarník等人在其研究中指出，TES的初值估算方法(NEM模块)导致其反演的地表温度与发射率之间不具有耦合性，表现为独立误差，且其算法自身并不适用于低发射率地表。针对这一问题，Pivovarník等人基于最大最小通道发射率的经验关系，提出了一种改进的TES算法，虽然其改进后的算法具有一定的地表参数耦合特性，且可应用于低发射率地表，但TES算法的本质问题仍没有得到有效解决，即对经验关系和大气校正的依赖。其他诸如NDVI算法，日/夜算法，单通道算法等，虽然都有其自身的算法优势，但其最大的局限仍是对经验关系，先验知识或大气校正的依赖，导致其算法的
稳定性和应用性都受到限制。
[0004] 如果说多光谱受限于其有效的观测信息，导致其不得不采用一种折衷的反演策略，那么随着高光谱数据的出现，为研究学者提供了更多的可能。近年来，越来越多的高光谱反演算法陆续被提出，由于大部分的算法都是建立在合理数学假设的
基础上，因此其在理想条件下均可实现高精度的地表参数反演，并且相关研究也表明，其算法普遍具有较好的抗白噪声特性，然而，与绝大部分的多光谱算法相同，现有高光谱算法仍主要依赖于大气校正水平，其对大气廓线误差，尤其是水汽含量误差极为敏感，导致目前高光谱算法较难应用在产品生产之中。此外，针对高光谱遥感数据，科学团队也提出了一体化的反演方法，但目前该类算法被更多的用于大气廓线的相关研究，并不适合地表参数的反演。
[0005] 针对前面所说的大气校正方法，目前主要存在两种主要的技术手段，一种是直接利用大气廓线产品，结合辐射传输方程，直接进行大气参数的反演，显然，此类方法最大的问题在于其反演精度直接取决于大气廓线产品精度，然而，现有的大气廓线产品显然很难满足这一要求，极大的限制了此类方法的应用。区别与直接的大气廓线产品，科研学者利用图像的多像元特征，相继提出了AAC、ISAC、WVS等经典算法。其中，ISAC与AAC相类似，均是从图像中选取出特征像元(近
黑体或灰体)，并通过一定的假设与约束，实现大气上行辐射与透过率的近似拟合，其算法并不直接反演大气下行辐射，而是通过大气参数间的相互经验关系进行间接反演。此类方法最大的局限性在于，其对大气下行的忽略在高发射率地表下的影响可近似忽略，但对于低发射率地表而言，此类假设的不合理性将十分突出，受限于此，ISAC与AAC的应用目前也较为有限。
[0006] 与上述方法均不同的是，Tonooka等人提出了一种基于既有大气廓线产品的大气水汽修正方法，因在热红外遥感领域，水汽是影响地表温度反演的主要因素之一，因此有效地对大气廓线进行水汽修正可显著提升反演效果。WVS理论上来讲是一种适用于多光谱和高光谱的大气廓线修正方法，目前也是ASTER和MODIS地表温度产品反演的辅助生产手段。然而，WVS算法仍没有逃离对经验关系的依赖，如其算法首先就是根据通道观测亮温，水汽含量以及经验系数实现离地辐射的反演。文中也指出其算法的实现需要首先提取影像范围内一定数量的类灰体像元，理想情况下可实现约0.6K的反演精度。显然，WVS算法不可避免的会存在误差传递问题，反演精度的稳定性无法保证。尽管如此，WVS方法仍为研究学者提出了一种针对热红外遥感反演领域中大气廓线修正的新思路，即针对某特定敏感系数的修正仍可显著提升反演效果。
发明内容
[0007] 本文发明的目的是提供一种基于大气尺度因子的地表温度估算方法。
[0008] 具体地，本发明提供一种基于大气尺度因子的地表温度估算方法，包括如下步骤；
[0009] 步骤100，分别基于普朗克定律和辐射传输方程实现双通道辐射
能量间的函数表示以及双通道地表温度的表达式，然后将两式联立方程，对所有观测通道的地表发射率取光谱的均值处理后，得到修正大气廓线误差的代价函数；
[0010] 步骤200，获取卫星观测辐射数据以及大气廓线产品；
[0011] 步骤300，从现有发射率产品中获取目标像元的平均地表发射率；
[0012] 步骤400，
选定尺度因子代入代价函数进行大气廓线的
迭代修正，并认为代价函数取最小值时所对应的地表温度最接近真值，确定大气廓线的尺度修正与地表温度的估算结果。
[0013] 在本发明的一个实施方式中，所述步骤100中代价函数的具体得到过程如下：
[0014] 根据普朗克公式
[0015]
[0016] 得到地表-大气间的辐射传输方程
[0017]
[0018] 其中，i表示观测通道；Bi(T)为黑体温度T下的观测通道i处的热辐射能量；λi为观测通道i的
波长，单位：μm；c1和c2为普朗克热常数，取值分别为：1.19104×108W·μm·m-2·sr-1和1.43877×104μm·K；Ts为地表温度；Li表示通道i处的观测辐射；Bi(Ts)表示通道i在地表温度Ts下的普朗克函数；BiTa↓和BiTa↑分别表示通道i在大气等效下行辐射亮温Ta↓以及大气等效上行辐射亮温Ta↑下的普朗克函数；τi为通道i处的大气透过率；
[0019] 根据式(1)，对于给定的温度T，Bi(T)可表示为通道波长λ和温度T的函数，因此，当观测通道i和j具有相同的温度T时，得到下式：
[0020]
[0021] 其中，i和j分别表示第i和第j个观测通道；Bi和Bj表示温度T下的普朗克函数；
[0022] 当Bi作为变量使用时，隐去了其函数中的原变量T，对于相同的温度T，通过化简式(3)，并将Bj表示为Bi的函数，得到
[0023]
[0024] 由式(2)可知，位于通道j处的观测辐射表示如下：
[0025]
[0026] 化简后为
[0027]
[0028] 然后联立式(4)和式(6)化简得到下式：
[0029]
[0030] 根据式(2)将Bi表示如下：
[0031]
[0032] 其中，式(2)计算得到的Bi记为Bi,RTE，式(3)得到的Bi记为Bi,PLK；
[0033] 将所有通道的地表发射率取光谱的均值处理的公式如下
[0034]
[0035] 其中，i表表示第i个观测通道；n代表观测通道数；εm表示平均发射率；
[0036] 由于地表温度的误差成分应包含大气廓线和地表发射率取均值误差两部分，即：
[0037]
[0038]
[0039] 其中，i为第i个观测通道；Tis,RTE和Tis,PLK分别表示基于Bi,RTE和Bi,PLK计算得到的地表温度Ts；B-1表示普朗克逆函数；Δε2和Δatm2分别表示地表发射率误差和大气廓线误差成分；
[0040] 结合式(7)-(11)，即得到代价函数公式：
[0041]
[0042] 其中，c为代价函数；m表示共m次迭代；k表示第k次迭代。
[0043] 在本发明的一个实施方式中，所述通道i，j为任意选定的通道
位置。
[0044] 在本发明的一个实施方式中，所述通道j为所述通道i的相邻通道。
[0045] 在本发明的一个实施方式中，将所述大气廓线以及地表发射率误差以综合方式体现为反演通道的地表温度误差，然后利用不同尺度因子下的地表温度差值统计特征，找出最优参数组合，即可实现大气廓线的误差修正与地表温度的估算。
[0046] 在本发明的一个实施方式中，针对多光谱和高光谱通道可提供的的多通道对组合，通过所述代价函数建立多组通道对的代价函数，可消除仪器观测产生的白噪声，进而提升反演精度。
[0047] 在本发明的一个实施方式中，所述步骤400中，选定尺度因子代入代价函数进行大气廓线迭代修正的过程如下：
[0048] 步骤401，首先选定尺度因子S的变化范围，默认选择0.5～1.5之间，以0.1为步长进行迭代；
[0049] 步骤402，利用辐射传输模型，计算在尺度因子S(k)作用下的大气参数，记为{B(Ta,loop↓),B(Ta,loop↑),τloop}k,，其中k表示第k次迭代；
[0050] 步骤403，基于观测通道设置，构建n组通道对(i,j)组合，并利用所述代价函数计算得到此时的代价函数c(k)m,loop，其中m表示共进行m此迭代；
[0051] 步骤404，重复步骤401至403，直至迭代全部完成。
[0052] 在本发明的一个实施方式中，所述直至迭代全部完成是指：以cm,loop的最小值时所对应的尺度因子S(k)c_min作为与地表温度(Ts)c_min最接近的真实值，完成大气廓线修正与地表温度的反演。
[0053] 在本发明的一个实施方式中，所述步骤401中，设可获取的大气廓线的水汽含量误差在±50％的区间内。
[0054] 本发明建立了一种基于双(多)通道地表-大气-卫星观测参数间的约束关系，并通过迭代修正的方式，实现多种观测条件下的大气水汽修正与地表温度的估算；由于算法的构建基于普朗克函数和辐射传输方程，不依赖任何经验关系，因此具有更高的反演稳定性和普适性；本发明不但适用极干-极湿大气环境，且能够在多观测角度下保证精度；此外，本发明还可适用于多光谱和高光谱
传感器，且修正后的大气廓线，可结合现有算法实现地表参数的二次反演，具有更高的应用价值。
附图说明
[0055] 图1是本发明一个实施方式的地表温度估算方法的步骤示意图；
[0056] 图2是本发明一个实施方式的地表温度估算方法的流程示意图。
具体实施方式
[0057] 本发明针对地表温度反演中普遍存在的技术问题，从普朗克黑体定律和辐射传输方程自身出发，提出了一种基于大气廓线产品的双(多)通道约束大气修正与地表温度估算方法，即通过迭代修正的方式，实现多种观测条件下的高精度大气水汽修正与地表温度的初步估算。本算法不依赖任何经验关系，具有更高的反演稳定性和普适性；本算法可适用极干/极湿大气环境，且能够在多观测角度下保证精度；同时，算法还可用于人造地表发射率的反演，并适用于多类型传感器，显著优于现有算法，具有更高的应用价值。
[0058] 如图1、2所示，在一个实施方式中，本发明公开一种基于大气尺度因子的地表温度估算方法，包括如下步骤；
[0059] 步骤100，分别基于普朗克定律和辐射传输方程实现双通道辐射能量间的函数表示以及双通道地表温度的表达式，然后将两式联立方程，对所有观测通道的地表发射率取光谱的均值处理后，得到修正大气廓线误差的代价函数；
[0060] 在物理学中，普朗克黑体辐射定律，也称为普朗克定律或黑体辐射定律，是指在任意温度T下，从一个黑体中发射出的
电磁辐射的辐射率与
频率彼此之间的关系，普朗克公式可表达如下：
[0061]
[0062] 其中，i表示观测通道；Bi(T)为黑体温度T下的观测通道i处的热辐射能量；λi为观测通道i的波长，单位：μm；c1和c2为普朗克热常数，取值分别为：1.19104×108W·μm·m-2·sr-1和1.43877×104μm·K。
[0063] 在热红外定量遥感领域，普朗克函数作为定量反演的基础定律，被用于描述地表自身的辐射特性，在此基础上，可得到地表-大气间的辐射传输方程如下：
[0064]
[0065] 其中，i表示观测通道；Ts为地表温度；Li表示通道i处的观测辐射；Bi(Ts)表示通道i在地表温度Ts下的普朗克函数；BiTa↓和BiTa↑分别表示通道i在大气等效下行辐射亮温Ta↓以及大气等效上行辐射亮温Ta↑下的普朗克函数；τi为通道i处的大气透过率。
[0066] 根据式(1)，对于给定的温度T，Bi(T)可表示为通道波长λ和温度T的函数，因此，当观测通道i和j具有相同的温度T时，得到下式：
[0067]
[0068] 其中，i和j分别表示第i和第j个观测通道；Bi和Bj表示温度T下的普朗克函数。本发明中,为了便于表示，当Bi作为变量使用时，隐去了其函数中的原变量T，以更好的
定位其在方程中的角色。
[0069] 对于相同的温度T，通过化简式(3)，并将Bj表示为Bi的函数，得到
[0070]
[0071] 由式(2)可知，位于通道j处的观测辐射表示如下：
[0072]
[0073] 化简后为
[0074]
[0075] 在式(3)-(6)中，关于通道j的选择原理上并不会显著影响算法精度，默认情况下，我们采用的均为相邻的通道对，即对于任意选定的通道位置i，通道j默认选择其相邻的通道。
[0076] 然后联立式(4)和式(6)化简得到下式：
[0077]
[0078] 根据式(2)将Bi表示如下：
[0079]
[0080] 其中，为了便于区分式(7)和式(8)所得到的Bi，将基于辐射传输方程，即公式(2)，计算得到的Bi记为Bi,RTE，而将基于普朗克函数，即公式(3)计算得到的Bi记为Bi,PLK。
[0081] 由式(7)和式(8)可知，Bi,RTE与Bi,PLK是分别基于不同通道参数(i和j)实现的地表发射辐射能量的计算，彼此独立，且当式(7)和式(8)中的大气参数以及地表发射率可以准确获取时，式(7)与式(8)的辐射能量应近似相等。
[0082] 由公式(2)可知，即便在可实现精准大气校正的前提下，地表温度与地表发射率彼此仍是显著耦合的，也就是说，在不考虑地表参数反演算法自身误差的基础上，当大气廓线存在误差时，其对反演结果的误差成分将会不同程度的同时以地表温度与发射率误差形式存在，正如ISSTES算法所基于的光谱平滑原理。
[0083] 将所有通道的地表发射率做取光谱的均值处理，即
[0084]
[0085] 其中，i表表示第i个观测通道；n代表观测通道数；εm表示平均发射率；
[0086] 采用式(9)的发射率处理方法的优点在于，由于将地表发射率作为常数，大气廓线误差将直接呈现在地表温度之上，即此时得到的地表温度不再是唯一单值，而实以误差的光谱形式存在，需要指出的是，针对低发射率地表或者光谱
波动明显的地表，此类处理也会产生附加误差，因此针对此类地表仍需提供完整通道发射率，但是，在基于卫星影像的反演生产而言，大部分地表均为自然地表，因此此类处理的影响有限，且仍可进一步地通过影像插值平滑等处理方法，实现对个别异常值的剔除，进一步降低均值化发射率的影响。
[0087] 在采用均值发射率后，此时地表温度的误差成分应包含大气廓线和地表发射率取均值误差两部分，即：
[0088]
[0089] 其中，
[0090] 其中，i为第i个观测通道；Tis,RTE和Tis,PLK分别表示基于Bi,RTE和Bi,PLK计算得到的地表温度Ts；B-1表示普朗克逆函数；Δε2和Δatm2分别表示地表发射率误差和大气廓线误差成分；
[0091] 需要注意的是，采用平均发射率的处理方法在具有较平滑光谱的下垫面上的反演精度较为稳定，而对于人造地表，其低发射率的特点导致大气误差的比重上升，为了保证精度，仍需使用其发射率光谱的近似估算。
[0092] 结合式(7)-(11)，即得到代价函数公式：
[0093]
[0094] 其中，i为第i个观测通道；n为观测通道数；c为代价函数；m表示共m次迭代；k表示第k次迭代。
[0095] 虽然，地表温度与地表发射率的耦合一方面为地表参数的反演提供了契机，但是在大气廓线无法准确提供时，误差的耦合方式反而对大气廓线的校正增加了难度。因此，在本方法中，第一步的目标首先是完成对大气廓线的修正，然后进一步的实现地表温度的反演，最终才能更好的建立大气廓线误差的评价参数。
[0096] 通过给定初值的方式，可不断迭代进行大气参数的修正，并计算公式(12)，认为当c取最小值cmin时，其所对应的地表温度和大气参数最为准确，即可完成修正与反演过程。
[0097] 需要注意的是，本方法考虑到大气廓线中的水汽含量是影响热红外定量遥感反演的主要参数，且相关研究也表明有效的水汽修正可显著提升反演精度，因此，虽然本方法虽然只给出对大气水汽含量的修正说明，但算法理论上仍适用于温度廓线的误差修正。此外，根据公式(1)-(12)，本代价函数还可适用于多光谱和高光谱通道特征，尤其是在多通道下，可通过式(12)建立多组通道对的代价函数以进一步消除仪器观测白噪声，从而提升反演精度。
[0098] 步骤200，获取卫星观测辐射数据以及大气廓线；
[0099] 步骤300，从现有发射率产品中获取目标像元的平均地表发射率；
[0100] 这里的现有发射率产品也可以通过NDVI等算法直接计算得到。
[0101] 步骤400，选定尺度因子代入代价函数进行大气廓线的迭代修正，并认为代价函数取最小值时所对应的地表温度最接近真值，确定大气廓线的尺度修正与地表温度的估算结果。
[0102] 选定尺度因子代入代价函数进行大气迭代修正的过程如下：
[0103] 步骤401，首先选定尺度因子S的变化范围，默认选择0.5～1.5之间，以0.1为步长进行迭代；
[0104] 即认为该步骤中获取的大气廓线的水汽含量误差在±50％的区间内。
[0105] 步骤402，利用辐射传输模型，计算在尺度因子S(k)作用下的大气参数，记为{B(Ta,loop↓),B(Ta,loop↑),τloop}k,，其中k表示第k次迭代；
[0106] 步骤403，基于观测通道设置，构建n组通道对(i,j)组合，并利用所述代价函数计算得到此时的代价函数c(k)m,loop，其中，m表示共进行m此迭代；
[0107] 步骤404，重复步骤401至403，直至迭代全部完成。
[0108] 其中，迭代全部完成是指：以cm,loop的最小值所对应的尺度因子S(k)c_min与地表温度(Ts)c_min最为接近真实值，完成大气廓线修正与地表温度的反演。
[0109] 本发明建立了一种基于双(多)通道地表-大气-卫星观测参数间的约束关系，并通过迭代修正的方式，实现多种观测条件下的大气水汽修正与地表温度的估算；由于算法的构建基于普朗克函数和辐射传输方程，不依赖任何经验关系，因此具有更高的反演稳定性和普适性；本发明不但适用极干-极湿大气环境，且能够在多观测角度下保证精度；此外，本发明还可适用于多光谱和高光谱传感器，且修正后的大气廓线，可结合现有算法实现地表参数的二次反演，具有更高的应用价值。
[0110] 本发明能够克服现有大气廓线产品中水汽含量普遍高估的问题且算法具有较高的抗噪特性，即可在同时存在仪器噪声和大气偏差的情况下，仍能够有效完成水汽修正与地表温度的反演。本发明可与现有TES算法结合，构建二步反演策略，具有显著的普适性。
[0111] 至此，本领域技术人员应认识到，虽然本文已详尽示出和描述了本发明的多个示例性
实施例，但是，在不脱离本发明精神和范围的情况下，仍可根据本发明公开的内容直接确定或推导出符合本发明原理的许多其他变型或
修改。因此，本发明的范围应被理解和认定为
覆盖了所有这些其他变型或修改。
