CN106909722B

CN106909722B - 一种近地面气温的大面积精准反演方法

Info

Publication number: CN106909722B
Application number: CN201710073515.8A
Authority: CN
Inventors: 秦江林; 杨秀好; 罗基同; 符合; 雷秀峰; 韦军; 秦元睿; 杨自鑫
Original assignee: Guangxi Zhuang Autonomous Region Harmful Control Quarantine Station; GUANGXI INSTITUTE OF METEOROLOGICAL DISASTER MITIGATION
Current assignee: Guangxi Zhuang Autonomous Region Harmful Control Quarantine Station; GUANGXI INSTITUTE OF METEOROLOGICAL DISASTER MITIGATION
Priority date: 2017-02-10
Filing date: 2017-02-10
Publication date: 2019-07-26
Anticipated expiration: 2037-02-10
Also published as: US10831949B2; US20190057171A1; WO2018145229A1; CN106909722A

Abstract

一种近地面气温的大面积精确反演方法，包括如下步骤：首先获取第1子模式modelVEC³至第f子模式modelVEC^f；构造栅格点原始数据记录集sBlkVEC并归一化处理；计算各个子模式的近地面气温反演全图；在近地面气温反演全图surfT^f中标记误差异常点及其邻近区域；利用选择性求和算法求出近地面气温精准反演影像图surfT；对surfT^f中像素点对所有f＝1,2,3,...都成立的区域的进一步处理。本发明通过引入一类超级非线性算法，解决了现有技术中存在的信息融合来反演近地气温问题，获得了反演误差与云量、地形、地势无关的满意结果。特别是融合了空间和地面信息，反演重建多云地区近地面气温，克服了云层的干扰而精确、可靠地获得大面积地区反演近地面气温影像图。

Description

一种近地面气温的大面积精准反演方法

技术领域

本发明涉及气象数据处理领域，特别是一种近地面气温的大面积精确反演方法。

背景技术

近地面气温，指离地面1.5-2.0米高的百叶箱气温[1-2]，是地表能量平衡中的一个重要参数，是地面与大气层相互作用过程中的一个重要物理量[3]，它可以反映土壤及地下可见的和潜在的能量变化[2-3]，是影响人类生产活动重要因素[4]。农作物、森林生长发育状况、病虫害(注：昆虫是温敏动物)的发展演变、低温冷冻灾害和高温危害无不与近地面温度有关[1-2]。另外，完备的近地面温度还是准确的预报各类气象要素的重要因子[5]；其它气象要素，例如，湿度、蒸发/蒸腾量、干旱、地表温度和湿度等等，都与近地面气温紧密相关。近地面气温可以在气象、农业[4]、林业、地理、地质等多个领域进行分析研究和应用[1-5]，特别是近几年来，利用热红外手段研究近地面气温受到大家广泛的重视。但是，利用以往的测量手段很难获取大面积的近地面气温数据，阻碍了近地面气温研究的发展。

MODIS和FY2/FY3影像的热红外波段可以获得大面积、动态变化的温度信息，为与近地面气温相关的研究提供了良好的数据基础。然而，由于地表热辐射的复杂性及大气等因素的影响，导致了MODIS和FY2/FY3影像所表征的温度与近地面气温之间还存在着一定的差异。而且卫星热红外波段表征的亮度温度并不能代表近地面的真实温度，它们之间的误差会给研究带来许多不确定的因素。因此，如何获取一种能够克服地面特征的复杂性和大气层结构的多变性，尤其广西，的近地气温的大面积反演方法，是许多学者多年来力图解决的难题。

传统温度反演多数是针对地表温度，本发明根据实际需要和我们能够获得的温度数据，针对的是近地面气温，概念上有差别，但是方法可以通用。

罗智勇等[6]提出单窗算法。李天宏等[7]发展了组合窗算法，其做法是根据下垫面的水分状况和土地覆盖类型进行分区，在整个流域不分区、水分状况分区和土地覆盖分区双因子分区2种情形下，对比了7种常用的地表温度遥感反演裂窗算法的结果，在分析每种反演算法适用性的基础上，针对不同分区单元分别选择效果最好的算法组合来进行流域地表温度的反演。杨虎[8]等人应用Becker等人的裂窗算法，但是仍然是一类线性模式。王春林等[9]借助气候学模型GIS技术把劈窗LST算法发展成全天候地表温度反演方法，其基本思路是晴天和多云天气分别建模式，但也只是部分克服了云的困难。陈少辉[10]等人给出了不同空间分辨率图像的融合算法，但是对不同时间分辨率图像的融合仍然无能为力。许多模式[8][10-17]要求在局部或者特定条件才能达到预期效果，例如柳钦火[14]的方法要求水汽廓线误差在±20％之内，大气模式误差一个模式时，反演的温度均方根误差达到0.85K。还有些模式要求晴天条件。这些缺点成为开展温度反演业务化的障碍性因素。尤其我国广大南方地区多云、地形地势和下垫面特征特别复杂，使现有方法对近地面温度遥感反演以及大面积的业务化变得极其困难。这些方法都需要特定条件才能达到预期效果。另外，这些模式的输入参数有些要通过各种订正，例如大气水汽廓线订正，否则会产生大误差；模式本身计算也有误差累积误差变得不确定了。还有这些模式考虑的变化因素太少，不能反演地形地势、地表植被特征的变化对气温的影响等等。更重要的是没有讨论过模式泛化能力。

本发明采用的SVM，即Support Vector Machines，以结构风险最小化为原则，在最小化样本点误差的同时缩小模型泛化误差的上界，从而提高模型的泛化能力。Kearns和Vzliant[18]证明了只要有足够的数据，弱学习算法就能通过集成的方法生成任意高精度的估计，弱学习等价于强学习。集成学习的有效性取决于子学习精确度和子学习之间的差异度。子学习精确度越高，差异度越大，集成学习精确度就越高，泛化能力也越大。1995年，Krogh和Vedelsby[19]还给出了集成学习的泛化误差的计算公式。

桌面高性能计算机和超大量核并行计算技术，近些年来发展突飞猛进，迅速漫延到各行各业并且取得了许多突破性成果。最新的Fermi架构，它是CUDA[Compute UnifiedDevice Architecture]础架构的最新一代产品，由于价格底、体积小、功耗低、速度快，堪称当今同类顶尖技术，其CPU+GPU，即Graphics Processing Unit，架构：CPU核管理功能和GPU计算功能的分离，特别适合图像处理，其基本原因在于我们可以让GPU的一核负责一像元或者数个像元的计算，从而简化了并行计算过程。基于Fermi架构组织数以千计的核和线程，开发超大量核并行反演算法、计算算法软件，可以满足本大面积的实时反演需求。我们在CPU+C2075上实现了近地面气温反演算法，为以后应用提供了核心技术之一。

我们在研究[20][专利号201510103565.7]基于GIS,应用非线性数学理论、超级计算技术反演近地面气温取得了突破性进步成果,通过不同季节、不同条件下的53万对样本的试验验证,该项专利克服了云层、地形、植被等的干扰，同时计算精度高稳定性还有大幅提高，计算精度缩小60％以上，获得比较高的时间和空间分辨率,可以全天候监测近地面气温,通过与国外同类方法比对，大幅高于国际上的七大算法[21-27]，引入超级计算技术使计算速度提高1000倍以上,而且速度扩展容易,获得国际同行专家的一致认可,有重大意义。但是，[20][专利号201510103565.7]的各个子模式结果误差分布不太均匀,存在小部分、小概率的点其误差没有规律、没有道理，因而本发明去掉它；同时去掉的点及其邻近区域的点由另外的子模式结果填补。即把原来构成模式结果的方法从对子模式算术平均改成有选择的求和平均。也就是丢弃不合理误差，保留合理误差，测试表明结果更加精确、稳定，各项指标值都有提高。

广西近几年来从92个人工气象站的8小时一次的人工测量温度迅速发展到1400多个无人气象站1小时一次测量温度，我国其他大部省人工气象站的布置也接近尾声。但目前的地面气象站大多处于地势平坦的地方，不能充分反映由于复杂地形因素导致的气温在空间分布上的差异性。在卫星技术方面，FY2E/C&FY2D/FY3双静止卫星可以提供1小时一次的、双保险的空中监测服务；MODIS卫星可以提供多光谱36频250×250米，90米DEM模型可以提供空间分辨率高达90×90米的地形特征。虽然采用MODIS和FY2/FY3等极轨卫星热红外遥感资料，可以较好地反演近地面气温，但是无法克服云的干扰，不能实现全天候监测。

在过去的几年中，我国FY2/FY3系列、MODIS卫星技术提高比较快，各省无人气象站布置也接近尾声，但是分析应用技术和软件分析方法相对滞后。当今我们迫切需要解决的关键科技问题是如何把我国投入巨资发展起来的无人气象站、国内外卫星信息资源、地理信息资源，即地形地势植被状况等，融合起来，形成高质量，即高时间和高空间分辨率、高精确度、高可靠性等，的近地面温度影像图，来满足益日提高的国民经济和社会发展对气象信息服质量的需求。

发明内容

本发明针对现有技术的不足，提供一种利用多种遥感图像、地形图地势DEM和无人气象站数据，在多云地区，获得大面积近地面气温的时空分布的精确反演方法。

本发明是通过如下技术方案实现的：

