劉 瓊，羅 沖，孟祥添，張新樂(lè)，唐海濤，馬士耐，劉煥軍

典型黑土區(qū)耕作土壤質(zhì)地遙感時(shí)間窗口及影響因素分析

劉 瓊1，羅 沖2，孟祥添2，張新樂(lè)1※，唐海濤1，馬士耐1，劉煥軍2

（1. 東北農(nóng)業(yè)大學(xué)公共管理與法學(xué)院，哈爾濱 150030；2. 中國(guó)科學(xué)院東北地理與農(nóng)業(yè)生態(tài)研究所，長(zhǎng)春 130012）

了解黑土區(qū)耕作土壤質(zhì)地的空間分布對(duì)于黑土區(qū)農(nóng)業(yè)精準(zhǔn)管理以及耕地保護(hù)至關(guān)重要。遙感技術(shù)是快速獲取土壤質(zhì)地空間分布的有效方法。該研究以黑龍江省友誼農(nóng)場(chǎng)耕地為研究對(duì)象，評(píng)估研究區(qū)土壤質(zhì)地遙感反演的最佳時(shí)間窗口并分析其影響因素。篩選覆蓋研究區(qū)的2019—2021年25幅Sentinel-2影像，將每幅影像的波段和構(gòu)建的光譜指數(shù)輸入隨機(jī)森林模型，建立土壤質(zhì)地遙感反演模型，比較不同時(shí)期影像反演土壤質(zhì)地的模型精度，確定土壤質(zhì)地遙感反演的最適宜影像，并分析造成反演土壤質(zhì)地精度變化的原因，獲取友誼農(nóng)場(chǎng)土壤質(zhì)地空間分布。結(jié)果表明：1）友誼農(nóng)場(chǎng)反演土壤質(zhì)地的最佳時(shí)間窗口為4月下旬至5月中旬；2）在25幅Sentinel-2影像中，2020年5月7日反演粉粒和砂粒的模型精度最高（粉粒的2為0.785，均方根誤差為6.697%；砂粒的2為0.776，均方根誤差為8.296%）；2019年5月3日反演黏粒的模型精度最高（2為0.776，均方根誤差為1.600%）；3）不同時(shí)期的Sentinel-2影像對(duì)土壤質(zhì)地反演的準(zhǔn)確性有很大的影響，而土壤含水率和秸稈覆蓋是造成不同時(shí)期土壤質(zhì)地預(yù)測(cè)精度差異的重要原因。研究為確定土壤質(zhì)地遙感反演的最佳時(shí)間窗口、實(shí)現(xiàn)區(qū)域尺度土壤質(zhì)地制圖提供關(guān)鍵技術(shù)。

土壤；隨機(jī)森林；遙感；黑土；時(shí)間窗口；Sentinel-2；影響因素；質(zhì)地

0 引 言

土壤質(zhì)地是影響土壤理化性質(zhì)形成的重要因素，特別是對(duì)土壤結(jié)構(gòu)、孔隙狀況、保肥性、保水性、耕性等均有重要影響[1-2]。對(duì)黑土區(qū)耕作土壤質(zhì)地進(jìn)行測(cè)定和分析是國(guó)家黑土地保護(hù)工程實(shí)施方案的重要內(nèi)容[3]。研究黑土區(qū)耕層土壤質(zhì)地的空間分布不但有利于分析黑土區(qū)耕地產(chǎn)出能力，還可為黑土區(qū)耕地保護(hù)提供準(zhǔn)確、定量的數(shù)據(jù)支持。

地統(tǒng)計(jì)分析方法是土壤質(zhì)地空間預(yù)測(cè)與制圖的重要方法。地統(tǒng)計(jì)學(xué)是基于區(qū)域化變量理論和變異函數(shù)模型，以具有空間相關(guān)性和依賴性且在空間分布上既有隨機(jī)性又有結(jié)構(gòu)性的自然現(xiàn)象為研究對(duì)象的科學(xué)[4]。但是其反演精度往往受到樣本密度和樣本變異程度的限制，只有具有足夠數(shù)量和密集的樣本才能獲得準(zhǔn)確的土壤質(zhì)地空間分布圖[5]。近年來(lái)，遙感數(shù)據(jù)以更直接、更具成本效益和更快速的方式被用來(lái)估算土壤的重要指標(biāo)[6]。多光譜和高光譜衛(wèi)星影像已被證明是反演土壤質(zhì)地的一種有效手段[7]。這為土壤質(zhì)地遙感反演的發(fā)展奠定了基礎(chǔ)。土壤質(zhì)地遙感反演的方法是在研究土壤反射光譜的基礎(chǔ)上，通過(guò)地面實(shí)測(cè)數(shù)據(jù)，分析采樣點(diǎn)遙感影像光譜反射率的變化，建立土壤質(zhì)地反演模型，實(shí)現(xiàn)土壤質(zhì)地的空間預(yù)測(cè)和制圖。目前，已有很多研究運(yùn)用普通克里金[8-9]、協(xié)同克里金[10]、回歸克里金[11]、偏最小二乘回歸模型[12-13]、多元線性回歸[14]、回歸樹[15]、人工神經(jīng)網(wǎng)絡(luò)[16]進(jìn)行土壤質(zhì)地的空間預(yù)測(cè)分析和制圖。隨機(jī)森林（Random Forest，RF）算法是通過(guò)集成學(xué)習(xí)的思想將多棵樹集成的一種算法。RF模型具有較高的預(yù)測(cè)精度，隨機(jī)性的引入使得RF不容易陷入過(guò)擬合，并且具有很好的抗噪聲能力。

