何旭剛， 買買提·沙吾提， 盛艷芳， 李榮鵬

（1.新疆大學地理與遙感科學學院，新疆 烏魯木齊 830046；2.新疆綠洲生態(tài)重點實驗室，新疆 烏魯木齊 830046；3.智慧城市與環(huán)境建模自治區(qū)普通高校重點實驗室，新疆 烏魯木齊 830046）

棉花作為最重要的經濟作物，在我國新疆廣泛種植，2022 年新疆棉花種植面積為249.69×104hm2，總產達539.1×104t，棉花種植面積、單產、總產連續(xù)28 a 位居全國首位［1］，在新疆經濟發(fā)展中占有舉足輕重的地位。然而長期以來干旱少雨和缺水極大的限制了該地區(qū)棉花健康生長，尤其近年來，區(qū)域內生活、工業(yè)和生態(tài)用水占比逐漸增加，農業(yè)用水占比及可用水量在不斷減?。?］，導致新疆棉花生產面臨巨大挑戰(zhàn)。因此，解決水資源短缺與提高作物產量，以及提高農業(yè)用水效率，已成為新疆尤其是南疆農業(yè)發(fā)展中面臨的關鍵難題所在［3］。

水分生產率是指單位水資源量在一定的作物品種和耕作栽培條件下所獲得的產量，是作物灌溉管理和水資源優(yōu)化配置的主要依據和核心［4］。定量評價區(qū)域作物水分生產率是緩解干旱區(qū)水資源短缺和生態(tài)環(huán)境危機的重要前提。以往水分生產率研究基于田間實驗觀測和采用行政單元統(tǒng)計，不能準確表達空間分異和區(qū)域特征的需求。隨著3S 技術和Google Earth Engine（GEE）等平臺的發(fā)展，遙感在農業(yè)領域的應用日益廣泛［5］，遙感估算作物耗水量理論和方法日益完善，為區(qū)域作物識別、產量及耗水量估算等提供了可能，其中MODIS遙感數據具有高時間分辨率，周柳萍［6］、蔣磊等［7］和楊建瑩［8］利用MODIS數據對小麥、玉米和大麥等水分生產率分析了時空變化，并探究了灌溉系數與影響因素。然而，MODIS 數據的空間分辨率有限，對于破碎化耕地難以達到精細要求［9］。Landsat 影像有更高的空間分辨率，可以滿足長時間序列的耕地作物識別精度，Zhou 等［10］和Ghorbanpour 等［11］利用遙感數據估算了冬小麥和夏玉米等水分生產率，進一步對管理尺度評價，但以上研究中缺少對棉花作物水分生產率的分析，雖然Mauget［12］、Chen 等［13］、Thorp 等［14］和Himanshu 等［15］使用田間控制灌溉的方法對小型棉田水分生產率估算，但研究集中于小區(qū)域，缺乏較大空間尺度、長時間序列的研究，并無法滿足棉花水分生產率時空變化的分析視角。

鑒于此，本文利用GEE平臺的優(yōu)勢構建棉花識別模型、估產模型和蒸散發(fā)模型，對渭干河-庫車河綠洲（簡稱渭庫綠洲）棉花水分生產率進行探究，揭示棉花種植和產量時空分布、棉花生育期內耗水量的時空變化規(guī)律，在此基礎上對棉花水分生產率進行定量評價，可為干旱區(qū)棉花生產提供科學管理和技術參考。

1 研究區(qū)概況

渭庫綠洲（82°00′～83°60′E，41°00′～41°85′N）位于新疆南部，地處塔里木盆地和塔克拉瑪干沙漠北緣，是典型而完整的山前沖洪積扇平原，隸屬阿克蘇地區(qū)的庫車、沙雅和新和3 個縣級行政區(qū)（圖1）。光熱資源豐富，降水稀少，蒸發(fā)強烈，氣溫差異顯著，年均氣溫9.7～12.9 ℃，年降水量28.9～194.7 mm，年均蒸發(fā)量2000.7～2092.0 mm［16］。渭干河成為流域內生態(tài)、農業(yè)、人民生產生活用水的主要來源［17］。渭庫綠洲是干旱區(qū)綠洲農業(yè)的典型代表，以種植棉花、小麥作物為主。

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