一种近地面气温的大面积精确反演方法，包括如下步骤：

步骤Ⅰ.建立原始数据记录集并归一化处理，得到归一化数据记录集AutoVEC'

步骤Ⅱ.构造第1子模式学习集supVEC¹和第1子模式验证集exVEC¹

步骤Ⅲ.获取第1子模式modelVEC¹

步骤Ⅳ.获取第2子模式modelVEC²

步骤Ⅴ.获取第3子模式modelVEC³至第f子模式modelVEC^f

(1)将归一化数据记录集AutoVEC'和第f-1子模式modelVEC^f代入下式，计算AutoVEC'中每条数据记录的第f-1子模式计算气温

式中，j＝1,2,...,L；f≥3,且f∈Z⁺；

(2)求出第f-1子模式计算气温与实测气温MT_j的误差绝对值AE_j ^f-1，

(3)求出第f子模式在第j个样本的误差绝对值估计AE_i ^f，AE_j ^f＝Min(AE_j ¹,…,AE_j ^f ^-2,AE_j ^f-1)，

(4)按照误差绝对值估计AE_j ^f由大到小对AutoVEC'中的数据记录进行排序，

(5)取误差绝对值估计AE_j ^f排名靠前的supNum^f条数据记录作为第f子模式学习集supVEC^f，其余的作为第f子模式验证集exVEC^f，

(6)按照步骤Ⅲ的方法获得第f子模式modelVEC^f；

步骤Ⅵ.构造栅格点原始数据记录集sBlkVEC并归一化处理，得到栅格点归一化数据记录集sBlkVEC'

步骤Ⅶ.各个子模式的近地面气温反演全图

(1)将一个小区域sBlk的栅格点归一化数据记录集sBlkVEC'依次与第1子模式modelVEC¹至第f子模式modelVEC^f代入下式，依次算出各个栅格点的第1子模式栅格点计算气温至第f子模式栅格点计算气温

式中，q＝1,2,...,sBlkSamples，f＝1,2,3,...

(2)第1子模式栅格点计算气温就是该小区域sBlk的第1子模式近地面气温反演影像图，

同理，依次得到该小区域sBlk的第2子模式至第f子模式的近地面气温反演影像图，

(3)选取另外一个小区域sBlk，重复上述步骤(1)和(2)，直到所有的BlkNum个小区域sBlk做完为止，

(4)将BlkNum个小区域sBlk的第1子模式近地面气温反演影像图合并形成目标区域的第1子模式近地面气温反演全图surfT¹，

同理，依次得到目标区域的第2子模式至第f子模式的近地面气温反演全图surfT^f；

步骤Ⅷ.在近地面气温反演全图surfT^f中标记误差异常点及其邻近区域

(1)在第1子模式的近地面气温反演全图surfT¹中，读取无人气象站k所在栅格点的第1子模式计算气温k＝1,2,...,K_t，K_t为目标时间点上有气温数据的无人气象站的个数；

(2)求出与无人气象站k的实测气温MT_k的误差绝对值RealAE_k ¹，并设定误差异常点的阈值δ，若RealAE_k ¹≥δ，则无人气象站k所在栅格点为误差异常点；

(3)划定误差异常点的邻近区域，将误差异常点的邻近区域定义为无人气象站k所在栅格点到所有的邻近无人气象站所在栅格点的距离的σ倍的点构成的多边形区域，其中，σ为邻近区域界限比率阈值，σ根据实际需要设定；如果某个方向上没有其他的无人气象站，则邻近区域的要抵达目标区域的边界；

(4)相对于阈值δ，若无人气象站k所在栅格点为误差异常点，则将无人气象站k的邻近区域标记为badPolygon_δk ¹；

(5)记

(6)重复步骤(1)～(5)，依次求出第2子模式至第f子模式相对于阈值δ的

(7)标记

步骤Ⅸ.利用选择性求和算法求出近地面气温精准反演影像图surfT

针对阈值δ＝δ₀(0＜δ₀)

(1)记为切割surfT^f所得的区域，

(2)定义surfT^f中像素点pixl^f的求和权值surfParameter_δ ^f，

如果则surfParameter_δ ^f的值为0，

如果则surfParameter_δ ^f的值为1，

如果对所有f＝1,2,3,...都成立，即pixl^f∈mostbadPolygon_δ，则surfParameter_δ ^f的值为1，

(3)选择性求和法计算近地面气温初步反演影像图，公式如下：

式中，f＝1,2,3,...；

步骤Ⅹ.对surfT^f中像素点对所有f＝1,2,3,...都成立的区域的进一步处理

(1)如果mostbadPolygon_δ为空集，则步骤Ⅸ得到的近地面气温初步反演影像图即为近地面气温精准反演影像图surfT；

(2)如果mostbadPolygon_δ非空，则在目标区域中单独将这个区域截取出来，设定下一个级别的阈值δ＝δ₁，重复步骤Ⅷ和Ⅸ，直到mostbadPolygon_δ为空或者阈值δ＝δ_s，这里多个阈值δ的临界值分级遵从0＜δ₀＜δ₁＜δ₂＜...δ_i...＜δ_s，得到这个区域的近地面气温初步反演影像图，最后将这个区域的近地面气温初步反演影像图代替步骤Ⅸ得到的近地面气温初步反演影像图中的相同区域，得到近地面气温精准反演影像图surfT。

误差异常点的阈值δ＝2～3℃。

邻近区域界限比率阈值σ＝0.75～1.0。

步骤Ⅰ.建立原始数据记录集的具体步骤如下：

(1)确定需要进行近地面气温反演的目标区域以及目标区域上的无人气象站所在位置，

(2)确定一个目标时间点，并以该目标时间点为中心的左右2～72个小时为采样时间段，以采样时间段内无人气象站的实测气温的时刻为采样时间点，

(3)设原始数据记录的格式为j＝1,2,…,L，i＝1,2,…,N，MT_j为第j条原始数据记录的实测气温，为第j条原始数据记录的第i个近地面气温影响因子的原始值，N为近地面气温影响因子的个数，L为原始数据记录的条数，最后Label标记该条记录的站号和时间点，若MT_j缺测则不予收录，

将一个无人气象站在一个采样时间点的实测气温MT_j输入原始数据记录，再根据气象卫星资料、DEM资料以及天文历法规律，将该无人气象站所在位置的该采样时间点对应的第1个至第i个近地面气温影响因子的数值依次输入获得一条原始数据记录，

(4)选取另一个无人气象站或者选取同一个无人气象站另一个采样时间点，重复步骤(3)，获得下一条原始数据记录，

(5)由步骤(2)和(3)获得的原始数据记录组成原始数据记录集AutoVEC，

(6)针对该目标区域，确定另一个目标时间点，重复步骤(2)到(4)，然后将获得的原始数据记录并入原始数据记录集AutoVEC，

(7)对原始数据记录集AutoVEC中每个近地面气温影响因子的原始值进行归一化处理，得到归一化数据记录集AutoVEC’，所述原始数据记录指无人气象站位置获得的原始数据记录；

步骤Ⅱ.构造第1子模式学习集supVEC¹和第1子模式验证集exVEC¹的具体步骤如下：

(1)将归一化数据记录集AutoVEC’中的L条数据记录按照实测气温MT_j进行排序，排名前c名的数据记录和排名倒数d名的数据记录选入第1子模式学习集supVEC¹，

(2)将实测气温MT_j排名前3～5％的归一化数据记录和排名倒数3～5％的数据记录交替地选入第1子模式学习集supVEC¹和第1子模式验证集exVEC¹，

(3)随机选取其余10～12％的数据记录进入第1子模式学习集supVEC¹，剩下的数据记录进入第1子模式验证集exVEC¹；

步骤Ⅲ.获取第1子模式modelVEC¹的具体步骤如下：

利用SVM的学习算法，选用高斯核函数，根据经验确定e个初始参数γ的数值，用第1子模式学习集supVEC¹参加训练，获得e个候选模式；根据候选模式求出第1子模式验证集exVEC¹中各条数据记录的计算气温FT_j，计算FT_j与MT_j的误差，将误差最小的一个作为第1子模式modelVEC¹；

步骤Ⅳ.获取第2子模式modelVEC²的具体步骤如下：

(1)将归一化数据记录集AutoVEC'和第1子模式modelVEC¹代入下式，计算AutoVEC'中每条数据记录的第1子模式计算气温

(2)求出第1子模式计算气温与实测气温MT_j的误差绝对值AE_j ¹，按照误差绝对值由大到小对AutoVEC'中的数据记录进行排序，

(3)取误差绝对值AE_j ¹排名靠前的supNum²条数据记录作为第2子模式学习集supVEC²，其余的作为第2子模式验证集exVEC²，

(4)按照步骤Ⅲ的方法获得第2子模式modelVEC²；

步骤Ⅵ.构造栅格点原始数据记录集sBlkVEC并归一化处理的具体步骤如下：

针对目标时间点的采样时间段内的指定的一个时间点，

(1)将目标区域划分成BlkNum个正方形的小区域sBlk，再将每个小区域sBlk划分成M×M＝sBlkSamples个栅格点，其中M小区域sBlk的行或列的个数；

