m第4章 表观热惯量遥感信息模型

第4章表观热惯量遥感信息模型4.1热惯量的概念4.2表观热惯量(ATIM感信息模型4.3ATI遥感图像的计算4.4小结物体的热惯量在热力学中是一个不变的物理量。因此，研究地物的特性可从不变的热惯量入手，例如不同的岩石、土壤、水体、植被等等。当热惯量有差别时，则可用热惯量作为识别地物的条件之一。研究热惯量一直是遥感中的难题之一，世界上只有少数几个国家的科学工作者推演了热惯量的遥感信息模型［1］，且至今尚没有发表过推演的过程。推演的假设有几种条件，本章也只能从一种假设出发进行推演。目的是使读者了解推导过程，以培养严谨的思维方法。4.1热惯量的概念热传导方程物体在不同空间点有不同的温度，且温度还随时间变化，因此可以用温度场T(x,y,z,t)来描述温度的时空分布。温度场中存在热量流动，热量由较热部分向较冷部分流动，且温度差越大，热量流动得也越快。在场中任取一曲面单元dS，且确定法线方向n。根据热传导理论，在时间dt内沿法线n流过dS的热量dQ与dt，dS以及温度沿曲面法线的方向导数8*成正比，以k(>0)为比例系数(称为热传导系数，单位是J/(m・s・K)，即焦/(米•秒•开))。at…，于是，有dQ=-k dSdt8n这样，当；8T-an为负时，即沿法线”的方向T变小时，热量dQ则为正。若以q表示热流向量，由热传导理论有q=-kgradT。于是，即有dQ=qndSdtat其中qn为在法线n上的投影，即q广-k言在物体内部任取一由闭曲面S所围成的区域V，可以看到：在时间dt内，从区域V流过曲面S向外流出的热量为Q=-dtjjk-gradT-ds这里认为gradT(=8T8n)>0，即沿外法线方向的温度是上升的。应用奥氏公式R］，并将上式变号，可得时间dt内向区域V内部流入的热量为dtjjjdiv(k-gradT)dV (4—1)(V)式中grad表示梯度，div表示散度。另一方面，设以p表示物体的密度(单位是kg/m3)，yh表示物体的比热容(单位是J/(kg-K))，在时间山内，体积元素dV的温度增加dT=8Tdt8n

所以在这段时间内，这个体积元素流入热量为 YpdVdT=ypdV^LdtH Han从而在时间dt以内，全区域就需要吸收热量 dJJLHp江dV (4—2)(V)根据能量平衡原理：流进和流出区域的热量相等，故(4—1)和(4—2)所表达的值相等，即JJLp°dV=ff\div(k-gradT)dV(V) (V)(4—(4—3)Yp-=diVk-gradTat这就是热传导方程。对于均匀物体来说，k，yh，p均为常数，如果展开(4—3)式，则有atYHpata(at)kIaxJa-+fat)k{ayJa-+fat)k=k'a2ta2t、ax2 ay2a2t}+az2JIazJaxayaz所以，可得ata2ta2ta2t=D(++)(4—4)atHax2ay2a乙2式中 Dh=-^ (4—5)YHPDh称为热扩散系数，也称为温度传导率(单位m2/s)。(4—4)式是均匀物体的热传导方程。热惯量若对方程(4—3)在一维情况下给予一个简单的温度周期函数为边值条件，求解方程，可以得到一个简单的公式(具体推导见4.2节)TQ=-" (4—6)式中Q为热量；To为温度变化的振幅；①为周日角频率；P是物体固有的常量，称为热惯量，其表达式为 P=\：kY丑「 (4—7)显然，热惯量P的单位是J/(m2•s1/2-K)。从(4—6)式中可见，当地物吸收或释放的热量相等时，地物温度变化的振幅与地物的热惯量成反比，即：P值大的地物，日温差小；P值小的地物，日温差大。在力学牛顿第二定律F=ma中，当作用力相等时，加速度与质量成反比，人们把物体不变的质量看成是运动的惯性。P与惯性有类同之处，因此称为热惯量。热惯量的实质是地物阻止其温度变化幅度的一种特性。从微观的角度看，热惯量是阻止物体内部分子运动速度变化的阻力。在地物的周民温度变化中，热惯量起着决定性的作用。热惯量是地物的体内特性，不是表面特性，因此不能用一般的遥感方法直接测量。热惯量P中的热传导系数k、密度p、比热容yh都难以直接与各不同波段的遥感信息相对应。因此，必须建立热惯量P的遥感信息模型。从20世纪40年代以来，有许多科学工作者研究热惯量，直到1976年，KahleA等人作出了一张加里福尼亚Pisgah火山口熔岩湖的热惯量图⑶。它引起了各国遥感工作者的重视。之后，我国也有人开始研究这一问题；1986年，我们利用航空遥感数据制作出表观热惯量图，建立了表观热惯量遥感信息模型囹。1990—1992年，又制作出真实热惯量图；建立了真实热惯量遥感信息模型，并在此基础上完成了土壤含水量的遥感信息模型⑸。本章特讨论表观热惯量遥感信息模型，下一章中将讨论真实热惯量和土壤含水量的遥感信息模型。4.2表观热惯量(ATI)遥感信息模型(一)一维热传导方程到目前为止，所有的热模型都是求解一维热传导微分方程。若只考虑地面向空中辐射的垂直分量，并且假定地面物质为均匀介质，则热传导方程(4—4)成为一维形式，方向指向地面，直至均匀土壤的下界面。方程的表达式为(4—8)式中x为一维空间坐标；t为时间坐标；T(x,t)为地面以下x深处、t时刻的温度；Dh是热扩散系数。由(4—5)式和(4—7)式，可以推出(4—9)D=(k/P)2(4—9)(二)边值条件及其简化求解物体表面温度的热传导微分方程，需要一个满足方程和边值条件的形式解。如果不给定边值条件，热传导微分方程就会有无穷多个解。因此，在找到热传导方程可能的形式解之前，我们要讨论边值条件，它必须满足周日循环的能量守衡。在正常情况下，控制地表温度的决定性能量是太阳向地球辐射能量，以及地球本身的辐射能量和地表与大气之间的热传输。这些辐射能共同作用，导致地表温度的周期性变化(有日周期，也有年周期)，我们主要讨论日周期。地球表面的辐射能量只。n(即地表净辐射通量)具有如下平衡关系式Rsn=H+E+G (4—10)式中的度量单位是W/m2(瓦/米2)，或者J/(m2・s)［焦/(米2-秒)］。其中只为由地表到大气的显热通量；E为由地表到大气的潜热通量；G为地面的热通量；Rsn为地表净辐射通量，是来自太阳的辐射Rsn、天空的辐射Rsk以及地球发射的辐射虬/所共同作用的结果，且有Rsu=Rsn+Rsk+Rea (4—11)对于不同形式的边值条件，有三种不同的解法。这三种解法是：Jaeger1958年所采用的拉普拉斯变换法⑹，Kahle1976年所采用的有限差分法［3］和PrattD1978年所采用的傅里叶级数法［7］。PrattD所采用的傅里叶级数法最符合温度周日变化的地学模型，且其物理过程清楚，参数少，数学方法简便，因此，我们以此为基础进行推导。其边值条件为-k°T(*")|=(1-ABE)SCcos0+£R -R-H-E (4—12)Q* x=0 0T sk,dea式中，ABE为反照度；S0为太阳常数；Cn是太阳短波辐射的大气透过率；0为太阳天顶角；

cos0=cosScos0^cos3t+sinSsin。/其中：8为太阳赤纬；①是地球自转的角速度，即周日的角频率；0a为被测点纬度)；£为地物的发射比；Rskd是天空向下的长波辐射；Rea是地面向上的长波辐射；H为由地表到大气的显热通量；E为由地表到大气的潜热通量。等式右边的第一项代表地表吸收的净太阳辐射，后边几项都与大气和地面有关。(4—12)式本身符合(4—10)式所描述的地表辐射能量的平衡关系式，实际上是地球表面辐射能量的热平衡方程。(4-12)式没有考虑地形影响，因此，下面的所有推导都只适用于平原地区。PrattD对(4—12)式进行了简化，得到线性化的边值条件-k°T(*")|=(1-ABE)SCcos9-[A+BT(0,t)] (4—13)故*=0 0T C其中，AC和B与气象条件和地表状况有关。在本章中，我们所在的实验区是平原地区，其地表平坦，面积较小，气象条件可视为是均一的，所以，AC和B均为常数。一维热传导方程的求解我们的目标是以(4—13)式为边值条件，求解一维热传导方程(4—8)。图4.1所示地表温度0 6 12 18 24t图4.1地表温度周日的变化图4.1 地表温度周日的变化周日的变化情况，(4—13)式的边值条件是符合这一规律的。从图4.1中我们可以看出，地表温度周日的变化并不是严格的正(余)弦波形。但是，为了给后面的求解作准备，作为一种最简单的近似，我们可以假设在x=0处的边值条件为T(0,t)=T0cos(st) ①=2n/t (4一14)式中T0为地表温度变化的振幅；tp为温度变化的时间周期。方程(4—8)满足边值条件(4—14)的形式解并不难求，其表达式为*T(*,t)=T0exp-J①/2Dh)cosc^t-J①/2气)令H=加/2DH则 T(*,t)=Te-H*cosSt-Hx) (4—15)在固体岩石、土壤等内部，热量的传导符合热传导定律，于是Q=aQ(*,t)|=—kaT(*,t)|=.2TkHcosSt+")=ET?cosSt+")f at *=0 a* *=0 0 4 04H通过x=0的总能量只要对上式积分即可。由于假设温度是周期性的，在一个周期内积分值，2k亳T°7TH为零，即在0〜tp/2与t，2k亳T°7THQ=ftp/2Qdt=ftp"2rsT— cos(st+—)dt=-0f0 。弋D 4H由于热惯量P-JkYHp及热扩散系数Dh，于是P=^-——，从而"hQ=-'2：rTP此即(4—6)式。现在，我们以前面的(4—13)式为边值条件，求解热传导方程(4—8)由于热惯量P-JkYHp及热扩散系数Dh，于是P=^-——，从而"hQ=-'2：rTP此即(4—6)式。现在，我们以前面的(4—13)式为边值条件，求解热传导方程(4—8)，变化的物理过程和各参数的意义。以便了解周日温度(4—16)、-k。T(xG| -(1-abe)sccos0- 0T-[Ac+BT(0,t)](4—17)(4—17)式可以写成a (1-abe)S令4(x,t)-T(x,t)-b at翌官，代入(4—16)式(4—18)(4—18)式在x=0处的边值条件为NOt) 史+(1-abe)S0CTcos0,B(4—19)由于cosO是周期函数，显然①(x,t)也为周期函数，①(x,t)可以展为傅里叶级数，即尖。(0,t)=Zacosn①t-£)其中，An和£n分别是cosO的傅里叶展开式各简谐分量的振幅与相位。注意到cosO是偶函数，故冒0。对各简谐分量用(4—15)式，可得。(x,t)=-—+ 0~~TZAexp(-H7nx)cos(n①t-Hnx)BB n(4—20)n-0利用Carlaw等人的积分法，由(4—20)式就可以求出表面温度圆T0A1abescAcosn(①t-8),B 0T；(4—21)n\B2+n①P2+<2nBPn—0(4—21)式建立了表面温度T(0,t)(从而也就有昼夜温差AT)和反照度A、热惯量P的函数关系式。式中初位相8-arctan■v2B+Ptn①，％B，JS°是可以测量的常数，而卜是cosz的傅里叶展开式各简谐分量的振幅，从而2rt/22rt/2八A=一jpcoS)cosn①t)dtnn=1,2,3,(4—22)t-tp/22〃 1 .由(4—22)式，可得A=—sin6simpsi珅+cos5cosp(sii20+2。)A=2s*sinpasinn。)+2'°'°°pa[nsinn。)co。-cosn。)sin。]n n兀 兀(n2-1)其中，。=arccos(tan8tanp)。不难看出，An随着n的增加很快地趋于0，对于昼夜温差AT来说，气(基频的振幅)起着决定性的作用。从图4.2中可以清楚地看到，在傅里叶级数解法中，用基频的振幅求AT=(Tmax-Tmin)是一种很好的近似。图4图4.2地表温度的傅里叶解法第一项和所有项的比较在同一个日周期内，地面昼夜的温差为AT(0,t)=T(0,t)-T(0,t) (4—23)dadanini式中Tda，Tni分别为白天、夜间的温度；tda，tni分别为卫星白天、夜间的过境时间。将(4--22)式代入(4—23)式，可得1-ABESc1-ABESc玉A[cosn(①t-8)-cosn①t0Tn=0-8)(4—24)式中，B，C所决定。i dan …• i \B2+n①P2+v2n①BPS0在短时期内均为常数，因此(1-ABE)/AT、的值唯一地由地物的热惯量P对NOAA卫星而言，选择两次过境时间分别为周日温度的最大和最小时刻，根据上面的分析，当AT=(Tmax-Tmin)时，可以用基频振幅来近似，(4--24)式可以简化为1-ABE1-ABE\B2+①P2+"2①BP(4—25)式中，A=—sin8sinpsin^+—~cos8cosp(sin2。+2。)，它是太阳赤纬和当地纬度1兀 a 2兀的因数，在同一地区短期内是常数。由(4—25)式可见，(1-ABE)/AT值是随P的增大单调上升的。令(1-ABE)/AT=ATI，我们将它称为表观热惯量。由此得出两点结论：(1-ABE)/AT值的大小反映地物热惯量P的相对大小。(1-ABE)/AT值大，P也大;(1-ABE)/AT值小，P也小。在日周期内，如果不同地物所吸收的太阳能量相同，即(1-ABE)相同，那么热惯量大的地物昼夜温差AT小，热惯量小的地物昼夜温差AT大。如果不要求测量热惯量的绝对值，只要相对大小，则只需要测出各种地物的ABE和AT，计算(1-ABE)/AT值，并编码显示，就可以得到一张表观热惯量图。4.3ATI遥感图像的计算根据上一节的推导，表观热惯量ATI=(1-ABE)/AT。这样，求表观热惯量即为如何求ABE和AT。下面分别加以讨论。(一)反照度ABE的计算反照度的定义是地物全波段的反射比，反射度为各个波段的反射系数。因此，反照度为地物波长从0到8的反射比。假定地物具有朗伯体特性，反照度可以表示为Pp(人)Q^(人)d人J"Q^(人)d人式中p(X)是地物分光反射比，Q(X)是太阳分光辐照度。由于太阳能量主要集中在0.2一1.5gm内很窄的波段，因此。可以用可见光和近红外波段的反照度近似代替波长从0~8的反照度；对于分波段情况，则有Qen(X)Ep(X)ABE= n (4—27)EQen(X)n式中n为波段序号，Q.抖(X)为各个波段上地表的太阳分光辐照度，pQ)为地物各个波段的分光反射比。若令加权系数Wn=Qse匕Q，于是se,nnABE=EWp (4—28)n此式表明反照度ABE是地物各个波段的分光反射比的加权平均，权重Wn代表了入射到地表的第n波段的太阳能量与所有波段的太阳总能量之比。对于NOAA卫星AVHRR传感器，可见光和近红外波段为ch.1(0.58〜0.68gm)和ch.2(0.725〜1.0gm)。热红外波段内太阳能量很弱，地物反射的能量非常弱，可以略而不计，因此，只能用ch.1和ch.2两个波段推求反照度，即n=2。由Wn的定义，Qsen是地表各个波段内太阳的分光辐射能。对于大气外，Qsen可以从大气外太阳分光辐照表按波段累加得到。如果考虑了大气对辐射能量的衰减作用，，把大气外太阳分光辐照度乘以相应波长的晴空大气分光透过率，然后按波段累加，得到Qsen和Wn的值，见表4.1。 口表4.1 Qsn和可牛的值分段地 ~8e,n 1表i 大气外Q…/(W/cm2)WQ…/(W/cm2)WCh.11.59700.42302.15040.4214Ch.22.19300.57702.93320.5786由表4.1可以看出，太阳在各个波段上的辐射能，大气外和地球表面差别很大，而Wn的值变化不是很大(<0.00016)。这就是说，在计算Wn时是否考虑大气的衰减问题，对反照度计算的影响不是很显著。我们使用的一组wn是地球表面的值。这样，用NOAA气象卫星的AVHRRch.1和1ch.2求反照度A8E，可以采用下述公式ABE=0.423Gch1+0.577Gch2 (4—29)(二)昼夜温差AT的计算昼夜温差是在白天与黑夜的实际温度图配准的基础上，相减而得。由卫星影像直接得到的是辐射亮度温度，如何将辐射温度反演成实际温度，近年来很多学者开展了大量的研究工作。各种方法的目的之一是避开比发射比£(入)的直接求得。对于一个窄的红外通道，如NOAA气象卫星的AVHRRch.4(10.3〜11.3um)，可以认为s(X)与入无关。根据维恩位移定律入maxT=2.8978x10-3mK，以及地面温度约为300K，可以知道地物发射的辐射能的中心波长"maxR10.3〜11.3pm也就是说，ch.4集中了主要的地物辐射能量。设£(")与入无关，根据灰体辐射的斯式蕃一玻耳兹曼定律[9] M=£oT4bb式中，Mbb为地物辐射的通量密度；z为斯式善一玻耳兹曼常量(5.6697x10-16W-m-2-K-4)；T为地表实际温度；£是地物的发射比。根据辐射亮度温度的定义，则有ST4=zT4.b式中，Tb为地物的辐射温度。由上式可得(TV£=f (4—30)ITJ另一方面，地物的辐射温度与实际温度之间，存在下述关系(4—31)T=(4—31)人ln[£exp(C/人T)—8+1]式中C2是第二辐射常量(1.4388x10-2m・K)。将(4—30)式代入(4—31)式，得(4—32)exp(C2/人T)—1(4—32)(4—33)

