褶皺山區(qū)NDVI與地形因子的多尺度耦合關(guān)系研究

王慧玲，熊亞蘭，陳小鳳，劉 穎

（四川農(nóng)業(yè)大學(xué)資源學(xué)院，成都 611130）

歸一化植被指數(shù)，即NDVI（normalized difference vegetation index）數(shù)值，其動(dòng)態(tài)主要受地形特征和人為因素的影響[1]。NDVI可以削弱大氣層和地形陰影的影響，對(duì)地表植被的覆蓋程度非常敏感，是目前檢測(cè)和指示植被覆蓋狀況和動(dòng)態(tài)的常用指標(biāo)之一[1]。在山區(qū)，地形因素顯得尤為重要。其通過(guò)影響太陽(yáng)輻射與降水的空間再分配過(guò)程，對(duì)植被分布格局的塑造起著不可或缺的作用[2]。利用地形差異來(lái)揭示植被空間分異規(guī)律，是自然地理、生態(tài)學(xué)等領(lǐng)域的重點(diǎn)研究方向之一。如仇寬彪等研究了基于NDVI太湖流域片植被覆蓋變化[3]。鄧元杰等研究了長(zhǎng)江流域中上游植被NDVI時(shí)空變化及其地形分異效應(yīng)[4]。劉梁美子等研究了黔桂喀斯特山區(qū)NDVI與地形之間存在的相關(guān)性[5]。而三峽庫(kù)區(qū)重慶段的渝東北地區(qū)是長(zhǎng)江流域重要生態(tài)屏障和長(zhǎng)江上游特色經(jīng)濟(jì)走廊[6]，但該區(qū)域內(nèi)褶皺山脈廣布、山高坡陡、土層淺薄，農(nóng)業(yè)生產(chǎn)較為落后，因此生態(tài)環(huán)境保護(hù)任務(wù)重，發(fā)展農(nóng)業(yè)生產(chǎn)的壓力大，有必要對(duì)該區(qū)植被覆蓋情況與地形因子的關(guān)系進(jìn)行全面深入和多尺度的研究。

地形要素作用于一系列的地理過(guò)程，從而對(duì)區(qū)域植被的格局產(chǎn)生影響，而簡(jiǎn)單的相關(guān)性分析不足以揭示這其中的影響關(guān)系。目前對(duì)于植被覆蓋指數(shù)和地形因素的相關(guān)影響程度伴隨著小波變換分析法的發(fā)展，已由靜態(tài)、單一尺度轉(zhuǎn)向動(dòng)態(tài)、多尺度劃分[7-8]。小波分析被稱為“數(shù)學(xué)顯微鏡”，可以揭示自然或生態(tài)因子的多尺度格局，并且不要求數(shù)據(jù)符合二階平穩(wěn)假設(shè)[9-11]。過(guò)去對(duì)NDVI與地形因子多尺度關(guān)系的研究主要是采樣一維小波變換方法對(duì)研究區(qū)多條樣帶上NDVI與地形因子的多尺度相關(guān)性或小波相干性進(jìn)行研究[12-15]。如李雙成等通過(guò)小波變換揭示了西藏高原生態(tài)系統(tǒng)的NDVI與地形因子的多尺度空間格局和相關(guān)性[14-15]。李小梅等利用小波多尺度分析工具分析了福州市NDVI和地形因子等生態(tài)環(huán)境參數(shù)的尺度—格局特征以及生態(tài)環(huán)境參數(shù)相關(guān)關(guān)系的尺度響應(yīng)[9]。邱炳文、鄒偉成等利用離散小波變換對(duì)武夷山自然保護(hù)區(qū)的NDVI與地形因子進(jìn)行了多尺度空間相關(guān)分析[16-17]。

