謝俊峰, 莫 凡, 奚紹禮,唐洪釗,褚 存

1 自然資源部國(guó)土衛(wèi)星遙感應(yīng)用中心, 北京 100048

2 遼寧技術(shù)工程大學(xué)測(cè)繪與地理科學(xué)學(xué)院, 阜新 123000

3 重慶交通大學(xué)土木工程學(xué)院 重慶 400074

4 遼寧科技大學(xué)土木工程學(xué)院, 鞍山 114000

玉龍雪山位于我國(guó)的云南省麗江市,呈南北走向,海拔最高為5500余米,其溫度氣候呈現(xiàn)垂直差異,山頂終年積雪,區(qū)域溫差較大。玉龍雪山自然保護(hù)區(qū)留有我國(guó)珍貴的古冰川遺跡,育有完整自然生態(tài)系統(tǒng),具有極其重要的科研及旅游價(jià)值[1]。由于其特殊的地理位置,玉龍雪山對(duì)全球氣候變化趨勢(shì)敏感性較強(qiáng),因此玉龍雪山的生態(tài)環(huán)境變化可以很好地反映全球氣候變化[2]。

陸地表面溫度是生態(tài)環(huán)境分析等研究應(yīng)用的關(guān)鍵指標(biāo),在全球環(huán)境變化和生態(tài)資源監(jiān)測(cè)等領(lǐng)域具有重要的研究?jī)r(jià)值[3]。隨著衛(wèi)星遙感技術(shù)的不斷發(fā)展,衛(wèi)星平臺(tái)的承載力逐漸增大、穩(wěn)定度不斷提高,為大區(qū)域陸地表面溫度定量化地反演提供了一種有效的技術(shù)手段。目前,利用衛(wèi)星熱紅外遙感影像反演陸地表面溫度在地礦資源的探測(cè)、城市熱島效應(yīng)分析、植被病蟲害狀態(tài)監(jiān)測(cè)等領(lǐng)域都取得了重要的研究成果。黃秀華等利用不同時(shí)相下熱紅外遙感數(shù)據(jù)提取熱異常標(biāo)志,指出了找油的最佳時(shí)機(jī)[4]。張小飛通過(guò)提取相關(guān)的下墊面類型、地表溫度和植被覆蓋等信息,證實(shí)城市區(qū)域植被覆蓋狀況可直接影響城市地表溫度[5]。毛克彪提出利用MODIS數(shù)據(jù)基于輻射傳輸方程來(lái)進(jìn)行地表溫度反演[6]。周紀(jì)構(gòu)建立了城市熱島容量計(jì)算模型,綜合城市熱島強(qiáng)度、足跡面積等多方面信息, 在最大程度上定量地反映和刻畫了熱島效應(yīng)發(fā)生的顯著程度, 能夠客觀描述城市熱島對(duì)局地氣候、人居環(huán)境質(zhì)量的影響[7]。徐永明利用MOIDS數(shù)據(jù)提出了改進(jìn)的溫度-植被指數(shù),在農(nóng)田區(qū)域及農(nóng)作物生長(zhǎng)期內(nèi)具有很好的適用性和精度,為有效獲取大范圍農(nóng)田氣溫提供了新的思路[8]。Brabyn Lars從Landsat 7 ETM +影像提取地表溫度,指出利用Landsat 7 ETM +得出的地表溫度可廣泛地應(yīng)用于生態(tài)和氣象研究[9]。盧顯等基于中國(guó)HJ-1B衛(wèi)星熱紅外數(shù)據(jù)驗(yàn)證了2010年4月14日青海玉樹M_S7.1地震LST臨震地表溫度異常的可靠性[10]。馮海英等人利用MODIS反演的地表溫度和土地覆蓋產(chǎn)品提出一種評(píng)估森林調(diào)節(jié)溫度生態(tài)服務(wù)價(jià)值的新方法,為正確認(rèn)識(shí)和定量評(píng)估森林生態(tài)服務(wù)價(jià)值及進(jìn)行生態(tài)系統(tǒng)管理提供參考[11]。吳文淵等利用Landsat 8 OLI/TIRS遙感影像反演地表溫度,結(jié)果表明地表熱環(huán)境分布與斷裂帶等自然因素存在一定相關(guān)性[12]。目前,基于兩線陣、多線陣傳感器獲取的立體像對(duì)重建高精度數(shù)字表面模型的研究也取得了很大的進(jìn)步[13]。林行剛基于DSM討論了精確地形對(duì)飛行器跟蹤制導(dǎo)技術(shù)的重要工程應(yīng)用前景[14]。王明華采用DSM完成了工程巖體三維地質(zhì)建模與可視化[15]。穆超等基于高分辨率遙感影像,實(shí)現(xiàn)了DSM建筑物點(diǎn)的提取,對(duì)建筑物的精確三維重建提供了有效幫助。Ola Friman等采用數(shù)字表面模型來(lái)估計(jì)入射光的不同分量,進(jìn)而預(yù)測(cè)不同光照條件下的被測(cè)光譜[16]。Mathieu Brédifd等提出一種從數(shù)字表面模型提取建筑物占地面積的全自動(dòng)框架,給出了利用DSM構(gòu)建 3D城市模型的全自動(dòng)過(guò)程[17]。Lasse Sander等利用SRTM數(shù)字地面模型獲取到寬灘脊系統(tǒng)的海拔趨勢(shì)[18]。凌成星等提出了一種基于GeoEye-1立體像對(duì)提取平均樹高的方法,可快速獲得研究區(qū)大范圍森林平均樹高[19]。解金衛(wèi)等人提出一種高精度DSM估計(jì)方法,能有效提高植被區(qū)DSM反演精度[20]。

