一种遥感地表温度产品角度效应校正方法与流程

本发明涉及热红外遥感地表温度反演领域，尤其涉及一种基于多种modis地表温度发射率产品的地表方向发射率确定方法，及其在校正modis温度产品角度效应中的应用。

背景技术：

地表温度常用来估计地球大气界面能量收支，也是众多生态、气候、水文等模型中不可缺少的参量。随着人类科技的进步，生产力的提高，大气中碳排放量大大增加，加剧了温室效应，全球变暖已经从一个概念变成了可以觉察得到的事实，全球变暖致使两极冰川融化，海平面上升，极端气候频发。地表温度是人类研究全球变暖规律，预测全球升温趋势重要依据，研究结果是制定碳排放政策，达成国际气候协定的重要理论基础，因此精确测量地表温度十分重要。

遥感技术的出现使得获取全球或区域尺度的地表温度成为了可能，常见的温度反演算法除了温度发射率分离算法，单通道算法，劈窗算法都需要地表发射率的先验知识。地表发射率确定的方法有很多，例如，根据ndvi确定地表发射率(valore,casellesv.mappinglandsurfaceemissivityfromndvi:applicationtoeuropean,african,andsouthamericanareas[j].remotesensingofenvironment,1996,57(3):167-184.)，根据tisi指数确定发射率(beckerf.temperature-independentspectralindicesinthermalinfraredbands[j].remotesensingofenvironment,1990,32(1):17-33.)，使用α残余法确定地表发射率(kealyps,hooksj.separatingtemperatureandemissivityinthermalinfraredmultispectralscannerdata:implicationsforrecoveringlandsurfacetemperatures[j].ieeetransactionsongeoscience&remotesensing,2002,31(6):1155-1164.)。地表发射率方向性和地表的强异质性的存在增加了准确获得像元级发射率的难度。中等分辨率成像光谱辐射计(modis：moderateresolutionimagingpectroradiometer)11_l2系列温度产品未考虑发射率的方向性，基于使用基于分类发射率的静态发射率确定方法也存在普适性弱的缺点。而发射率每变化0.01，对温度反演产生的误差为0.5k以上。

为了提高modis11_l2温度产品的精度，需要获取地表像元尺度上的发射率方向性变化规律。目前，建立发射率角度效应模型主要有纯经验法、半经验的核驱动模型法和物理模型法。纯经验模型有如使用一个余弦函数来拟合地表的发射率变化规律，不足之处是地表的强异质性使得该模型有一定的局限性(a.j.landsurfacetemperaturesderivedfromtheadvancedveryhighresolutionradiometerandthealong-trackscanningradiometer.ii:experimentalresultsandvalidationofavhrralgorithms[j].journalofgeophysicalresearchatmospheres,1993,99(d6):13025-13058.)；半经验的方法代表性的有经验核驱动模型，vinnikov等提出了地表温度的各向异性使用发射率核和太阳核来表示，不过，其模型是利用六个surfrad站点数据训练得到，surfrad站点代表了有限的地表覆盖状况，使得该模型有了一定的局限性(vinnikovky,yuy,goldbergmd,etal.angularanisotropyofsatelliteobservationsoflandsurfacetemperature[j].geophysicalresearchletters,2012,39(23):23802.)；除此之外还有一些几何光学模型和辐射传输模型，这些模型为了准确模拟地表的辐射传输规律，需要准确确定地表的众多参数，而大多数参数的确定都十分困难，所以这些模型一般在地表或者实验室选择代表性样本来训练模型，普适性也不佳。

技术实现要素：

本发明所要解决的技术问题是针对现有技术校正温度角度效应具有空间的局限性以及地形对热辐射的影响未受到重视的问题，结合已有的研究，数据和方法，提出一种普适性强的像元级发射率方向性模型，并将其应用于modis温度产品角度效应校正中。

为了解决上述技术问题，本发明提供的遥感地表温度产品角度效应校正方法，包括：

步骤1：计算像元modis31,32波段的亮度温度(t31、t32)；