樣帶的選擇是一維小波變換研究中的關(guān)鍵，而克服樣帶代表性等的限制則需要運(yùn)用到二維小波變換的方法。二維小波變換能將柵格數(shù)據(jù)進(jìn)一步分解為低頻和高頻部分，高頻部分又可劃分為水平、垂直和對(duì)角3個(gè)部分。低頻系數(shù)為實(shí)測(cè)值經(jīng)過(guò)多尺度分解后的近似信息，可顯示出實(shí)測(cè)值通過(guò)小波變換剩余的基本輪廓信息，若進(jìn)一步分解，即分解尺度擴(kuò)大。其更能濾去細(xì)節(jié)信息，趨向清晰化和簡(jiǎn)單化，從而較為直觀地反映出研究對(duì)象的總體變化規(guī)律[9-11]。高頻系數(shù)即實(shí)測(cè)值經(jīng)過(guò)多尺度分解后的細(xì)節(jié)信息，反映出該變量在不同空間位置的波動(dòng)變化情況，某個(gè)區(qū)域的波動(dòng)變化越大，則反映了該區(qū)域研究對(duì)象變化越復(fù)雜，受到其他因素干擾的強(qiáng)度越大。故二維小波變換較一維小波變換更能從多方位、多角度、深入貼合地反映出植被與地形因子在大空間和小空間的耦合關(guān)系。因此本研究擬采用二維小波變換對(duì)褶皺山區(qū)NDVI與地形因子的多尺度耦合關(guān)系進(jìn)行研究，研究結(jié)果在一定程度上能夠?yàn)樯絽^(qū)植被保護(hù)、生態(tài)環(huán)境恢復(fù)、水土流失治理和提高農(nóng)業(yè)生產(chǎn)規(guī)劃提供相關(guān)數(shù)據(jù)支撐。

1 材料和方法

1.1 研究區(qū)概況

研究區(qū)（108°12′13″～108°54′27″E，30°48′33″～31°10′43″N）位于長(zhǎng)江上游，川東平行嶺谷東北部，面積約2 700 km2（圖1）。該區(qū)域?qū)賮啛釒Ъ撅L(fēng)氣候，立體氣候顯著，四季分明，光照充足，雨量豐沛，降雨量1 200 mm左右，年平均氣溫18℃左右。全區(qū)為兩列褶皺山脈由北東—南西走向轉(zhuǎn)為西東走向的轉(zhuǎn)折地帶，地勢(shì)西高東低，海拔高度為78～1 350 m；出露地層的地質(zhì)年代多見(jiàn)于中生代三疊紀(jì)和侏羅紀(jì)，侏羅紀(jì)地層分布最廣，三疊紀(jì)次之；山地和丘陵廣布；土壤類型以黃壤、紫色土、石灰土和水稻土等為主[18]。海拔500 m以下的丘陵區(qū)主要為農(nóng)耕區(qū)，海拔500～1 350 m的低山和中山下部，自然林被主要為常綠闊葉與落葉闊葉混交林[2]。

圖1 NDVI和地形因子的多尺度空間分布圖Figure 1 Multi-scale spatial patterns of NDVI and topographic factors

圖1 研究區(qū)位圖Figure 1 The location of study area

1.2 研究方法

1.2.1 數(shù)據(jù)來(lái)源與預(yù)處理

研究區(qū)年降雨充沛，植被整體覆蓋情況較好，在7月份和8月份植被覆蓋度達(dá)到最大值，但其空間分布的差異較小，因此為了更好地反映地形對(duì)植被覆蓋度的影響，本研究選擇的是一年中植被覆蓋度最低、空間差異最大的1月份（2016年）的植被覆蓋度數(shù)據(jù)。MODIS合成的標(biāo)準(zhǔn)NDVI產(chǎn)品MODND1M為月內(nèi)每天最大值，空間分辨率為500 m，坐標(biāo)系為WGS84。

在山地生態(tài)學(xué)研究中，高程、坡度和坡向是地形的3個(gè)主要屬性，也是決定植被生境其他要素分異（如土壤、小氣候、水文等）的主導(dǎo)因子。本研究利用數(shù)字高程模型（digital elevation model，DEM）來(lái)計(jì)算并用ArcGIS提取得到高程、坡度和坡向。數(shù)據(jù)空間分辨率為30 m，坐標(biāo)系為WGS84。

NDVI和DEM數(shù)據(jù)均在地理空間數(shù)據(jù)云平臺(tái)下載，并在ArcGIS中將DEM重采樣為500 m。考慮到研究區(qū)地形以山地和丘陵為主，且重采樣前后空間分辨率差異較大，因此采用三次卷積插值法重采樣。

