劉少軍 易雪 張京紅 蔡大鑫 田光輝 李偉光

摘 要 在橡膠風害評估中，橡膠林零平面位移和粗糙度是建立定量評估模型的重要動力學參數(shù)。以MODIS遙感數(shù)據(jù)為信息源，計算橡膠林的葉面積指數(shù)，并利用形態(tài)學Raupach方法，估算橡膠林歸一化零平面位移和粗糙度。該方法充分發(fā)揮了遙感大面積獲取數(shù)據(jù)的優(yōu)勢，同時解決了傳統(tǒng)氣象學方法獲取橡膠林零平面位移和粗糙度需要大量長期的野外觀測資料的限制。結(jié)果表明：橡膠林葉面積指數(shù)與歸一化零平面位移和粗糙度分別存在正相關(guān)和負相關(guān)；估算結(jié)果精度能滿足橡膠風害評估模型參數(shù)的需求。研究結(jié)果可為開展大范圍內(nèi)橡膠林風害評估提供技術(shù)支撐。

關(guān)鍵詞 橡膠林；零平面位移；粗糙度；海南島

中圖分類號 P49 文獻標識碼 A

Abstract Zero-plane displacement and roughness length of rubber are crucial parameters in the typhoon impact assessment model on rubber plant. Leaf area index was derived based on the MODIS remote sensing data. Zero-plane displacement and roughness length of rubber were derived by the morphology model. This method not only made full use of the advantage of remote sensing， but also solved the traditional meteorology method to get zero-plane displacement and roughness depending on a long-term field observation data. The results showed that leaf area indices were positively and negatively correlated with d/h and z0/h. The accuracy of result could meet the need of rubber wind damage assessment model parameters. The results provided technical support for carrying out large scale rubber wind damage assessment.

Key words Rubber plant；Zero-Plane Displacement；Roughness Length；Hainan Island

doi 10.3969/j.issn.1000-2561.2016.10.029

零平面位移和粗糙度是描述下墊面空氣動力學特征的重要參數(shù)[1]。零平面位移為下墊面開始吸收動量的高度，地表粗糙度是下墊面風速為零的高度[2]，但在不同天氣背景下， 零平面位移和粗糙度長度并不是常數(shù)[3]。在海南橡膠林臺風災害動力學評估模型中，橡膠樹斷、倒條件的判識，必須解決大范圍、不同季節(jié)、不同地形條件下橡膠林零平面位移和粗糙度。由于早期的研究對零平面位移和粗糙度的理解不確切，一般通過冠層高度和葉面積指數(shù)等分別建立零平面位移和粗糙度的經(jīng)驗公式進行估算[4-7]。零平面位移和粗糙度是風速和大氣層結(jié)的函數(shù)，而不僅僅是由冠層幾何高度來確定[3，8]。目前，零平面位移和粗糙度的計算一般存在氣象學和形態(tài)學兩種方法[9]。氣象學方法常采用三維超聲風溫儀，觀測多層的風速觀測資料，然后利用風速廓線公式通過迭代求解[1，10]。如張學仕等[1]利用應用溫度方差方法計算了江蘇句容下蜀次生櫟林空氣動力學參數(shù)；鐘中等[3]根據(jù)空氣質(zhì)量守恒原理， 對長白山天然闊葉紅松林冠層空氣動力學參數(shù)進行了數(shù)值計算；趙曉松等[10]利用廓線法中的牛頓迭代法計算了長白山闊葉紅松林的零平面位移和粗糙度；Lo[11]根據(jù)冠層上下多層風、溫觀測資料計算高大森林冠層零平面位移和粗糙度的方法。而形體學方法主要是采用粗糙元高度、粗糙元頂部面積、粗糙元迎風面積等計算零平面位移和粗糙度[12]。如楊阿強等[9]采用粗糙元形態(tài)學方法，計算中國東部不同季節(jié)空氣動力學零平面位移和粗糙度；Macdonald等[13]根據(jù)粗糙元迎風面積和頂部面積的方法計算零平面位移和粗糙度。

