苗旺元，王加勝，許人偉，王 剛，陳 波

（1.云南師范大學(xué)信息學(xué)院；2.西部資源環(huán)境地理信息技術(shù)教育工程研究中心，云南昆明 650500）

0 引言

三七是云南省的名貴中藥材和特色經(jīng)濟(jì)作物，藥用歷史悠久［1］。三七對產(chǎn)地的氣候、海拔、土壤和光照等生長環(huán)境條件有特殊要求，其主要分布在中國西南部海拔1 300-1 800m，北緯23.5°地區(qū)，云南文山州為原產(chǎn)地，而丘北縣為其核心主產(chǎn)區(qū)。利用遙感技術(shù)［2-4］監(jiān)測三七種植變化趨勢是農(nóng)業(yè)監(jiān)測的重要任務(wù)之一，可快速準(zhǔn)確分析其分布變化情況，為指導(dǎo)丘北縣農(nóng)業(yè)發(fā)展提供重要依據(jù)。

傳統(tǒng)的三七種植監(jiān)測及分析變化主要通過普查統(tǒng)計和實(shí)地調(diào)查，不僅耗時耗力，且難以從宏觀尺度全面反映三七種植面積及其變化［5］。遙感技術(shù)具有探測范圍廣、信息豐富、數(shù)據(jù)周期短、處理快等特點(diǎn)，已成為實(shí)時監(jiān)測土地利用變化的主要手段［6］。周應(yīng)群等［7］將5m 分辨率的SPOT5 影像和30m 分辨率的Landsat 5 影像融合，采用目視解譯的方法提取馬關(guān)縣的三七種植面積，融合后影像的三七面積提取精度為92.7%，取得了相對理想的效果；戴晨曦等［5］利用Landsat 數(shù)據(jù)通過支持向量機(jī)提取2010、2012、2014 和2015 年文山州和紅河州的三七種植區(qū)，并分析了三七種植區(qū)時空變化與地形、政策和三七價格因素的關(guān)系，但是沒有考慮山體陰影會被誤分為三七種植區(qū)的問題，其結(jié)果有待進(jìn)一步改進(jìn)；朱贊等［8］以文山州4 個縣為研究區(qū)域，采用GF-1 影像，通過改進(jìn)的決策樹模型提取三七種植區(qū)，其提取精度達(dá)到87%，證明該方法可行。但以上研究只針對三七種植信息提取，沒有長時間系列的監(jiān)測分析。目前，關(guān)于三七的研究大多聚焦在提取方法上，作物監(jiān)測研究主要集中在水稻、小麥和玉米等常見農(nóng)作物，名貴中草藥的時空變化監(jiān)測是農(nóng)業(yè)遙感一個新的應(yīng)用方向，而三七最大的生長特點(diǎn)是需要較好的遮陽條件，因此種植三七的地塊都覆蓋有黑色遮陽網(wǎng)，這種特殊的黑色遮陽網(wǎng)與山體陰影的光譜相似，造成易將山體陰影誤判為三七種植區(qū)現(xiàn)象，影響三七種植的時空變化分析，而顧及山體陰影的三七種植監(jiān)測研究較少。

本文針對丘北縣三七種植遙感監(jiān)測問題，通過GEE（Google Earth Engine）云平臺，基于Java-Script 語言，選取云南省丘北縣為研究區(qū)，基于2010 年、2014 年、2019 年Landsat 和DEM 數(shù)據(jù)，采用效果較好的回歸校正法恢復(fù)山體陰影亮度，在此基礎(chǔ)上構(gòu)建特征集，并通過J-M 距離選擇最優(yōu)特征集，利用隨機(jī)森林分類算法［9-11］得到3 年丘北縣的土地利用分類結(jié)果，分別提取3 年三七種植區(qū)面積，統(tǒng)計各鄉(xiāng)鎮(zhèn)三七種植分布情況，并分析三七時空變化趨勢及其與河流、道路、海拔、坡度和坡向的關(guān)系。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