1.2.2 小波多尺度分析的基本原理

設(shè)任意一個(gè)采樣點(diǎn)坐標(biāo)為（x，y），點(diǎn)坐標(biāo)（x，y）變化時(shí)各要素的變化就形成了一個(gè)離散變化的二維信號(hào)（函數(shù)）f（x，y）。小波變換就是利用母小波（mother wavelet）與原函數(shù)內(nèi)積，從而達(dá)到分解原函數(shù)的目的。二維小波變換的公式如下[19]：

式中，WTf為小波系數(shù)，a為尺度因子，b1和b2分別為x和y方向的平移因子，ψ（x）為母小波，f（x，y）為待分析函數(shù)（網(wǎng)格實(shí)測(cè)值）。將母小波與待分析函數(shù)內(nèi)積后，可得到不同尺度下的小波系數(shù)并對(duì)它們進(jìn)行分解。本研究選取小波分析中常用的db6算法作為小波母函數(shù)，采用Mallat塔式算法。根據(jù)塔式算法，實(shí)測(cè)值與分解后小波系數(shù)的關(guān)系為[20-22]：

式中，j為分解的級(jí)數(shù)；A為低頻系數(shù)；H、V和D分別為水平、垂直和對(duì)角方向的高頻系數(shù)；G為多級(jí)分解尺度下3個(gè)方向高頻系數(shù)之和。在MATLAB中完成了小波低頻系數(shù)和高頻系數(shù)的提取。

根據(jù)小波變換的原理，最大分解尺度與柵格行（列）的數(shù)量密切相關(guān)，若行（列）的數(shù)量為2n，那最佳分解尺度為n-2[11]。本研究采樣點(diǎn)為128行、64列，因此可對(duì)NDVI和地形因子進(jìn)行四級(jí)分解，分解后每一級(jí)所對(duì)應(yīng)的網(wǎng)格大小分別為1 km×1 km、2 km×2 km、4 km×4 km和8 km×8 km。每一級(jí)小波分解后低頻系數(shù)和高頻系數(shù)對(duì)應(yīng)的空間尺度單元（樣點(diǎn)間距）會(huì)放大為上一級(jí)的2倍，即每級(jí)小波系數(shù)的個(gè)數(shù)會(huì)減少為上一級(jí)的1/2。通過(guò)二維小波重構(gòu)可以確保不同分解尺度下小波系數(shù)的個(gè)數(shù)一致。小波分解后的小波系數(shù)是沒(méi)有量綱的，但重構(gòu)后的小波系數(shù)是有量綱的。

2 結(jié)果與分析

2.1 NDVI與地形因子的基本統(tǒng)計(jì)特征分析

小波低頻系數(shù)反應(yīng)全局變化態(tài)勢(shì)，高頻系數(shù)表征出細(xì)節(jié)信息，由于本研究中柵格尺度較大，因此只對(duì)低頻系數(shù)進(jìn)行分析。NDVI和地形因子的基本統(tǒng)計(jì)特征見(jiàn)表1。0.5 km尺度，NDVI的取值為0.09～1.00，均值為0.50，變異系數(shù)為24%，說(shuō)明研究區(qū)植被覆蓋情況良好，但具有中等程度的變異。高程、坡度和坡向的變異系數(shù)分別為49%、53%和58%，均為中等程度的變異。從0.5～8 km尺度，各要素最大值與最小值的距離縮小，變異系數(shù)、峰度和偏度均逐漸降低，但均值保持不變，這表明在小波分解的過(guò)程中隨著尺度的增大，各要素低頻信息得到平滑。

表1 NDVI及地形因子的基本統(tǒng)計(jì)特征Table 1 Basic statistical characteristics of NDVI and terrain factors

2.2 NDVI與地形因子的多尺度空間分布特征

基于小波低頻系數(shù)的NDVI與地形因子空間分布見(jiàn)圖2。根據(jù)國(guó)家《土壤侵蝕分級(jí)標(biāo)準(zhǔn)》[23]中的植被覆蓋度分級(jí)標(biāo)準(zhǔn)，并結(jié)合研究區(qū)的實(shí)際植被覆蓋情況，將植被覆蓋度劃分5個(gè)等級(jí)：＜0.1（裸地）、0.10～0.35（低植被覆蓋）、0.36～0.45（中低植被覆蓋）、0.46～0.60（中等植被覆蓋）、＞0.6（高植被覆蓋）。