已有研究將RF模型應(yīng)用到土壤質(zhì)地的空間預(yù)測(cè)中[17]。Lie等[15]比較了RF模型和回歸樹模型反演土壤質(zhì)地的精度，并表明RF模型具有更高的反演精度。da Silva等[14]利用RF模型和多元線性回歸模型分別對(duì)半干旱地區(qū)土壤質(zhì)地的空間分布進(jìn)行預(yù)測(cè)，發(fā)現(xiàn)RF模型有更高的精度，Landsat 5的波段2、波段5、波段3和波段7的比值以及歸一化差值植被指數(shù)（Normalized Difference Vegetation Index，NDVI）為重要的環(huán)境變量。Hengl等[18]使用多元線性回歸和RF模型預(yù)測(cè)了非洲大陸的土壤質(zhì)地空間分布，并認(rèn)為RF模型能得出最優(yōu)的結(jié)果。Pahlavan-Rad等[19]使用RF模型分析了Zahak County地區(qū)土壤質(zhì)地和pH值的空間分布格局，發(fā)現(xiàn)該地區(qū)土壤質(zhì)地空間變異大，RF模型反演的土壤質(zhì)地和pH值結(jié)果精度高。

現(xiàn)有的土壤質(zhì)地遙感反演模型研究一般利用遙感影像，結(jié)合不同的景觀環(huán)境因子作為輸入量，探究其對(duì)土壤質(zhì)地反演精度的影響。但即使在同一個(gè)研究區(qū)獲得的不同遙感影像會(huì)產(chǎn)生不同的土壤質(zhì)地反演精度。這主要是由于影像可能受到降水、秸稈覆蓋、表面形態(tài)等因素的影響，影像特征反映的是部分地區(qū)的異常情況，導(dǎo)致遙感影像部分地區(qū)地物反射率異常，降低土壤質(zhì)地反演模型的穩(wěn)定性與精度[20]。為了更好地探究耕地土壤質(zhì)地的反演精度，有必要確定一個(gè)特定的時(shí)間窗口，在該窗口上，可以獲得相對(duì)較好的光譜數(shù)據(jù)，而基于此建立的光譜模型可以更準(zhǔn)確地預(yù)測(cè)土壤質(zhì)地。因此，本研究以友誼農(nóng)場(chǎng)耕地為研究對(duì)象，采用了高空間分辨率的Sentinel-2衛(wèi)星遙感影像作為數(shù)據(jù)源，基于RF算法，將波段和構(gòu)建的光譜指數(shù)作為輸入量，根據(jù)反演的精度確定友誼農(nóng)場(chǎng)土壤質(zhì)地反演的最佳時(shí)間窗口，分析其造成反演精度變化的原因，實(shí)現(xiàn)土壤質(zhì)地空間預(yù)測(cè)與制圖，以期為確定土壤質(zhì)地遙感反演的最佳時(shí)間窗口和實(shí)現(xiàn)區(qū)域尺度土壤質(zhì)地制圖提供關(guān)鍵技術(shù)。

1 材料與方法

1.1 研究區(qū)概況

友誼農(nóng)場(chǎng)是國(guó)家重要的商品糧基地之一和現(xiàn)代化農(nóng)業(yè)示范窗口，位于雙鴨山市境內(nèi)，地理坐標(biāo)為131°27′～132°15′E，46°28′～46°59′N，地處黑龍江、松花江、烏蘇里江三江平原腹地，內(nèi)地地勢(shì)相對(duì)平坦，西南為丘陵，東北為低洼地，由西南向東北傾斜。研究區(qū)耕地范圍面積為1 800 km2，此區(qū)域土壤肥沃，物理結(jié)構(gòu)好，主要土壤類型有草甸土、黑土和沼澤土。年平均氣溫2.5 ℃，年平均降雨量500 mm，適宜各類農(nóng)作物生長(zhǎng)，以種植水稻、大豆、玉米和麥類為主，作物一年一季。圖1為友誼農(nóng)場(chǎng)裸土影像及采樣點(diǎn)空間位置分布。

圖1 研究區(qū)土壤采樣點(diǎn)分布圖

1.2 數(shù)據(jù)獲取與處理

1.2.1 土壤數(shù)據(jù)獲取

由于土壤質(zhì)地在1～10 a的較小時(shí)間尺度內(nèi)不會(huì)發(fā)生顯著變化[21]，所以土壤樣本的野外采集、影像數(shù)據(jù)的獲取和土壤質(zhì)地分析可以隨著時(shí)間的推移而單獨(dú)進(jìn)行。本研究所使用的耕地土壤樣本于2021年3月22日—4月5日進(jìn)行野外采集，此時(shí)，土壤直接暴露于表面，既無(wú)農(nóng)作物殘留、又無(wú)積雪覆蓋。在研究區(qū)范圍內(nèi)進(jìn)行GPS定位，樣本點(diǎn)的位置選擇視野開(kāi)闊、地勢(shì)相對(duì)平坦、地質(zhì)條件均一的耕作土壤，且兼顧了地形特征和不同土類等信息，共計(jì)獲得188個(gè)野外實(shí)測(cè)樣本，記錄其坐標(biāo)和地形等相關(guān)信息，其分布圖如圖1所示。將獲取的土壤樣本置于布制土袋中存放，帶回實(shí)驗(yàn)室進(jìn)行風(fēng)干處理，并去除結(jié)石、雜草根和其他雜質(zhì)，室內(nèi)將土壤稱質(zhì)量后研磨、風(fēng)干，過(guò)2 mm篩，使用Malvern MS-2000激光粒度分析儀（英國(guó)馬爾文儀器有限公司生產(chǎn)，殘差小于1%）測(cè)定樣本的土壤顆粒組成（質(zhì)量分?jǐn)?shù)）。表1為188個(gè)采樣點(diǎn)土壤顆粒組成統(tǒng)計(jì)結(jié)果。