2 數據與方法

2.1 研究數據與處理

為保證研究區(qū)該年有充足的影像估算棉花全育期的蒸散發(fā)，結合2014年國家政策對新疆棉花產量和價格的調整，選用前后年份的數據，即2009—2020年Landsat影像逐年6—7月數據用于棉花提取和產量估算；其次篩選2009年（含2008年5月8日、2008 年6 月9 日）和2020 年（含2019 年4 月5 日）棉花生長期（4—10 月）影像14 景，結合使用氣候再分析數據ERA5、數字高程數據SRTM等進行蒸散發(fā)遙感反演；同時使用蒸散發(fā)數據MOD16A2 和氣象數據驗證蒸散發(fā)的精度，Landsat數據基本信息見表1。

表1 研究數據Tab.1 List of the data used in this study

實地調研數據來自研究區(qū)歷年農作物調查數據，2009年至今每年一次，獲取30個樣方，樣方大小為500 m×500 m，樣方包括作物類型、面積、長勢、水分、溫度等；另歷年采集各類作物隨機樣本點若干。

氣象數據來源于國家氣象信息中心（http://data.cma.cn）發(fā)布的渭庫綠洲內3 個氣象站（沙雅站、庫車站、新和站）的逐日平均氣溫、最低、最高氣溫、降雨、風速、日照時間等；2009—2020 年《阿克蘇統(tǒng)計年鑒》作為棉花產量估算參考和驗證數據。

遙感數據均來自于GEE 云平臺（https://earthengine.google.com），利用GEE云平臺對影像數據去云、裁剪和拼接等預處理，將空間分辨率統(tǒng)一重采樣為30 m，然后使用該平臺提供的JavaScript語言和函數進行數據處理和模型構建。

2.2 研究方法

本研究根據作物水分生產率原理，通過遙感數據等數據，在GEE平臺建立研究區(qū)域棉花種植區(qū)分布識別模型、估產模型和SEBAL 遙感蒸散發(fā)模型，進而構建棉花水分生產率估算模型。同時結合局部空間自相關、標準差橢圓與重心遷移分布等地理分析方法，并準確和定量化評價研究區(qū)棉花水分生產率時空分布特征。

2.2.1 棉花水分生產率估算水分生產率反映了農業(yè)水量投入和產出的效率，是衡量農業(yè)科學和合理用水的綜合指標，即作物經濟產量與實際蒸散發(fā)的比值［18］，計算公式如下：

式中：WP 為作物水分生產率（kg·m-3）；Y為作物產量（kg·hm-2），通過建立遙感估產模型估算得到；ET為作物的蒸散發(fā)（mm·d-1），通過遙感定量反演獲得，并通過氣象站和蒸散發(fā)產品數據驗證。

2.2.2 棉花產量估算根據預處理后的Landsat系列影像，選取各地物樣本點的純凈像元，利用GEE 平臺構建主要作物生育期歸一化植被指數（NDVI）［19-22］時間序列曲線，根據棉花NDVI特征，獲取棉花出苗期和花鈴期的影像，通過時間特征，結合光譜特征、紋理特征、地形特征等構建分類特征。將提取的棉花種植面積與統(tǒng)計年鑒數據進行逐年對比，其平均相對誤差均小于5%。最終，選取分類精度較高的隨機森林方法進行分類，總體精度為83%，Kappa系數為0.78。因此，本文利用隨機森林分類方法提取研究區(qū)棉花種植面積。

進一步使用NDVI指數進行產量估算［23-26］，以往研究中棉花的產量與NDVI 具有較強的相關性，并在實際得到充分的驗證［23-25］，因此本研究使用柵格圖像推算，將縣級產量圖擴展到像元級空間分布，結合統(tǒng)計年鑒中的棉花經濟產量數據，得到研究區(qū)3 個縣級區(qū)劃的棉花產量像元級的均值，然后利用棉花的NDVI 指數擴展到單個像元的棉花產量［27］，計算公式如下：