丘北縣（見圖1）位于云南省東南部（彩圖掃OSID 碼可見，下同），位于東經(jīng)103°34′-104°45′、北緯23°45′～24°28′之間，總面積5 038km2，海拔732-2 504m。丘北地處滇東南巖溶山原丘陵地帶，地勢西南高、東北低，由于地處低緯季風(fēng)區(qū)域，氣候總體上屬中亞熱帶高原季風(fēng)氣候，年平均氣溫13.2～19.7 ℃，年平均降雨量1 000-1 270mm，適合三七生長［12］。由于三七的社會需求量不斷增加，加上其種植需要較長時間的輪作，丘北縣作為三七種植的主產(chǎn)區(qū)，當(dāng)?shù)夭块T制定了一系列政策扶植三七產(chǎn)業(yè)，使得近年來三七產(chǎn)業(yè)在丘北得到了飛速發(fā)展。本文選取丘北縣為研究區(qū)監(jiān)測三七種植的分布變化，為指導(dǎo)當(dāng)?shù)剞r(nóng)業(yè)發(fā)展提供重要依據(jù)。

1.2 研究數(shù)據(jù)

1.2.1 遙感數(shù)據(jù)

本文數(shù)據(jù)包括Landsat 8 表觀反射率數(shù)據(jù)、Landsat 5 表觀反射率數(shù)據(jù)、DEM 數(shù)據(jù)，這些數(shù)據(jù)來源于Google 地球引擎（Google Earth Engine，GEE），其中Landsat 數(shù)據(jù)分辨率都為30m 且經(jīng)過了大氣校正和幾何校正。分別挑選出云量較少且與被采樣日期相近的影像，如表1 所示。對這些數(shù)據(jù)運(yùn)用CFMask［13］算法進(jìn)行去云處理，得到無云影像。

Fig.1 Overview of the study area圖1 研究區(qū)概況

Table 1 Landsat image data表1 Landsat 影像數(shù)據(jù)

1.2.2 樣本數(shù)據(jù)

本文將研究區(qū)地物分為水體、裸地、林地、建筑用地、草地、耕地、三七種植區(qū)7 類，并通過Google Earth 分別采集2014、2010 和2019 年年內(nèi)的7 類樣本，其中2010 年樣本數(shù)366 個，2014 年樣本數(shù)630 個，2019 年樣本數(shù)384 個，采樣過程中盡可能使樣本均勻分布在整個研究區(qū)。針對不同時期的樣本點(diǎn)選用被采日期相近的影像來訓(xùn)練，最后將采集好的樣本數(shù)據(jù)導(dǎo)入GEE 云平臺。

2 研究方法

技術(shù)路線如圖2 所示。首先對利用太陽入射角與DEM數(shù)據(jù)檢測出每幅影像的山體陰影區(qū)域采用回歸校正法修復(fù)山體陰影，在修復(fù)山體陰影影像基礎(chǔ)上進(jìn)行去云處理，構(gòu)建光譜特征、指數(shù)特征、地形特征、纓帽變換特征和紋理特征，通過J-M 距離得到最優(yōu)特征集，并通過隨機(jī)森林算法利用樣本數(shù)據(jù)和最優(yōu)特征集對每年的合成影像分類處理，得到三七種植區(qū)提取結(jié)果，并應(yīng)用驗(yàn)證數(shù)據(jù)開展精度評估；分析三七種植重心變化趨勢，統(tǒng)計各鄉(xiāng)鎮(zhèn)三七種植面積，最后分析三七種植面積變化及其與河流、道路、海拔、坡度和坡向的關(guān)系。

Fig.2 Technical route圖2 技術(shù)路線

2.1 山體陰影恢復(fù)

利用DEM 數(shù)據(jù)與Landsat 影像的太陽高度角和太陽方位角檢測出該幅影像的山體陰影區(qū)域，檢測方法如式（1）所示：

式中，z為太陽天頂角，s為坡度角，az為太陽方位角，as為坡向，cosi為計算出的太陽入射角余弦值，若cosi小于0則此像元為山體陰影區(qū)域。

利用回歸校正法［14］對檢測出的山體陰影區(qū)域進(jìn)行恢復(fù)處理，回歸校正法是通過統(tǒng)計整幅影像亮度值，利用各波段亮度值與太陽入射角的強(qiáng)相關(guān)性進(jìn)行回歸統(tǒng)計，找到合適的回歸方程對山體陰影進(jìn)行修復(fù)［15］，其表達(dá)式如式（2）、式（3）所示。

式中，DNm為修復(fù)后的像元亮度值；DN為修復(fù)前的像元亮度值；DNa為各波段像元亮度均值；i為太陽入射角；b和m通過回歸計算得到。