步骤2：获得像元大气下界温度(tair)、大气总水蒸汽含量(cwv)、传感器观测天顶角(θv)；

步骤3：使用像元大气下界温度(tair)、大气总水蒸汽含量(cwv)、传感器天顶观测角(θv)从劈窗算法系数查找表中获得像元传感器观测方向通用劈窗算法系数c，a1，a2，a3，b1，b2，b3；

步骤4：确定像元的天顶发射率(ε0)；

步骤5：计算像元的坡度(α)和坡向

步骤6：计算本地反射角(θ0)；

步骤7：获得地表不同覆盖类型纯像元在传感器观测角度﹣65°到﹢65°之间的131个观测角度每个角度的平均发射率数据集；

步骤8：将步骤7中结果都减去获得到不同观测天顶角相对天顶方向的发射率差数据集；

步骤9：将像元在观测方向相对天顶方向发射率变异的表达式假设为：

vε(θ)＝a1θ+a2θ²+a3θ³+...+anθⁿ

确定多项式中的系数a1，a2，a3，...an以及n的值；

步骤10：确定像元的在观测方向的发射率表达式；

步骤11：确定观测方向的发射率，结合亮度温度，劈窗算法系数计算地表温度(ts)。

在上述方法中，在步骤1中使用mod021km数据和普朗克函数近似公式计算像元modis31,32波段的亮度温度(t31、t32)，

t31＝1304.413871/ln(1+729.541636/b31)

t32＝1196.978785/ln(1+474.684780/b32)。

在上述方法中，在步骤2中，使用mod07_l2数据获得像元大气下界温度(tair)、大气总水蒸汽含量(cwv)，以及从mod021km数据中获得传感器天顶观测角θv。

在上述方法中，在步骤4中，使用mod21a22001年至2010年共十年数据计算像元十年历史月平均发射率，将结果作为地表天顶发射率数据集，根据像元获取月份确定像元的天顶发射率(ε0)。

在上述方法中，在步骤5中，使用nasaglobe1kmdem数据和arcmap中aspect和slope功能计算像元的坡度(α)和坡向

在上述方法中，在步骤6中，使用如下公式计算本地反射角

cosθ0＝acos(cosαcosθv+sinαsinθvcos(φs-φ))，

其中，式中φs，θv分别是传感器观测方位角和观测天顶角。

在上述方法中，在步骤7中，使用2006年至2010年的modis地表覆盖类型数据mod12q1和modis地表发射率数据mod11b1，地表覆盖数据经过聚合升尺度为和发射率数据相同的分辨率，获得地表不同覆盖类型纯像元在传感器观测角度﹣65°到﹢65°之间的131个观测角度每个角度的平均发射率数据集。

在上述方法中，在步骤9中，根据像元的地表覆盖类型，使用相应的数据集拟合确定多项式中的系数a1，a2，a3，...an，使用留一法确定n的值。

在上述方法中，步骤10中，结合步骤4的结果确定像元的在观测方向的发射率的表达式为：

ε(θ)＝ε(0)+vε(θ)＝ε(0)+a1θ+a2θ²+a3θ³+...+anθⁿ。

在上述方法中，在步骤11中，根据步骤6中结果θ0代入步骤10中的发射率表达式从而确定观测方向的发射率ε(θ0)，结合亮度温度，劈窗算法系数使用下式计算地表温度(ts)：

本发明的遥感地表温度产品角度效应校正方法结合已有的研究、数据和方法，提供了一种普适性强的像元级发射率方向性模型，并且可以应用于modis温度产品角度效应校正中。通过本发明的方法能够有效地对角度效应进行校正之后，提高了地表温度的测量精度。

附图说明

图1为本发明的遥感地表温度产品角度效应校正方法的技术流程图。

图2为使用modis021km数据，应用本发明的方法温度角度效应校正前后的对比图。

图3为步骤7中mod12q1和mod11b1数据的覆盖范围。

具体实施方式

根据本发明的实施例，本发明提供了一种遥感地表温度产品角度效应校正方法，其技术流程图如图1所示。具体地，本发明的遥感地表温度产品角度效应校正方法，包括以下步骤：