從0.5～8 km尺度，隨著分解尺度的增大，小波變換方法的平滑效應(yīng)逐步強(qiáng)化，輪廓信息凸顯出來(lái)而內(nèi)部細(xì)節(jié)信息逐漸被弱化，使得NDVI空間分布規(guī)律更加簡(jiǎn)單明了。NDVI＞0.6的區(qū)域主要為兩列北東-南西走向的條帶狀區(qū)域，邊緣呈齒狀分布。利用多尺度分解后的NDVI空間分布圖和高清遙感影像進(jìn)行參照比較分析，可得知，該區(qū)植被主要為常綠闊葉與針闊葉混交林，NDVI＞0.6的區(qū)域與平行嶺谷的分布和走向極為吻合，齒狀邊緣主要是受山脈兩側(cè)嶺谷相間的地貌特征影響。另外，NDVI＞0.6的區(qū)域在長(zhǎng)江兩岸也有局部分布，這些區(qū)域的海拔在500米以下。NDVI為0.10～0.35的區(qū)域主要集中在平行嶺谷兩側(cè)的丘陵地帶，地表植被稀少，多為農(nóng)耕區(qū)。NDVI＜0.1的區(qū)域?yàn)槁愕?，裸地主要分布在農(nóng)田附近。

綜合5個(gè)尺度看，在高程、坡度和坡向3個(gè)地形因子中，高程和NDVI的空間分布格局最為相似，坡度和坡向與NDVI的空間分布相似度較低。坡度和坡向的多尺度空間變化特征不如NDVI和高程明顯，但仍可識(shí)別出其在0.5～4 km尺度上呈北東-南西走向。

2.3 NDVI與地形因子的多尺度小波相關(guān)分析

多尺度相關(guān)分析可以更好地區(qū)分不同要素，將其耦合程度與單一尺度相關(guān)分析的耦合程度比較，顯得更加合理。基于小波低頻系數(shù)的NDVI與地形因子的相關(guān)性見(jiàn)表2。地形因子與NDVI相關(guān)性的大小的順序?yàn)楦叱蹋酒露龋酒孪颍㈦S小波分解尺度的增加而增大。從0.5～8 km尺度，NDVI與高程均呈極顯著正相關(guān)，相關(guān)性隨分解尺度的增加而增大，由0.256逐漸增加到0.507。主要原因在于該區(qū)域受海拔影響，在較低處地形為丘陵和平壩，受人類活動(dòng)影響程度大，農(nóng)業(yè)生產(chǎn)率高，伴隨著海拔的增加，人類活動(dòng)逐漸減少，植被生長(zhǎng)茂盛。同時(shí)，由于研究區(qū)地形起伏空間尺度較大，因此隨著空間尺度的增加，高程對(duì)NDVI的影響作用越來(lái)越明顯。坡度與NDVI在各尺度下均呈正相關(guān)，相關(guān)性隨分解尺度的增加而增大。呈正相關(guān)的原因在于該區(qū)坡度小的區(qū)域往往開(kāi)墾為農(nóng)田，坡度大于25°的區(qū)域受國(guó)家退耕還林工程影響植被覆蓋較高。坡向與NDVI在5個(gè)尺度上均無(wú)相關(guān)性，其原因是該區(qū)山脈呈北東-南西走向，因此受太陽(yáng)輻射影響的差異較小。同時(shí)，本研究中坡向的最小空間分辨率為0.5 km，而坡向是在較小尺度上對(duì)NDVI產(chǎn)生影響的因子，因此在較大尺度上對(duì)NDVI的影響不明顯。

表2 不同尺度下NDVI及地形因子的相關(guān)系數(shù)Table 2 Correlation coefficients between NDVI and topographic factors at multi-scale