2.2 特征構(gòu)建

過多的特征參與分類會導(dǎo)致信息冗余，降低分類精度，本文利用J-M（Jeffries-Matusita distance）距離［16］計算在某特征下三七種植區(qū)與其他地類的可分度，剔除可區(qū)分度小的特征，選擇可區(qū)分度大的特征構(gòu)建最優(yōu)特征集。J-M距離表達(dá)式如式（4）、式（5）所示。

式中，Bij為在某特征下類別i和j的巴氏距離；Vi和Vj分別為類別i和j的樣本協(xié)方差矩陣；Mi和Mj分別為類別i和j的樣本均值向量。J-M 距離取值范圍為0～2，值越高，表示兩類別的可區(qū)分度越高。

2.2.1 特征選取

采用Landsat 數(shù)據(jù)常用的B1-B7 波段作為起始光譜特征集，通過計算基于每個波段的J-M 距離，在Landsat 8 數(shù)據(jù)中剔除J-M 值最小的B1 波段，使用B2-B7 波段作為最優(yōu)光譜特征集，Landsat 5 數(shù)據(jù)中剔除J-M 值最小的B6 波段，使用其余波段作為最優(yōu)光譜特征集。

采用歸一化植被指數(shù)（NDVI）和比值植被指數(shù)（RVI）增強(qiáng)三七與林地、草地、耕地的區(qū)分度，采用歸一化建筑指數(shù)（NDBI）增強(qiáng)三七種植區(qū)與建筑用地的區(qū)分度，采用改進(jìn)型水體指數(shù)（MNDWI）增強(qiáng)三七種植區(qū)與水體的區(qū)分度。

通過NDVI 與灰度共生矩陣［17］得到18 個紋理特征，通過J-M 距離選取角二階矩（Angular Second Moment，ASM）、對比度（Contrast，CON）和方差（Correlation，CORR）作為最優(yōu)紋理特征集。

由于三七生長環(huán)境與地形緊密相關(guān)，因此通過DEM 數(shù)據(jù)選取海拔（Elevation）、坡度（Slope）和坡向（Aspect）3 個特征參與特征構(gòu)建。

選擇纓帽變換得到的前3 個分量作為纓帽變換特征，其中綠度（Greenness）反映了植被覆蓋信息，亮度（Bright?ness）反映了裸露土地信息，濕度（Wetness）反映了水分信息［18］。將優(yōu)選后的光譜特征、指數(shù)特征、紋理特征、地形特征和纓帽變換特征構(gòu)建成最優(yōu)特征集。

3 結(jié)果與分析

將以上訓(xùn)練樣本與最優(yōu)特征集投入隨機(jī)森林分類器中，并應(yīng)用到研究區(qū)，得到2019 年、2014 年和2010 年三七提取結(jié)果，提取的3 年三七種植區(qū)分布圖如圖3 所示（彩圖掃OSID 碼可見，下同），對其結(jié)果進(jìn)行分析，其分類精度與統(tǒng)計面積如表2 所示。

Table 2 Classification results and accuracy evaluation表2 分類結(jié)果與精度評價

Fig.3 Distribution of Panax notoginseng planting area圖3 三七種植區(qū)分布

3.1 三七種植的時空變化

由表2 可知，丘北縣三七種植面積從2010 年的35.73 km2增長到2014 年的56.1 km2；2014-2019 年，丘北縣三七種植面積由56.1 km2擴(kuò)張到147.07 km2。從圖3 可以看出三七種植的特點(diǎn)：三七分布呈現(xiàn)零散種植的特點(diǎn)，由于三七一般需要種植3 年才能采收，其種植過的區(qū)域需要經(jīng)過最長20 年的土地休養(yǎng)，因此需要更換區(qū)域種植［19］。