表1 土壤顆粒組成統(tǒng)計(jì)結(jié)果

1.2.2 影像數(shù)據(jù)介紹與處理

如表2所示，本研究選取了2019—2021年共25幅可用的Sentinel-2的L2A產(chǎn)品（https://code.earthengine. google.com/）進(jìn)行土壤質(zhì)地反演研究。Sentinel-2的L2A產(chǎn)品為經(jīng)過(guò)幾何校正和大氣處理的地表反射率數(shù)據(jù)，其中，波段1為氣溶膠波段，波段9為水蒸氣波段，波段10為反應(yīng)大氣波段[22]，故在后續(xù)的遙感反演過(guò)程中去除這3個(gè)波段。然后在ENVI5.3軟件中根據(jù)GlobeLand30- 2020（http://www.globallandcover.com/）中的耕地范圍對(duì)遙感影像進(jìn)行裁剪和鑲嵌，并將Sentinel-2影像空間重采樣為10 m分辨率。

表2 2019—2021年Sentinel-2影像獲取時(shí)間

1.3 光譜指數(shù)的構(gòu)建

由于多光譜影像中的波段很少，光譜指數(shù)的構(gòu)建可以充分利用多光譜影像中的有用信息[23]。差值指數(shù)（Difference Index，DI）、比值指數(shù)（Ratio Index，RI）、歸一化差值指數(shù)（Normalized Difference Index，NDI）和其他光譜指數(shù)的構(gòu)建能減少水分、粗糙度和大氣的影響[24]。因此，對(duì)Sentinel-2影像數(shù)據(jù)波段反射率進(jìn)行數(shù)學(xué)變換，計(jì)算式如下

式中ρ和ρ為和波段的反射率，其中>。本研究共生成了135個(gè)光譜指數(shù)（45個(gè)差值指數(shù)，45個(gè)比值指數(shù)，45個(gè)歸一化差值指數(shù)）。

1.4 隨機(jī)森林反演模型

隨機(jī)森林模型是基于決策樹的非線性機(jī)器學(xué)習(xí)模型。RF建模通過(guò)R語(yǔ)言中的Random Forest包實(shí)現(xiàn)。Random Forest包在運(yùn)算過(guò)程中對(duì)影響土壤質(zhì)地的自變量重要性進(jìn)行排序，其數(shù)值越大，代表自變量對(duì)土壤質(zhì)地影響越大，相關(guān)性越強(qiáng)。此外，RF模型具有準(zhǔn)確性高、無(wú)需剪枝、較少出現(xiàn)過(guò)擬合現(xiàn)象、能容忍一定的干擾和異常值、訓(xùn)練速度較快等優(yōu)點(diǎn)[25]。本研究在RF實(shí)際建模操作中，通過(guò)觀察袋外誤差選擇最佳回歸樹的數(shù)量（ntree）和分裂節(jié)點(diǎn)數(shù)（mtry），從而建立最優(yōu)RF預(yù)測(cè)模型。經(jīng)多次訓(xùn)練，確定ntree為500，ntry為輸入量個(gè)數(shù)的1/3，此時(shí)，可以產(chǎn)生較穩(wěn)定的袋外誤差率，模型較為穩(wěn)定。

1.5 模型建立與驗(yàn)證

本研究采用RF算法模型，使用不同時(shí)期Sentinel-2影像數(shù)據(jù)的10個(gè)波段與135個(gè)光譜指數(shù)作為輸入量，分別進(jìn)行土壤質(zhì)地反演建模。188個(gè)樣本中按照3∶1的比例劃分建模集和驗(yàn)證集，其中建模樣本141個(gè)，驗(yàn)證樣本47個(gè)。決定系數(shù)（2）用于檢驗(yàn)?zāi)Ｐ头€(wěn)定性。均方根誤差（Root Mean Square Error，RMSE）用于衡量模型的精度[26]。通常，2值在0到1的范圍內(nèi)，數(shù)值的絕對(duì)值越接近1，預(yù)測(cè)值越接近實(shí)際測(cè)量值；RMSE是用來(lái)檢驗(yàn)土壤質(zhì)地實(shí)際測(cè)量值和預(yù)測(cè)值之間差異程度的依據(jù)，其數(shù)值越小，則表明兩者的差異程度越小，模型的反演精度越好。

1.6 土壤質(zhì)地遙感反演的影響因素

為分析土壤含水率和秸稈覆蓋對(duì)土壤質(zhì)地反演造成的影響，本研究采用地表水分指數(shù)（Land Surface Water Index，LSWI）和歸一化耕作指數(shù)（Normalized Difference Tillage Index，NDTI）分別表示土壤含水率和秸稈覆蓋。計(jì)算式為：

LSWI=(nir-swir1)/(nir+swir1) （4）

式中nir為近紅外波段的地表反射率；swir1為短波紅外1波段（11）的地表反射率。

NDTI=(swir1-swir2)/(swir1+swir2) （5）

式中swir2為短波紅外2波段（12）的地表反射率。

地表水分指數(shù)包含對(duì)土壤濕度和植被水分敏感的短波紅外波段反射率，可以用于土壤濕度變化的監(jiān)測(cè)，濕度越大，LSWI值越大[27]。NDTI越大秸稈覆蓋程度越高，當(dāng)NDTI小于0.05時(shí)，基本無(wú)秸稈覆蓋；當(dāng)NDTI為0.05～0.1時(shí)，秸稈覆蓋度在10%～20%之間；當(dāng)NDTI大于0.1時(shí)，秸稈覆蓋度大于20%[28]。