式中：YP為單個像素級的棉花產量（kg·hm-2）；Yavg為縣級棉花平均產量，從統(tǒng)計年鑒獲得；NDVIP為花鈴期NDVI 值；NDVIavg為縣級花鈴期平均NDVI值。然后根據數值建立棉花產量線性關系。

2.2.3 棉花蒸散發(fā)遙感反演SEBAL（Surface energy balance algorithm for land）模型是以陸面能量平衡為理論基礎，通過遙感影像和氣象參數對各地表分量進行估算，得到區(qū)域的實際蒸散發(fā)。地表能量平衡方程為［28］：

式中：Rn為地表的凈輻射通量（W·m-2）；G為土壤熱通量（W·m-2）；H為下墊面到大氣的感熱通量（W·m-2）；λET 為潛熱通量（W·m-2），其中λ為潛熱蒸發(fā)系數，ET 為蒸散發(fā)（W·m-2），在計算水分生產率時，將ET單位轉換為mm·d-1。

因為衛(wèi)星過境時所觀測的是地面瞬時數據，因此把瞬時蒸散發(fā)擴展到日蒸散發(fā)，蒸散比在白天變化較小，可以假定白天蒸發(fā)比保持不變，根據瞬時蒸散比和日有效輻射能量計算出日蒸散發(fā)，利用日蒸散比插值法，得到衛(wèi)星過境日前后一段時間的逐日蒸發(fā)比［29］，然后根據各日的有效凈輻射量，計算出逐日蒸散發(fā)：

式中：ETd為第d天蒸散發(fā)值；ETi為衛(wèi)星過境時i的瞬時蒸散發(fā)值。

2009—2020年間，由于缺乏中間年份棉花生長期蒸散發(fā)估算的連續(xù)遙感影像，本文根據棉花生長物候歷［30］，選取起始和最終2 a的棉花生長期影像估算蒸散發(fā)。以研究區(qū)新和、沙雅、庫車3個氣象站點作為樣點，使用基于氣象站的P-M方法與SEBAL模型對3 個站點的平均蒸散發(fā)精度驗證，得到2009年蒸散發(fā)平均絕對誤差為0.56 mm·d-1、2020 年為0.27 mm·d-1，同時利用MODIS16A2 數據8 d 產品棉花生長期累積蒸散發(fā)精度檢驗，與SEBAL模型反演的2009年蒸散發(fā)平均絕對誤差為0.25 mm·d-1、2020年為0.27 mm·d-1，兩種數據都具有較高一致性，并取得較為滿意的估算結果。

3 結果與分析

3.1 棉花產量時空分析

棉花產量的時空變化是研究水分生產率的前提，將遙感估算的縣級總產量與統(tǒng)計年鑒數據逐年精度對比，其平均誤差小于8%。如圖2所示，2009—2020 年渭庫綠洲棉花產量呈增加趨勢，從2009 年的1610.10 kg·hm-2增長到2020年的1855.05 kg·hm-2，產量增長率為13.20%。從各縣來看，新和縣增長6.09%，沙雅縣增長12.72%，庫車市增長20.74%，新和縣的產量最高，其次為沙雅縣和庫車市。渭庫綠洲產量在2010年最低（1520.04 kg·hm-2），2019年達到最高（1984.80 kg·hm-2），2014 年出現低谷（1645.05 kg·hm-2）。由于2014 年棉花價格調整，開始大面積種植，水土和管理措施不成熟等原因使得棉花產量較低，而在2015 年之后趨于穩(wěn)定上升，主要由于棉花種植面積增加導致產量逐年上升［31］，加之政策補貼、水肥精細管理和機械化程度持續(xù)提高，對棉花產量起了推動作用。

圖2 2009—2020年棉花產量統(tǒng)計Fig.2 Statistics on cotton yield from 2009 to 2020

