基于遙感的海南島橡膠林零平面位移和粗糙度估算

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

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

關鍵詞 橡膠林；零平面位移；粗糙度；海南島

中圖分類號 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

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

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

1 材料與方法

1.1 材料

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

1.2 方法

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

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

LAI=7.82NDVI+0.781 （1）

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

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

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

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

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

模型中其他參數參考相關文獻獲取[4-8，10-13]。

2 結果與分析

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

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

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

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

3 討論與結論

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

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

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

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

參考文獻

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

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

[3]鐘 中， 韓士杰. 長白山闊葉紅松林冠層空氣動力學參數的計算[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] 趙曉松， 關德新， 吳家兵， 等.長白山闊葉紅松林的零平面位移和粗糙度[J]. 生態學雜志， 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]劉少軍， 張京紅， 何政偉. 基于面向對象的橡膠分布面積估算研究[J]. 廣東農業科學， 2010， 37（1）： 168-170.

[16] 唐群鋒， 郭澎濤， 劉志崴， 等. 基于FMT-AHP的海南農墾花崗巖類多雨區橡膠園地力評價[J]. 生態學報， 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] 陳小敏， 陳匯林， 李偉光， 等. 海南島天然橡膠林春季物候期的遙感監測[J]. 中國農業氣象， 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] 劉和平， 劉樹華， 朱廷曜， 等.森林冠層空氣動力學參數的確定[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.