王修信,湯谷云,羅漣玲,孫 濤,朱啟疆
1 廣西師范大學計算機科學與信息工程學院,桂林 541004 2 廣西師范大學廣西多源信息挖掘與安全重點實驗室,桂林 541004 3 北京師范大學遙感科學國家重點實驗室,北京 100875 4 廣西師范大學生命科學學院, 桂林 541004
城市化進程所引發的熱環境問題是當前“城市氣候與環境”的研究熱點。在我國西南地區的廣西、貴州、云南等存在大面積喀斯特地貌,構成世界面積最大的喀斯特地區之一[1]。桂林自古享有“山水甲天下”的美譽,是該地區典型的喀斯特城市,特點是在城區中鑲嵌著一些喀斯特自然植被覆蓋的石灰巖孤峰或峰叢,巖溶山峰接收的太陽輻射受復雜地形因素的影響。近20多年來桂林城市隨著向郊區、向高空的快速擴展呈現越來越明顯的熱環境問題,中高層建筑逐漸增多,建筑物、道路等不透水地表逐漸代替了近郊有植被覆蓋的農田、林地以及水塘等自然地表,一些喀斯特山峰也逐漸進入城區,但少數城市周邊的喀斯特山峰被開辟為采石場,成為獲取城市建設工程石料的便捷途徑,導致喀斯特植被遭到毀滅性破壞,山體的基巖大面積完全裸露。2014年桂林獲準列入“中國南方喀斯特第二期”世界遺產名錄,如何應對城市化進程中的生態環境問題、實現可持續發展成為其城市建設急需研究的內容。
地表水熱通量(潛熱與顯熱通量)是描述近地層大氣和下墊面間水分、能量交換的參數,城市地物覆蓋變化使得地表物理特性參數發生變化,影響地表能量平衡和水熱通量,從而影響城市熱環境[2]。Ando和Ueyama[3]利用渦度相關系統對日本Sakai城市建筑區連續兩年的觀測數據分析潛熱、顯熱等能量通量的變化規律。中科院地理所匡文慧[4]研發了城市地表結構組分與熱環境生態調控模型(EcoCity),將潛熱通量和顯熱通量列為模型的10個輸入參數之中。研究城市地表水熱通量的時空變化,可從能量角度分析熱環境問題形成的物理機理,為減緩城市熱島提供科學數據。
水熱通量的地面觀測方法主要有空氣動力學法、波文比能量平衡法、渦度相關法[5-6]和大孔徑閃爍儀法[6-7],然而其數值的尺度為局地觀測點上風向一定范圍內,或大孔徑閃爍儀幾百米到數公里的較大尺度。對于城市生態系統,水熱通量受地表覆蓋類型和環境因素等的影響而高度時空變化[3,8-9],以有限的觀測點的局地實驗值來描述整個城市將帶來較大誤差,而遙感反演方法提供了準確獲取整個城市區域水熱通量信息的唯一經濟可行方法。
遙感反演水熱通量的方法有經驗方法(包括Penman-Monteith法[10]、波文比法、水量平衡法、Dalton法等)和物理方法(梯度法、熱慣量法等)。經驗方法通過擬合通量地面觀測值與相關遙感參數的數學方程,較大程度受局地環境因素的影響。熱慣量法不再依賴非遙感數據,熱通量的計算利用地表土壤對吸收輻射的響應方程,但目前只適用于裸土或沙地情形。梯度法近年來常用有單層模型(SEBAL模型、SEBS模型、METRIC模型等)和雙層模型(TSEB模型等),其首先遙感反演地表溫度,利用阻抗公式結合氣溫來計算顯熱通量,然后潛熱通量的計算利用地表能量平衡方程。雙層模型的計算過程較復雜,其將地表分為植被和土壤,增加了模型相應的輸入參數,需計算大量的阻抗。SEBS模型提出參數化計算kB-1方法,當在沒有研究區先驗了解時存在kB-1不確定性,需更多下墊面信息才能合理估計kB-1。SEBAL模型的物理概念清晰,運行時只需輸入少量氣象數據,研究對象主要為自然生態系統[10],近年來模型被嘗試應用于中國南京、印度德里等地形平坦的城市區域[11-12],由于模型的輸入參數未考慮地形因素,而喀斯特城市中存在喀斯特山峰使得研究區地形復雜,使用SEBAL模型將導致較大誤差。針對SEBAL模型存在的問題,METRIC模型綜合考慮了高程、坡度和坡向對輻射的影響[13],已成功應用于綠洲、農田、流域等[14- 21],但未見研究城市區域,原因主要是我國西南喀斯特位于經濟欠發達的西部地區,缺少項目資金支持對喀斯特城市復雜地形地表水熱通量的研究。因此,針對我國西南地區桂林喀斯特城市近20多年來快速擴展所引發的熱環境問題,改進METRIC模型使其適用于喀斯特城市實際狀況,利用模型和Landsat遙感時序圖像反演地表水熱通量,定量分析水熱通量的時空變化規律。
研究區位于中國西南巖溶地區的廣西東北部桂林主城區,石灰巖孤峰、峰林廣布地面,是典型的喀斯特巖溶地貌城市。屬亞熱帶濕潤季風氣候,氣候溫和,年平均氣溫18—19℃,最熱的7、8月份平均氣溫約28℃,最冷的1、2月份平均氣溫約9℃,最低氣溫偶爾降到0℃以下,年平均相對濕度73%—79%。四季分明,夏長冬短,無霜期長,年平均無霜期309 d。光照充足,年日照時數1447.1 h。雨量充足,年平均降雨量1887.6 mm,降雨量主要集中于4—7月,年分配不均,雨熱基本同季。全年風向以偏北風為主,平均風速2.2—2.7 m/s。
由于桂林位于亞熱帶氣候區,春季、夏季和冬季多云霧和降雨,遙感圖像中云量太大,無法準確提供地物信息,因此遙感數據選取覆蓋桂林主城區1994年10月22日TM、2000年10月30日ETM+、2006年9月21日TM、2010年11月11日TM和2015年10月16日OLI/TIRS秋季5景Landsat衛星圖像,經輻射校正、配準、幾何精糾正和大氣校正。數字高程DEM數據為ASTER GDEM,空間分辨率為30 m。遙感模型運行所需的衛星過境氣象數據來自桂林國家基本氣象站。
地表能量平衡方程是模型的原理[11]:
LE=Rn+A-HS-G-S
(1)
其中LE為潛熱通量,包括植被蒸騰和地面蒸發的能量,LE=ET,ET為蒸散量,為水的汽化潛熱;Rn為凈輻射,A為人為熱源釋放的能量,HS為顯熱通量;G為土壤熱通量,ΔS為下墊面熱能儲存。
Rn=(1-α)Sin+(Lin-Lout)-(1-ε)Lin
(2)
(3)
其中Sin為太陽入射短波輻射,Lin為大氣入射長波輻射,Lout為地表發射長波輻射,ε、α、TS分別為地表的比輻射率、反照率、溫度,Stefan-Boltzman常數σ=5.67×10-8W/m2K4。
HS=ρaCP(T1-T2)/rah
(4)
其中T1、T2分別為高度z1、z2空氣溫度,ρa、rah分別為空氣的密度、動力學阻抗,空氣定壓比熱CP=1004 J/kg K。
METRIC模型的特點為:(1)假設地-氣溫差與地表溫度TS呈線性相關方程,通過選取極端干濕“熱點”和“冷點”來計算回歸系數,避免像元空氣溫度的確定。(2)利用迭代法計算像元顯熱,確保熱量傳輸粗糙度、溫度梯度和顯熱通量之間正確的物理耦合關系。最后瞬時潛熱通量LE由能量平衡方程計算。
研究區“熱點”選取桂林城區中心廣場,為較大面積的裸露水泥地表,潛熱近似為零;研究區“冷點”選取漓江水體中心位置,顯熱近似為零,所有能量都轉變為潛熱。利用“熱點”和“冷點”的地表參數,經多次循環迭代運算獲得地-氣溫差和地表溫度TS之間的線性相關方程系數[13]。
模型輸入參數改動部分如下,其他輸入參數按Allen等[13]文獻計算。
1.4.1地表反照率的計算
地表反照率的計算采用Liang算法[22]:
αtoa=0.356ρB+0.130ρR+0.373ρNIR+0.085ρSWIR1+0.072ρSWIR2-0.0018
(5)
α=(αtoa-αpath)/2
(6)
其中,ρB、ρR、ρNIR、ρSWIR1、ρSWIR2分別為藍、紅、近紅外、短波紅外1和短波紅外2等波段反射率,α為地表反照率,αtoa為大氣外反照率,αpath為程輻射(取值0.025—0.04之間),為大氣單向透射率。
1.4.2地表比輻射率計算的改進
根據喀斯特城市的實際情況,提出估算喀斯特山峰混合像元比輻射率的方法。首先從ASTER光譜實驗數據庫(http://speclib.jpl.nasa.gov)確定典型地物的比輻射率。對TIRS波段10數據,比輻射率水體εW取0.997,植被εV取0.987,水泥建筑εB取0.965,干燥土壤εS取0.968,石灰巖εR取0.958;對TM波段6數據,比輻射率水體εW取0.995,植被εV取0.985,水泥建筑εB取0.968,干燥土壤εS取0.973,石灰巖εR取0.960。Landsat熱紅外波段的地面空間分辨率TIRS為100 m,TM為120 m,城市的遙感像元實際上是多種地物構成的混合像元。由于亞熱帶地區較好的水熱條件利于植被生長,植被基本完全覆蓋裸巖間的土壤,喀斯特山峰在遙感圖像中是由植被與裸巖構成混合像元,比輻射率計算公式為:
ε=PVRVεV+(1-PV)RRεR+dε
(7)
dε=(1-εR)(1-F)εV
(8)
建筑、道路與綠化植被構成的混合像元的比輻射率計算公式為:
ε=PVRVεV+(1-PV)RBεB
(9)
綠化植被與土壤構成的混合像元的比輻射率計算公式為:
ε=PVRVεV+(1-PV)RSεS
(10)
其中,PV為植被覆蓋度;RV、RB、RS、RR分別為植被、建筑/道路、土壤、裸巖的溫度比率,可由PV估算[23];dε為混合像元中由于地形起伏引起的比輻射率修正項;地形因子F根據不同幾何分布取值[24]。
1.4.3地表溫度反演的改進
在經典METRIC模型中使用熱紅外波段通過對比輻射率的校正計算地表溫度,而本文地表溫度的反演利用包含大氣與地表影響的單窗算法[25]。對TM/ETM+數據使用熱紅外波段6,而對TIRS 數據有兩個熱紅外波段,鑒于波段11的定標存在較大的不確定性,故使用波段10[26]。
TB= k2/ln(1+k1/L)
(11)
TS={a(1-C-D)+[(b-1)(1-C-D)+1]TB-DTa}/C
(12)
C=ε,D=(1-ε)[1+(1-ε)]
(13)
其中,L、TB分別為熱紅外波段的輻射亮度、亮度溫度; k1、k2為常數;a、b為常量[25],TS為地表溫度,Ta為大氣平均作用溫度,可由地面附近氣溫計算。
1.4.4土壤熱通量和熱能儲存、地表粗糙度的計算
參考已有的研究,利用經驗公式對土壤熱通量G進行估算[27]:
對植被覆蓋下墊面
G=(TS-273.16)(0.0038α+0.0074α2)(1-0.98NDVI4)Rn/α
(14)
對非植被下墊面
G=CgRn
(15)
其中,NDVI為歸一化植被指數。系數Cg取值,建筑和道路、水體取0.40,裸土、裸巖取0.30。城市下墊面建筑和道路等的熱能儲存,參考以往的研究用其占凈輻射比例合并到土壤熱通量中進行估算。
地表粗糙度參考已有文獻的研究結果[27]。建筑取1.0 m,以林地為主的植被取0.8 m,水體取0.000034 m,裸土取0.001 m。
由于桂林是我國西南地區的旅游城市,主城區基本沒有大型工業廠房;又屬于我國西部地區經濟欠發達的三線城市,汽車數量較少,汽車交通流密度相對不太高;遙感圖像獲取時間在氣溫舒適的秋季,一般不使用空調,因此人為熱源釋放的熱量較少,遠小于太陽輻射,在能量平衡方程中可以近似忽略。
使用機器學習方法支持向量機法對遙感圖像進行地物分類,然后根據分類結果估算地表比輻射率[26]。在遙感圖像中根據地面GPS定位數據提取典型地物樣本,分類特征向量的元素由歸一化植被指數NDVI、歸一化建筑指數NDBI、水體指數MNDWI 、光譜特征和數字高程DEM組成。支持向量機SVM的內積函數取徑向基RBF函數,使用交叉驗證法確定最優懲罰系數和間隔松弛因子。SVM的學習使用樣本訓練集,分類精度使用測試集進行檢測,將地物分為建筑(道路)、植被、水體、裸土、裸巖。2015年、2010年、2006年、2000年、1994年的總體分類精度分別為89.6%、88.7%、89.3%、87.5%、88.2%,接近90%,Kappa系數分別為0.881、0.856、0.877、0.827、0.843。
使用改進的METRIC模型反演桂林城區地表顯熱和潛熱通量,結果見圖1、圖2。反演結果的驗證,利用BREB法對2015年衛星過境時廣西師范大學校園內草地上架設的HOBO小型自動氣象站、四分量凈輻射傳感器CNR4的空氣溫濕度梯度、輻射等氣象數據估算地表顯熱與潛熱通量,根據GPS數據定位,模型值與觀測值的相對誤差見表1,模型改進前顯熱、潛熱通量反演值與觀測值的相對誤差分別為16.5%、17.8%,模型改進后顯熱、潛熱通量反演值與觀測值的相對誤差分別為14.4%、13.7%,比模型改進前精度有所提高,反演誤差與已有的研究相近[14,17],在該模型正常誤差范圍之內,說明模型估算結果是合理的。

圖1 研究區顯熱通量空間分布 /(W/m2)Fig.1 Spatial distribution of sensible heat flux over the study area

圖2 研究區潛熱通量空間分布 /(W/m2)Fig.2 Spatial distribution of latent heat flux over the study area

表1 通量模型估算值與實測值比較
統計研究區典型地物水熱通量均值和波文比(顯熱通量與潛熱通量比值)均值,結果見表2、表3。
由圖1、圖2、表2、表3明顯看出,水體和植被覆蓋地表的潛熱遠高于顯熱,水體的潛熱高于植被覆蓋地表的潛熱而顯熱低于植被覆蓋地表的顯熱,水體的顯熱低于100 W/m2,波文比數值范圍:水體在0.10—0.25,地面植被(以林地為主)和喀斯特山峰陽坡植被在0.35—0.50,喀斯特山峰陰坡植被在0.30—0.45;在太陽輻射下受地形影響,潛熱和顯熱通量均值以喀斯特山峰陽坡植被、地面植被明顯高于喀斯特山峰陰坡植被。喀斯特山峰植被像元實際上是包含少量比例裸巖的混合像元。水體(穿城而過的漓江和城區的桂湖、榕湖、杉湖等)的潛熱最高而顯熱最低;而城區茂密植被覆蓋的較大型喀斯特山峰(普陀山、疊彩山、西山、穿山、騮馬山、老人山、南溪山等)的潛熱也相對較高而顯熱相對較低。
建筑/道路、地面裸土和喀斯特山峰裸巖的顯熱高于潛熱,波文比數值范圍:建筑/道路在1.10—1.65,地面裸土在1.10—1.80,喀斯特山峰裸巖在1.40—2.50。顯熱通量均值以喀斯特山峰裸巖和建筑/道路高于地面裸土,潛熱通量均值以建筑/道路和地面裸土高于喀斯特山峰裸巖,原因是建筑/道路周邊常種植樹木等綠化植被,建筑/道路實際上是包含少量比例綠化植被的混合像元,植被蒸騰產生潛熱;研究區中地面裸土面積不多,主要是建筑工地開挖土方,常用澆水降塵措施減少揚塵對周邊環境的影響,裸土中具有一定的含水量,土壤蒸發產生潛熱;喀斯特山峰裸巖間土壤土層瘠薄、含水量低,植被的生長速度十分緩慢,故實際上是包含少量比例草本植被、低矮灌木的混合像元,植被蒸騰和土壤蒸發量較少。

表2 典型地物水熱通量均值 /(W/m2)

表3 典型地物波文比均值
研究區水熱通量隨時間的變化受政府城市新開發區建設和綠化建設導致的地表覆蓋類型變化的綜合影響。由表2可以看出研究區波文比在1994年最高,達到1.62,原因是1991年國務院批準建設桂林國家級高新技術產業開發區,1994年廣西政府批準建設桂林西城經濟開發區,帶動了桂林城區的快速擴展。
然而桂林作為旅游城市,新開發區的綠化建設也同時受到政府的重視,常使用大樹移植方式快速提高新開發區的植被覆蓋率,加上水熱條件較好,樟樹、榕樹等本土樹種生長速度較快,1999年桂林被命名“全國園林綠化先進城市”,2000年研究區波文比下降到20年來的最低點,為1.24。
2003年建設高新區鐵山工業園和信息產業園;2005年建設高新區英才科技園;為保護漓江水環境,2007年8月臨桂新區獲批建設計劃到2015年完成13.96 km2建城區目標;2010年開始建設臨桂新區秧塘工業園(占地10 km2),新開發區的顯熱較高,而潛熱較低,2000—2015年研究區波文比呈現逐漸升高的趨勢,2015年上升到1.51,接近1994年數值。
分別將顯熱、潛熱通量的數值范圍從低到高的40%、30%、40%定義為低值區、中值區、高值區,統計各區間像元比例,結果見表4、表5。研究區顯熱高值區和潛熱低值區主要為植被覆蓋率非常低的較大面積的新建高密度建筑區和喀斯特山峰裸巖,比例低于10.0%;顯熱低值區和潛熱高值區主要為水體(漓江、桂湖、榕湖、杉湖等)、茂密植被覆蓋的喀斯特山峰、城市公園林地等,受到人為保護,比例較高,顯熱低值區比例在35.0%—50.0%,而潛熱高值區比例在35.0%—55.0%。喀斯特城市擴展過程出現的顯熱高值區比例以1994年最高,為10.0%,2000年下降到5.4%,之后呈逐漸上升到2010年的9.4%,但2015年下降到7.1%,而潛熱低值區比例的變化趨勢與顯熱高值區比例基本相同。顯熱高值區和潛熱低值區比例的變化導致顯熱中、低值區比例和潛熱中、高值區比例的變化,城市擴展將顯熱和潛熱中值區的低密度建筑區,顯熱低值區和潛熱高值區的農田、林地、池塘等,轉變為顯熱高值區和潛熱低值區的新建高密度建筑區,部分新建高密度建筑區經過若干年的綠化建設后植被覆蓋率增加可以轉變為顯熱和潛熱中值區。

表4 地表顯熱通量分布比例/%

表5 地表潛熱通量分布比例/%
城市地表水熱通量受地物覆蓋變化、氣候條件等影響,植被覆蓋度變化是最重要的因素之一,將其劃分為10等份,統計各等份水熱通量均值,結果見圖3。統計植被覆蓋度增加0.1時水熱通量變化的均值,結果見表6、表7。植被覆蓋度與顯熱通量呈明顯負相關關系,而與潛熱通量呈明顯正相關關系。在不同植被覆蓋度范圍,水熱通量隨植被覆蓋度的變化程度存在差異,植被覆蓋度在0.0—0.1范圍極稀疏和在0.8—1.0范圍極茂密時其變化對水熱通量的影響相對較弱,植被覆蓋度增加0.1,顯熱降低、潛熱升高4—10 W/m2;而植被覆蓋度在0.1—0.8范圍時其變化對水熱通量的影響相對較強,水熱通量與植被覆蓋度近似呈線性關系,植被覆蓋度增加0.1,顯熱降低8—27 W/m2,而潛熱升高8—24 W/m2。

圖3 研究區顯熱與潛熱通量均值隨植被覆蓋度變化Fig.3 Average sensible and latent heat fluxes over the study area against vegetation fraction

表6 地表顯熱通量隨植被覆蓋度增加0.1的變化 /(W/m2)
使用METRIC模型反演桂林喀斯特城市地表潛熱、顯熱通量時,由于實驗條件限制,未能對不同類型下墊面的土壤熱通量、熱能儲存、地表粗糙度、零平面位移等進行實驗觀測,而參考已有的研究進行估算,將導致一定的反演誤差;模型驗證時盡管潛熱、顯熱通量遙感模型像元反演值與點上實驗觀測值由于所代表的空間尺度不同,存在驗證誤差。但由于反演誤差主要影響地表通量數值,對地表通量空間整體分布趨勢、顯熱通量與潛熱通量比值的波文比、水熱通量與植被覆蓋度關系的影響較小。

表7 地表潛熱通量隨植被覆蓋度增加0.1的變化 /(W/m2)
桂林喀斯特城市城區中鑲嵌著一些拔地而起的喀斯特孤峰,常被設置為公園,喀斯特山峰上的植被受人為保護而生長狀況較好,植被覆蓋度普遍相對較高,植被覆蓋阻擋、反射太陽對石灰巖山體的部分直接輻射,茂密的喀斯特植被葉面整體的蒸騰作用較大而吸收較多熱量,潛熱較高而顯熱較低,影響其上方及周邊空氣,形成局地小氣候環境,對調節喀斯特城市微氣候環境具有重要作用。
喀斯特山峰裸巖一般是采石場,其上以灌叢和闊葉林為主的喀斯特植被遭受人為破壞,形成大面積的喀斯特山峰基巖裸露。由于裸巖間的土壤土層瘠薄,植被生長緩慢,很難在短時間內恢復喀斯特山峰植被覆蓋,容易造成土壤侵蝕,形成石漠化現象。而熱容量較大的大面積裸露石灰巖在太陽直接輻射下吸收、貯存大量的熱量,顯熱較高而潛熱較低,對城市熱環境產生不透水地表相似的影響效應。因此鑲嵌于城區中喀斯特山峰的植被保護對喀斯特城市熱環境的改善至關重要。
植被覆蓋率非常低的較大面積的新建高密度建筑區種植幼樹等少量的植被,植被覆蓋率非常高的較大面積的林地繼續增加植被覆蓋,對降低顯熱和提高潛熱的效果不顯著,對環境狀況的改變影響不大。
在太陽輻射下喀斯特山峰的水熱通量受地形因素的影響,但桂林主城區的喀斯特山峰距地面的相對高度一般不超過200 m,高程的影響較小。
論文的創新點在于根據我國西南地區桂林喀斯特城市城區存在喀斯特山峰復雜地形的實際情況改進METRIC模型,使模型適用于我國西南地區喀斯特城市地表水熱通量的反演,使用地面觀測數據進行驗證,分析喀斯特城市地表水熱通量的時空變化規律,迄今未見相關研究論文發表。
(1)水體和植被覆蓋地表的潛熱通量遠高于顯熱通量,而建筑/道路、裸土和喀斯特山峰裸巖的顯熱高于潛熱;潛熱通量從高到低依次為水體、喀斯特山峰陽坡植被、地面植被、喀斯特山峰陰坡植被、建筑/道路和裸土、喀斯特山峰裸巖,顯熱通量從高到低依次為喀斯特山峰裸巖和建筑/道路、裸土、喀斯特山峰陽坡植被、地面植被、喀斯特山峰陰坡植被、水體。
(2)喀斯特城市水熱通量隨時間的變化受政府城市新開發區建設和綠化建設導致的地表覆蓋類型變化的影響。研究區波文比在1994年最高,達到1.62,2000年下降到20年來的最低點,為1.24,之后逐漸升高到2015年的1.51,接近1994年數值。
(3)喀斯特城市擴展過程出現的顯熱高值區和潛熱低值區比例低于10.0%,其變化引發顯熱中低值區和潛熱中高值區比例的變化。顯熱高值區比例在1994年最高,為10.0%,2000年下降到5.4%,之后至2010年逐漸上升到9.4%,但2015年下降到7.1%,潛熱低值區比例的變化趨勢與顯熱高值區比例基本相同。
(4)植被覆蓋度在不同范圍對水熱通量的影響程度存在差異,植被覆蓋度在0.0—0.1范圍很低和在0.8—1.0范圍很高時其增加0.1,顯熱通量降低、潛熱通量升高4—10 W/m2,影響相對較弱;而植被覆蓋度在0.1—0.8范圍時其增加0.1,顯熱通量降低8—27 W/m2,而潛熱通量升高8—24 W/m2,影響相對更顯著。
(5)METRIC遙感模型提供了定量分析喀斯特城市水熱通量時空變化的快速、經濟可行的方法,并且取得較好的結果。
致謝:本研究得到了“廣西八桂學者創新研究團隊”、“廣西區域多源信息集成與智能處理協同創新中心”的支持,特此致謝。