通過三七提取結(jié)果得到2010、2014 和2019 年丘北縣各個鄉(xiāng)鎮(zhèn)的三七種植面積，如表3 所示。2010 年丘北縣三七主要種植在八道哨彝族鄉(xiāng)和雙龍營鎮(zhèn)，其面積分別為7.61km2和7.24km2；2014 年三七主要種植在樹皮彝族鄉(xiāng)和雙龍營鎮(zhèn)，相對于2010 年，樹皮彝族鄉(xiāng)三七種植面積增加了14.86km2；2019 年三七在各個鄉(xiāng)鎮(zhèn)都有廣泛種植，其中種植面積最小的是5.96km2的溫瀏鄉(xiāng)，種植面積最大的是樹皮彝族鄉(xiāng)，達(dá)24.67km2。相對于2014 年，各個鄉(xiāng)鎮(zhèn)種植面積都有增長，其中平寨鄉(xiāng)和舍得彝族鄉(xiāng)增長面積最為明顯，分別增長了12.28km2和16.47km2。

Table 3 Annual changes of Panax notoginseng planting area in each township of the study area表3 研究區(qū)各鄉(xiāng)鎮(zhèn)三七種植面積年際變化（km2）

由圖4 可以看出，2010-2014 年，三七種植重心由丘北縣雙龍營鎮(zhèn)向西南方向轉(zhuǎn)移7.35km，到了八道哨彝族鄉(xiāng)，2014-2019 年三七種植重心在八道哨彝族鄉(xiāng)內(nèi)向東北方向轉(zhuǎn)移了2.7km，這是由于八道哨彝族鄉(xiāng)位于丘北縣中心區(qū)域，緊鄰丘北縣城，交通方便，地勢較為平坦，水源充足，適宜耕作，是三七產(chǎn)業(yè)的主產(chǎn)區(qū)和功能區(qū)。

Fig.4 Shift of center of gravity of Panax notoginseng planting圖4 三七種植重心轉(zhuǎn)移情況

3.2 三七種植與河流、道路關(guān)系

統(tǒng)計分析2010、2014 和2019 年內(nèi)各三七種植區(qū)與其最近的河流、道路的最小距離，并計算出三七種植區(qū)距離河流、道路的平均最小距離，如圖5 所示。其中河流包括常年河、時令河和溝渠，道路包括4 級以上的公路?？梢钥闯鱿鄬τ?010 年，2014 年的三七種植區(qū)與河流的平均最小距離減少到了2.73km，但在2019 年，該距離為2.81km；2010-2014 年，三七種植區(qū)與道路的平均最小距離由1.78km 減少到了1.28km，2014-2019 年，該距離又增加到了1.57km。2010-2019 年間，三七種植區(qū)與河流、道路的平均最小距離呈現(xiàn)先大幅度減小后小幅度增大的趨勢，距離大幅度減小是由于種植戶考慮耕作和運(yùn)輸距離成本，將三七種植在距離河流、道路較近的區(qū)域，而距離大幅度增大是因種植面積飛速增長，需要開荒更多的適宜種植區(qū)，其中有些種植區(qū)距離河流、道路相對較遠(yuǎn)，造成三七種植區(qū)域距河流、道路的平均最小距離小幅度變大，但沒有超過2010 年的距離。

Fig.5 Average minimum distance between Panax notoginseng planting area rivers and roads圖5 三七種植區(qū)與河流、道路的平均最小距離

3.3 三七種植與海拔的關(guān)系

通過2010、2014 和2019 年三七種植區(qū)域提取結(jié)果，結(jié)合DEM 數(shù)據(jù)，以300m 為海拔間隔，統(tǒng)計分析三七種植區(qū)分布的時空變化與海拔關(guān)系，結(jié)果表明丘北縣三七主要種植在海拔1 300-2 200m 之間，其中海拔1 300-1 900m 最適宜三七種植。如表4 所示，可以看出2010 年到2014 年、2014年到2019 年，三七種植分布向中低海拔擴(kuò)張，2010 年分布在海拔1 900m 以下的三七種植面積為25.3km2，2014 年分布在海拔1 900m 以下的三七種植面積為49.38km2，而2019年同海拔范圍的三七種植面積達(dá)到了127.63km2。

Table 4 Area of Panax notoginseng planting area at different altitudes表4 不同海拔區(qū)間三七種植區(qū)面積（km2）

3.4 三七種植與坡度的關(guān)系