采用500 m×500 m 格網，對2009—2020 年逐年棉花產量按縣級行政進行重心遷移規(guī)律分析，新和縣產量重心在玉奇喀特鄉(xiāng)和渭干鄉(xiāng)之間移動（圖3a）。沙雅縣棉花產量重心穩(wěn)居于古勒巴格鄉(xiāng)（圖3b），呈現出先西后東移動趨勢。庫車市產量重心主要位于墩闊坦鎮(zhèn)、阿克吾斯塘鄉(xiāng)、齊滿鎮(zhèn)（圖3c）。渭庫綠洲產量重心在沙雅縣紅旗鎮(zhèn)和努爾巴格鄉(xiāng)之間波動遷移（圖3d），12 a 間自西向東移動2485 m，表明渭庫綠洲棉花產量隨著種植面積增加向東西方向擴展。總體上，棉花產量各時間點的標準差橢圓長軸各自沿綠洲下游方向逐漸移動，庫車市和沙雅縣橢圓面積大于新和縣，新和縣總體趨勢自東北向西南方向擴展，沙雅縣自西向東趨勢移動，庫車市則自西南向東北方向延伸（圖3e）。

3.2 棉花生長期蒸散發(fā)量時空分析

棉花生長期蒸散發(fā)的準確估算是研究水分生產率的重要環(huán)節(jié)，渭庫綠洲棉花生長期蒸散發(fā)2009年為237.86～843.95 mm（均值為686.80 mm，日均蒸散發(fā)為3.71 mm），蒸散發(fā)最大值為花鈴期，其次為吐絮期、蕾期、苗期和出苗（圖4k）。蒸散發(fā)高值主要分布在渭庫綠洲內部河流周圍及棉田小規(guī)模種植區(qū)域（圖4a～e）。2020 年棉花生長期蒸散發(fā)為493.13～996.3 mm（均值為738.66 mm，日均蒸散發(fā)為3.99 mm），其中蒸散發(fā)最大值為花鈴期，其次為吐絮期、苗期、蕾期、出苗（圖4l）。蒸散發(fā)高值主要在綠洲內部及塔里木河北岸邊緣（圖4f～j），由于地勢低洼土壤含水量較高又毗鄰沙漠邊緣，蒸散發(fā)偏高，低值主要分布在渭庫綠洲東北和西南邊緣，由于棉花集中連片種植，使用滴管技術降低水分蒸發(fā)。整體上，12 a 間棉花生長期蒸散發(fā)2020 年高于2009 年，其增長率為7.02%，尤其在花鈴期和吐絮期蒸散發(fā)均值明顯高于2009 年，但苗期蒸散發(fā)均值低于2009 年，然而取值范圍變大，表明該階段渭庫綠洲氣候差異性較為明顯，其余生長階段蒸散發(fā)均值較為穩(wěn)定（圖4k～l）。隨著全球氣溫上升，表明該區(qū)域蒸散發(fā)也呈上升趨勢，棉花蒸騰作用加劇。

圖4 棉花生長期蒸散發(fā)Fig.4 Evapotranspiration during cotton growing

3.3 棉花水分生產率時空分析

3.3.1 棉花水分生產率時空變化分析基于像元統(tǒng)計水分生產率頻率統(tǒng)計，可以精準量化逐像元棉花水分生產率情況，由圖5可知，2009年棉花水分生產率變化范圍為0～0.40 kg·m-3（均值為0.21 kg·m-3），2020 年棉花水分生產率變化范圍為0～0.45 kg·m-3（均值為0.25 kg·m-3），12 a 間增長率為16%。2009年棉花水分生產率在0.18～0.25 kg·m-3之間，有明顯的一個峰值期，在0.21 kg·m-3對應的像元數最多，達到14632個。2020年棉花水分生產率在0.17～0.33 kg·m-3之間，有明顯的一個峰值期，在0.28 kg·m-3對應的像元數最多，達到26478 個，2009 年峰值比較集中，2020 年向橫軸右側移動，高值區(qū)間增大且較為分散。12 a 間水分生產率較高的區(qū)域增多，棉花水分生產率整體提高，布局更加合理。將水分生產率累積頻率為95%的像元值統(tǒng)計作為臨界值，2009年對應的值為0.37 kg·m-3，2020年為0.41 kg·m-3，大于該臨界值的為水分生產率較高區(qū)域，該區(qū)域在2009年主要集中分布在渭干鄉(xiāng)，2020 年分布于玉奇喀特鄉(xiāng)，2個鄉(xiāng)均位于新和縣且兩鄉(xiāng)毗鄰，土壤以殘余內陸鹽土為主，灌渠設施密集，其中棉花呈現連片化和規(guī)?；N植。