首先，本文構(gòu)建了基于小波低頻系數(shù)的NDVI回歸方程（表3）。是為了更好地分析整理NDVI與地形因子，這兩者之間的定量關(guān)系。本研究又在SPSS軟件中采用逐步回歸法對(duì)NDVI進(jìn)行擬合，解決了建模因子共線性的問(wèn)題。在0.5、1和2 km尺度，最優(yōu)擬合方程的建模因子僅有高程，說(shuō)明在該尺度范圍雖然坡度與NDVI也有一定的相關(guān)性，但NDVI的高低主要由高程決定。在4 km和8 km尺度，擬合方程的建模因子為高程和坡度，說(shuō)明在較大尺度上NDVI高低受控于高程和坡度。

表3 不同尺度下NDVI的擬合方程Table 3 Regression prediction model and R2for NDVI at multi-scale

由回歸方程的決定系數(shù)分析可知，伴隨分解尺度的增加，擬合方程的R2由0.066提高到0.313，說(shuō)明伴隨分解尺度的增加，圖片中微小的信息被過(guò)濾，小波低頻系數(shù)反映的空間總體規(guī)律趨于清晰，NDVI和地形因子的關(guān)系更加明顯。

3 結(jié)論與討論

本文選擇二維小波變換的實(shí)驗(yàn)方法，對(duì)褶皺山區(qū)NDVI與地形因子之間的多尺度空間耦合性進(jìn)行分析研究，研究結(jié)果如下：①多尺度下NDVI和高程的空間分布相似度最高，高植被覆蓋（NDVI＞0.6）區(qū)域主要集中于褶皺山脈和長(zhǎng)江兩岸，低植被覆蓋區(qū)域（NDVI＜0.35）主要集中于丘陵農(nóng)地。②綜合5個(gè)尺度看，NDVI與地形因子相關(guān)性的大小的順序?yàn)楦叱蹋酒露龋酒孪?，其中，NDVI與高程和坡度的值都呈正相關(guān)，與坡向幾乎無(wú)相關(guān)性，且相關(guān)性隨小波分解尺度的增加而增大。③在0.5、1和2 km尺度，NDVI的最佳擬合方程中建模因子僅需高程，而在4 km和8 km尺度，最佳擬合方程的建模因子為高程和坡度，擬合方程的決定系數(shù)隨著分解尺度的增加而增加。

高程、坡度與NDVI的相關(guān)性較高，且隨空間尺度的增加而增大，說(shuō)明兩者傾向于在宏觀尺度上影響NDVI的空間分布，這與胡云鋒等[8，14-16]在蒙古高原、青藏高原、浙江省和武夷山的研究結(jié)果一致。NDVI回歸方程的決定系數(shù)在5個(gè)尺度下均偏低，其原因在于NDVI與地形因子并非簡(jiǎn)單的線性關(guān)系，并且一些重要的影響因素如土地利用類型、距河流的距離等沒(méi)有考慮，在今后的研究中可以嘗試增加這些因子。

小波變換能夠反映出NDVI和地形因子的多尺度特征，對(duì)于小波變換的方法來(lái)說(shuō)仍然要考慮以下幾個(gè)問(wèn)題，首先因?yàn)楸狙芯恐兴捎玫腗ODIS影像的空間分辨率為0.5 km，所以不宜對(duì)小波分解后的細(xì)節(jié)信息進(jìn)行分析，今后的研究中可以考慮采用分辨率較高的數(shù)據(jù)源，以便對(duì)NDVI和地形因子小波分解后的多尺度細(xì)節(jié)信息進(jìn)行研究，以進(jìn)一步揭示兩者的多尺度耦合關(guān)系；其次，數(shù)學(xué)領(lǐng)域的研究已經(jīng)表明，選擇不同的小波基函數(shù)，會(huì)導(dǎo)致小波分解后結(jié)果發(fā)生變化[7，20]，在研究中仍需進(jìn)一步尋找有關(guān)小波基的篩分方法；最后，目前所采用的小波分析方法僅能給出斷續(xù)尺度（即尺度必須為2n）上的小波系數(shù)，無(wú)法在連續(xù)尺度上開(kāi)展計(jì)算和分析，此問(wèn)題在今后的研究過(guò)程中仍需得到進(jìn)一步解決[8，24]。