熱紅外遙感能夠獲取地表溫度信息,可見光遙感具有出色的地面分辨率,在對(duì)地進(jìn)行精確定位及立體成像方面具有無(wú)可比擬的優(yōu)勢(shì)。現(xiàn)階段圍繞數(shù)字表面模型以及熱紅外遙感地表溫度反演兩個(gè)領(lǐng)域的應(yīng)用研究都有了很多重要成果,但是通過(guò)兩者的融合應(yīng)用開展多維遙感信息探測(cè)還處于初始階段,尚有大量空白空間需要國(guó)內(nèi)外學(xué)者進(jìn)行探索。將星載熱紅外數(shù)據(jù)與可見光數(shù)據(jù)融合使用,突破單一數(shù)據(jù)的桎梏,可進(jìn)一步拓寬并提升遙感應(yīng)用的邊界與高度,構(gòu)建地表三維溫度場(chǎng)。相比較單一的熱紅外遙感應(yīng)用,更真實(shí)地反應(yīng)研究區(qū)的實(shí)際地表情況,結(jié)合溫度信息可以實(shí)現(xiàn)立體式熱紅外探測(cè),通過(guò)利用多維信息主被動(dòng)復(fù)合觀測(cè),顯著提升目標(biāo)多維度信息綜合監(jiān)測(cè)能力。

本文從星載熱紅外與可見光立體遙感數(shù)據(jù)融合應(yīng)用出發(fā),以研究玉龍雪山的溫度變化為目的,基于Landsat系列衛(wèi)星的熱紅外數(shù)據(jù)反演二維地表溫度,基于資源三號(hào)衛(wèi)星的三線陣可見光影像構(gòu)建數(shù)字表面模型,通過(guò)兩者嵌套完成地表三維溫度場(chǎng)構(gòu)建,定量分析研究區(qū)域的地表三維溫度場(chǎng)時(shí)空變化。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)

研究數(shù)據(jù)選取區(qū)域?yàn)橛颀堁┥絽^(qū)域,大致位于我國(guó)麗江市內(nèi)100°4′2″E至100°16′30″E、27°3′2″N至27°18′57″N處,地理位置如圖1所示。研究數(shù)據(jù)選用陸地衛(wèi)星(Landsat)系列的熱紅外影像以及立體測(cè)繪衛(wèi)星資源三號(hào)(ZY3)的三線陣影像為基礎(chǔ)研究數(shù)據(jù)。Landsat TM以及Landsat 8 TIRS提供了從1984年至今的熱紅外數(shù)據(jù),本文選取了從1987年至2018年的時(shí)間跨度中,以兩到三年為時(shí)間間隔的13個(gè)年份中每年的12月熱紅外影像數(shù)據(jù)進(jìn)行處理,并且通過(guò)對(duì)比選擇該月份中質(zhì)量較高的日期的影像數(shù)據(jù)(由于數(shù)據(jù)質(zhì)量原因,在2004年至2008年的數(shù)據(jù)為四年的時(shí)間間隔)。數(shù)據(jù)在美國(guó)地質(zhì)勘探局網(wǎng)站(USGS)下載得到。資源三號(hào)衛(wèi)星的三線陣影像數(shù)據(jù)包括前視、后視及全色影像,軌道行列號(hào)Path為33,Row為157,數(shù)據(jù)通過(guò)自然資源部國(guó)土衛(wèi)星遙感應(yīng)用中心獲取得到。由于較早時(shí)期的熱紅外影像數(shù)據(jù)無(wú)相對(duì)應(yīng)的資源三號(hào)數(shù)據(jù),以相近時(shí)相下的Landsat ETM+ 全色影像模擬為高分辨率可見光影像對(duì)其進(jìn)行影像配準(zhǔn)。以上數(shù)據(jù)均以Geotiff格式儲(chǔ)存,采用了UTM-WGS84坐標(biāo)系統(tǒng)進(jìn)行地圖投影。