圖5 2009年與2020年水分生產率頻率分布Fig.5 Frequency distribution of water productivity in 2009 and 2020

為準確分析水分生產率在空間的聚集程度和演變趨勢，以縣級行政和綠洲為單位，建立棉花種植區(qū)標準差橢圓和水分生產率重心點。2009 年和2020 年新和縣水分生產率重心點分布位于渭干鄉(xiāng)和玉奇喀特鄉(xiāng)（圖6a），標準差橢圓的長軸持續(xù)增長且短軸持續(xù)縮短（圖6e），自2009年（82.25°E，41.27°N）移動至2020年（82.21°E，41.25°N），自東北向西南方向平均移動速度為888 m·a-1，由于棉花種植面積的擴張和產量的增加，水分生產率重心趨勢逐漸向綠洲外緣延伸。2009 年和2020 年沙雅縣水分生產率重心點均位于古勒巴格鄉(xiāng)（圖6b），自2009年（82.53°E，41.11°N）移動至2020年（82.52°E，41.12°N），重心點自東南向西北平均移動速度為320 m·a-1，12 a間變化趨勢較小，由于靠近塔里木河岸和沙雅縣南部沙漠邊緣棉花生長期耗水量較大，使得水分生產率重心向綠洲內部移動。2009 年和2020 年庫車市水分生產率重心均位于墩闊坦鎮(zhèn)（圖6c），生產率標準差橢圓的長軸縮短而短軸增長，其水分生產率重心點從2009年（83.05°E，41.28°N）移動至2020年（83.12°E，41.29°N），自西向東平均移動速度為1071 m·a-1。12 a間水分生產率重心移動速度依次為：庫車市＞新和縣＞沙雅縣，表明庫車市棉花水分生產率以最快速度自西向綠洲東部邊緣移動，棉花產量和蒸散發(fā)的變化直接影響水分生產率重心的遷移，同時說明庫車市棉花種植面積和產量重心向綠洲邊緣延伸面積最大。

圖6 棉花水分生產率空間分布及重心遷移Fig.6 Spatial distribution and gravity movement of cotton water productivity

總體上，渭庫綠洲水分生產率重心由2009 年（82.51°E，41.22°N）移動至2020年（82.52°E，41.21°N），自東北向西南移動距離為1832 m（圖6d），年均移動速度為166.55 m·a-1，表明綠洲棉花水分生產率呈現擴張趨勢且東西方向大于南北方向。渭庫綠洲水分生產率標準差橢圓的長軸持續(xù)增長同時短軸縮短（圖6f），長軸短軸的共同作用下離心率逐漸增長，同時標準差橢圓的面積由2009 年的3487.72 km2逐漸增加到2020年的4122.78 km2，縣級尺度（圖6e）和流域尺度（圖6f）標準差橢圓面積均增加，表明12 a間渭庫綠洲棉花水分生產率空間分布方向趨勢增強，空間格局趨向集聚化。

3.3.2 棉花水分生產率局部空間自相關分析以500 m×500 m 格網數據為基礎，對棉花水分生產率進行LISA 聚集性和Moran'sI指數分析，其中2009 年和2020 年棉花水分生產率的Moran'sI指數分別為0.349 和0.496，表明水分生產率具有顯著的空間自相關。LISA 聚類圖是對Moran'sI指數中通過了顯著性檢驗的區(qū)域單元的地理表達，渭庫綠洲正相關聚集（高-高聚集或低-低聚集）較為明顯，聚集程度呈現上升趨勢，異常區(qū)（高-低聚集或低-高聚集）呈零散分布，且波動較?。▓D7）。2009 年高值聚集區(qū)主要位于新和縣西側綠洲邊緣和沙雅縣古勒巴鄉(xiāng)，其余高值聚集區(qū)呈零散分布，該區(qū)域棉花種植較為密集，且產量較高，同時耗水量低，使得水分生產率高而聚集（圖7a）。低值聚集區(qū)零散條帶狀分布于綠洲中部且在沙雅縣境內最多，說明棉花產量較低且耗水量大，導致水分生產率低。異常區(qū)零散均勻分布于整個綠洲。2020年正相關程度逐漸上升且向綠洲邊緣擴展，空間分布基本上與2009年一致，但低值在綠洲邊緣和塔里木河北岸增加，表明該地水分生產率較低，由于新開墾的土地土壤貧瘠，沒有形成大規(guī)?；N植，產量較低（圖7b）。以上分析得到，高值區(qū)主要分布于沙雅縣新墾農場和新和縣桑塔木農場，所以在一定程度上說明，農場規(guī)?；N植增產和科學化灌溉，促使棉花產量和水分生產率提高。

