韓阿茹汗 包斯琴 馬雅婷 那日嘎拉 張佳華

(內(nèi)蒙古農(nóng)業(yè)大學(xué)沙漠治理學(xué)院，內(nèi)蒙古 呼和浩特 010018)

土地沙化是指因風(fēng)力、重力等外營力作用，導(dǎo)致地表組成物質(zhì)中的細(xì)顆粒損失，或者原有地表因?yàn)橥鈦砩车[覆蓋等原因而導(dǎo)致土地退化的過程[1，2]。科左中旗境內(nèi)有著中國4大沙地之一的科爾沁沙地，是防風(fēng)治沙的重要區(qū)域，曾是水草豐沛的科爾沁大草原，進(jìn)入21世紀(jì)以來，受人口劇增、濫墾濫牧等多種因素的影響，導(dǎo)致土地退化嚴(yán)重、水土流失加劇，生態(tài)環(huán)境急劇惡化[3]。科爾沁沙地殲滅戰(zhàn)是習(xí)近平總書記親自部署的“三北”攻堅(jiān)戰(zhàn)3大標(biāo)志性戰(zhàn)役之一，2023年7月31日，內(nèi)蒙古自治區(qū)建設(shè)我國北方重要生態(tài)安全屏障促進(jìn)條例正式發(fā)布，明確表示了土地沙化荒漠化綜合防治應(yīng)因地制宜，突出重點(diǎn)治理，全面提升生態(tài)系統(tǒng)的質(zhì)量和穩(wěn)定性，因此分析科左中旗土地沙化時(shí)空動(dòng)態(tài)變化，了解土地沙化的空間分布具有現(xiàn)實(shí)意義。

隨著遙感技術(shù)的快速發(fā)展，衛(wèi)星圖像成為監(jiān)測(cè)和分析土地利用時(shí)空變化的重要數(shù)據(jù)支撐，谷歌地球引擎遙感大數(shù)據(jù)平臺(tái)因其獲取數(shù)據(jù)便捷，處理計(jì)算數(shù)據(jù)準(zhǔn)確，已廣泛成為遙感監(jiān)測(cè)的重要工具[4]。目前對(duì)于土地沙化的監(jiān)測(cè)方法，野外實(shí)地調(diào)查與遙感目視解譯結(jié)合，這種方法耗時(shí)耗力，不適用于長(zhǎng)時(shí)間序列的監(jiān)測(cè)；多指標(biāo)綜合監(jiān)測(cè)法，單一的地表參量指標(biāo)不能夠全面反映區(qū)域沙化，因此想要監(jiān)測(cè)結(jié)果更加合理，應(yīng)選取多個(gè)地表參量進(jìn)行共同監(jiān)測(cè)。如，曾永年等[5]利用歸一化植被指數(shù)(Normalized Difference Vegetation Index，NDVI)和地表反照率指數(shù)(Surface Albedo)構(gòu)建了特征空間，選取沙漠化遙感監(jiān)測(cè)差值指數(shù)(Desertification Difference Index，DDI)作為監(jiān)測(cè)指標(biāo)。劉曉茜等[6]通過建立Albedo-NDVI特征空間沙化監(jiān)測(cè)模型法對(duì)騰格里沙漠東南緣格林灘綠洲20年間沙化土地進(jìn)行動(dòng)態(tài)監(jiān)測(cè)及分析。該模型充分利用了多維遙感信息，可以自動(dòng)提取土地沙化信息，可以避免單一地表參量的誤差，有利于土地沙化的監(jiān)測(cè)與分析。

科爾沁左翼中旗作為科爾沁沙地的重要部分，亟需對(duì)土地沙化監(jiān)測(cè)與驅(qū)動(dòng)力分析?；诖?，本文以科爾沁左翼中旗為研究區(qū)，通過提取歸一化植被指數(shù)與地表反照率指數(shù)分析2000—2021年土地沙化的時(shí)空演變特征及轉(zhuǎn)移特征，并利用地理空間探測(cè)器模型探究自然驅(qū)動(dòng)因子和社會(huì)經(jīng)濟(jì)驅(qū)動(dòng)因子對(duì)土地沙化的影響，進(jìn)而對(duì)土地沙化提供新的治理思路，為防沙治沙提供科學(xué)依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

科爾沁左翼中旗位于內(nèi)蒙古自治區(qū)通遼市(E121°08′～123°32′，N43°32′～44°32′)，地處大興安嶺東南邊緣，西遼河北岸，地勢(shì)西北高東南低。旗內(nèi)氣候?qū)僦袦貛Т箨懶约撅L(fēng)氣候，氣候干燥多風(fēng)沙。土壤類型主要包括草甸土、風(fēng)沙土、栗鈣土、鹽土、沼澤土，植被類型有沙生植被、草甸植被、鹽生植被等。