(2)对小区域sBlk的每个栅格点建立一条栅格点原始数据记录，栅格点原始数据记录的格式为q＝1,2,…,sBlkSamples，i＝1,2,…,N；为第q栅格点的第f子模式计算气温，为第q栅格点的第i个近地气温影响因子，根据气象卫星资料、DEM资料以及天文历法规律，将各个栅格点所在位置的目标时间点对应的第1个至第i个近地面气温影响因子的数值输入Label为缺省的标记，构造栅格点原始数据记录集sBlkVEC；

(3)对栅格点原始数据记录集sBlkVEC中每个近地面气温影响因子的数值进行归一化处理，得到栅格点归一化数据记录集sBlkVEC'。

进一步的，近地面气温的大面积精确反演方法，

(1)从MODIS卫星中读取的近地面气温影响因子如下：

①无人气象站所在位置或者栅格点代表位置的植被叶绿素吸收1C、云和植被覆盖变换2C、陆地和云的性质7C、土壤植被差异3C以及绿色植被4C，

②无人气象站所在位置或者栅格点代表位置的云和植被变换2C、云_大气性质19C以及大气层含水量，

③无人气象站所在位置或者栅格点代表位置的对流层中层湿度27C、对流层中层湿度28C、云_地面气温22C、云_地面气温23C、云_地面气温31C以及云高和地面气温32C；

(2)从风云二号卫星中读取的近地面气温影响因子如下：

无人气象站所在位置或者栅格点代表位置的FY2C_IR1亮温、FY2C_IR2亮温、FY2C_IR3亮温、FY2C_IR4亮温以及FY2C_VIS反射率；

(3)从DEM中读取的近地面气温影响因子如下：

无人气象站所在位置或者栅格点代表位置的经度、纬度、海拔高度、坡度、坡向、太阳阴影、太阳高度、太阳高度sin、水平面太阳方位、坡面太阳高度、坡面太阳高度sin、坡面小时太阳辐射能以及之前累积辐射能；

(4)从天文历法规律中读取的近地面气温影响因子如下：

无人气象站所在位置或者栅格点代表位置的当日坡地总太阳辐射能、日落时间、日出时间以及日照时数。

进一步的，近地面气温的大面积精确反演方法，先把动力因素和影响因素分别对待，引入理想条件下的变量，所述动力因素指通过假设无大气条件下地面单位面积获得的太阳辐射能量，夜间则考虑地表单位面积辐射能量的变化趋势，动力因素驱动近地面气温的变化，所述影响因素指应用近期晴天植被遥感图获得无云条件下的植被状态，还考虑了周边地势对计算点造成的阴影，实际状况的影响由临界时间的遥感技术图片提供，还特别计算了大气廓线的水汽含量，他们最终对结果的影响和贡献则由非线性模式的优化处理后决定，影响因素对近地面气温的变化有限制作用。

进一步的，近地面气温的大面积精确反演方法，所述归一化处理的方法如下：确定原始数据记录集AutoVEC中第i个近地面气温影响因子的最大值F(i)_max和最小值F(i)_min，然后根据依次求出每个近地面气温影响因子的原始值的归一化值得到归一化数据记录集AutoVEC’。

进一步的，所述的近地面气温的大面积精确反演方法，步骤Ⅱ(1)所述的c＝3～20，d＝3～20。

进一步的，所述的近地面气温的大面积精确反演方法，步骤Ⅲ所述的e＝3～10。

进一步的，所述的近地面气温的大面积精确反演方法，所述第f子模式modelVEC^f中包含下列参数：支持向量的个数supNum^f，γ^f以及b^f，其中f＝1,2,3,…。

进一步的，所述的近地面气温的大面积精确反演方法，其特征在于，所述步骤Ⅴ.近地面气温反演在超算平台上完成，采用CPU+GPU异构协同并行技术，由CPU的多个线程负责构造栅格点原始数据记录集sBlkVEC、栅格点原始数据记录集sBlkVEC的归一化处理、电脑内存数据与硬盘数据的交换、电脑内存数据与GPU内置内存数据的交换并根据GPU内置内存大小来调节每次数据的交换量；由GPU负责依次求算各个栅格点的第1子模式栅格点计算气温至第f子模式栅格点计算气温

所述植被叶绿素吸收1C为MODIS卫星1通道数据，所述云和植被覆盖变换2C为MODIS卫星2通道数据，所述陆地和云的性质7C为MODIS卫星4通道数据，所述土壤植被差异3C为MODIS卫星3通道数据，所述绿色植被4C为MODIS卫星4通道数据，所述云和植被变换2C为MODIS卫星2通道数据，所述云_大气性质19C为MODIS卫星19通道数据，所述对流层中层湿度27C为MODIS卫星27通道数据，所述对流层中层湿度28C为MODIS卫星28通道数据，所述云_地面温度22C为MODIS卫星22通道数据，所述云_地面温度23C为MODIS卫星23通道数据，所述云_地面温度31C为MODIS卫星31通道数据，所述云高和地面气温32C为MODIS卫星32通道数据。另外，气层含水量是计算量，由MODIS的2、19通道的数据根据Kaufman&Gao公式计算得到。

本发明所述太阳高度sin是指太阳高度与太阳高度角的正弦值之积。

本发明所述坡面太阳高度sin是指坡面太阳高度与坡面太阳高度角的正弦值之积。

本发明所述当日坡地总太阳辐射能是指目标时间点所在日坡地的坡面的单位面积全天太阳辐射能。

本发明所述之前累积辐射能是指当日采样时间点之前正的单位面积辐射能之和。

本发明所述坡面小时太阳辐射能是指目标时间点所在时次的单位面积太阳辐射能。

对于步骤Ⅰ(5)所述的无人气象站原始数据记录集AutoVEC中每个近地气温影响因子的数值进行的归一化处理和步骤Ⅴ(3)所述的栅格点原始数据记录集sBlkVEC中每个近地气温影响因子的数值进行的归一化处理，由于气温数值T_j不属于近地气温影响因子，因而各条数据记录对应的气温数值T_j不需进行归一化处理而要保持原来数值。

本发明获得广西自治区科技厅“科技攻关计划”项目(桂科攻1355010-5)、广西自治区林业有害生物防治检疫站自选项目《林业有害生物监测预警与应急防控关键技术研究与应用》(桂科鉴字[2014]323号)、广西区气象局气象科研计划(桂气科2016M21)的资助。在此表示致谢。

与现有技术相比较，本发明具备的有益效果：

本发明通过引入一类超级非线性算法，比较满意的解决了现有技术中存在的信息融合来反演近地气温问题，获得了反演误差与云量、地形、地势无关的满意结果。特别是融合了空间和地面信息，反演重建多云地区近地面气温，同时利用非线性技术和超级计算技术，克服了云层的干扰而精确、可靠地获得大面积地区反演近地面气温影像图。实现了多元信息的结合应用，提高了反演精度、扩大了适应的范围，达到大范围应用；与GIS结合，实现了应用对接。比许多传统做法更有意义，与我们的上个发明[20][专利号201510103565.7]相比也有较大幅度的进步、提高。

在我们的实验中，预测误差大于1℃的约25％的样本中，我们发现这些点的误差差异度很大。依据SVM理论、Kearns和Vzliant[18]关于子学习与集成学习的理论，以及Krogh和Vedelsby[19]集成学习的泛化误差的论述，我们通过恰当地设计样本分类算法构造学习样本集，来提高集成学习的精确度和泛化性能；通过设计若干子模式集成选择性算术平均合成算法，使精确度和泛化性能进一步提高。依据这些理论和我们的实验结果，结合应用GIS的空间分析技术，充分挖掘空间、地面等与近地面气温相关的信息，使模式消除了云层、地形、地势、天气背景和大气层垂直结构对近地面气温遥感反演的影响。

近地面气温遥感反演与CPU+GPU异构协同并行计算技术结合，在该研究领域内罕见。CPU+GPU异构协同并行计算构造是当今性能价格比最高的超算技术，该发明通过巧妙的分配CPU和GPU在整个问题中的处理功能。即：CPU的多个线程负责数据记录的采样、归一化处理、电脑内存数据与硬盘数据的交换、电脑内存数据与GPU内置内存的交换；GPU负责大量核心计算。本发明依据GPU内置内存大小来调节设计每次数据的交换量，达到了提高计算速度的目的。总的计算速度比传统计算的单个CPU提高1000倍以上。

发明所述算法能够得到全天候的、空间和时间高分辨率的近地面气温反演信息，为作物冷冻害、寒害、高温危害、森林有害生物等监测和评估提供比较精确的近地面气温实时监测信息。在近地面气温反演技术和结果基础上，我们可以轻易地获得蒸发蒸腾量、湿度，应用于作物滴灌溉调控等方面。项目还可以为精细区划：农业、气候、林业、水果、林业有害生物区划提供更加精细的技术支持。还有可以为精细化气象预报提供实时初始场。本项目还可以为我国正在建设的数字地球、智慧地球、数字林业、数字城市等工程提供近地面气温的数字化影像图和定点、定量数字化信息产品，意义重大。

本发明的反演模式的因子构造比传统模式的有创新意义。传统温度反演模式主要可以分为两类：一类是实验方法，主要是在实际工作中利用地面定标，根据实测处在卫星传感器过境时的地表温度，来建立图像灰度值和地表温度的回归方程，求出地表温度图像。另一类是理论方法，该方法是通过求解辐射传输方程，来消除大气影响，从而求出地表温度，如“分裂窗法”。传统温度反演多数是针对地表温度，本发明根据实际需要和我们能够获得的温度数据，针对的是近地面气温，概念上有差别，但是方法可以通用。