2 結(jié)果與分析

2.1 星載土壤反射光譜特征分析

土壤光譜反射率與土壤質(zhì)地的表土粒徑有著密切的關(guān)系。圖2為同一時(shí)期土壤顆粒組成的反射光譜特征。結(jié)果表明，不同顆粒下的Sentinel-2光譜曲線特征相似，但光譜反射率存在顯著的差異。從可見(jiàn)光到短波紅外1（1 565～1 655 nm）光譜反射率隨著波長(zhǎng)的增加而增加，從短波紅外1到短波紅外2（2 100～2 280 nm）光譜反射率隨著波長(zhǎng)的增加而降低。其中，光譜反射率隨著粉粒和黏粒質(zhì)量分?jǐn)?shù)的增加而減少，隨著砂粒質(zhì)量分?jǐn)?shù)的增加而增加。這是由于不同粒徑顆粒的混合效應(yīng)[29]。此外，不僅土壤顆粒的大小本身對(duì)土壤光譜產(chǎn)生影響，而且不同比例的礦物也會(huì)產(chǎn)生光譜反射差異。黏粒中含有的層狀硅酸鹽會(huì)降低土壤反射率；而砂粒中含有的大量石英會(huì)增加光譜反射率[30]。

a. 粉粒a. Siltb. 黏粒b. Clayc. 砂粒c. Sand

注：波段序號(hào)1～10分別為B2、B3、B4、B5、B6、B7、B8、B8A、B11、B12。

Note: No. of band of 1-10 is B2, B3, B4, B5, B6, B7, B8, B8A, B11, B12, respectively.

圖2 不同質(zhì)地的土壤反射光譜特征

Fig.2 Characteristics of soil reflectance spectrum with different textures

2.2 反演土壤質(zhì)地的最佳時(shí)間窗口

盡管在使用遙感方法反演土壤質(zhì)地時(shí)，通常會(huì)排除NDVI高于0.4的光譜數(shù)據(jù)，但是NDVI高于0.4的光譜數(shù)據(jù)也可能反映土壤質(zhì)地[31]。因此，可以適當(dāng)?shù)睦眠@些光譜數(shù)據(jù)反演土壤質(zhì)地。將2019—2021年獲取的單期Sentinel-2影像10個(gè)波段反射率和135個(gè)光譜指數(shù)作為輸入量，輸入RF模型得到不同時(shí)期影像反演土壤質(zhì)地的精度。表3為2019—2021年不同時(shí)期Sentinel-2影像反演土壤質(zhì)地的驗(yàn)證模型精度。從圖中可以看出，在2019—2021年中，每年的土壤質(zhì)地反演精度的變化呈現(xiàn)了大致相同的趨勢(shì)：增加-降低-再增加。反演精度最低的時(shí)期均在7—8月。2020年5月7日反演粉粒（2=0.785）和砂粒（2=0.776）的精度最高，2019年5月3日反演黏粒（2=0.776）的精度最高。因?yàn)榉囱菥容^高的4個(gè)時(shí)期（2019年5月3日、2020年4月29日和5月7日、2021年5月17日）在4月下旬至5月中旬這個(gè)時(shí)間范圍內(nèi)，所以本研究認(rèn)為該時(shí)間窗口為友誼農(nóng)場(chǎng)反演土壤質(zhì)地的最佳時(shí)間窗口。

2.3 最佳輸入量的分析

表4為在土壤質(zhì)地反演中，基于RF模型的前10個(gè)輸入量的重要性排序。如表所示，Sentinel-2影像的波段12（短波紅外2）為遙感反演粉粒、黏粒和砂粒的最佳波段。這與Shahriari等[17]的研究結(jié)果相似。該波段不僅包含許多土壤礦物的信息，而且與土壤含水率緊密相關(guān)。此外，光譜指數(shù)DI12_2、DI11_2、DI12_3、NDI4_2和NDI8_2都對(duì)土壤質(zhì)地的遙感反演具有重要意義。

表3 2019—2021年不同時(shí)期土壤質(zhì)地反演精度對(duì)比

2.4 影響土壤質(zhì)地反演精度的因素

通過(guò)對(duì)2019—2021年不同日期的Sentinel-2影像反演土壤質(zhì)地的精度分析，發(fā)現(xiàn)不同日期的耕地環(huán)境差異較大。這就造成了不同時(shí)期的Sentinel-2影像對(duì)土壤質(zhì)地反演的準(zhǔn)確性有很大的影響，秸稈覆蓋或土壤含水率有可能是造成不同時(shí)期土壤質(zhì)地反演精度差異的重要原因。

表4 土壤質(zhì)地隨機(jī)森林反演中輸入量的重要性排序

注：DI12_2為波段12和波段2的差值指數(shù)。NDI4_2為波段4和波段2的歸一化差值指數(shù)。RI4_2為波段4和波段2的比值指數(shù)。12為波段12的差值指數(shù)。其他類似解讀。

Note: DI12_is the difference index between band 12 and 2. NDI4_2is normalized difference index between band 4 and 2. RI4_2is ratio index between band 4 and 2.12is band 12. The others are explained as the same way.

如圖3所示，通過(guò)計(jì)算2019—2021年25幅影像的LSWI值，發(fā)現(xiàn)粉粒、黏粒和砂粒反演過(guò)程中不同程度的受到了土壤含水率的影響，隨著LSWI的增加，RMSE整體呈現(xiàn)一個(gè)增加的趨勢(shì)。在圖4中，2019年的8幅影像中，LSWI值最低的3幅影像分別為11月4日、5月3日和3月29日；2020年的7幅影像中，LSWI值最低的3幅影像為4月29日、10月26日和5月7日；2021年的10幅影像中，LSWI值最低的3幅影像為10月29日、5月17日和4月19日。反演土壤質(zhì)地的最佳的Sentinel-2影像具有較低的土壤含水率。