統(tǒng)計分析丘北縣2010、2014 和2019 年三七種植區(qū)在不同坡度的分布情況。根據(jù)水利部頒發(fā)的中華人民共和國行業(yè)執(zhí)行標(biāo)準(zhǔn)SL190-96《土壤侵蝕分類分級標(biāo)準(zhǔn)》［20］，坡度分為微坡（＜5°）、較緩坡［5°，8°）、緩坡［8°，15°）、較陡坡［15°，25°）、陡坡［25°，35°），急陡坡（≥35°）6 級。如表5 所示，2010 年和2014 年三七主要種植在微坡、較微坡、緩坡和較緩坡區(qū)域，位于陡坡和急陡坡的三七種植面積不到5km2，2019 年三七種植在陡坡和急陡坡的面積達(dá)到了46.94km2，這說明三七種植有向坡度較大的區(qū)域擴(kuò)張的趨勢。

Table 5 Area of Panax notoginseng planting area with different slopes表5 不同坡度的三七種植區(qū)面積（km2）

3.5 三七種植與坡向的關(guān)系

統(tǒng)計分析丘北縣在2010、2014 和2019 年三七種植區(qū)在不同坡向的分布變化。將坡向分為5 種：坡向角在［0，45 °）和［315 °，360 °）間為北坡；坡向角在［45 °，135 °）間為東坡；坡向角在［135 °，225 °）間為南坡；坡向角在［225 °，315 °）為西坡，還有平坡地，結(jié)果如表6 所示。丘北縣三七很少種植在平坡地，北坡與西坡的種植面積高于其他坡向，這是因?yàn)樵颇系貐^(qū)受東南季風(fēng)影響，迎風(fēng)坡受到的水熱條件優(yōu)于背風(fēng)坡，而三七適宜種植在陰涼區(qū)域，背風(fēng)坡更適宜種植三七。

Table 6 Planting area of Panax notoginseng in different directions表6 不同坡向三七種植區(qū)面積（km2）

4 結(jié)語

本文基于山體陰影修復(fù)的三七種植區(qū)提取方法研究了云南省丘北縣2010 年、2014 年和2019 年三七種植區(qū)情況，探討三七種植區(qū)時空變化，并分析了其變化與河流、道路、海拔、坡度和坡向因素的相互關(guān)系。主要研究結(jié)論如下：

（1）結(jié)合回歸校正法進(jìn)行山體陰影修復(fù)，在處理后的影像基礎(chǔ)上提取三七種植區(qū)，其總體精度都在88%以上，Kappa 系數(shù)都在85%以上，能夠滿足遙感監(jiān)測需求，證明遙感技術(shù)對三七種植區(qū)進(jìn)行動態(tài)監(jiān)測是可行的。

（2）對丘北縣三七種植區(qū)面積變化分析表明，2010-2014 年，丘北縣三七種植區(qū)面積從35.73km2增加到了56.1km2，而在2019 年，三七種植面積達(dá)到了147.07km2。總體來看，三七種植范圍在2010-2019 年期間呈現(xiàn)逐漸擴(kuò)張趨勢。2010-2014 年，丘北縣三七種植重心向西南方向轉(zhuǎn)移7.35km，到了縣域中心附近的八道哨彝族鄉(xiāng)，2014-2019年三七種植重心在八道哨彝族鄉(xiāng)內(nèi)向東北方向轉(zhuǎn)移了2.7km。

（3）相對于2010 年，2014 年的三七種植區(qū)與河流的平均最小距離減少到2.73km，但2019 年該距離為2.81km；2010-2014 年，三七種植區(qū)與道路的平均最小距離由1.78km 減少到1.28km，2014-2019 年，該距離又增加到1.57km。由于種植的距離成本和開墾新的種植區(qū)，平均最小距離呈現(xiàn)先大幅減小后小幅度增加的趨勢。

（4）從地形特征來看，三七主要種植在海拔［1 300，2 200）m 之間，其中在海拔［1 300，1 900）m 種植的最多，2010-2019 年，三七種植區(qū)有向海拔［1 300，1 900）m 擴(kuò)張的趨勢；三七主要分布在坡度＜25°區(qū)域，并且有向≥35°坡度區(qū)域擴(kuò)張的趨勢；就坡向因素而言，三七主要種植在北坡與西坡，東坡與南坡也有少量三七種植。

后續(xù)要探討更多的山體陰影修復(fù)方法，找到對山體陰影修復(fù)效果最好的方案。另外，本研究所用數(shù)據(jù)源為分辨率30m 的Landsat 數(shù)據(jù)，后續(xù)可用分辨率更高的數(shù)據(jù)，以滿足更精確的三七種植檢測研究。