圖1 玉龍雪山地理位置示意圖

1.2 方法

1.2.1普適性單通道算法反演地表溫度

普適性單通道算法(generalized single-channel method, GSC)顧及了地表比輻射率和大氣輻射參數(shù),是一種可靠性較高的溫度反演算法。其可適用于Landsat 4 TM熱紅外數(shù)據(jù)、Landsat 7 ETM+熱紅外數(shù)據(jù)和Landsat 8 TIRS數(shù)據(jù),算法公式如下[21-22]：

(1)

(2)

(3)

式中,TS為地表溫度,單位為K。Ti為波段通道的亮度溫度,單位為K。Li為大氣頂層衛(wèi)星傳感器所接收到的輻射亮度,單位為W m-2sr-1μm-1；bγ=c2/λ,c2為普朗克函數(shù)中的常量,c2=1.43877×104μm·K,λ是有效波長(zhǎng)。對(duì)于Landsat 4 TM第6波段,bγ為1290K。對(duì)于Landsat 5 TM第6波段,bγ為1256K。對(duì)于Landsat 7 ETM+第6波段,bγ為1277K。對(duì)于Landsat8 TIRS第10波段,bγ為1324K。Ψ1、Ψ2及Ψ3是關(guān)于大氣水汽含量ω的大氣參數(shù),由下式及表中的計(jì)算系數(shù)利用大氣水汽含量計(jì)算得到：

(4)

1.2.2鉛垂線軌跡法制作數(shù)字表面模型

圖2 鉛垂線軌跡法匹配制作數(shù)字表面模型

鉛垂線軌跡法(VLL——Vertical Line Locus)是制作數(shù)字表面模型過(guò)程中常用的技術(shù)策略之一。其基本思想是：假設(shè)在物方有一條鉛垂線軌跡,則它在畫幅式像片上的投影也是一條直線,鉛垂線與地面交點(diǎn)A在像片上的構(gòu)像必定位于相應(yīng)的投影輻射線上,原理示意圖如下[13,23]。

利用鉛垂線軌跡法搜索其相應(yīng)的像點(diǎn)a1和a2,確定A點(diǎn)高程的基本計(jì)算步驟如下：給定地面點(diǎn)的平面坐標(biāo)(X0,Y0)與可能的最低高程Zmin,高程搜索步距dZ可以按照所要求的高程精度決定。由地面點(diǎn)的平面坐標(biāo)(X0,Y0)與可能的高程：Zi=Zmin+i·dZ,(i=0,1,…n),計(jì)算地面點(diǎn)分別在左右像片上的像點(diǎn)坐標(biāo)(x1,y1),(x2,y2)。分別以(x1,y1)、(x2,y2)為中心在各自影像上選取匹配窗口,計(jì)算其匹配測(cè)度,如計(jì)算相關(guān)系數(shù)ρi(也可以利用其他測(cè)度)。令i=i+1,重復(fù)上述步驟,得到一組匹配測(cè)度{ρ0,ρ1,ρ2,…,ρk},選取其中的最大值ρk,其對(duì)應(yīng)的高程為Zk=Zmin+k·dZ,則認(rèn)為地面點(diǎn)P的高程值為Zk。還可以利用ρk及其相鄰的幾個(gè)相關(guān)系數(shù)擬合一條拋物線,以其極值對(duì)應(yīng)的高程作為地面點(diǎn)的高程,以進(jìn)一步提高精度,或以更小的高程步距dZ,在一小范圍內(nèi)重復(fù)以上過(guò)程。

1.2.3三維溫度場(chǎng)構(gòu)建