圖7 水分生產率LISA聚集圖Fig.7 LISA aggregation diagram of water productivity

4 討論

本文充分利用GEE云平臺強大的計算能力，建立研究區(qū)遙感模型，分析棉花水分生產率的時空分布，為渭庫綠洲提高農業(yè)用水效率和緩解水資源供需矛盾提供新思路。同時，為今后新疆小麥、玉米等主要作物水分生產率遙感估算提供基礎。

在棉花產量遙感估算中，本文使用NDVI 指數法和統(tǒng)計數據空間化估算棉花產量，該方法比較通用且簡單，得到的結果與統(tǒng)計產量數據一致，能夠滿足像元級產量推算且效果較好。以往較多研究［32-33］使用增強植被指數（EVI）建立估產模型，今后應與該指數結合對棉花面積提取和估產，進一步提高估產精度。從渭庫綠洲棉花生長期蒸散發(fā)來看，本文使用的SEBAL 模型遙感估算棉花生長期蒸散發(fā)較好的結果，在精度分析中，除了使用氣象站點數據外，還使用了MOD16A2數據進一步驗證，為估算結果提供有利的數據支持和可靠性。針對干旱區(qū)作物蒸散發(fā)遙感估算的模型較多，如周柳萍［6］使用SSEB 模型、Zhou 等［10］使用SEBS 模型、于兵［34］使用HTEM-ABL 模型，今后研究中應結合多種模型對作物生長期的蒸散發(fā)估算。充分定量化利用空間自相關、標準差橢圓與重心轉移等空間分析方法進行棉花水分生產率時空分布。本研究遙感數據估算了棉花水分生產率起始和最終2 a，今后可以補充并加強連續(xù)性數據，如使用高空間分辨率Landsat數據和高時間分辨率MODIS影像時空融合的方法［10］提高水分生產率時間序列的分析。此外，以往研究中缺乏新疆棉花水分生產率的研究，本研究拓展并利用GEE平臺優(yōu)勢和遙感數據，大面積對棉花水分生產率進行縣域尺度和流域尺度時空分析與研究。影響水分生產率的主要因素為產量和蒸散發(fā)，12 a 間產量的增長速度超過了蒸散發(fā)的上升速度，使得水分生產率提高。

5 結論

（1）2008—2020年棉花種植面積整體呈增加趨勢，棉花種植面積逐年向綠洲邊緣擴展。同時，棉花產量也呈增加趨勢，從2009年的1610.10 kg·hm-2增長到2020 年的1855.05 kg·hm-2，產量增長率為13.20%，平均每年增長1.1%，新和縣的棉花產量最高，庫車市產量最低。產量重心整體自西向東移動2485 m。

（2）渭庫綠洲棉花生長期2009 年蒸散發(fā)均值為686.80 mm，日均值為3.71 mm，2020年為738.66 mm，日均值為3.99 mm，整體呈上升趨勢，增長率為7.02%，2009 年和2020 年蒸散發(fā)最大值主要分布在綠洲內部與塔里木河北部邊緣。

（3）2009 年水分生產率均值為0.21 kg·m-3，2020 年為0.25 kg·m-3，12 a 間水分生產率均值增長率為16%。在空間上，渭庫綠洲水分生產率重心由2009年（82.51°E，41.22°N）移動至2020年（82.52°E，41.21°N），在紅旗鎮(zhèn)自東北向西南移動了1832 m，年均移動速度為152.67 m·a-1，棉花水分生產率呈現擴張趨勢且東西方向大于南北方向，棉花水分生產率空間分布方向趨勢增強，空間格局趨向集聚化。