本发明的新思维是先把动力因素和影响因素分别对待处理，引入理想条件下的变量，例如：应用近期晴天植被遥感图获得无云条件下的植被状态；通过假设无大气条件下地面获得的太阳辐射能量，夜间则考虑地表辐射变化趋势，来构造近地面气温的动力因子；还考虑了周边地势对研究点造成的阴影。实际状况的影响由临界时间的遥感技术图片提供，还特别计算了大气廓线的水汽含量。他们最终对结果的影响和贡献则由非线性模式的优化处理后决定，以求获得非线性纠正干扰的效果，中间没有计算过程，减少了累积误差，本发明还用几个差别大的子模式合成新模式来提高模式的计算误差和泛化能力。实施例1的大面积反演绝对误差平均的MAE约为0.75-0.79℃，同时RMSEs稳定在约0.60-0.86℃。

另外，本发明除了传承了上个发明[20][专利号201510103565.7]的所有优点之外，还具有以下优点：

1.精确度比较高，结果可靠度高，泛化能力强。传统方法没有讨论过模式的泛化能力。

2.克服了云层、地形地势、地面植被等的干扰或影响。这点可以从影像图，及其误差分析得知。

3.针对该过程计算量巨大的问题，应用CPU+GPU异构协同并行计算技术，使得计算速度提高1000倍以上，可以大面积应用，而且可以通过扩展GPU数量轻易地扩展计算能力。

本发明的独特创新在于针对合成子模式的调试，提出了获得精确近地面气温的优质反演模式的4个基本原则：

(1)确保各个单独的子模式具有比较强的泛化特征，比较简单的做法是调试时尽可能使各个单个模式用于模拟的数据记录集合的MAE和RMSEs数值与用于验证的数据记录集合的MAE和RMSEs数值比较接近。

(2)δ₀趋向于临界值时，合成模式的效果最好。

δ₀的临界值定义为：使各个子模式AE≤δ的点数约为75％，本例的临界值约为δ₀＝2.0，见表3和4的第二到六行的对比。δ₀从远大于2.0而接近2.0时，各项指标快速收敛到合理最小。

(3)多个差异大的子模式的合成，可以明显改进效果。

(4)多个临界值分级0＜δ₀＜δ₁＜δ₂＜…＜δ_s可以使各个子模式AE≥δ₀的交集部分的MAE迅速收敛到比较小的值，尤其在邻近δ_i处存在一个深度凹槽，但是需要成倍增加的内存空间。

附图说明

图1为本发明实施例1得到的2015年8月21日16时广西全境的近地面气温精准反演影像图surfT，空间分辨率为90米×90米。在实际实施中，得到的surfT是能够更加容易识别的彩色图。

图2为所述误差异常点及其邻近区域的划定方法示意图。

图3为本发明实施例1所得的不同大气水分含量下的预期误差分析图。

图4为本发明实施例1所得的不同植被指数下的预期误差分析图。

图5为本发明实施例1所得的不同海拔高度下的预期误差分析图。

图6为本发明实施例1所得的不同空间区域下的预期误差分析图。

图7为本发明所述近地面气温的大面积精确反演方法的流程图。

图8为CPU+GPU异构协同并行计算技术反演近地面气温的处理流程图。CPU0触发多个CPU的多个线程并行计算，各个线程并行地处理数据，然后触发GPU大规模核做超级并行计算。

图9为实施例2各个子模式近地面气温反演全图中标记的误差异常点及其邻近区域的示意简图

图3～6中，frequency distribution是指分布频率；standard deviation是指标准误，statistics of WaterContent analysis是指针对大气水含量的统计分析，statistics of ENV analysis是指针对植被指数的统计分析，statistics of Heightanalysis是指针对海拔高度的统计分析，statistics of Region analysis是指针对空间区域的统计分析。

具体实施方式

下面结合附图和实施例详细描述本发明的具体实施方式。

实施例1

一种近地面气温的大面积精确反演方法，包括如下步骤：

步骤Ⅰ.建立原始数据记录集并归一化处理

(1)确定需要进行近地面气温反演的目标区域和目标时间点，确定目标区域上的无人气象站所在位置，并以目标时间点为中心的左右2-72个小时为采样时间段，以采样时间段内无人气象站记录实测气温的时刻为采样时间点，

本实施例以广西全境为目标区域，并搜集广西全境内1359个无人气象站每小时一次的实测气温，搜集2015年8月21日0时到2015年8月31日23时的MODIS卫星和FY2C/D卫星逐个小时的卫星资料；以2015年8月21日16时为目标时间点，采样时间段为2015年8月21日0时到2015年8月31日23时(见表1)，这里的基本依据是2015年8月21日03时(世界时UT)的MODIS卫星代表的大气、云的状态会在03时(世界时UT)左右几个小时到几天内维持不变这样的事实，持续不变的时间的长短取决于大气的稳定性，冬季长些而夏季会短些，本例的采样时段内共计配了9张MODIS图片。另外同时应用2015年1月15日06时(世界时UT)的MODIS卫星1、2、7、3、4通道的数据代表植被状态，该时间点为晴天，用来获得实际的地面植被状况，这里也是基于一个事实，即地面植被状况相对于空中云的变化是稳定的，因而可以取近期晴天遥感图来代替计算时间点的地面植被状况。

表1：本例的MODIS卫星图片分配时段表(世界时UT)

起点时间(UT)	终点时间(UT)	MODIS图片(UT)
			08/20/2015 16:00:00PM	08/21/2015 15:00:00PM	terra_2015_08_21_03_18_sy.MOD0250M_PRJ.tif
08/21/2015 16:00:00PM	08/22/2015 15:00:00PM	TERRA_2015_08_22_03_59_GZ.MOD0250M_PRJ.tif
			08/22/2015 16:00:00PM	08/23/2015 15:00:00PM	terra_2015_08_23_03_05_sy.MOD0250M_PRJ.tif
08/24/2015 16:00:00PM	08/25/2015 15:00:00PM	terra_2015_08_25_02_53_sy.MOD021KM_PRJ.tif
			08/25/2015 16:00:00PM	08/27/2015 15:00:00PM	aqua_2015_08_27_05_46_sy.MOD021KM_PRJ.tif
08/27/2015 16:00:00PM	08/28/2015 15:00:00PM	TERRA_2015_08_28_03_22_GZ.MOD021KM_PRJ.tif
			08/28/2015 16:00:00PM	08/30/2015 15:00:00PM	terra_2015_08_30_03_11_sy.MOD021KM_PRJ.tif
08/30/2015 16:00:00PM	09/01/2015 05:00:00AM	TERRA_2015_09_01_02_57_GZ.MOD021KM_PRJ.tif
			09/01/2015 06:00:00AM	09/01/2015 08:00:00AM	aqua_2015_09_01_06_04_sy.MOD021KM_PRJ.tif

(2)确定一个时间点2015年8月21日16时；

(3)现在举一个例子说明第1条原始数据记录的建立方法：

首先，柳城无人气象站，站号59041，于2015年08月22日14时记录的无人气象站实测气温MT₁＝34.0℃；然后，从MODIS卫星中读取的近地面气温影响因子如下：无人气象站所在位置的植被叶绿素吸收1C为94.0、云和植被覆盖变换2C为251.0、陆地和云的性质7C为87.0、土壤植被差异3C为115.0以及绿色植被4C为99.0，无人气象站所在位置或者栅格点代表位置的云和植被变换2C为196.0、云_大气性质19C为128.0以及大气层含水量为0.45496，无人气象站所在位置或者栅格点代表位置的对流层中层湿度27C为2719.0、对流层中层湿度28C为2.0、云_地面气温22C为3051.0、云_地面气温23C为3030.0、云_地面气温31C为2936.0以及云高和地面气温32C为2699.0；从风云二号卫星中读取的近地气温影响因子如下：无人气象站所在位置的FY2C_IR1亮温为0.000135、FY2C_IR2亮温为0.305916、FY2C_IR3亮温为0.248436、FY2C_IR4亮温为0.296953、FY2C_VIS反射率为0.298367；从DEM中读取的近地气温影响因子如下：无人气象站所在位置的经度为109.25、纬度为24.65、海拔高度为136.0、坡度为12.19449、坡向为357.13576、太阳阴影为222.0、太阳高度为67.931、太阳高度sin为0.92673、水平面太阳方位为236.1212、坡面太阳高度为59.87402、坡面太阳高度sin为0.864924、坡面小时太阳辐射能为4.15、之前累积辐射能为23.65；从天文历法中读取的近地气温影响因子如下：无人气象站所在位置的当日坡地总太阳辐射能为24.13、日落时间(北京时)为19.11、日出时间(北京时)为6.41、日照时数为12.70。

将上述气温影响因子的数值读入原始数据记录完毕后，得到的一条原始数据记录如下：

<34.0，109.25，24.65，136.0，12.19449，357.1376，94.0，251.0，87.0，115.0，99.0，196.0，128.0，0.45496，2719.0，2.0，3051.0，3030.0，2936.0，2699.0，0.000135，0.305916，0.248436，0.296953，0.298367，222.0，67.931，0.92673，236.1212，59.87402，0.864924，4.15，23.65，24.13，19.11，6.41，12.70,B59041D2015082214>