圖1 研究區(qū)地理位置

1.2 數(shù)據(jù)來源

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

科爾沁左翼中旗2000年、2005年、2010年、2015年、2021年的Landsat影像數(shù)據(jù)均下載自GEE(Google Earth Engine)云平臺(tái)。其中，2000年、2005年、2010年為TM影像，2015年和2021年為OLI影像，分辨率均為30m。

1.2.2 驅(qū)動(dòng)因子數(shù)據(jù)

本文選擇自然驅(qū)動(dòng)因素和社會(huì)經(jīng)濟(jì)驅(qū)動(dòng)因素的10個(gè)驅(qū)動(dòng)因子?xùn)鸥駭?shù)據(jù)進(jìn)行土地沙化驅(qū)動(dòng)研究。其中，年平均氣溫和年降水量數(shù)據(jù)通過國家氣象科學(xué)中心(http：//www.geodata.cn)途徑獲得，人口空間分布、GDP、土壤類型、土壤侵蝕、平均風(fēng)速、日照時(shí)數(shù)均來源于中國科學(xué)院資源環(huán)境科學(xué)與數(shù)據(jù)中心(https：//www.resdc.cn)，坡度數(shù)據(jù)使用Arcgis 10.8坡度功能計(jì)算，DEM數(shù)據(jù)來源于地理空間數(shù)據(jù)云(http：//www.gscloud.cn)。

1.3 研究方法

1.3.1 Albedo-NDVI特征空間

地表反照率(Albedo)是重要的遙感反演參數(shù)，反映了地球表面對(duì)太陽短波輻射反射特效。Albedo受到地表粗糙度、地表含水量因素影響而發(fā)生變化。參照梁順林[7]建立的反演模型來估算地表反照率，計(jì)算公式：

Albedo=0.356B+0.130R+0.373NIR+0.085SWIR1+0.072SWIR2-0.0018

(1)

式中，B為藍(lán)光波段的反射值；R為紅光波段的反射值；SWIR為短波紅外反射值，利用Arcgis柵格計(jì)算器計(jì)算出各年份的地表反照率。

歸一化植被指數(shù)(NDVI)是一種使用較為廣泛的植被指數(shù)，是現(xiàn)階段充分體現(xiàn)植被具體生長(zhǎng)情況的一種指數(shù)，依據(jù)該指數(shù)可知植被的覆蓋量、植物量以及生長(zhǎng)狀況等數(shù)據(jù)[8]?？赏ㄟ^遙感影像近紅外波段和紅光波段的反射值計(jì)算，計(jì)算公式：

NDVI=(NIR-R)/(NIR+R)

(2)

式中，NIR為近紅外波段的反射值；R為紅光波段的反射值。

為了后續(xù)數(shù)據(jù)處理的方便，對(duì)二者進(jìn)行歸一化處理：

N=[(NDVI-NDVImin)/(NDVImax-NDVImin)]×100%

(3)

A=[(Albedo-Albedomin)/(Albedomax-Albedomin)]×100%

(4)

式中，N為歸一化后的NDVI；A為歸一化后的Albedo；NDVImin、Albedomin分別為NDVI和Albedo的最小值；NDVImax、Albedomax分別為NDVI和Albedo的最大值。

對(duì)歸一化后的Albedo和NDVI數(shù)據(jù)，利用Origin軟件進(jìn)行線性回歸分析并繪制特征空間散點(diǎn)圖，可以得出，在Albedo-NDVI特征空間中，Albedo和NDVI有著明顯的負(fù)線性關(guān)系，計(jì)算得到擬合方程Albedo=a×NDVI+b，a為方程斜率，b為方程截距，即隨著土地沙化程度的加劇，NDVI的值不斷變小，而Albedo的值逐漸變大。高中低等級(jí)的植被覆蓋區(qū)域可以明顯的劃分，反映了沙化土地的地表信息和水熱組合狀況。其線性擬合方程：

Albedo2000=-0.6243×NDVI+0.8697

(5)

Albedo2005=-0.6490×NDVI+0.8552

(6)

Albedo2010=-0.7785×NDVI+0.9322

(7)

Albedo2015=-0.6413×NDVI+0.8742

(8)

Albedo2021=-0.4038×NDVI+0.6445

(9)

圖2 Albedo-NDVI特征空間散點(diǎn)圖

1.3.2 沙漠化遙感監(jiān)測(cè)差值指數(shù)計(jì)算

在以上構(gòu)建的Albedo-NDVI特征空間中不難發(fā)現(xiàn)，兩者的關(guān)系為線性負(fù)相關(guān)，若將在代表土地沙化變化趨勢(shì)的垂直方向劃分，即在線性擬合的趨勢(shì)線上劃分二維空間，由此反演出沙漠化差值指數(shù)模型(DDI)[9]，通過此模型可以區(qū)分不同等級(jí)的沙化土地，進(jìn)而可以對(duì)不同程度的沙化土地進(jìn)行監(jiān)測(cè)和分析，計(jì)算公式：

DDI=(-1/a)×NDVI-Albedo

(10)

式中，DDI為沙漠化差值指數(shù)；a為特征空間線性擬合的斜率。計(jì)算結(jié)果表達(dá)式：

DDI2000=1.6015×NDVI-Albedo

(11)

DDI2005=1.5408×NDVI-Albedo

(12)

DDI2010=1.2844×NDVI-Albedo

(13)

DDI2015=1.5593×NDVI-Albedo

(14)

DDI2021=2.4765×NDVI-Albedo

(15)

結(jié)合《第五次全國土地荒漠化和沙化監(jiān)測(cè)技術(shù)規(guī)定》和Jenk自然間斷點(diǎn)法[10]，將DDI分為極重度沙化、重度沙化、中度沙化、輕度沙化和非沙化。自然間斷點(diǎn)法具有分級(jí)的良好界限，可以把同一類別中的數(shù)據(jù)值差異值差異最小化，把類別之間的數(shù)據(jù)值差異最大化。

1.3.3 Theil-Sen Median趨勢(shì)分析和Mann-Kendall檢驗(yàn)分析法

Theil-Sen median斜率估計(jì)被廣泛運(yùn)用在分析長(zhǎng)時(shí)間序列數(shù)據(jù)的趨勢(shì)分析中，是一種使用較普遍計(jì)算非參數(shù)統(tǒng)計(jì)的方法。作為一種檢驗(yàn)非參數(shù)統(tǒng)計(jì)方法，Mann-Kendall優(yōu)點(diǎn)在于對(duì)趨勢(shì)是否線性不做要求，也無需測(cè)量值服從正態(tài)分布[11]。Sen趨勢(shì)分析和M-K顯著性檢驗(yàn)已成為長(zhǎng)時(shí)間序列數(shù)據(jù)趨勢(shì)判斷的重要方法，與線性回歸分析相比，具有很強(qiáng)的抵抗數(shù)據(jù)誤差的能力，不易受異常值干擾，在一定程度上提高了準(zhǔn)確性。對(duì)于土地沙化的變化趨勢(shì)類型，根據(jù)實(shí)際情況并參考郭曉娜研究方法[12]劃分為嚴(yán)重退化、輕微退化、穩(wěn)定不變、輕微改善、明顯改善。

1.3.4 土地沙化轉(zhuǎn)移矩陣

土地沙化轉(zhuǎn)移矩陣可知其在不同階段，存在不同程度沙化類型相互轉(zhuǎn)化的特征，反映不同程度沙化之間的轉(zhuǎn)移速率和土地沙化的時(shí)空格局變化，其表達(dá)式：

(16)

式中，i與j分別為研究始末的土地沙化類型；Pij為土地沙化類型i轉(zhuǎn)為土地沙化類型j的數(shù)量。

1.3.5 地理探測(cè)器

地理探測(cè)器是王勁峰等[13]通過研究得出一種全新的統(tǒng)計(jì)學(xué)方法，用以詮釋背后驅(qū)動(dòng)因子，其核心理念是探討若某個(gè)自變量影響某個(gè)因變量，在空間分布上面，自變量與因變量存在相似性特征，可基于空間異質(zhì)性原理分析各驅(qū)動(dòng)因素對(duì)土地沙化的影響。根據(jù)計(jì)算得到的q值判斷各影響因子的影響力大小，依據(jù)q值可知，自變量影響因變量的程度，兩者呈現(xiàn)正相關(guān)關(guān)系；p值則代表統(tǒng)計(jì)學(xué)中的置信度檢驗(yàn)標(biāo)準(zhǔn)，p值<0.05則表示為差異顯著。地理探測(cè)器由因子、交互、風(fēng)險(xiǎn)以及生態(tài)等4種探測(cè)手段組成，本文使用因子探測(cè)進(jìn)行分析，其計(jì)算公式：

(17)

2 結(jié)果與分析

2.1 科爾沁左翼中旗土地沙化空間分布特征

利用Arcgis 10.8軟件，通過沙化差值指數(shù)法計(jì)算得到研究區(qū)2000年、2005年、2010年、2015年和2021年5期土地沙化遙感影像，空間分布圖如圖3所示。

圖3 科爾沁左翼中旗土地沙化空間分布圖

由圖4可知，2000—2021年土地分布情況出現(xiàn)較大變化，有一部分土地沙化程度好轉(zhuǎn)，但也有一部分土地沙化程度出現(xiàn)惡化。在西部和東北部出現(xiàn)極重度和重度土地沙化；西北部和東部存在中度沙化土地，其分布較為分散；中部和西南部主要分布的是輕度沙化土地；西南部則主要分布的是非沙化土地。總體而言，研究區(qū)土地沙化程度類型主要以中度沙化和輕度沙化土地分布為主。

圖4 2000—2021年科爾沁左翼中旗土地沙化等級(jí)轉(zhuǎn)移矩陣

2.2 科左中旗土地沙化面積變化特征

由表1可以看出，在研究區(qū)中，不同等級(jí)沙化土地面積均出現(xiàn)不同變化。2000—2021年，研究區(qū)極重度沙化土地面積出現(xiàn)先減后增加的情況，相較于2000年，2021年極重度沙化土地面積減少了0.02%；重度沙化土地面積整體上呈增長(zhǎng)趨勢(shì)，在2015年下降至19.72%，隨后在2021年又上漲到21.38%；中度沙化土地面積整體呈減少的趨勢(shì)，2021年占比26.15%；輕度沙化土地面積整體呈下降趨勢(shì)，在2015年增加到27.94%，隨后在2021年又下降到20.70%；未沙化土地面積在2021年增加到23.53%。

表1 2000—2021年科爾沁左翼中旗沙化土地面積及比例

2.3 土地沙化程度轉(zhuǎn)移分析

由圖5可知，2000—2005年，中度轉(zhuǎn)重度和輕度轉(zhuǎn)中度較為顯著，轉(zhuǎn)化面積分別為901.71km2和813.19km2，同時(shí)也有699.88km2的中度轉(zhuǎn)為輕度；2005—2010年，輕度轉(zhuǎn)中度和重度轉(zhuǎn)中度較為明顯，轉(zhuǎn)化面積分別為804.82km2和761.56km2；2010—2015年，面積轉(zhuǎn)化較大的是中度轉(zhuǎn)輕度和重度轉(zhuǎn)中度，轉(zhuǎn)化面積分別為888.78km2和751.10km2，說明此階段有較多面積的沙化土地好轉(zhuǎn)；2015—2021年主要是輕度和中度之間的轉(zhuǎn)化；2000—2021年，整體上也是輕度和中度之間的轉(zhuǎn)化，有704.61km2的中度轉(zhuǎn)成了輕度，說明總體上土地沙化有好轉(zhuǎn)趨勢(shì)。

圖5 2000—2021年科爾沁左翼中旗土地沙化趨勢(shì)空間分布圖

2.4 土地沙化變化趨勢(shì)分析

通過Theil-Sen Median趨勢(shì)分析和Mann-Kendall檢驗(yàn)結(jié)果對(duì)2000—2021年研究區(qū)土地沙化22年的變化趨勢(shì)進(jìn)行分析，從表2可以看出，變化趨勢(shì)以穩(wěn)定不變?yōu)橹?，占研究區(qū)總面積的46.71%。呈改善趨勢(shì)的土地占比為30.75%，主要分布于研究區(qū)中部；但也有22.52%的土地呈退化趨勢(shì)，主要分布于研究區(qū)西部和東北部。其中，明顯改善土地和嚴(yán)重退化土地均占0.01%?？偟膩砜?，沙化狀況趨于穩(wěn)定發(fā)展，說明應(yīng)在加大治理沙化土地的同時(shí)，堅(jiān)持邊治理邊鞏固治理成果的原則。

