专利名称：：从modis数据反演地表温度方法
技术领域：
：本发明涉及一种利用对地观测卫星上MODIS传感器获得的地面热辐射信息反演地表温度的方法，该方法克服了以往反演方法中参数获取不同步和复杂难以实用的缺点。能够应用在气象、农业、环境监测和旱情监测等遥感部门。
背景技术：
：地表温度是指地表表层的温度，英文是skintemperature。地表温度在区域资源环境研究中的重要性已经使热红外遥感成为遥感研究的一个重要领域，目前已经开发了很多实用的地表温度遥感反演方法，如热辐射传输方程法、劈窗算法、单窗算法和多通道算法。许反演算法是针对具体的传感器开发的，例如Qinetal[2001]针对NOAA-AVHRR开发了一个劈窗算法[QinZ.H.,OlmoGD.，andKarnieliA.,DerivationofsplitwindowalgorithmanditssensitivityanalysisforretrievinglandsurfacetemperaturefromNOAA-advancedveryhighresolutionradiometerdata,Gec^A".i饥，2001，22,655-22,670.]，由于AVHRR传感器上没有近红外波段反演大气水汽含量，只能通过气象站点或者其他传感器获得大气参数。1999、2002年搭载MODIS遥感器的对地观测卫星发射成功，为全球和区域资源环境动态监测开辟了又一新的途径。MODIS是一个拥有36个波段的中分辨率遥感系统(如图1)，每12天可获得一次全球观测数据，其飞行与太阳同步，每天同一区域至少可获得昼夜两景图像，并且是免费接收，因此非常适合于中大尺度的区域资源环境动态监测。在MODIS的36个波段中有8个是热红外波段(如表O，因而非常合适于区域尺度的地表热量空间差异分析。但是，目前针对MODIS遥感数据的地表温度反演方法不是很多，万正明等在1996和1997年在MODIS传感器还没有搭载卫星上天之前，针对这个传感器提出了两种反演算法[WanZ.andDozierJ.，Ageneralizedsplit-windowalgorithmforretrievinglandsurfacetemperaturemeasurementfromspace,rra肌Gmsc/,/emoMSe打s.，1996,34:892-905.;WanZ.M.,LiZ..L.,APhysics-BasedAlgorithmforRetrievingland-surfaceemissivityandtemperaturefromEOS/MODISdata,Gmyc/.及e附We&肌，1997,35:980-996.]。其中一个是劈窗算法，只反演地表温度；另外一个算法是能同时反演地表温度和发射率，该算法同时需要白天和晚上的数据，由于数据匹配原因和该算法需要的一些参数需要其它同时搭载的传感器或者其它反演产品获得，这两个算法被美国宇航局(NASA)采用作为产品算法。但算法比较复杂，一般科研人员难以实现，比较复杂。国内许多研究者针对MODIS传感器的地表温度反演做过一些研究，比如毛克彪和覃志豪等[毛克彪，覃志豪，施建成，宫鹏，针对MODIS数据的劈窗算法研究，武汉大学学报(信息科学版),2005(8):703-708.;毛克彪，覃志豪，施建成,用MODIS影像和劈窗算法反演山东半岛的地表温度,中国矿业大学学报(自然科学版)，2005(1):46-50.]。虽然MODIS卫星传感器搭载在美国卫星上，但我国MODIS数据地面接收站比较多，比如农业部资源遥感与数字农业重点室，中国气象局卫星气象中心，中国科学院地理所等都有接收站点。虽然美国宇航局向全球发布MODIS温度产品，其精度在美国本土比较准确，但在世界其它地区有些地方精度还是不够，需要进一步提高[WAN，Z.,ZHANG,Y.,ZHANG,Q.andLI,Z.-L.,2002，Validationoftheland-surfacetemperatureproductsretrievedfromTerraModerateResolutionImagingSpectroradiometerdata.'RemoteSensingofEnvironment,2002，83,163-180.]。需要我们根据世界各地本地的实际情况，进一步研究新算法或者做进一步的校正，以提高精度。表1MODIS遥感器技术参数<table>tableseeoriginaldocumentpage5</column></row><table><table>tableseeoriginaldocumentpage6</column></row><table><table>tableseeoriginaldocumentpage7</column></row><table>
发明内容本发明的目的在于提供一种从遥感数据MODIS反演地表温度的方法，以克服现有地表温度反演方法复杂，且难以满足科研人员个人幵发使用的需要，以及更加充实国内研究人员针对热红外遥感器开发的产品算法，为我国2020之前计划发射100颗卫星上热红外传感器提供地表温度反演方法参考，而且还能进一步提高目前农业部业务运行中地表温度反演精度，提高旱情监测和农作物估产精度[已经拟定作为农业部农情监测中的地表温度反演方法选择的方法之一]。为实现上述目的，本发明提供的从遥感数据MODIS反演地表温度方法步骤为第一步，通过大气水汽含量，计算MODIS第31和32波段的透过率1-1)利用MODIS第19波段和第2波段计算比值T，然后利用这个比值计算大气水汽含量W[KaufmanY.J.,GaoBo-Cai.，RemoteSensingofWaterVaporintheNearIRfromEOS/MODIS,压五五Traw.Geosc/.Wewo化Se朋.，1992,30:871掘.〗<formula>formulaseeoriginaldocumentpage8</formula>1-2)利用大气辐射传输模型MODTRAN4模拟计算得到大气水汽含量与MODIS第31和32波段透过率的关系，将l-l)中计算得到的大气水汽含量代入下面公式，计算得到第31波段和32波段的透过率(r31，)。<formula>formulaseeoriginaldocumentpage8</formula>第二步、通过植被指数NDVI(NormalizedDifferenceVegetationIndex)计算MODIS传感器第31和32波段的发射率。2-1)读入MODIS第l波段(近红外波段)和第2波段(红光波段)，计算<formula>formulaseeoriginaldocumentpage8</formula>其中A表示第1波段，A表示第二波段。2-2)利用NDVI来判断1平方公里分辨率像元的地物类型，然后根据美国喷汽推进实验室(JPL)测试的ASTER波谱库计算MODIS传感器第31波段和32波段的发射率(f31，￡32分别表示第31和32波段发射率)。计算如下当NDVK0时，是水体和雪，f31=0.992，￡32=0.988当0<NDVI<0.05,是裸土，s31=0.986,s32=0.991当0.05<NDVI<0.65,为植被和裸土的混合像元。利用PV指数混合像元中植被的覆盖的比率[KERR,Y.H.,LAGOUARDE,J.P.andIMBERNON,J.，1992，AccuratelandsurfacetemperatureretrievalfromAVHRRdatawithuseofanimprovedsplitwindowalgorithm.RemoteSensingofEnvironment,1992，41,197-209.]。计算公式PV=(NDVI-0.05)/0.6(式5)s31=0.986*(l-PV)+0.972*PV,s32=0.991*(1-PV)+0.976*PV;当NDVI〉0.65，为植被f31=0.972，f32=0.976第三步简化辐射传输方程3-1)简化普朗克(Planck)函数。在能量平衡方程中，每一项都包含普朗克(Planck)函数，A,r=，(e^-1)-1(式6)5"是分谱辐射通量密度，单位是『.m—2.，、A是波长，单位辉；A是普朗克常数(6.6256xl0-34J*s);c是光速(3xl08m/s);/fc是玻耳兹曼常数(1.38x10-23J/K);T是绝对温度(K)。普朗克(Planck)函数是一个非线性函数，为了简化方程组，需要对普朗克(Planck)函数进行线性简化。分别对MODIS的第31波段C10.78011.280戶)和32波段(11.77-12.27，))的热辐射与温度在273K-322K区间内的变化关系进行计算，得到如图所示的计算结果。从图中可以看出，热辐射强度随温度的变化接近于线性关系。因此，对散点图建立线性回归方程，得到对第31波段331(7)=0.13787731-31.65677，i2=0.9971对第32波段532(7>0.11849r32-26.50036，及2=0.99783-2)简化辐射传输方程组，9对于MODIS的第31和32波段，简化的辐射传输方程组可以写成如下[毛克彪，针对M0DIS数据的地表温度反演方法研究，硕士学位论文，凌二貧乂学，2004.5.]:a31(t31)=r31的&(6)万31(t;)+[i-r31的][i+(i-f31⑨^⑨]^(t;)(式7A)A2(4)=&(e)s32(0)532(rs)+[i_r32的][i+(i-f32(争32(-，32(r。)(式7B)上式中r,.(P)是第31和32波段在角度0方向的透过率，s.(0是第31和32波段在角度0方向的发射率，r。是大气平均作用温度，t;是地表温度，左边第一项s,(;.)是第31和32波段的星上亮度温度，等式右边第一项是地表辐射，等式右边第二项是大气影响。把3-1)中第31和32波段的简化方程S31(7)=0.13787r31-31.65677和532(rH).11849r32-26.50036分别代入辐射传输方程组，得到0.13787￡"31r317;=0.13787:r31+31.65677￡"31r31-(l-r31)[l+(a-f31A31](0.13787ra-31.65677)-31.65677(式8A)0.11849s32T327;=0.11849:r32+26.50036￡"32z"32-(l-ir32)[l+(l-s32)r32](0.118497>26.50036)-26.50036(式8B)为了便于计算，将上面方程组中的系数分别记为A尸O.13787*r3/+31.65677*r31*￡"31墨31.65677&尸(1〈31)*(1+(1<31)*r31)*0.13787332=0.11849*r32+26.50036*T"32*f32-26.50036C32=(l-r32)*(l+(l-s32)*r32)*0.11849""=(l-r32)*(l+(l-Q2)*r32)*26.50036辐射传输方程组可以写成r"(式9A)&(式9B)解方程组，陆地表面温度的计算公式为K=(C"^+i^)-C3//(C^;;-C"》(式10)第四步反演地表温度，4-1)利用普朗克函数(Planck)反函数计算第31和32波段的星上亮度温度T31和T32,下式中B31和B32分别是MODIS第31和32波段的星上亮度。T31=14380/(11.03*ln(2*59500000/(B31*11.03A5)+1))T32=14380/(12.02*ln(2*59500000/(B32*12.02A5)+l))4-2)将第一步和第二步计算得到的透过率和发射率，以及4-l)计算得到的星上亮度温度分别代入系数爿w，,C^，C，利用地表温度计算公式计算得到地表温度。本发明的有益效果是，考虑MODIS像元l公里分辨率的特征，利用MODI近红外和红波段计算的植被指数来计算热红外波段31和32的发射率，保证了发射率的同步获取，而且简单并保证了一定的精度。利用近红外波段反演大气水汽含量，并利用MODTRAN4建立大气透过率与水汽的关系，从而实时的获得大气透过率参数。这样解决了地表温度反演中的两个关键参数。提高了反演精度和计算时间，克服以往美国宇航局算法复杂，一般科研人员难以实现的缺点。为气候变化研究，气象预报，蒸散发，农情监测以及灾害监测等提供了有效手段和技术支撑，这个方法已经确定被农业部农情监测系统采用。下面结合附图和实施例对本发明进一步说明。图1MODIS遥感器。图2是本发明的主流程示意图。图3是环渤海地区的大气水汽含量遥感反演结果。图4是MODIS的第31波段大气透过率图。ii图5是MODIS的第32波段大气透过率图。图6是NDVI分布图。图7是MODIS第31波段的地表比辐射率分布图。图8是MODIS第32波段的地表比辐射率分布图。图9是反演而得的环渤海地区地表温度分布图。图10对比验证结果。具体实施方式这里用2003年8月上旬中国环渤海地区的一景MODIS影像数据做具体反演，本方法主要包括四个步骤，如图2。第一步，通过大气水汽含量，计算MODIS第31和32波段的透过率1-1)利用MODIS第19波段和第2波段计算比值T，然后利用这个比值计算大气水汽含量W，如图3:A0.02-lnr20.6511-2)利用大气水汽含量与MODIS第31和32波段的关系f3l=2.89798—1.88366e2'.22704t32=—3.59289+4.60414e32639将l-l)中计算得到的大气水汽含量代入上面公式，计算得到第31波段和32波段的透过率(r31，r32)，得到第31和32波段的透过率，如图4和5。第二步、通过植被指数NDVI(NormalizedDifferenceVegetationIndex)计算MODIS传感器第31和32波段的发射率。2-1)读入MODIS第l波段(近红外波段)和第2波段(红光波段)，计算NDVI;5，+及计算得到的结果如图6。2-2)利用NDVI来判断1平方公里分辨率像元的地物类型，然后根据美国喷汽推进实验室(JPL)测试的ASTER波谱库计算M0DIS传感器第31波段和32波段的发射率(&，^分别表示第31和32波段发射率)。计算如下当NDVK0时，是水体和雪，&=0.992,&=0.988当0<NDVI<0.05,是裸土，=0.986，&2=0.991当0.05<NDVI<0.65,为植被和裸土的混合像元。PV=(NDVI-0.05)/0.6f31=0.986*(l-PV)+0.972*PV,f32=0.991*(1-PV)+0.976*PV;当NDVI〉0.65，为植被s31=0.972，s32=0.976计算后得到MODIS第31和32波段的发射率分别如图7和图8。第三步建立MODIS第31波段和32波段简化的辐射传输方程组第四步利用下面公式计算MODIS第31波段和32波段的星上亮度温度T31=14380/(11.03*ln(2*59500000/(B31*11.03A5)+1》T32=14380/(12.02*ln(2*59500000/(B32*12.02A5)+l》将第一步和第二步计算得到的发射率和透过率，计算得到的星上亮度温度分别代入系数A"A;,A〃A尸O.13787*7^+31.65677*r31-31.65677G;=(l-z"31)*(l+(l-f31)*r31)*0.13787仏尸("31)*(1+(1,)、)*31.65677^32=0.11849*￡32*r32532=0.11849*7^+26.50036*r32*￡"32-26.50036C^=(l-r32)*(l+(l-s32)*r32)*0.11849A2=(l-r32)*(l+(l^2)*r32)*26.50036利用地表温度计算公式计算得到地表温度，如图9。7>(C32(^+D3;)-/(c"广c"》在应用中，对方法反演精度评价是非常重要的。因为MODIS的像元分辨率是1公里*1公里)，地表测量通常是点测量。我们很难用几个点测量的数据的平均值来代表一个大尺度的像元，而且是要在卫星过境时同时获取。在这里我们用美国通量网数据来对我们的方法进行评价。通量数据库0lttp:〃public.ornl.gov/ameriflux/datahandler.cfrn)是一个全球微气象观测网络，主要是用来监测二氧化碳交换、水蒸汽、陆地生态系统和大气能量交换。也被用来对MODIS地表温度产品进行验证[Wang,W.,S.Liang,ValidatingMODISlandsurfacetemperatureproduct,ispmsrs05,17-19，October,Beijing,China,2005.〗。这里我们也选择了3个地表比较比较平坦且均一站点(Brookings,Audubon，FortPeck)的地表温度作为实测量数据。具体介绍请参考(http:〃public.oml.gov/ameriflux/site-selectcfin)。从MODIS反演地表温度和通量站点数据比较如图10所示。平均精度大约是1.2K，能满足我们目前的应用需求。权利要求1、从MODIS数据反演地表温度方法，其步骤为第一步，通过大气水汽含量，计算MODIS第31和32波段的透过率1-1)利用MODIS第19波段和第2波段计算比值T，然后利用这个比值计算大气水汽含量W。利用大气辐射传输模型MODTRAN4模拟计算得到大气水汽含量与MODIS第31和32波段透过率的关系，计算得到的大气水汽含量代入下面公式，计算得到第31波段和32波段的透过率。<mathsid="math0001"num="0001"><math><![CDATA[<mrow><msub><mi>τ</mi><mn>31</mn></msub><mo>=</mo><mn>2.89798</mn><mo>-</mo><mn>1.88366</mn><msup><mi>e</mi><mfrac><mi>w</mi><mn>21.22704</mn></mfrac></msup></mrow>]]></math></maths><mathsid="math0002"num="0002"><math><![CDATA[<mrow><msub><mi>τ</mi><mn>32</mn></msub><mo>=</mo><mo>-</mo><mn>3.59289</mn><mo>+</mo><mn>4.60414</mn><msup><mi>e</mi><mfrac><mrow><mo>-</mo><mi>w</mi></mrow><mn>32.70639</mn></mfrac></msup></mrow>]]></math></maths>第二步、通过植被指数NDVI计算MODIS传感器第31和32波段的发射率。2-1)读入MODIS第1波段(近红外波段)和第2波段(红光波段)，计算NDVI。利用NDVI来判断1平方公里分辨率像元的地物类型，然后根据美国喷汽推进实验室(JPL)测试的ASTER波谱库计算MODIS传感器第31波段和32波段的发射率(ε31，ε32分别表示第31和32波段发射率)。计算如下当NDVI＜0时，是水体和雪，ε31＝0.992，ε32＝0.988当0＜NDVI＜0.05，是裸土，ε31＝0.986，ε32＝0.991当0.05＜NDVI＜0.65，为植被和裸土的混合像元。利用PV指数混合像元中植被的覆盖的比率。计算公式PV＝(NDVI-0.05)/0.6(式5)ε31＝0.986*(1-PV)+0.972*PV，ε32＝0.991*(1-PV)+0.976*PV；当NDVI＞0.65，为植被ε31＝0.972，ε32＝0.976第三步简化辐射传输方程3-1)简化普朗克(Planck)普朗克(Planck)函数是一个非线性函数，为了简化方程组，需要对普朗克(Planck)函数进行线性简化。分别对MODIS的第31波段(10.780～11.280μm)和32波段(11.77-12.27μm))的热辐射与温度在273K-322K区间内的变化关系进行计算，得到如图所示的计算结果。从图中可以看出，热辐射强度随温度的变化接近于线性关系。因此，对散点图建立线性回归方程，得到对第31波段B31(T)＝0.13787T31-31.65677，R2＝0.9971对第32波段B32(T)＝0.11849T32-26.50036，R2＝0.99783-2)简化辐射传输方程组，对于MODIS的第31和32波段，简化的辐射传输方程组可以写成如下B31(T31)＝τ31(θ)ε31(θ)B31(Ts)+[1-τ31(θ)][1+(1-ε31(θ)τ31(θ)]B31(Ta)B32(T32)＝τ32(θ)ε32(θ)B32(Ts)+[1-τ32(θ)][1+(1-ε32(θ)τ32(θ)]B32(Ta)把3-1)中第31和32波段的简化方程B31(T)＝0.13787T31-31.65677和B32(T)＝0.11849T32-26.50036分别代入辐射传输方程组，得到0.13787ε31τ31Ts＝0.13787T31+31.65677ε31τ31-(1-τ31)[1+(1-ε31)τ31](0.13787Ta-31.65677)-31.656770.11849ε32τ32Ts＝0.11849T32+26.50036ε32τ32-(1-τ32)[1+(1-ε32)τ32](0.11849Ta-26.50036)-26.50036为了便于计算，将上面方程组中的系数分别记为A31＝0.13787*ε31*τ31B31＝0.13787*T31+31.65677*τ31*ε31-31.65677C31＝(1-τ31)*(1+(1-ε31)*τ31)*0.13787D31＝(1-τ31)*(1+(1-ε31)*τ31)*31.65677A32＝0.11849*ε32*τ32B32＝0.11849*T32+26.50036*τ32*ε32-26.50036C32＝(1-τ32)*(1+(1-ε32)*τ32)*0.11849D32＝(1-τ32)*(1+(1-ε32)*τ32)*26.50036辐射传输方程组可以写成A31TS＝B31-C31Ta+D31A32TS＝B32-C32TA+D32解方程组，陆地表面温度的计算公式为Ts＝(C32(B31+D31)-C31(D32+B32))/(C32A31-C31A32)。全文摘要本发明涉及一种从MODIS数据反演地表温度的方法，能够应用在气象、环境监测、土地管理、农情监测、以及灾害监测等遥感应用部门。该方法，包含四个步骤第一步骤是通过大气水汽含量，计算MODIS第31和32波段的透过率。第二个步骤是通过植被指数NDVI计算MODIS传感器第31和32波段的发射率。第三步骤是简化辐射传输方程。第四步是对MODIS数据进行反演计算，得到地表目标对象温度和发射率分布情况，可以用于气象预报、环境监测、农情监测和灾情监测等部门。文档编号G01J5/00GK101629850SQ200910091030公开日2010年1月20日申请日期2009年8月24日优先权日2009年8月24日发明者周清波,张立新,李三妹,毛克彪,王道龙,覃志豪申请人:中国农业科学院农业资源与农业区划研究所;国家卫星气象中心