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

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

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

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

中圖分類號 P49 文獻(xiàn)標(biāo)識碼 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

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

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

1 材料與方法

1.1 材料

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

1.2 方法

在形態(tài)學(xué)方法中，一般認(rèn)為Raupach[17]方法精度較高[18]，因此選擇Raupach方法進(jìn)行海南島橡膠林零平面位移和粗糙度估算。通過提取橡膠林葉面積指數(shù)、相關(guān)參數(shù)計(jì)算，代入Raupach 的形態(tài)學(xué)模型，計(jì)算海南島橡膠林零平面位移和粗糙度。具體技術(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)力模型和全球氣候、水文、生物地球化學(xué)和生態(tài)學(xué)模型的關(guān)鍵輸入?yún)?shù)。葉面積指數(shù)的計(jì)算方法有：直接方法（葉面積的測定、描形稱重法、儀器測定法）；間接方法（點(diǎn)接觸法、消光系數(shù)法、經(jīng)驗(yàn)公式法、遙感方法）。本文通過橡膠林葉面積指數(shù)（LAI）和歸一化植被指數(shù)（NDVI）統(tǒng)計(jì)關(guān)系，計(jì)算得到如下回歸方程，進(jìn)行天然橡膠葉面積指數(shù)的反演。

LAI=7.82NDVI+0.781 （1）

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

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

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

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

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

模型中其他參數(shù)參考相關(guān)文獻(xià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)了此時(shí)橡膠冠層物候的差異。其主要原因是該時(shí)段海南島橡膠普遍進(jìn)入第一蓬葉展葉期，但由于不同區(qū)域的溫度、降水、光照、風(fēng)速等的差異，對此時(shí)橡膠葉的生長有著不同的影響，一般以東南和南部地區(qū)最早出現(xiàn)，最后是西北部和中部地區(qū)[20]。

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

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

3 討論與結(jié)論

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

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

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

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

參考文獻(xiàn)

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

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

[3]鐘 中， 韓士杰. 長白山闊葉紅松林冠層空氣動力學(xué)參數(shù)的計(jì)算[J]. 南京大學(xué)學(xué)報(bào)（自然科學(xué)）， 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] 楊阿強(qiáng)， 孫國清， 盧立新， 等.基于MODIS資料的中國東部時(shí)間序列空氣動力學(xué)粗糙度和零平面位移高度估算[J]. 氣象科學(xué)， 2011， 31（4）： 516-524.

[10] 趙曉松， 關(guān)德新， 吳家兵， 等.長白山闊葉紅松林的零平面位移和粗糙度[J]. 生態(tài)學(xué)雜志， 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è)科學(xué)， 2010， 37（1）： 168-170.

[16] 唐群鋒， 郭澎濤， 劉志崴， 等. 基于FMT-AHP的海南農(nóng)墾花崗巖類多雨區(qū)橡膠園地力評價(jià)[J]. 生態(tài)學(xué)報(bào)， 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] 連士華. 橡膠樹風(fēng)害成因間題的探討[J]. 熱帶作物學(xué)報(bào)， 1984， 5（1）： 59-72.

[23] 劉和平， 劉樹華， 朱廷曜， 等.森林冠層空氣動力學(xué)參數(shù)的確定[J]. 北京大學(xué)學(xué)報(bào)（自然科學(xué)版）， 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.