圖3 土壤質(zhì)地反演的均方根誤差與地表水分指數(shù)、歸一化耕作指數(shù)的關(guān)系

計(jì)算25幅Sentinel-2影像的NDTI值后，發(fā)現(xiàn)秸稈覆蓋度影響土壤質(zhì)地的遙感反演，如圖3d～3f，NDTI與RMSE呈現(xiàn)正向相關(guān)關(guān)系。在圖4中，2019年NDTI最低的3幅影像日期為3月29日、3月6日和11月4日，其中5月3日的秸稈覆蓋度仍在0.2以下；2020年NDTI最低的3幅影像為4月29日、10月26日和5月7日；2021年NDTI最低的3幅影像為5月17日、10月29日和4月19日。反演土壤質(zhì)地最佳的Sentinel-2影像秸稈覆蓋度低。

a. 2019b. 2020c. 2021

綜上所述，在最佳時(shí)間窗口中的4幅影像中，其土壤含水率和秸稈覆蓋度都比較低，且隨著土壤含水率和秸稈覆蓋的增加，土壤質(zhì)地遙感反演精度逐漸降低，現(xiàn)認(rèn)為土壤含水率和秸稈覆蓋度是影響土壤質(zhì)地反演精度的重要原因。土壤質(zhì)地與可見(jiàn)光近紅外波段-短波紅外波段（400～2 500 nm）區(qū)域的土壤反射率顯著相關(guān)。而在這個(gè)光譜范圍內(nèi)，土壤水分和秸稈覆蓋會(huì)影響土壤質(zhì)地遙感反演精度。Gomez等[31]的研究發(fā)現(xiàn)，由于土壤濕度的變化，2017年4月24和2017年4月4日的Sentinel-2數(shù)據(jù)在對(duì)土壤質(zhì)地進(jìn)行反演時(shí)產(chǎn)生了最高和最低的精度。Swain等[32]在使用來(lái)自不同采集日期的63幅Sentinel-2進(jìn)行土壤質(zhì)地空間預(yù)測(cè)時(shí)，發(fā)現(xiàn)土壤含水率或植被覆蓋影響了土壤質(zhì)地反演的精度。Luo等[33]在使用Landsat 8 遙感影像對(duì)土壤有機(jī)質(zhì)進(jìn)行反演的研究中也發(fā)現(xiàn)，秸稈覆蓋度和土壤含水率是造成不同時(shí)期土壤有機(jī)質(zhì)預(yù)測(cè)精度差異的原因。

2.5 土壤質(zhì)地的空間分布圖

采用2019年5月3日的Sentinel-2反演黏粒的空間分布，2020年5月7日的Sentinel-2反演粉粒和砂粒的空間分布，得到了最高的預(yù)測(cè)精度，如圖5a～5c所示。由于選取的采樣點(diǎn)均為耕地土壤樣本，因此利用耕地邊界范圍對(duì)影像進(jìn)行裁剪，并反演土壤質(zhì)地空間分布圖。從本研究的土壤質(zhì)地遙感反演結(jié)果可以看出友誼農(nóng)場(chǎng)粉粒和黏粒在東北部、北部和南部分布較多，中部和西南分布較少；砂粒則呈現(xiàn)出相反的趨勢(shì)，尤其是友誼農(nóng)場(chǎng)中部的砂粒含量整體較高。

圖5 10與250 m分辨率的土壤質(zhì)地反演結(jié)果的比較

3 討 論

本文分析了土壤含水率和秸稈覆蓋影響土壤質(zhì)地的遙感反演。根據(jù)圖4發(fā)現(xiàn)2019年3月29日和2019年3月6日的秸稈覆蓋低而土壤含水率相對(duì)較高；2021年4月4日的土壤含水率低而秸稈覆蓋度高。上述兩種情況的土壤質(zhì)地遙感反演的精度不佳。而在反演精度較高的4幅影像中，土壤含水率和秸稈覆蓋都比較低?？梢?jiàn)，只有同時(shí)具備較低的土壤含水率和秸稈覆蓋度的影像才能產(chǎn)生較優(yōu)土壤質(zhì)地反演精度。2019年11月4日、2020年10月26日和2021年10月29日的土壤含水率和秸稈覆蓋都比較低，但是反演的精度仍然不佳，這可能是由于土壤表面粗糙度影響了反射率的變化。粗糙度越高，預(yù)測(cè)的性能越低。冬季由于深耕作業(yè)造成地表具有明顯的粗糙度；與冬季不同，由于春季作物播種的淺耕作業(yè)，土壤表面粗糙度則比較低[34]。在最佳時(shí)間窗口中的2020年4月29日的土壤含水率和秸稈覆蓋是最低的，但是反演精度卻不是最高的，這可能是受大氣干擾等其他的因素影響。

在土壤質(zhì)地遙感反演的結(jié)果中，粉粒的2范圍為0.341～0.785；黏粒的2范圍為0.200～0.776；砂粒的2范圍為0.308～0.776。這表明可用于土壤質(zhì)地反演的遙感影像隨著不同采集日期的地表形態(tài)的變化而變化。干擾反射光譜的任何因素都可能會(huì)影響土壤質(zhì)地的遙感反演。光譜干擾可能來(lái)自于以下方面：1）土壤水分、植被覆蓋以及地表的作物殘留等；2）產(chǎn)生不同光效應(yīng)的土壤表面形態(tài)（粗糙度）；3）隨季節(jié)、云層、太陽(yáng)方位角和海拔變化的大氣干擾會(huì)對(duì)野外反射信號(hào)產(chǎn)生影響。