原始数据记录中第1个数值为无人气象站实测气温，第2至第37个数值按顺序依次为下列气温影响因子的数值：无人气象站所在位置的经度、纬度、海拔高度、坡度Slope、坡向Aspect、植被叶绿素吸收1C、云和植被覆盖变换2C、陆地和云的性质7C、土壤植被差异3C、绿色植被4C、云和植被变换2C、云_大气性质19C、大气层含水量、对流层中层湿度27C、对流层中层湿度28C、云_地面气温22C、云_地面气温23C、云_地面气温31C以及云高和地面气温32C、FY2C_IR1亮温、FY2C_IR2亮温、FY2C_IR3亮温、FY2C_IR4亮温、FY2C_VIS反射率、太阳阴影、太阳高度、太阳高度sin、水平面太阳方位、坡面太阳高度、坡面太阳高度sin、坡面小时太阳辐射能、之前累积辐射能、当日坡地总太阳辐射能、日落时间、日出时间以及日照时数。

(4)选取另一个无人气象站或者选取同一个无人气象站在2015年8月21日0时内的数据，重复步骤(2)和步骤(3)，获得下一条原始数据记录，

(5)由步骤(2)和步骤(3)获得的原始数据记录组成原始数据记录集AutoVEC，本实施例原始数据记录集AutoVEC共含有324668条原始数据记录，

(6)本实例只设置了1个目标时间点,

(7)对原始数据记录集中每个近地面气温影响因子的原始值进行归一化处理，得到归一化数据记录集AutoVEC’；

步骤Ⅱ.构造第1子模式学习集supVEC¹和第1子模式验证集exVEC¹

(1)将归一化数据记录集AutoVEC’中的所有数据记录按照实测气温MT_j的高低进行排序，

(2)将实测气温MT_j排名前5名的归一化数据记录和排名倒数5名的归一化数据记录选入第1子模式学习集supVEC¹，

(3)将实测气温MT_j排名前5％的归一化数据记录和实测气温MT_j排名倒数5％的归一化数据记录交替地选入第1子模式学习集supVEC¹和第1子模式验证集exVEC¹，

(4)随机选取其余15％的归一化数据记录进入第1子模式学习集supVEC¹，余下的第1子模式验证集exVEC¹；

按照上述方法，本实施例选取了30000条归一化数据记录进入第1子模式学习集supVEC¹，其余294668条归一化数据记录进入第1子模式验证集exVEC¹；

步骤Ⅲ.获取第1子模式modelVEC¹

利用SVM的学习算法，选用高斯核函数，根据经验确定5个初始参数γ的数值，γ分别为0.5、1.0、1.2、1.5、2，用第1子模式学习集supVEC¹参加训练，获得5个候选模式；根据候选模式求出第1子模式验证集exVEC¹中各条归一化数据记录的计算气温FT_j，计算FT_j与MT_j的误差，将误差最小的一个作为第1子模式modelVEC¹；

步骤Ⅳ.获取第2子模式modelVEC²

(3)取误差绝对值AE_j ¹排名靠前的supNum条数据记录作为第2子模式学习集supVEC²，其余的作为第2子模式验证集exVEC²，

(4)按照步骤Ⅲ的方法获得第2子模式modelVEC²；

步骤Ⅴ.获取第3子模式modelVEC³

(1)将归一化数据记录集AutoVEC'和第2子模式modelVEC²代入下式，计算AutoVEC'中每条数据记录的第2子模式的计算气温

(2)求出第2子模式计算气温与实测气温MT_j的误差绝对值AE_j ²，

(3)求出AE_j ¹和AE_j ²的最小作为第3子模式在第j个样本的误差绝对值估计AE_j ³，即AE_j ³＝Min(AE_j ¹,AE_j ²)，按照误差绝对值估计AE_j ³由大到小对AutoVEC'中的数据记录进行排序，

(4)取误差绝对值估计AE_j ³排名靠前的supNum条数据记录作为第3子模式学习集supVEC³，其余的作为第3子模式验证集exVEC³，

(5)按照步骤Ⅲ的方法获得第3子模式modelVEC³；

步骤Ⅵ.构造栅格点原始数据记录集sBlkVEC并归一化处理

(1)将目标区域划分成37个正方形的小区域sBlk，再将每个小区域sBlk划分成1202×1202＝sBlkSamples个栅格点，

(2)对小区域sBlk的每个栅格点建立一条栅格点原始数据记录，栅格点原始数据记录的格式为q＝1,2,…,sBlkSamples，i＝1,2,…,N；为第q栅格点的第f子模式计算气温，为第q栅格点的第i个近地面气温影响因子，根据气象卫星资料、DEM资料以及天文历法规律，将各个栅格点所在位置的目标时间点对应的第1个至第i个近地面气温影响因子的数值输入构造栅格点原始数据记录集sBlkVEC，Label为缺省的标记，

(3)对栅格点原始数据记录集sBlkVEC中每个近地面气温影响因子的数值进行归一化处理，得到栅格点归一化数据记录集sBlkVEC’，

步骤Ⅶ.各个子模式的近地面气温反演全图

(1)将一个小区域sBlk的栅格点归一化数据记录集sBlkVEC'依次与第1子模式modelVEC¹至第3子模式modelVEC³代入下式，依次算出各个栅格点的第1子模式栅格点计算气温至第3子模式栅格点计算气温

式中，q＝1,2,...,sBlkSamples，f＝1,2,3

(2)第1子模式栅格点计算气温就是该小区域sBlk的第1子模式近地面气温反演影像图，同理，依次得到该小区域sBlk的第2子模式至第3子模式的近地面气温反演影像图，

(3)选取另外一个小区域sBlk，重复上述步骤(1)和(2)，直到37个小区域sBlk做完为止，

(4)将37个小区域sBlk的第1子模式近地面气温反演影像图合并形成目标区域的第1子模式近地面气温反演全图surfT¹，

同理，依次得到目标区域的第2子模式至第3子模式的近地面气温反演全图surfT³；

(5)可以在目标区域、目标时间段内选择其他时间点重复上述(1)到(4)步，获得连续多个近地面气温反演影像图，本实例只针对指定时间点2015年8月21日0时，作反演图；

(2)求出与无人气象站k的实测气温MT_k的误差绝对值RealAE_k ¹，并设定误差异常点的阈值δ，若RealAE_k ¹≥δ，则无人气象站k所在栅格点为误差异常点；误差异常点的阈值δ通常可以取2～3℃，本实施例δ＝2.0℃；

(3)划定误差异常点的邻近区域，将误差异常点的邻近区域定义为无人气象站k所在栅格点到所有的邻近无人气象站所在栅格点的距离的σ倍的点构成的多边形区域，其中，σ为邻近区域界限比率阈值，σ根据实际需要设定；如果某个方向上没有其他的无人气象站，则邻近区域的要抵达目标区域的边界；邻近区域界限比率阈值σ通常可以取0.75～1.0，本实施例σ＝0.75；

(5)记

(6)重复步骤(1)～(5)，依次求出的第2子模式至第f子模式的badPolygon^f；

(7)标记

针对阈值δ＝δ₀＝2.0℃

(1)记为切割surfT^f所得的区域，

(2)定义surfT^f的像素pixl^f的求和权值surfParameter_δ ^f，

如果surfParameter_δ ^f的值为0，

如果surfParameter_δ ^f的值为1，

如果对所有f＝1,2,3,...都成立，即pixl^f属于mostbadPolygon_δ，surfParameter_δ ^f的值为1，

(3)选择性求和法计算近地面气温精准反演影像图surfT，计算公式如下：

式中，f＝1,2,3,...。

最终，本实施例得到的2015年8月21日16时广西全境的近地面气温精准反演影像图surfT，如图1所示，图中深色的区域为高温区域，浅色区域为低温区域。

以上实施例在下列硬件和软件的支持下完成：

硬件：AMAX超算平台，含有16个CPU；四个GPU，即四个C2075，每个C2075有448个核。软件：①Windows7操作系统。②SVM-lightVersionV6.01，用于完成步骤Ⅲ.。③ArcGIS9.3、Arc Engine9.3、GDAL1.9，用于GIS空间数据的分析计算和GIS数据的相关关键软件开发。④VS2010C++、VS2008C++/C#、MicrosoftOffice Excel2007和MicrosoftOfficeAccess 2007、Arrafire1.9用于完成步骤Ⅰ.中数据的读取、原始数据记录集AutoVEC的管理以及异构协同并行计算的相关技术的开发。

实施例2

为了说明步骤Ⅷ误差异常点及其邻近区域，我们给出如图2所示的简图。

图中每个小方格都是一个栅格点，图中包括无人气象站k₁以其邻近无人气象站k₂、k₃、k₄和k₅，表示无人气象站所在的栅格点，无人气象站k₁的邻近区域即为四边形ABCD。

为了说明步骤Ⅷ和Ⅸ，我们给出如图9所示的简图。

图9中，黑点表示无人气象站所在位置，方框表示无人气象站的邻近区域，细线条方框中无人气象站是误差正常的所以没有用黑点表示，粗线条方框表示误差异常点的邻近区域。