步骤1：使用mod021km数据和普朗克函数近似公式计算像元modis31,32波段的亮度温度t31和t32，

t31＝1304.413871/ln(1+729.541636/b31)；

t32＝1196.978785/ln(1+474.684780/b32)；(式1)

步骤2：使用mod07_l2数据获得像元大气下界温度tair，大气总水蒸汽含量cwv；从mod021km数据中获得传感器天顶观测角θv；

步骤3：使用tair，cwv，θv从劈窗算法系数查找表中获得该像元传感器观测方向通用劈窗算法系数c，a1，a2，a3，b1，b2，b3；

步骤4：使用mod21a2发射率产品，2001-2010年共十年数据，通过数据质量控制选择计算该像元十年历史月平均发射率，将该结果作为地表天顶发射率数据集，根据像元获取月份确定其天顶发射率ε0；

步骤5：使用noaaglobe1km数字高程(dem)数据和arcmap中aspect和slope功能计算该像元的坡度α和坡向

步骤6：使用下式计算本地反射角θ0；

cosθ0＝acos(cosαcosθv+sinαsinθvcos(φs-φ))；(式2)

φs，θv分别是传感器观测方位角和天顶角；

步骤7：使用modis地表覆盖类型数据mod12q1(2006-2010)和modis地表发射率数据mod11b1(2006-2010)，图3为该步骤中的mod12q1和mod11b1数据的覆盖范围，地表覆盖数据经过聚合升尺度为和发射率数据相同的分辨率，获得地表不同覆盖类型纯像元在传感器观测角度﹣65°到﹢65°之间的131个观测角度每个角度的平均发射率数据集；

步骤8：将步骤7中结果都减去得到不同观测天顶角相对天顶方向的发射率差数据集；

步骤9：将该像元在观测方向相对天顶方向发射率变异的表达式假设为：

vε(θ)＝a1θ+a2θ²+a3θ³+...+anθⁿ；(式3)

根据像元的地表覆盖类型，使用相应的数据集拟合确定多项式中的系数a1，a2，a3，...an；为了避免过度拟合，使用留一法(leave-one-outcrossvalidation，loo-cv)确定n的值；

步骤10：结合步骤4的结果确定该像元的在观测方向的发射率表达式为：

ε(θ)＝ε(0)+vε(θ)＝ε(0)+a1θ+a2θ²+a3θ³+...+anθⁿ；(式4)

步骤11：根据步骤6中结果θ0代入(式4)从而确定该观测方向的地表发射率ε(θ0)，结合步骤1中的亮温，步骤3中的劈窗算法系数计算地表温度ts

本发明选取了2010年7月26日16：45(utc时)包含了四个surfrad观测站点的一幅影像反演温度来检验方法的有效性。如图2，其中(a)为角度效应校正之前的modis11_l2温度，(b)校正之后的温度。已有结论证明modis温度偏低，尤其是在干旱和半干旱地区，偏低可达到3-5k。从总体上看，角度效应校正后温度有提高，这说明了角度效应校正之后温度精度有提高。从细节上来看，角度效应校正之后，地表的细节显现了出来，这与地表温度的实际分布是相符的。

角度效应消除后的温度反演结果，选取位于surfrad站点处的像元，使用surdrad站点处上行和下行长波辐射计算地表真实温度，并检验modis角度效应消除后的温度，具体结果请见表1：。

表1

从偏差(bias)和均方根误差(rmse)两个指标检验温度验证的结果，在四个站点处选取了1555个样本用来验证，发现偏差(bias)从-1.47k减小到-0.82k，均方根误差(rmse)从2.13k减小到了1.55k，这都反映了方法的有效性。

本发明的遥感地表温度产品角度效应校正方法结合已有的研究、数据和方法，提供了一种普适性强的像元级发射率方向性模型，并且可以应用于modis温度产品角度效应校正中。通过本发明的方法能够有效地对角度效应进行校正之后，提高了地表温度的测量精度。

以上实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。