在圖5中，將本研究中的10 m分辨率的土壤質(zhì)地空間分布圖與Tomislav等[35-37]發(fā)布的基于機(jī)器學(xué)習(xí)的250 m分辨率的全球土壤質(zhì)地分布結(jié)果進(jìn)行對(duì)比。根據(jù)粉粒含量的空間分布結(jié)果（圖5d），本研究發(fā)現(xiàn)友誼農(nóng)場(chǎng)粉粒主要分布在東北部、北部和南部小部分地區(qū)，而Tomislav等研究的粉粒分布結(jié)果則表現(xiàn)出相反的趨勢(shì)，但根據(jù)采樣點(diǎn)數(shù)據(jù)，友誼農(nóng)場(chǎng)的東北部、北部和南部的粉粒含量都比較高。根據(jù)圖5e可以看出，友誼農(nóng)場(chǎng)黏粒含量至西南向東北逐漸升高，和Tomislav等的研究一樣，本研究發(fā)現(xiàn)友誼農(nóng)場(chǎng)西南部黏粒的含量低，但是由于中部地區(qū)的土壤質(zhì)地組成以砂粒為主，所以黏粒含量會(huì)相對(duì)比較低。在圖5f中，友誼農(nóng)場(chǎng)西南部砂粒的含量高，中部地區(qū)的砂粒含量低。而在本研究砂粒反演結(jié)果中（圖5c），除了西南部砂粒的含量較高以外，中部地區(qū)存在一個(gè)沙地圈。根據(jù)采樣點(diǎn)數(shù)據(jù)，中部地區(qū)土壤的含沙量高，本研究的砂粒反演結(jié)果也更準(zhǔn)確。

4 結(jié) 論

本研究基于Sentinel-2影像數(shù)據(jù)，運(yùn)用隨機(jī)森林算法，將波段和構(gòu)建的光譜指數(shù)作為輸入量，根據(jù)隨機(jī)森林反演模型精度逐步確定土壤質(zhì)地空間預(yù)測(cè)的最佳時(shí)間窗口，對(duì)區(qū)域土壤質(zhì)地的空間分布進(jìn)行預(yù)測(cè)，得到如下結(jié)論：

1）在25幅Sentinel-2影像中，2020年5月7日反演粉粒和砂粒的模型精度最高（2分別為0.785和0.776）；2019年5月3日反演黏粒的模型精度最高（2=0.776）。

2）選擇合適時(shí)間的衛(wèi)星影像對(duì)土壤質(zhì)地反演有著重要的影響，友誼農(nóng)場(chǎng)反演土壤質(zhì)地的最佳時(shí)間窗口為4月下旬—5月中旬。在此時(shí)間窗口不僅可以獲得質(zhì)量較好的光譜數(shù)據(jù)，而且能開(kāi)發(fā)穩(wěn)定的光譜模型來(lái)反演土壤質(zhì)地的空間分布。

3）根據(jù)2019—2021年的25幅Sentinel-2影像，得出了不同的土壤質(zhì)地反演精度結(jié)果，而土壤含水率和秸稈覆蓋是造成這種結(jié)果的重要原因。

本研究的土壤質(zhì)地遙感反演模型精度高，可為區(qū)域尺度土壤質(zhì)地制圖以及黑土區(qū)耕地保護(hù)提供關(guān)鍵技術(shù)。

[1] Liao K, Xu S, Wu J, et al. Spatial estimation of surface soil texture using remote sensing data[J]. Soil Science and Plant Nutrition, 2013, 59(4): 488-500.

[2] 王德彩，鄔登巍，趙明松，等. 平原區(qū)土壤質(zhì)地的反射光譜預(yù)測(cè)與地統(tǒng)計(jì)制圖[J]. 土壤通報(bào)，2012，43(2)：257-262.

Wang Decai, Wu Dengwei, Zhao Mingsong, et al. Prediction and mapping of soil texture of a plain area using reflectance spectra and geo-statistics[J]. Chinese Journal of Soil Science, 2012, 43(2): 257-262. (in Chinese with English abstract)

[3] 中華人民共和國(guó)農(nóng)業(yè)農(nóng)村部. 國(guó)家黑土地保護(hù)工程實(shí)施方案(2021—2025年)[J]. 中國(guó)農(nóng)業(yè)綜合開(kāi)發(fā)，2021(8)：4-11.

[4] Sullivan D G, Shaw J N, Rickman D. IKONOS imagery to estimate surface soil property variability in two Alabama physiographies[J]. Soil Science Society of America Journal, 2005, 69(6): 1789-1798.