(a)为第1子模式的近地面气温反演全图surfT¹中针对阈值δ＝δ₀＝2.0℃标记的误差异常点及其邻近区域，(b)为第2子模式的近地面气温反演全图surfT¹中针对阈值δ＝δ₀＝2.0℃标记的误差异常点及其邻近区域，(c)为第3子模式的近地面气温反演全图surfT¹中针对阈值δ＝δ₀＝2.0℃标记的误差异常点及其邻近区域。

经计算，(a)中，无人气象站k_E3和k_A6为误差异常点，

无人气象站k_E3的邻近区域标记为即为粗线条方框E3；

无人气象站k_A6的邻近区域标记为即为粗线条方框A6；

记

经计算，(b)中，无人气象站k_E3和k_C4为误差异常点，

无人气象站k_E3的邻近区域标记为即为粗线条方框E3；

无人气象站k_C4的邻近区域标记为即为粗线条方框C4；

记

经计算，(c)中，无人气象站k_E3、k_C4和k_A6为误差异常点，

无人气象站k_E3的邻近区域标记为即为粗线条方框E3；

无人气象站k_C4的邻近区域标记为即为粗线条方框C4；

无人气象站k_A6的邻近区域标记为即为粗线条方框A6；

记

标记

在surfT¹中，粗线条方框E3、粗线条方框EA6中像素点pixl¹的求和权值surfParameter_δ ¹等于0，其余部分的像素点pixl¹的求和权值surfParameter_δ ¹等于1；

在surfT²中，粗线条方框E3、粗线条方框C4中像素点pixl²的求和权值surfParameter_δ ²等于0，其余部分的像素点pixl²的求和权值surfParameter_δ ²等于1；

在surfT³中，粗线条方框E3、粗线条方框C4、粗线条方框A6中像素点pixl³的求和权值surfParameter_δ ³等于0，其余部分的像素点pixl³的求和权值surfParameter_δ ³等于1。

经过以上surfT¹～surfT³的求和权值判断后，发现粗线条方框E3中像素点对所有f＝1,2,3,...都成立，即pixl^f∈mostbadPolygon_δ，因此，将surfT¹～surfT³中粗线条方框E3中像素点pixl^f的求和权值surfParameter_δ ^f赋值为1。

根据选择性求和法计算得到近地面气温初步反演影像图，公式如下：

式中，f＝1,2,3,...；

由于本实施例中，存在粗线条方框E3中像素点对所有f＝1,2,3,...都成立，即因此mostbadPolygon_δ非空；所以需要对粗线条方框E3中像素点进一步处理。

在目标区域中单独将粗线条方框E3这个区域截取出来，根据下一个级别的阈值δ＝δ₁＝3.0℃在surfT¹、surfT²和surfT³中标记误差异常点及其邻近区域，结果如(e)、(f)和(g)所示。此时，无人气象站k_E3及其邻近区域粗线条方框E3仅仅在surfT¹中属于误差异常点及其邻近区域，在surfT²和surfT³中属于误差正常点。从而，在surfT¹中，粗线条方框E3中像素点pixl¹的求和权值surfParameter_δ ¹等于0，在surfT²中，粗线条方框E3中像素点pixl²的求和权值surfParameter_δ ²等于1，在surfT³中，粗线条方框E3中像素点pixl³的求和权值surfParameter_δ ³等于1。再次根据选择性求和法计算得到粗线条方框E3的近地面气温反演影像图。最后，将粗线条方框E3的近地面气温反演影像图代替步骤Ⅸ得到的近地面气温初步反演影像图中的相同区域得到近地面气温精准反演影像图surfT。

对比例1

将本发明所述的选择性算术平均合成算法和算术平均值算法[20][专利号201510103565.7]的计算精度进行了对比分析，结果如表2-4所示，表2是表3和表4的对照，表2-4使用相同的数据集，只是处理方法不同。

表2：算术平均合成算法的收敛效果示意表

注：MAE是误差绝对值的平均数值，单位为℃；RMSEs是标准差，单位为℃；N＝1359，为样本总数，即无人气象站的数量。

表2中，第二列是对2015年8月21日16时广西全境的第一子模式近地面气温反演全图surfT¹的误差统计结果，第三列是对2015年8月21日16时广西全境的第二子模式近地面气温反演全图surfT²的误差统计结果，第四列是对2015年8月21日16时广西全境的由专利[20]方法获得的第三子模式近地面气温反演全图surfT³的误差统计结果，第五列是对2015年8月21日16时广西全境的由本发明获得的第三子模式近地面气温反演全图surfT^3’的误差统计结果，第六列所述算术平均方法合成结果是指子模式1(γ＝3.0)、子模式2(γ＝5.0)、由专利[20]方法获得的子模式3(γ＝2.0)以及由本发明方法获得的子模式3(γ＝5.0)四者通过算术平均算法进行合成得到的近地面气温反演影像图。

根据以上第二列到第六列所述的近地面气温反演影像图，首先在各图中读取各个无人气象站所在栅格点的计算气温，共得到1359个样本，然后求出计算气温与实测气温的误差的绝对值AE，再计算MAE(结果在第二行)和RMSEs(结果在第二行)，最后计算AE≥1℃的样本数占样本总数的百分比(结果在第四行)，计算AE≥2℃的样本数占样本总数的百分比(结果在第五行)，计算AE≥3℃的样本数占样本总数的百分比(结果在第六行)。

表3和表4中给出了两对(δ₀＝2.0；δ₁＝3.0)和(δ₀＝2.5；δ₁＝3.0)参数选择下，选择性算术平均合成算法的各个合成阶段的结果对比。

表3：两个δ(δ₀＝2.0、δ₁＝3.0)的合成模式下选择性算术平均及其每步合成结果

表3是基于表2所得结果的进一步分析，第二列是指将surfT¹与surfT²通过本发明所述步骤Ⅷ和Ⅸ的方法得到的2015年8月21日16时广西全境的近地面气温精准反演影像图surfT；第三列是指将surfT¹、surfT²以及surfT³通过本发明所述步骤Ⅷ和Ⅸ的方法得到的2015年8月21日16时广西全境的近地面气温精准反演影像图surfT；第四列是指将surfT¹、surfT²、surfT³以及surfT^3’通过本发明所述步骤Ⅷ和Ⅸ的方法得到的2015年8月21日16时广西全境的近地面气温精准反演影像图surfT。在所述步骤Ⅷ和Ⅸ的方法中，设置两级误差异常点的阈值δ，分别为δ₀＝2.0和δ₁＝3.0。另外，表3中第二行至第六行的对比指标计算方法与表2相同。

表4：两个δ(δ₀＝2.5、δ₁＝3.0)的合成模式下选择性算术平均及其每步合成结果

表4也是基于表2所得结果的进一步分析，在第二列至第四列的计算中，设置两级误差异常点的阈值δ，分别为δ₀＝2.5和δ₁＝3.0，其他步骤与表3的一致。但是每步合成结果与最终结果(表3和表4的第一至第四列)的效果没有表3的好，原因是δ₀＝2.5不靠近临界值。另外，表3中第二行至第六行的对比指标计算方法与表2相同。当算术平均算法改进成选择性算术平均算法后，通过分析2015年8月21日16时广西全境的近地面气温精准反演影像图surfT可以看出，反演效果有比较大幅提高，具体见表3第五列。

表3和表4中可以总结得到，实验分析及合成模式调试的4个基本原则：

(2)δ₀趋向于临界值时，合成模式的效果最好。

(3)多个差异大的子模式的合成，可以明显改进效果。

对比例2

针对总样本(324668条数据记录)的三个算法，即算术平均(Arithmetic MeanMethod，论文[20]和[专利号201510103565.7])、选择性算术平均(SelectedArithmeticMean Method，本发明)、选择最小(Selected MinimumMethod，为了说明问题而设计的算法)的对比分析结果参见表5。

利用实施例1的方法依次计算08/20/201516:00:00PM～09/01/201508:00:00AM之间逐个小时的广西全境的近地面气温精准反演影像图surfT；再分别利用专利[20]方法和选择最小方法依次计算08/20/201516:00:00PM～09/01/201508:00:00AM之间逐个小时的广西全境的近地面气温反演影像图。分别在广西全境的近地面气温精准反演影像图surfT和广西全境的近地面气温反演影像图中读取无人气象站的计算气温，分别获得与上述三种方法对应的三组数据样本，每组有324668条总数据样本。

表5：总样本的三个算法的分析结果

注：MAE是误差绝对值的平均数值，单位为℃；RMSEs标准差；N＝324668，为测试样本总数。其中，选择最小方法(Selected Minimum Method)是为了说明问题而设计的算法。

由表5可以看出，选择最小算法过于贪婪而失去了泛化特征，不能保证控制点外也具有同等的精准度；算术平均算法比较保守，效果改进比较慢；只有选择性算术平均算法，去掉了不具有合理误差的点，而且结果是多个模式合成结果，确保了效果好泛化特征能力强。针对总样本也有同等提高，见表5第五列；而且总体(表5第二列)和局部(表3第四列)的MAE、RMSEs、AE≥1℃的点数百分率的具体数值接近，分别为0.775℃和0.762℃、1.091℃和0.8428℃、26.78％和25.08％，总体效果和局部效果比较一致，说明结果是比较稳定的。

对比例3