地表三維溫度場(chǎng)作為星載熱紅外數(shù)據(jù)與可見光數(shù)據(jù)融合產(chǎn)品,具有幾何高精度、信息多維度的優(yōu)勢(shì)。在構(gòu)建地表三維溫度場(chǎng)時(shí),根據(jù)應(yīng)用需求,在數(shù)據(jù)選擇方面需要考慮熱紅外數(shù)據(jù)的時(shí)間分辨率和溫度分辨率,例如Landsat系列衛(wèi)星在理想狀態(tài)下獲取從1984年至今的16天時(shí)間分辨率平均誤差值在1K左右的地表溫度數(shù)據(jù),可廣泛應(yīng)用于各種熱紅外遙感的研究應(yīng)用。資源三號(hào)衛(wèi)星三線陣影像具有高精度的幾何定位的能力,采用同軌立體像對(duì)生產(chǎn)數(shù)字表面模型,為地表三維溫度場(chǎng)提供精確可靠的地面三維坐標(biāo)數(shù)據(jù)。

三維溫度場(chǎng)的構(gòu)建基本流程如下：首先對(duì)原始可見光影像數(shù)據(jù)進(jìn)行預(yù)處理,消除影像的幾何與輻射畸變等。采用熱紅外數(shù)據(jù)反演地表溫度,得到表征地表溫度信息如圖3所示的影像數(shù)據(jù)。使用可見光數(shù)據(jù)的前視、后視及正視影像制作如圖3的數(shù)字表面模型。最后將熱紅外溫度影像與數(shù)字表面模型在統(tǒng)一地理基準(zhǔn)的情況下進(jìn)行套合,完成熱紅外與可見光遙感數(shù)據(jù)地融合,得到既包含了溫度信息又具有高精度三維坐標(biāo)如圖4所示的地表溫度場(chǎng)數(shù)據(jù),技術(shù)流程如圖5所示。

圖3 溫度反演與立體影像產(chǎn)品

圖4 地表三維溫度場(chǎng)

圖5 三維溫度場(chǎng)構(gòu)建流程圖

3 結(jié)果與分析

利用普適性單通道算法進(jìn)行溫度反演,在獲得研究區(qū)域的一個(gè)時(shí)間序列的地表溫度反演結(jié)果后,為了更好的研究分析該區(qū)域地表三維溫度場(chǎng)的時(shí)間演變規(guī)律及空間分布特征,需要對(duì)地表溫度產(chǎn)品進(jìn)行溫度等級(jí)的劃分。通過(guò)對(duì)不同溫度區(qū)間進(jìn)行分層賦色,可以直觀的表現(xiàn)出該研究區(qū)域中不同溫度區(qū)間的覆蓋范圍及相對(duì)分布情況。本文通過(guò)密度分割并基于等間距分割法來(lái)劃分地表溫度等級(jí),劃分低溫等級(jí)后的溫度產(chǎn)品結(jié)果如圖6所示：

圖6 玉龍雪山區(qū)域歷年溫度平面分布圖

將以上反演得到的地表溫度產(chǎn)品,與DSM進(jìn)行融合疊加處理,獲得所有溫度產(chǎn)品對(duì)應(yīng)的地表三維溫度場(chǎng)。如圖7所示,不僅可以看出最終生成的地表三維溫度場(chǎng)直觀明了地展現(xiàn)出玉龍雪山區(qū)域在1987年到2018年時(shí)間序列中具體地表的溫度情況,同時(shí),還可以獲取到不同地形地物的具體溫度信息及三維坐標(biāo)位置信息。

圖7 玉龍雪山區(qū)域歷年地表三維溫度場(chǎng)

為了方便對(duì)玉龍雪山區(qū)域的地表三維溫度場(chǎng)進(jìn)行更直觀的定量分析,根據(jù)影像獲取時(shí)刻的地表溫度,將研究區(qū)域簡(jiǎn)單劃分為四個(gè)溫度區(qū)間：低溫區(qū)(T<-10℃)、次低溫區(qū)(-10℃10℃)。由于地表三維溫度場(chǎng)具有地面及高程的高精度坐標(biāo)信息,可以計(jì)算得到不同于平面投影面積的真實(shí)地表面積,通過(guò)溫度屬性進(jìn)行分割可以計(jì)算得到不同溫區(qū)的真實(shí)地表面積。對(duì)每個(gè)溫區(qū)的地表投影及實(shí)際面積進(jìn)行統(tǒng)計(jì)分析對(duì)比(表1、表2),并以1987年至2018年中的12月為時(shí)間序列繪制不同溫區(qū)地表面積趨勢(shì)折線圖。

表1 研究區(qū)歷年不同溫區(qū)地表投影面積統(tǒng)計(jì)表

根據(jù)上述圖表給出的結(jié)果,對(duì)玉龍雪山區(qū)域從1987—2018年的地表三維溫度場(chǎng)的時(shí)空變化特征進(jìn)行研究分析：