形態(tài)學方法比氣象學方法更具有優(yōu)勢。氣象學方法物理意義更為確切，不同參數(shù)具有不同的物理含義，但計算過程中需要在穩(wěn)定中性大氣條件下準確測量氣象數(shù)據(jù)，推導橡膠林零平面位移和粗糙度對于大范圍長時間序列的計算需要大量長期的野外觀測，成本較高；而形態(tài)學方法是根據(jù)物理基礎(chǔ)經(jīng)過簡化綜合經(jīng)驗系數(shù)得到的，精度需要不斷提高，但其方法可以采用遙感數(shù)據(jù)作為數(shù)據(jù)源，而且可以進行大范圍內(nèi)參數(shù)的反演[9]。目前，關(guān)于采用遙感和形態(tài)學方法估算橡膠林零平面位移和粗糙度的研究未見報道。因此，本研究采用遙感數(shù)據(jù)源，嘗試利用形態(tài)學模型的方法，估算海南島橡膠林零平面位移和粗糙度，為快速大范圍獲取橡膠林動力參數(shù)提供了新的技術(shù)途徑。

1 材料與方法

1.1 材料

海南島MODIS遙感數(shù)據(jù)、TM遙感數(shù)據(jù)提取海南天然橡膠空間分布圖（具體方法參考文獻[14-15]）、2011年海南天然橡膠普查數(shù)據(jù)、橡膠生產(chǎn)分布以農(nóng)場為單位進行統(tǒng)計，數(shù)據(jù)來源于文獻[16]。省界和縣界行政區(qū)劃圖來源于國家基礎(chǔ)地理信息網(wǎng)站提供的1 ∶ 400萬基礎(chǔ)地理信息數(shù)據(jù)（http：//ngcc.sbsm.gov.cn/）。

1.2 方法

在形態(tài)學方法中，一般認為Raupach[17]方法精度較高[18]，因此選擇Raupach方法進行海南島橡膠林零平面位移和粗糙度估算。通過提取橡膠林葉面積指數(shù)、相關(guān)參數(shù)計算，代入Raupach 的形態(tài)學模型，計算海南島橡膠林零平面位移和粗糙度。具體技術(shù)流程見圖1。

1.2.1 橡膠葉面積指數(shù)提取 葉面積指數(shù)（Leaf Area Index，LAI）又叫葉面積系數(shù)，是指單位土地面積上植物葉片總面積占土地面積的倍數(shù)。即：葉面積指數(shù)=葉片總面積/土地面積。葉面積指數(shù)（LAI）作為植被冠層的重要結(jié)構(gòu)參數(shù)之一，是許多生態(tài)系統(tǒng)生產(chǎn)力模型和全球氣候、水文、生物地球化學和生態(tài)學模型的關(guān)鍵輸入?yún)?shù)。葉面積指數(shù)的計算方法有：直接方法（葉面積的測定、描形稱重法、儀器測定法）；間接方法（點接觸法、消光系數(shù)法、經(jīng)驗公式法、遙感方法）。本文通過橡膠林葉面積指數(shù)（LAI）和歸一化植被指數(shù)（NDVI）統(tǒng)計關(guān)系，計算得到如下回歸方程，進行天然橡膠葉面積指數(shù)的反演。

LAI=7.82NDVI+0.781 （1）

其中，歸一化植被指數(shù)的計算公式如下：

式中，NDVI表示歸一化植被指數(shù)，Rnir表示近紅外波段反照率，Rred表示可見光波段（紅光）反照率。

1.2.2 零平面位移和粗糙度估算 Raupach的形態(tài)學公式表示如下[17]：

d為零平面位移，ZO為粗糙度，h為粗糙元高度，λ為單位地面粗糙元迎風面積，ψh為粗糙層影響函數(shù)，b為冠層寬度，Cd為經(jīng)驗系數(shù)。β為粗糙元阻力系數(shù)與地表阻力系數(shù)之比，γ為距地面高度z 處的水平風速與摩擦速度之比。其中單位地面粗糙元迎風面積λ可以為：

A表示為單位面積上植被對風的阻擋面積，值與葉面積指數(shù)很接近[19]，因此采用遙感反演的葉面積指數(shù)LAI表示單位面積上植被對風的阻擋面積。

模型中其他參數(shù)參考相關(guān)文獻獲取[4-8，10-13]。

2 結(jié)果與分析

本次研究選擇了2015年3月中旬云量相對較少的MODIS遙感數(shù)據(jù)，并結(jié)合海南島橡膠分布[14-16]，分別提取了海南島橡膠林的葉面積指數(shù)、歸一化零平面位移、歸一化粗糙度分布圖。

結(jié)果顯示：海南島橡膠林的葉面積指數(shù)在1.0～8.6之間，高值區(qū)分布在海南島的中部和東南部區(qū)域，低值區(qū)域分布在海南島的北部和西部區(qū)域（圖2），遙感葉面積指數(shù)的差異體現(xiàn)了此時橡膠冠層物候的差異。其主要原因是該時段海南島橡膠普遍進入第一蓬葉展葉期，但由于不同區(qū)域的溫度、降水、光照、風速等的差異，對此時橡膠葉的生長有著不同的影響，一般以東南和南部地區(qū)最早出現(xiàn)，最后是西北部和中部地區(qū)[20]。