我们还发现学习过程和大面积反演过程各自的样本的采样程序，以及反演程序有差别，以2015年8月21日16时广西全境的近地面气温精准反演影像图surfT为例，会造成MAE约0.1925℃的差别，见表6。这个差别可能由于学习过程的采样和反演用的是原来的程序而大面积反演过程用的是移植改进后的超算程序。但是本发明由于计算误差异常数值时，直接从各个子模式的反演影像图取得气象站点的气温的计算数值，然后与该气象站点的实际测值来计算误差，因而克服了这个误差出差别，同时还省略掉了误差订正步骤。

表6：超级计算程序由于移植原来代码而带来的误差变化情况表

注：MAE是误差绝对值的平均数值，单位为℃；RMSEs是标准差，单位为℃；N＝1359，为测试样本总数，即无人气象站的数量。

从表7可以看出，选择性算术平均算法显著的优于国际上的七大算法。

表7：各种反演算法对比的结果

注：传统上，AE<3℃的点数百分率用于判断算法优劣的指标[7]。通常传统算法在50％左右。对于90×90米分辨率，如果一个气象站在这个像素内，气象站的气温测量数值可以很好的代表这个像素的实际气温；传统算法的分辨率多在1000×1000米或者更加低(例如，10km×10km米)，这样像素的实际气温的代表成了问题，计算误差的参照对象就用卫星产品了，而这些卫星产品也是反演计算出来的。

对比例4

根据对比例2所述的在广西全境的近地面气温精准反演影像图surfT中获得的324668条总数据样本，求出无人气象站的计算气温与对应时间点的气温实测值之间的误差绝对值AE，再对AE进行统计分析，结果如表8～9和图3～6所示。表8中，将324668条数据样本按照大气含水量进行排序并划分为第一列所述的区间，第二列N表示各个区间包括的样本数量；表8中第四列的第二行和第三行的意思是，选择性算术平均算法得到结果中误差绝对值AE在1.8～2.2℃范围内的样本数量为1条，而算术平均算法得到结果中误差绝对值AE在1.8～2.2℃范围内的样本数量为8条；第五列MAE是该区间样本的误差绝对值的平均数值，单位为℃；第六列RMSEs是指该区间样本的标准差，单位为℃；第七列是指AE<1℃的样本数占该区间样本总数的百分率；表8中第一行对应图3中第一行的两个小图，图3左列对应算术平均算法，右列对应选择性算术平均算法。表9到表11与表8相似，且分别对应图4、图5和图6。图3、图4、图5和图6中，横坐标为绝对误差AE，纵坐标为误差出现的频率(样本数)，左列小图为算术平均算法的预期误差分布图，而右列小图选择性算术平均算法的预期误差分布图。从表8、表9、图3和图4中可以看出，实施例1得到的近地面气温精准反演影像图surfT的误差与大气水分含量(云量)和植被是无关的。使用相同的数据、相同的数据数据处理、相同的子模式，仅仅子模式合成方式不同，但是各项指标值得到了明显改进，最为明显的特点是右列的每个小图在邻近δ₀处，在(δ₀±0.2)区间内有个深度凹槽，在这个误差区间的样本数从最大几千降到几十甚至是零，我们称之为δ₀凹槽，这是选择性算术平均算法形成的效应。其次，右列图的AE≥δ₀区间的预期误差显著减少，而AE分布的中心点被挤压向左移并且稳定在MAE约为0.75-0.79℃，同时RMSEs稳定在约0.60-0.86℃。另外，在AE接近0处还有个凹槽使AE接近0的样本数大幅减少，表面上看是选择性算术平均算法使AE接近0处的误差增大了，但是换来了这些点的泛化能力加强，而且AE接近0没有太多实际意义。同理，在δ₁处也可能存在一个凹槽。从表10、表11、图5和图6中可以看出，实施例1得到的近地面气温精准反演影像图surfT的误差与海拔高度和空间区域也是无关的，说明与上段的误差与云量、植被关系相同。在图5和图6中，右列的每个小图在邻近δ₀处，在(δ₀±0.2)区间内有个深度凹槽，右列图的AE≥δ₀区间的预期误差显著减少，而AE分布的中心点被挤压向左移并且稳定在MAE约为0.75-0.79℃，同时RMSEs稳定在约0.60-0.86℃。

表8：实施例1中不同大气水分含量下的预期误差分析表

表9：实施例1中不同植被指数EVI下的预期误差分析表

注：EVI指数的计算使用晴天MODIS图象见表1，

EVI＝(Channel2-Channel)/(Channel2+Channel)；

大气含水量(g cm^-2)的计算采用了[28]的B.C.Gao andY.J.Kaufman公式；

δ₀邻近处(n)表示邻近δ₀处，在(δ₀±0.2)区间内分布的样本数；

δ₀邻近处改进效果(n)用两种算法邻近δ₀处，在(δ₀±0.2)区间内分布的样本数减少的个数来表示。

表10：实施例1中不同海拔高度下的预期误差分析表

表11：实施例1中不同空间区域下的预期误差分析表

注：在广西境内按照如下规则划分区域

空间区域01，经度<＝106.5；

空间区域02，经度>111.0；

空间区域03，经度>106.5但是<108.5并且纬度>＝23.5；

空间区域04，经度>106.5但是<108.5并且纬度<23.5；

空间区域05，经度>＝108.5但是<111.0并且纬度>＝24.5；

空间区域06，经度>＝108.5但是<111.0并且纬度<＝24.5但是>23.0；

空间区域07，经度>＝108.5但是<111.0并且纬度<＝23.0。

空间区域01和空间区域03属于卡斯特区域；空间区域02、05、06属于山区，空间区域05属于高山区域，空间区域06属于中等高山区；空间区域04属于半卡斯特、半土山区域；大部分空间区域07几乎是平地农田；空间区域07和空间区域04的南部属于海岸区域。

δ₀邻近处(n)表示邻近δ₀处，在(δ₀±0.2)区间内分布的样本数；δ₀邻近处改进效果(n)用两种算法在邻近δ₀处，在(δ₀±0.2)区间内分布的样本数减少的个数来表示。

结论：

1.由表8到表11及其分别对应图3到图6可以看出，不同区间段的大气含水量、EVI指数、海拔高度、空间区域，选择性算术平均算法都一致的优于算术平均算法，各个比对指标都优越。

2.选择性算术平均算法优于算术平均算法最突出的表现在δ₀邻近处，即在(δ₀±0.2)区间。

3.选择性算术平均算法结果的稳定性也比较好，不同区间段的大气含水量、EVI指数、海拔高度、空间区域的AE在0.75-0.79℃之间，同时RMSEs稳定在约0.60-0.86℃，稳定性优于算术平均算法。局部结果与总体结果也比较一致。

4.该近地面气温反演方法可以在大面积应用，因为他有两个方面的特征：其一，反演模式适应于无云、有云、多云，模式适应范围广，可以应用于大面积反演；其二，模式计算应用了超算平台，使大面积反演变成可能。

参考文献

[1]T.R.Blandford et al.,“Seasonal and synoptic variations in near-surface air temperature lapse rates in a mountainous basin,”J.Appl.Meteor.Climatol.,Vol.47,No.1,pp.249-261,Jan.2008.

[2]C.Vancutsem,P.Ceccato,T.Dinku,and S.J.Connor,“Evaluation of MODISland surface temperature data to estimate air temperature in differentecosystems over Africa,”Remote Sens.Environ.,Vol.114,No.2,pp.449–465,Feb.2010.

[3]A.A.Hassaballa,and A.B.Matori,“The use and evaluation of split-window techniques for NOAA/AVHRR surface temperature extraction overdifferent surface covers:case study(Perak Tengah&Manjong)area,Malaysia,”IJCEE-IJENS，Vol.11,No.5,pp:22-27,Oct.2011.

[4]J.L.Qin,“Preliminary inquiry on technical support system for theprecision-farming with Chinese characteristics,”Transaction of the ChineseSociety of Agricultural Engineering,Vol.17,No.3,pp.1-6,Mar.2001(in Chinesewith English Abstract).

[5]赵高祥，汪宏七.由卫星测量确定地面温度和比辐射率的算法.科学通报1997，Vo142(18)1957～1960.

[6]罗智勇，刘汉湖，杨武年.单窗算法在成都市地面温度反演中的应用研究.热带气象学报，2007Vo123(4).

[7]李天宏，莫献坤，韩鹏.基于组合不同裂窗算法反演黄河流域地表温度研究.测试技术学报.2008，Vo122(04).

[8]杨虎，杨忠东.中国陆地区域陆表温度业务化遥感反演算法及产品运行系统.遥感学报.2006，Vo110(04).

[9]王春林，唐力生，陈水森，黄珍珠，何健.寒冷灾害监测中的全天候地表温度反演方法研究.会议：2006中日低温灾害及其防御对策研讨会.

[10]陈少辉，张秋文，王乘，周建中.基于亮度相关矩的MODIS和SPOT影像融合研究.遥感学报.2006,1，Vo12(1)：90～96.

[11]江东，王乃斌，杨小唤，刘红辉.地面温度的遥感反演：理论、推导及应用.甘肃科学学报，2001，Vo113(4).

[12]历华，曾永年，怡培东等.利用多源遥感数据反演城市地面温度.遥感学报.2007，11，Vo111(6)：891～898.

[13]历华，曾永年，怡培东等.利用多角度热图像提取冠层组分温度和方向量温.北京师范大学学报.2O07，O6，Vo143(3)：891～898.