a(ex-1)=x4则(4—32)式简化为 a(ex-1)=x4 (4—34)a(ex-1)=x4对于函数f(x)=ex-1,我们希望用多项式来进行逼近。根据f(x)=ex-1的图形(见图4.3),(4—35)可设f(x)=ex一1=ax4+ax3+a.x2(4—35)这样，方程(4—34)可整理成 (aa3-1)x2+aa2x+aa1=0求解这个一元二次方程，可得(舍去其中一个负根)—aa+\：a2a2—4aa(aa—1)(4—36)(4—37)x= 2 1 3 (4—36)(4—37)2(aa—1)于是 T=C2/Xx这样，对每一辐射温度Tb，有(C/入T)4exp(C/入T)—1由此求出a,再由(4—36)式求出2,进而由(4—37)式推出地表对应的实际温度。上面推导中的(4—35)式的系数，可以用待定系数法求得，a1=40.8274,a2=-17.1795,缶=2.04°利用上述方法，计算出昼夜温度图，相减就可以得到昼夜温差图。根据实测结果，一般情况下ATb=27K，而实际AT=27.71K,显然AT-ATb。用上述函数拟合法得到的昼夜温差与辐射温度的昼夜温差，最大绝对误差小于1K°用地物的辐射温度昼夜温差来近似代替实际温度，引起的绝对误差很小，并且可以把大气衰减和忽略地物发射比(£=1)造成的误差降到最小。这样，用天温度最大时刻的ch.4图像(记为Tch.4da)和夜间温度最小时刻的ch.4图像(记为Tch4ni)，就可以求出昼夜温差 °,a4=Tch.4,da-Lni (4一其中Tch4da是白天的辐射温度，Tch4ni是夜间的辐射温度。辐射温度也可以叫做表观温度。表观温度与实际温度并不相同，但是，我们在上面已经指出，表观温度的温差与实际温度的温差是十分接近的。因此，昼夜温差AT可以用上述公式计算。(三)表观热惯量ATI的计算利用(4—25)式的模型，对于NOAA的AVHRR图像，在配准的基础上，通过(4—29)式和(4—38)式分别计算出反照度ABE和昼夜温差』T,则可以得到表观热惯量ATI=(1-ABE)/AT对其编码表示，就可以获得表观热惯量图。本章中，我们从求解一维热传导方程出发，推导出表观热惯量的计算方法。综合上面的讨论，对于NOAA的AVHRR传感器，热惯量