[5] 劉煥軍，張美薇，楊昊軒，等. 多光譜遙感結(jié)合隨機(jī)森林算法反演耕作土壤有機(jī)質(zhì)含量[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2020，36(10)：134-140.

Liu Huanjun, Zhang Meiwei, Yang Haoxuan, et al. Invertion of cultivated soil organic matter content combining multi-spectral remote sensing and random forest algorithm[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(10): 134-140. (in Chinese with English abstract)

[6] Angelopoulou T, Tziolas N, Balafoutis A, et al. Remote sensing techniques for soil organic carbon estimation: A review[J]. Remote Sensing, 2019, 11(6): 676.

[7] Gomez C, Lagacherie P, Coulouma G. Regional predictions of eight common soil properties and their spatial structures from hyperspectral Vis-NIR data[J]. Geoderma, 2012, 189: 176-185.

[8] Hengl T, Heuvelink G, Rossiter D G. About regression-kriging: From equations to case studies[J]. Computers & Geosciences, 2007, 33(10): 1301-1315.

[9] Piccini C, Marchetti A, Francaviglia R. Estimation of soil organic matter by geostatistical methods: Use of auxiliary information in agricultural and environmental assessment[J]. Ecological Indicators, 2014, 36: 301-314.

[10] Delbari M, Afrasiab P, Loiskandl W. Geostatistical analysis of soil texture fractions on the field scale[J]. Soil & Water Research, 2011, 6(4): 173-189.

[11] Zhu Q, Lin H S. Comparing ordinary kriging and regression kriging for soil properties in contrasting landscapes[J]. Pedosphere, 2010, 20(5): 594-606.

[12] Vaudour E, Gomez C, Fouad Y, et al. Sentinel-2 image capacities to predict common topsoil properties of temperate and Mediterranean agroecosystems[J]. Remote Sensing of Environment, 2019, 223: 21-33.

[13] Gomez C, Oltra-Carrió R, Bacha S, et al. Evaluating the sensitivity of clay content prediction to atmospheric effects and degradation of image spatial resolution using Hyperspectral VNIR/SWIR imagery[J]. Remote Sensing of Environment, 2015, 164: 1-15.

[14] da Silva Chagas C, de Carvalho Junior W, Bhering S B, et al. Spatial prediction of soil surface texture in a semiarid region using random forest and multiple linear regressions[J]. Catena, 2016, 139: 232-240.

[15] Lie M, Glaser B, Huwe B. Uncertainty in the spatial prediction of soil texture: Comparison of regression tree and random forest models[J]. Geoderma, 2012, 170(3): 70-79.

[16] Moonjun R, Farshad A, Shrestha D P, et al. Artificial neural network and decision tree in predictive soil mapping of Hoi Num Rin sub-watershed, Thailand[M]. Digital Soil Mapping: Netherland; Springer, 2010: 151-164.

[17] Shahriari M, Delbari M, Afrasiab P, et al. Predicting regional spatial distribution of soil texture in floodplains using remote sensing data: A case of southeastern Iran[J]. Catena, 2019, 182: 104149.

[18] Hengl T, Heuvelink G B M, Kempen B, et al. Mapping soil properties of Africa at 250 m resolution: Random forests significantly improve current predictions[J]. PloS One, 2015, 10(6): e0125814.

[19] Pahlavan-Rad M R, Akbarimoghaddam A. Spatial variability of soil texture fractions and pH in a flood plain (case study from eastern Iran)[J]. Catena, 2018, 160: 275-281.

[20] Vaudour E, Baghdadi N, Gilliot J M. Mapping tillage operations over a peri-urban region using combined SPOT4 and ASAR/ENVISAT images[J]. International Journal of Applied Earth Observation and Geoinformation, 2014, 28: 43-59.

[21] Lehmann J, Bossio D A, K?gel-Knabner I, et al. The concept and future prospects of soil health[J]. Nature Reviews Earth & Environment, 2020, 1(10): 544-553.

[22] Gao B C, Goetz A, Wiscombe W J. Cirrus cloud detection from airborne imaging spectrometer data using the 1. 38 μm water vapor band[J]. Geophysical Research Letters, 1993, 20(4): 301-304.

[23] Dou X, Wang X, Liu H J, et al. Prediction of soil organic matter using multi-temporal satellite images in the Songnen Plain, China[J]. Geoderma, 2019, 356: 113896.

[24] Meng X T, Bao Y L, Liu H J, et al. Regional soil organic carbon prediction model based on a discrete wavelet analysis of hyperspectral satellite data[J]. International Journal of Applied Earth Observation and Geoinformation, 2020, 89: 102111.

[25] 黃翀，許照鑫，張晨晨，等. 基于Sentinel-1數(shù)據(jù)時(shí)序特征的熱帶地區(qū)水稻種植結(jié)構(gòu)提取方法[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2020，36(9)：177-184.

Huang Chong, Xu Zhaoxin, Zhang chenchen, et al. Extraction of rice planting structure in tropical region based on Sentinel-1 temporal features integration[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(9): 177-184. (in Chinese with English abstract)

[26] 唐海濤，孟祥添，蘇循新，等. 基于CARS算法的不同類型土壤有機(jī)質(zhì)高光譜預(yù)測(cè)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2021，37(2)：105-113.

Tang Haitao, Meng Xiangtian, Su Xunxin, et al. Hyperspectral prediction on soil organic matter of different types using CARS algorithm[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(2): 105-113. (in Chinese with English abstract)

[27] 劉丹，馮銳，于成龍，等. 基于植被指數(shù)的春玉米干旱響應(yīng)遙感監(jiān)測(cè)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2019，35(20)：152-161.

Liu Dan, Feng Rui, Yu Chenglong, et al. Remote sensing monitoring of drought response of spring maize based on vegetation indexes[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(20): 152-161. (in Chinese with English abstract)

[28] Jin X, Ma J, Wen Z, et al. Estimation of maize residue cover using Landsat-8 OLI image spectral information and textural features[J]. Remote Sensing, 2015, 7(11): 14559-14575.

[29] Xiao J, Shen Y, Tateishi R, et al. Development of topsoil grain size index for monitoring desertification in arid land using remote sensing[J]. International Journal of Remote Sensing, 2006, 27(12/14): 2411-2422.

[30] Conforti M, Matteucci G, Buttafuoco G. Using laboratory Vis-NIR spectroscopy for monitoring some forest soil properties[J]. Journal of Soils and Sediments, 2018, 18(3): 1009-1019.

[31] Gomez C, Dharumarajan S, Féret J B, et al. Use of Sentinel-2 time-series images for classification and uncertainty analysis of inherent biophysical property: Case of Soil Texture Mapping[J]. Remote Sensing, 2019, 11(5): 565.

[32] Swain S R, Chakraborty P, Panigrahi N, et al. Estimation of soil texture using Sentinel-2 multispectral imaging data: An ensemble modeling approach[J]. Soil and Tillage Research, 2021, 213: 105134.

[33] Luo C, Zhang X, Meng X, et al. Regional mapping of soil organic matter content using multitemporal synthetic Landsat 8 images in Google Earth Engine[J]. Catena, 2022, 209: 105842.

[34] Vaudour E, Gomez C, Loiseau T, et al. The impact of acquisition date on the prediction performance of topsoil organic carbon from Sentinel-2 for croplands[J]. Remote Sensing, 2019, 11(18): 2143.

[35] Tomislav H. Silt content in % (kg/kg) at 6 standard depths (0, 10, 30, 60, 100 and 200 cm) at 250 m resolution (v0. 2)[EB/OL]. Zenodo, 2018. (2018-12-24) [2022-05-10]. https://doi.org/10.5281/zenodo.2525676

[36] Tomislav H. Clay content in % (kg/kg) at 6 standard depths (0, 10, 30, 60, 100 and 200 cm) at 250 m resolution (v0. 2)[EB/OL]. Zenodo, 2018. (2018-11-02) [2022-05-10]. https://doi.org/10.5281/zenodo.2525663

[37] Tomislav H. Sand content in % (kg/kg) at 6 standard depths (0, 10, 30, 60, 100 and 200 cm) at 250 m resolution (v0. 2)[EB/OL]. Zenodo, 2018. (2018-12-24)[2022-05-10]. https://doi.org/10.5281/zenodo.2525662

Time window and influencing factors analysis of tillage soil texture remote sensing in the typical black soil region

Liu Qiong1, Luo Chong2, Meng Xiangtian2, Zhang Xinle1※, Tang Haitao1, Ma Shinai1, Liu Huanjun2

(1.,,150030,;2.,,130012,)

Spatial distribution of cultivated soil texture is generally crucial to the precision management of agriculture and farmland protection in black soil areas. Remote sensing technology can be an effective way to rapidly acquire the spatial distribution of soil texture. Taking the Youyi Farm in Heilongjiang Province of China as the study area, this study aims to evaluate the best time window for the remote sensing inversion of soil texture in the cultivated land. 25 Sentinel-2 images were selected in the study area from 2019 to 2021. The band and spectral index of each image were input into the random forest model, in order to establish a remote sensing inversion model of soil texture. A comparison was made on the model accuracy of soil texture inversion from the images in different periods. The most suitable image was determined for the remote sensing inversion of soil texture. The spatial distribution of soil texture was obtained to evaluate the accuracy of soil texture inversion from the Sentinel-2 multi-spectral images on different dates. The results showed that: 1) From the visible light to short wave infrared 1 (1 565-1 655nm) spectral reflectance increased with the increase of wavelength, from the short wave infrared 1 to short wave infrared 2 (2 100-2 280nm) spectral reflectance decreased significantly. The spectral reflectance decreased with the increase of silt and clay content, whereas there was an increase with the increase of sand content. 2) There was the maximum accuracy of inversion model in the silt and sand on May 7, 2020 (the coefficient of determination (2) of silt was 0.785 in 25 Sentinel-2 images, and Root Mean Square Error (RMSE) was 6.697%; the2of sand was 0.776, and RMSE was 8.296 %). The maximum accuracy was also achieved in the clay inversion model on May 3, 2019 (2=0.776, RMSE=1.600%). 3) The appropriate time of satellite images was selected as an important impact on soil texture inversion. The best time window was from late April to mid-May in the study area. The time window and good-quality spectral data were obtained to develop a stable spectral model for the spatial distribution of soil texture. 4) Different inversion accuracy of soil texture were obtained using the 25 Sentinel-2 images from 2019 to 2021. This data was attributed to the soil water content and straw mulching. 5) The silt and clay particles were distributed more in the northeast, north, and south, and less in the middle and southwest of the study area. There was the opposite trend of the sand, especially the generally high sand content in the middle of the study area. Therefore, the high-precision remote sensing inversion model was achieved in the soil texture. The finding can provide the key technologies for regional soil texture mapping and farmland protection in the black soil areas.

soils; random forest; remote sensing; black soil; time window; Sentinel-2; influencing factor; texture

10.11975/j.issn.1002-6819.2022.18.013

S157

A

1002-6819(2022)-18-0122-08

劉瓊，羅沖，孟祥添，等. 典型黑土區(qū)耕作土壤質(zhì)地遙感時(shí)間窗口及影響因素分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2022，38(18)：122-129.doi：10.11975/j.issn.1002-6819.2022.18.013 http://www.tcsae.org

Liu Qiong, Luo Chong, Meng Xiangtian, et al. Time window and influencing factors analysis of tillage soil texture remote sensing in the typical black soil region[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2022, 38(18): 122-129. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2022.18.013 http://www.tcsae.org

2022-04-26

2022-08-10

國(guó)家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目（2021YFD1500100）；王寬誠(chéng)教育基金

劉瓊，研究方向?yàn)檗r(nóng)業(yè)遙感。Email：liuqiong0422@163.com

張新樂(lè)，博士，副教授，研究方向?yàn)樯鷳B(tài)遙感。Email：zhangxinle@gmail.com