[14]柳钦火，徐希孺，陈家宜.遥测地面温度与比辐射率的迭代反演方法一一理论推倒与数值模拟.遥感学报.1998,2，Vo12(1)：1～9.

[15]闽文彬，罗秀陵，陈忠明等.应用LSF概念模型反演草冠层叶面温度的实验。遥感学报.2006，11，Vo11O(6)：941～948.

[16]张元生，郭晓，张小美等.应用静止卫星热红外遥感亮温资料反演地表温度的方法研究.西北地震学报.2004，6，Vo126(2)：113～117.

[17]赵高祥，汪宏七.由卫星测量确定地面温度和比辐射率的算法.科学通报1997，Vo142(18)1957～1960.

[18]KearnsM,LiM,Valiant L1994 J.ACM 41 1298.

[19]Krogh A,VedelsbyJ.Neural networkensembles,cross validation andactive learning.in:Tesauro G.Touretzky D.Leen Teds.Advances in NeuralInformation Processing Systems 7.Cambridge,MA:MIT Press,1996.231-238

[20]Jiang-Lin Qin,Xiu-Hao Yang,He Fu,and Xiu-Feng Lei.A NovelNonlinear Algorithm for Area-wide Near Surface Air Temperature Retrieval.IEEEJ.Sel.Top.Appl.Earth Obs.Remote Sens.VOL.9,NO.7,pp3283-3296,JULY2016.(SCI特刊,影响因子3.0以上，SCI和EI双收录).

[21]F.Becker andZ.L.Li,“Towards a local split window method over landsurface,”Int.J.Remote Sens.,vol.11,no.3,pp.369–393,Mar.1990.

[22]C.Coil,V.Caselles,and A.Sobrino,“On the atmospheric dependence ofthe split-window equation for land surface ternperature,”Int.J.Remote Sens.,vol.15,no.1,pp.105–122,Jan.1994.

[23]J.A.Sobrino,C.Coll,and V.Caselles,“Atmospheric correction forland surface temperature using NOAA-11 AVHRR channels 4 and 5,”RemoteSens.Environ.,vol.38,no.1,pp.19–34,Oct.1991.

[24]J.C.Price,“Land surface temperature measurements from the splitwindow channels of the NOAA 7 advanced very high resolution radiometer,”J.Geophys.Res.,vol.89,no.D5,pp.7231–7237,Aug.1984.

[25]A.J.Prata,“Land surface temperature derived from the advancedvery high resolution radiometer and the along-track scanning radiometertheory,”J.Geophys.Res.,vol.98,no.D9,pp.16689–16702,Sep.1993.

[26]G.B.Franca and A.P.Cracknell,“Retrieval of land and sea surfacetemperature using NOAA 11 AVHRR data in north-eastern Brazil,”Int.J.RemoteSens.,vol.15,no.1,pp.1695–1712,Jan.1994.

[27]Y.H.Kerr,J.P.Lagouarde,and J.Imbernon,“Accurate land surfacetemperature retrieval from AVHRR data with use of an improved splitwindowalgorithm,”Remote Sens.Environ.,vol.41,nos.2–3,pp.197–209,Mar.1992.

[28]B.C.Gao and Y.J.Kaufman,“Water vapor retrievals using ModerateResolution Imaging Spectroradiometer(MODIS)near-infrared channels,”J.Geophys.Res.,vol.108,no.D13,p.4389,Jul.2003,doi:10.1029/2002JD 003023.

Claims

1.一种近地面气温的大面积精确反演方法，其特征在于，包括如下步骤：

(4)选取另一个无人气象站或者选取同一个无人气象站另一个采样时间点，重复步骤Ⅰ.(3)，获得下一条原始数据记录，

(5)由步骤Ⅰ.(2)和步骤Ⅰ.(3)获得的原始数据记录组成原始数据记录集AutoVEC，

(6)针对该目标区域，确定另一个目标时间点，重复步骤Ⅰ.(2)到步骤Ⅰ.(4)，然后将获得的原始数据记录并入原始数据记录集AutoVEC，

步骤Ⅱ.构造第1子模式学习集supVEC¹和第1子模式验证集exVEC¹

步骤Ⅲ.获取第1子模式modelVEC¹

步骤Ⅳ.获取第2子模式modelVEC²

(4)按照步骤Ⅲ的方法获得第2子模式modelVEC²；

步骤Ⅴ.获取第3子模式modelVEC³至第f子模式modelVEC^f

(1)将归一化数据记录集AutoVEC'和第f-1子模式modelVEC^f-1代入下式，计算AutoVEC'中每条数据记录的第f-1子模式计算气温

式中，j＝1,2,...,L；f≥3,且f∈Z⁺；

(3)求出第f子模式在第j个样本的误差绝对值估计AE_j ^f，AE_j ^f＝Min(AE_j ¹,…,AE_j ^f-2,AE_j ^f ^-1)，

(6)利用SVM的学习算法，选用高斯核函数，根据经验确定e个初始参数γ的数值，用第f子模式学习集supVEC^f参加训练，获得e个候选模式；根据候选模式求出第f子模式验证集exVEC^f中各条数据记录的计算气温FT_j，计算FT_j与MT_j的误差，将误差最小的一个作为第f子模式modelVEC^f；

(7)令f-1＝2，则f＝3，执行步骤Ⅴ.(1)～步骤Ⅴ.(6)的步骤，得到第3子模式modelVEC³；依次令f-1＝3、4、…、f-1，则f＝4、5、…、f，执行步骤Ⅴ.(1)～步骤Ⅴ.(6)的步骤，得到第4子模式modelVEC⁴到第f子模式modelVEC^f；

针对目标时间点的采样时间段内的指定的一个时间点，

(3)对栅格点原始数据记录集sBlkVEC中每个近地面气温影响因子的数值进行归一化处理，得到栅格点归一化数据记录集sBlkVEC'；

步骤Ⅶ.各个子模式的近地面气温反演全图

式中，q＝1,2,...,sBlkSamples，f＝1,2,3,...

(3)选取另外一个小区域sBlk，重复上述步骤Ⅶ.(1)和步骤Ⅶ.(2)，直到所有的BlkNum个小区域sBlk做完为止，

(5)记

(6)重复步骤Ⅷ.(1)～步骤Ⅷ.(5)，依次求出第2子模式至第f子模式相对于阈值δ的

(7)标记f＝1,2,3,...；

针对阈值δ＝δ₀(0＜δ₀)

(1)记为切割surfT^f所得的区域，

(2)定义surfT^f中像素点pixl^f的求和权值surfParameter_δ ^f，

如果则surfParameter_δ ^f的值为0，

如果则surfParameter_δ ^f的值为1，

式中，f＝1,2,3,...；

2.如权利要求1所述的近地面气温的大面积精确反演方法，其特征在于，误差异常点的阈值δ＝2～3℃；邻近区域界限比率阈值σ＝0.75～1.0。

3.如权利要求1或2所述的近地面气温的大面积精确反演方法，其特征在于，

(1)从MODIS卫星中读取的近地面气温影响因子如下：

(2)从风云二号卫星中读取的近地面气温影响因子如下：

(3)从DEM中读取的近地面气温影响因子如下：

(4)从天文历法规律中读取的近地面气温影响因子如下：

4.如权利要求1或2所述的近地面气温的大面积精确反演方法，其特征在于，先把动力因素和影响因素分别对待，引入理想条件下的变量，所述动力因素指通过假设无大气条件下地面单位面积获得的太阳辐射能量，夜间则考虑地表单位面积辐射能量的变化趋势，动力因素驱动近地面气温的变化，所述影响因素指应用近期晴天植被遥感图获得无云条件下的植被状态，还考虑了周边地势对计算点造成的阴影，实际状况的影响由临界时间的遥感技术图片提供，还特别计算了大气廓线的水汽含量，他们最终对结果的影响和贡献则由非线性模式的优化处理后决定，影响因素对近地面气温的变化有限制作用。

5.如权利要求1或2所述的近地面气温的大面积精确反演方法，其特征在于，所述归一化处理的方法如下：确定原始数据记录集AutoVEC中第i个近地面气温影响因子的最大值F(i)_max和最小值F(i)_min，然后根据依次求出每个近地面气温影响因子的原始值的归一化值得到归一化数据记录集AutoVEC’。

6.如权利要求1或2所述的近地面气温的大面积精确反演方法，其特征在于，步骤Ⅱ(1)所述的c＝3～20，d＝3～20。

7.如权利要求1或2所述的近地面气温的大面积精确反演方法，其特征在于，步骤Ⅲ所述的e＝3～10。

8.如权利要求1或2所述的近地面气温的大面积精确反演方法，其特征在于，所述第f子模式modelVEC^f中包含下列参数：支持向量的个数supNum^f，γ^f以及b^f，其中f＝1,2,3,…。

9.如权利要求1或2所述的近地面气温的大面积精确反演方法，其特征在于，所述步骤Ⅴ.近地面气温反演在超算平台上完成，采用CPU+GPU异构协同并行技术，由CPU的多个线程负责构造栅格点原始数据记录集sBlkVEC、栅格点原始数据记录集sBlkVEC的归一化处理、电脑内存数据与硬盘数据的交换、电脑内存数据与GPU内置内存数据的交换并根据GPU内置内存大小来调节每次数据的交换量；由GPU负责依次求算各个栅格点的第1子模式栅格点计算气温至第f子模式栅格点计算气温