表2 2000—2021年科爾沁左翼中旗土地沙化趨勢(shì)分析

表3 科爾沁左翼中旗土地沙化驅(qū)動(dòng)因素解釋力統(tǒng)計(jì)

2.5 土地沙化驅(qū)動(dòng)因素分析

利用地理探測(cè)器模型對(duì)選取的驅(qū)動(dòng)因子進(jìn)行離散化處理和分類發(fā)現(xiàn)，研究區(qū)土地沙化受到多種驅(qū)動(dòng)因素的影響。由圖3可知，各個(gè)驅(qū)動(dòng)因子顯著性影響土地沙化的程度，且得出其影響程度的順序?yàn)橥寥李愋?DEM>平均氣溫>平均風(fēng)速>日照時(shí)數(shù)>年降水量>GDP>人口空間分布>土壤侵蝕>坡度。其中，土壤類型、DEM和年平均氣溫3種因子影響程度均超過10%，是主要影響因子。社會(huì)經(jīng)濟(jì)因子中GDP和人口空間分布q值排序靠后，這說明自然因素比社會(huì)經(jīng)濟(jì)因素對(duì)土地沙化影響程度更大。

3 討論

本文基于GEE大數(shù)據(jù)平臺(tái)，通過Landsat遙感影像提取了研究區(qū)NDVI和Albedo數(shù)據(jù)，從土地沙化時(shí)空演變、轉(zhuǎn)移矩陣、地理探測(cè)器模型等方面對(duì)科爾沁左翼中旗土地沙化進(jìn)行監(jiān)測(cè)與驅(qū)動(dòng)因素分析。綜合分析，2000—2021年研究區(qū)土地沙化主要分布在中部和東北部，整體呈穩(wěn)定趨勢(shì)發(fā)展，改善土地面積大于退化土地面積，這與眾多研究科爾沁沙地沙化的研究結(jié)果較相似[14]。土壤類型、DEM和年平均氣溫3種因子是土地沙化的主要影響因子?？茽柷咦笠碇衅焯幱谖覈狈睫r(nóng)牧交錯(cuò)帶，生態(tài)環(huán)境脆弱。從研究方法來看，選用歸一化植被指數(shù)和地表反照率2種地表參量相結(jié)合的方法，綜合了對(duì)植被敏感的可見光和近紅外波段的反射光譜信息，二者擬合的線性函數(shù)具有明顯的負(fù)相關(guān)關(guān)系，此研究結(jié)果與曾永年等[5]研究結(jié)果一致，曾永年等指出，隨著植被指數(shù)的增加，地表反照率減小，土地沙化過程得到了明顯的反映。李嘉豪等[15]利用該指數(shù)模型研究得到小灤河流域土地沙化動(dòng)態(tài)變化，本研究也認(rèn)為，沙漠化差值指數(shù)可作為土地沙化的監(jiān)測(cè)指標(biāo)。地理探測(cè)器模型現(xiàn)已成為主流地理分析模型，利用此模型可定量化地揭示不同驅(qū)動(dòng)因子對(duì)土地沙化的影響大小。本文選取了地理探測(cè)器模型中的因子探測(cè)模塊，后續(xù)應(yīng)結(jié)合其他3種模塊做進(jìn)一步分析。土地沙化是我國當(dāng)前最為嚴(yán)重的生態(tài)環(huán)境問題之一，近年來，隨著退耕還林還草、實(shí)施“山水林田湖草沙一體化保護(hù)和修復(fù)工程”等一系列治理措施，科爾沁沙地治理成效顯著，但仍有小部分區(qū)域改善不明顯。因此，應(yīng)加大治沙防沙治理力度，因地制宜，阻止土地沙化進(jìn)一步發(fā)展，優(yōu)化生態(tài)環(huán)境。

4 結(jié)論

2000—2021年，研究區(qū)沙化土地主要是輕中度沙化且集中在東北部和中部。對(duì)比2015年和2021年數(shù)據(jù)可知，雖然增加了2.80%的重度和極重度的沙化土地，但同時(shí)增加了8.52%的非沙化土地，說明沙化土地有好轉(zhuǎn)的跡象；2000—2021年，研究區(qū)沙化土地變化趨勢(shì)以穩(wěn)定不變?yōu)橹鳎娣e占比46.71%，輕微改善面積占比30.75%，輕微退化面積占比22.52%；通過地理探測(cè)器模型探測(cè)，在自然因素和社會(huì)經(jīng)濟(jì)因素中，研究區(qū)土地沙化主要受到自然因素影響，其中土壤類型是最主要的驅(qū)動(dòng)因子，q值為0.1279，其次是DEM和年平均氣溫，q值分別為0.1244、0.1055。