(1)從圖6玉龍雪山區(qū)域歷年溫度平面分布圖可以看出,由于研究區(qū)域中部的海拔較高,溫度一致呈現(xiàn)出由中間區(qū)域向外遞減的趨勢(shì),圖7玉龍雪山區(qū)域歷年地表三維溫度場(chǎng)則更直觀的呈現(xiàn)出了歷年溫度最低的區(qū)域一致集中在了海拔最高的山峰頂部,而在海拔較低及人類活動(dòng)區(qū)域溫度逐漸升高。

(2)由圖8和表3可以看出,針對(duì)玉龍雪山的時(shí)空變化分析,地表三維溫度場(chǎng)相較于二維數(shù)據(jù)分析優(yōu)勢(shì)明顯,對(duì)不同溫區(qū)真實(shí)地表面積變化情況更加敏感。其中次低溫區(qū)位于高程變化顯著區(qū)域,真實(shí)地表面積與投影面積差異最為明顯,中溫區(qū)次之?；诘乇砣S溫度場(chǎng),充分利用地形高程信息,可以表現(xiàn)出不同二維數(shù)據(jù)的真實(shí)地表情況。玉龍雪山從河谷到山頂高程差距較大,對(duì)于其垂直氣候的分層差異,表現(xiàn)出不同高程下育有不同的氣候帶。由于地表三維溫度場(chǎng)相比較平面溫度場(chǎng)在統(tǒng)計(jì)地表面積方面會(huì)因高程差異顯著更為精確,因此建立地表三維溫度場(chǎng)對(duì)于玉龍雪山等山地區(qū)域的生態(tài)研究分析意義重大。

表3 不同溫區(qū)地表面積變化斜率

(3)在1987—2018年間,研究區(qū)不同溫區(qū)的地表面積變化呈現(xiàn)出不穩(wěn)定的變化現(xiàn)象,整體變化趨于溫度升高。由表2及圖8可知,所選取的13年中,1989年玉龍雪山的整體溫度分布最低,低溫區(qū)(T<-10℃)所占面積大幅度提高,之后呈現(xiàn)出逐年遞減并在1994年后趨于穩(wěn)定的狀態(tài)。次低溫區(qū)(-10℃10℃)地表面積呈現(xiàn)出明顯上升,并在之后的較長(zhǎng)的之間序列中變化較小,但是在2000年后,不穩(wěn)定性增強(qiáng),年間溫度變化較大。從以上現(xiàn)象可以看出,隨著全球溫度上升的趨勢(shì)愈加顯著,玉龍雪山及類似山地或者冰川的自然生態(tài)所受影響也愈加明顯。同時(shí),根據(jù)相關(guān)研究表明,喜馬拉雅山冰川于1990—2015年也呈現(xiàn)出持續(xù)退縮狀態(tài)[24],并且其變化會(huì)隨著高程有所不同。全球變暖造成的影響是整體性的,同時(shí)與高程也密切相關(guān),本文方法對(duì)雪山或冰川相關(guān)的生態(tài)環(huán)境的科學(xué)研究有著一定的可行性及優(yōu)勢(shì)。

表2 研究區(qū)歷年不同溫區(qū)地表實(shí)際面積統(tǒng)計(jì)表

圖8 不同溫區(qū)地表面積變化總體趨勢(shì)圖

4 結(jié)論

本文利用陸地系列衛(wèi)星和資源三號(hào)衛(wèi)星的影像數(shù)據(jù),通過(guò)套合溫度反演產(chǎn)品和數(shù)字表面模型,構(gòu)建了云南省麗江市玉龍雪山區(qū)域的地表三維溫度場(chǎng),研究分析了該區(qū)域由1987—2018年的時(shí)空變化特征,并分析不同溫區(qū)的面積變化。根據(jù)本文地表三維溫度場(chǎng)與二維數(shù)據(jù)比較結(jié)果,針對(duì)高程起伏明顯區(qū)域,地表三維溫度場(chǎng)優(yōu)勢(shì)明顯,能夠更加有效地反映地表溫度情況,對(duì)于山地區(qū)域的研究具有重要意義。同時(shí),玉龍雪山區(qū)域研究時(shí)間序列內(nèi)的中溫區(qū)和暖溫區(qū)的地表面積有逐年增加的趨勢(shì),相較之下中溫區(qū)的地表面積升幅稍大；低溫區(qū)和次低溫區(qū)的地表面積大體呈現(xiàn)減少的趨勢(shì)；可以發(fā)現(xiàn)玉龍雪山地區(qū)的整體溫度變暖,驗(yàn)證了全球溫度變暖的趨勢(shì)。