海南橡膠林歸一化零平面位移在0.5～0.82之間，歸一化零平面位移最大值為0.82，最小值為0.5，平均值為0.75，標準差0.046。高值區(qū)分布在海南島的中部和東南部區(qū)域，低值區(qū)域分布在海南島的北部和西部區(qū)域（圖3）。由上可知，橡膠林的歸一化平均零平面位移與葉面積指數(shù)變化趨勢相一致，高值和低值區(qū)分布基本一致。通過ArcGIS軟件的Spatial analyst tools-Multivariate-Band collection statistics功能模塊，計算兩者的相關(guān)性為0.85。

海南橡膠林歸一化粗糙度在0.056～0.10之間，歸一化粗糙度最大值為0.10，最小值為0.056，平均值為0.073，標準差0.012。低值區(qū)分布在海南島的中部區(qū)域，高值區(qū)域分布在海南島的北部和西南部區(qū)域（圖4）。歸一化粗糙度的變化趨勢與橡膠林葉面積指數(shù)高值和低值區(qū)的分布相反。通過ArcGIS軟件計算兩者的相關(guān)性為-0.84。其主要原因可能是當橡膠林葉面積指數(shù)高值時，橡膠樹葉片相對繁茂，使氣流從作物頂部通過，從而減小了氣流同葉片的摩擦降低了空氣動力學粗糙度。

3 討論與結(jié)論

本文在借鑒前人研究方法的基礎(chǔ)上，結(jié)合資料可獲得性及滿足大范圍獲取橡膠林零平面位移和粗糙度的需求，嘗試利用MODIS遙感數(shù)據(jù)和Raupach的形態(tài)學模型，采用的計算方法是建立在嚴格的理論推導基礎(chǔ)上，同時兼顧經(jīng)濟性原則，希望突破以往有關(guān)零平面位移和粗糙度取值的經(jīng)驗性約束，為橡膠林動力學風害評估模型提供動力參數(shù)。

橡膠林零平面位移和粗糙度實際上是隨大氣穩(wěn)定情況發(fā)生改變的，可能同一天都會相差很大，不同大氣層結(jié)狀態(tài)，零平面位移和粗糙度均會有所差異，同時橡膠林平面位移和粗糙度隨季節(jié)會發(fā)生變化。遙感反演的海南橡膠林d/h和z0/h與橡膠林葉面積指數(shù)分別存在正相關(guān)和負相關(guān)的關(guān)系。這與Shaw等[21]、趙曉松等[10]利用氣象觀測塔實測數(shù)據(jù)計算的規(guī)律一致。

本研究由于缺乏與遙感數(shù)據(jù)相對應時間段的通量塔觀測資料與之進行對比驗證，導致無法進行遙感反演精度的直接驗證。但橡膠樹屬于闊葉林，相關(guān)參數(shù)的驗證僅與前人[9-10，22-24]的森林的結(jié)果進行對比分析，根據(jù)2015年3月15日MODIS數(shù)據(jù)反演的海南島橡膠林歸一化零平面位移和粗糙度均在同一量級范圍內(nèi)，反演結(jié)果整體上反映了橡膠林歸一化零平面位移和粗糙度在空間分布的差異，能滿足橡膠林風害模型的參數(shù)要求，但具體的誤差仍需要與后期建立橡膠林觀測數(shù)據(jù)進行對比分析，不斷改進Raupach的形態(tài)學模型相關(guān)參數(shù)，以提高遙感反演橡膠林零平面位移和粗糙度的精度。

由于本文僅停留在利用Raupach方法，探索性地估算橡膠林零平面位移和粗糙度，而且相關(guān)參數(shù)主要參考了經(jīng)驗值；同時由于橡膠林的面積的來源主要參考文獻[14-16]的成果，沒能及時更新，難免存在誤差和不足。

參考文獻

[1] 張學仕， 蔣 琰， 薛建輝， 等.江蘇句容下蜀次生櫟林的空氣動力學參數(shù)研究[J]. 南京林業(yè)大學學報：自然科學版， 2010， 34（6）： 61-65.

[2] 蔡 福， 周廣勝， 明惠青， 等. 玉米農(nóng)田空氣動力學參數(shù)動態(tài)及其與影響因子的關(guān)系[J]. 生態(tài)學報， 2013， 33（17）： 5 339-5 352.

[3]鐘 中， 韓士杰. 長白山闊葉紅松林冠層空氣動力學參數(shù)的計算[J]. 南京大學學報（自然科學）， 2002， 38（4）： 565-571.

[4] Szeicz G， Endrodi G， Tajchman S. Aerodynamic and surface factors in evaporation[J]. Water Resources Research， 1969， 5（2）： 380-394 .

[5] Stanhill G. A Simple instrument for field measurement of turbulent diffusion flux[J]. Journal of Applied Meteorology， 1969（8）： 509-513 .

[6] Zoumakis N M. Determination of the mean wind speed and momentum diffusivity profiles above tall vegetation and forest canopies using a mass conservation assumption[J]. Journal of Applied Meteorology， 1994， 33： 295-302.

[7] Raupach M R. Simplified expressions forvegetation roughness length and zeroplane displacement as functions of canopy height and area index[J]. Boundary-Layer Meteorology， 1994， 71（1-2）： 211-216.

[8] Lesnik G E.Estimation of displacement layer thickness and roughness parameter from measurements above a forest[J]. Izvestiya Atmospheric And Oceanic Physics， 1973， 9（1）： 88-90. [9] 楊阿強， 孫國清， 盧立新， 等.基于MODIS資料的中國東部時間序列空氣動力學粗糙度和零平面位移高度估算[J]. 氣象科學， 2011， 31（4）： 516-524.

[10] 趙曉松， 關(guān)德新， 吳家兵， 等.長白山闊葉紅松林的零平面位移和粗糙度[J]. 生態(tài)學雜志， 2004， 23（5）： 84-88.

[11] Lo A K. Determination of zero-plane displacement and roughness length of a forest canopy using profiles of limited height[J]. Boundary-Layer Meteorology， 1995， 75： 81-402.

[12] RaupachM R， Antonia R A， Rajagopalan S. Rough-wall turbulent boundary layers[J]. Applied Mechanics Reviews， 1991， 44（1）： 1-25.

[13] Macdonald R W， Griffiths R F， Hall D J. An improved method for estimation of surface roughness of obstacle arrays[J]. Atmospheric Environment， 1998， 32（11）： 1 857-1 864.

[14] 張京紅， 陶忠良， 劉少軍， 等.基于TM影像的海南島橡膠種植面積信息提取[J]. 2010， 31（4）： 661-665.

[15]劉少軍， 張京紅， 何政偉. 基于面向?qū)ο蟮南鹉z分布面積估算研究[J]. 廣東農(nóng)業(yè)科學， 2010， 37（1）： 168-170.

[16] 唐群鋒， 郭澎濤， 劉志崴， 等. 基于FMT-AHP的海南農(nóng)墾花崗巖類多雨區(qū)橡膠園地力評價[J]. 生態(tài)學報， 2014， 34（15）： 4 435-4 445.

[17] Jasinski M F， Borak J， Crago R. Bulk surface momentum parameters for satellite-derived vegetation elds[J]. Agricultural and Forest Meteorology， 2005， 133： 55-68.

[18] Grimmond C S B， Oke T R.Aerodynamic properties of urban areas derived from analysis of surface form[J]. Journal of Applied Meteorology， 1999， 38： 1 262-1 292.

[19] ZENG Xubin， Shaikh Muhammad， Dai Yongjiu， et al. Coupling of the common land model to the NCAR community climate model[J]. Journal of Climate， 2001， 15： 1 832-1 854.

[20] 陳小敏， 陳匯林， 李偉光， 等. 海南島天然橡膠林春季物候期的遙感監(jiān)測[J]. 中國農(nóng)業(yè)氣象， 2016， 37（1）： 111-116.

[21] Shaw R H， Pereira A R. Aerodynamic roughness of a plant canopy： a numerical experiment[J]. Agricultural Meteorology， 1982， 26： 51-65.

[22] 連士華. 橡膠樹風害成因間題的探討[J]. 熱帶作物學報， 1984， 5（1）： 59-72.

[23] 劉和平， 劉樹華， 朱廷曜， 等.森林冠層空氣動力學參數(shù)的確定[J]. 北京大學學報（自然科學版）， 1997， 33（4）： 522-528.

[24] Borak Jordan S， Jasinski Michael F， Crago Richard D. Time seriesvegetation aerodynamic roughness fields estimated from MODISobservations[J]. Agricultural and Forest Meteorology， 2005， 135： 252-268.