基于地理探測器的山區作物生產潛力指數修正研究

楊永俠 張 函 郭雅萍 張麗紅 桑 婧

(1.中國農業大學土地科學與技術學院， 北京 100083；2.自然資源部農用地質量與監控重點實驗室， 北京 100035)

0 引言

耕地質量等級調查評價成果是基本農田調整劃定、耕地差別化管理的首要依據[1-2]，是我國落實耕地數量、質量、生態三位一體保護的重要參考依據，是實現耕地數量-質量占補平衡，落實耕地質量管理的重要依據。根據文獻[3]中規定的方法，作物生產潛力指數是我國農用地分等工作的基礎參數之一[4]，其準確性在一定程度上決定了我國耕地質量等級調查評價成果的科學性。國內學者關于生產潛力指數的研究多集中在利用模型對不同尺度地區進行生產潛力指數估算或分析[5-8]，以及根據氣候變化對生產潛力指數造成的影響進行研究或預測[9-12]等方面。在耕地質量等級調查評價工作中，針對生產潛力指數劃分最小尺度(即縣級行政單位)下的地形或氣候對作物生產潛力指數的影響及修正研究尚不多見。

鑒于耕地質量等級調查評價工作對生產潛力指數科學性和精細化的需求，本文在原耕地質量評價體系的基礎上，提出一種基于地理探測器的山區作物生產潛力指數修正方法，以期提高耕地質量等級評價的準確性。

1 研究方法

作物生產潛力指數修正方法主要步驟為：分析生產潛力指數的測算過程，結合已有研究初步篩選影響因素；通過回歸分析確定顯著影響因素；不斷調整分層方式，比較地理探測器計算的不同分層方式下顯著影響因素對生產潛力指數的解釋力，得到每一顯著影響因素的最大解釋力并以該種分層方式作為這一顯著影響因素的最優分層方式；計算各顯著影響因素對生產潛力指數的解釋力歸一化值作為該顯著因素的權重；最后依據層內生產潛力指數和影響因素權重計算修正的生產潛力指數。作物生產潛力指數修正流程如圖1所示。

圖1 作物生產潛力指數修正流程Fig.1 Flow chart of correction for crop potential productivity

1.1 自然因素對作物生產潛力的影響

我國農用地分等體系中使用的作物生產潛力指數計算采用在作物光合生產潛力基礎上對溫度、含水率以及光照強度等因子逐級訂正的模型。數據基礎是國家氣象局在全國范圍內設立的氣象站點的觀測數據。決定作物生產潛力指數的各項“作物生產潛力測算參數”既有大尺度的輻射回歸系數，即華南地區、華中(西南)地區、華北(東北)地區、西北地區全國四大區所對應的輻射回歸系數；也有小尺度的各區縣氣象站的氣象數據如年均降雨量、溫度、水汽壓等。經過各種尺度參數綜合計算之后，所得生產潛力指數只代表各區縣氣象站點所在地形區的狀況，反映的是氣象站點所在地的農業氣候資源狀況[13]。

國內有研究揭示了作物生產潛力指數與某些自然因素的聯系[14-17]。本文選取海拔、坡度、坡向和年均降雨量與光溫(氣候)生產潛力指數進行相關性分析，篩選相關性顯著的因素對生產潛力指數進行修正。

通過回歸分析篩選影響因素，選擇線性回歸、對數回歸、二次回歸3種方式，分析中R2>0.1且通過顯著性水平p<0.05的假設檢驗的因素確定為顯著影響因素。

1.2 地理探測器應用

地理探測器是研究地理環境因素對疾病分布影響的方法[18-20]，其原理是基于因變量應與對其產生決定性影響的自變量在空間分布上具有顯著一致性的假設。其中，因變量可以是數值量也可以是類型量，自變量只能為類型量，自變量劃分出的不同類型即不同的層[21]。近幾年來，地理探測器被廣泛應用于研究土地利用變化[22-23]、景觀[24-25]、生態[26-27]以及農村經濟[28-29]等問題的機理。地理探測器從因子探測、風險探測、生態探測和交互探測4個方面對因變量和自變量之間的關系進行探測。本文主要應用地理探測器的因子探測功能。

因子探測：探測每個影響因子對生產潛力指數的解釋力，表達式為

(1)

式中h——根據作物生產潛力指數或其影響因子的差異劃分的層序號

nh——層h內研究對象個數

L——劃分層的個數

n——全區研究對象個數

σ2——全區的作物生產潛力指數的總方差

根據因子探測的原理，本研究對地理探測器進行逆向應用，即確定影響較大的因素后，不斷調整因素的分層方式，記錄地理探測器測算的解釋力，選擇解釋力最高的分層方式作為最優分層結果。

目前地理探測器常見分層方式有：①專家經驗法。組織相關領域的專家學者對因素進行分析，提出分層意見。②等間隔分層法。將影響因素屬性進行等間隔劃分類別，每個間隔內的對象劃分為同一層。③數據聚類法。使用K-means等聚類算法對影響因素進行分層。本文選用第2種方法，客觀且簡單易行。

1.3 生產潛力指數修正方法

在確定影響因素及其最優分層結果后，計算層間平均生產潛力指數作為該層生產潛力指數。當多因素同時影響，則利用歸一化后的各因素對生產潛力指數的解釋力作為該因素權重，解釋力q越大則權重越大，通過該因素修正的生產潛力指數可信度越高。生產潛力指數計算公式為

(2)

其中

(3)

式中εt——未知點的生產潛力指數

εi——未知點受第i種因素影響的層間平均生產潛力指數

qi——第i種因素在最優分層情況下對生產潛力指數的解釋力

ωi——第i種因素在最優分層情況下對生產潛力指數產生影響力的權重

2 研究區概況與數據來源

2.1 研究區概況

平安區坐落在青海省東部，海東市中心腹地，地處101.82°～102.17°E、36.25°～36.57°N之間，東臨海東市樂都區，南與化隆縣相鄰，西連省會西寧市和湟中縣，北隔湟水河與互助縣相望，南北相距33.6 km，東西相距33.69 km，轄區總面積769 km2。全區轄8個鄉鎮，共111個行政村、7個社區[30]。縣域內大部分地區為山地，地形復雜，溝壑縱橫，湟水河自西向東流經全境，地勢南高北低，由西南向東北傾斜，大部分地區海拔在2 066～3 000 m之間，海拔最高處為4 166.7 m，最低處為2 066 m，相對高差2 100.7 m，耕地分布在海拔2 066～3 300 m之間。

平安區屬于大陸性氣候，年平均日照時長為2 711 h，日均7.4 h。年平均氣溫7.3℃，最熱月平均氣溫18.7℃(7月)；最冷月平均氣溫-6.2℃(1月)。年內降雨分布不均，常年各地區降雨量240～550 mm，年均降雨量338.5 mm，年均蒸發量1 845.9 mm，為年均降雨量的5.5倍。

平安區內耕地面積為11 987.96 hm2，水澆地主要分布在平安區東北部海拔較低的溝谷地帶，僅占全區耕地總面積的26%；南部、西部大片耕地均為無法保證灌溉的旱地，約占全區耕地總面積的74%。最新一期耕地質量年度評價成果顯示，平安區內耕地的國家自然等為11～13等，耕地質量整體處于較低水平。

2.2 數據來源

本文基礎數據主要包括研究區域內各區縣氣象站點的地形、氣象數據。其中，海拔、坡度和坡向等地形數據是基于GIS平臺對GDEMDEM 30 m分辨率數字高程模型進行坡度分析和坡向分析后提取得到；70個氣象站點的降水數據來源于中國氣象局與中國氣象科學數據共享服務網提供的國家級氣象站點地面氣候資料日值數據集；平安區各耕地縣級分等單元的降水數據是利用平安區氣象站點觀測數據和測土配方施肥項目的降水監測數據獲取的；所有區縣的光溫(氣候)生產潛力指數摘錄自GB/T 28407—2012 農用地質量分等規程附錄中生產潛力指數速查表。

3 平安區作物生產潛力修正

3.1 數據預處理

平安區標準耕作制度為油菜、春小麥，熟制為一年一熟，基準作物和指定作物均為春小麥，作物生產潛力指數為春小麥的生產潛力指數。在生產潛力指數修正階段選取參考樣本時，同時考慮了標準耕作制度二級區和地形的影響，最終選取海拔1 500～3 500 m、基準作物為春小麥的區縣，在平安區周邊篩選出甘肅省、青海省、內蒙古自治區和寧夏回族自治區共70個區縣的春小麥生產潛力指數。

獲取的70個區縣數據中，指定作物生產潛力指數分布如圖2所示，其中光溫生產潛力指數最高值為甘肅省蘭州市市轄區1 744，最低為青海省門源回族自治縣780；氣候生產潛力指數最高值為青海省互助土族自治縣1 072，最低為甘肅省民勤縣150。所選取區縣的氣象站點坡度在1.43°～33°之間，海拔在1 225.7～3 663.6 m之間，年降雨量在58～700 mm之間。

圖2 鄰近區縣作物生產潛力指數分布Fig.2 Distribution of crop production potential index of surrounding counties

3.2 顯著影響因素確定

依次對研究區海拔、坡度、坡向、年均降雨量與作物生產潛力指數進行回歸分析，篩選相關性大的因素作為作物生產潛力指數修正因素。篩選標準為3種回歸模型的決定系數R2均大于0.1且p小于0.05。

篩選結果如表1所示，坡向與光溫生產潛力指數和氣候生產潛力指數回歸R2小于0.1，坡度與光溫生產潛力指數回歸p大于0.05，與氣候生產潛力指數回歸R2小于0.1。排除坡度和坡向為作物生產潛力指數的顯著影響因素；海拔和年均降雨量與光溫生產潛力指數和氣候生產潛力指數的回歸結果符合篩選標準，選擇海拔和年均降雨量為作物生產潛力指數修正的主要影響因素。

表1 影響因素-生產潛力指數回歸結果Tab.1 Regression result of influencing factor and crop production potential index

在農用地分等體系中，作物生產潛力指數的計算涉及作物生長期內光照條件、溫度、水汽壓、作物需水量等因素。在研究區內，海拔的變化對溫度的影響顯著；在旱地，降雨量直接影響作物的生長；但是在作物生產潛力指數計算模型中，光照條件的量化方式為利用作物生長季內各月平均大氣上界日輻射量(緯度每隔1°使用同一值)和當地氣象站點記錄的月平均日照百分率通過經驗公式計算當地作物生長期內各月的總輻射量，該計算方法得出的輻射量數據尺度較大且無法考慮地形對于光照的影響，因此回歸分析結果顯示坡度和坡向與作物生產潛力指數不存在顯著相關性。

3.3 最優分層方式

經過回歸分析可以確定，顯著影響光溫和氣候生產潛力指數的因素均為海拔和年均降雨量。利用地理探測器探測海拔、年均降雨量在不同分層方式下對指定作物光溫和氣候生產潛力指數的解釋力，選擇解釋力最大時的分層方式作為最優分層方式。為了提高本方法的普適性，便于分級，結合耕地質量調查評價分段劃等的等別確定方式，以固定步長調整分層基數和分層間隔。

通過測試可得海拔不同分層方式對生產潛力指數的解釋力影響，可以確定海拔對光溫生產潛力指數解釋力最高為0.837(圖3)，對氣候生產潛力指數解釋力最高為0.453(圖4)，海拔因素對光溫生產潛力指數和氣候生產潛力指數的最優分層方式均以1 100 m為基數，200 m為間隔。

圖3 海拔-光溫生產潛力指數分層探測結果Fig.3 Detection result of classes of altitude-light and temperature production potential index

圖4 海拔-氣候生產潛力指數分層探測結果Fig.4 Detection result of classes of altitude-climate production potential index

通過測試可得年均降雨量的不同分層方式對生產潛力指數的解釋力影響，可以確定年均降雨量對光溫生產潛力指數解釋力最高為0.526(圖5)，對氣候生產潛力指數解釋力最高為0.676(圖6)，年均降雨量因素對光溫生產潛力指數的最優分層方式為以0 mm為基數，100 mm為間隔，對氣候生產潛力指數的最優分層方式為以25 mm為基數，50 mm為間隔。

3.4 層間生產潛力指數確定

利用地理探測器計算得解釋力，確定影響因素最優分層方式后，計算該影響因素在最優分層方式下每一層的生產潛力指數。光溫生產潛力指數計算結果見表5、6，氣候生產潛力指數計算結果見表7、8。

圖5 年均降雨量-光溫生產潛力指數分層探測結果Fig.5 Detection result of classes of annual average precipitation-light and temperature production potential index

圖6 年均降雨量-氣候生產潛力指數分層探測結果Fig.6 Detection result of classes of annual average precipitation-climate production potential index

表5 海拔-光溫生產潛力指數Tab.5 Altitude-light and temperature production potential index

表6 年均降雨量-光溫生產潛力指數Tab.6 Annual average precipitation-light and temperature production potential index

表7 海拔-氣候生產潛力指數Tab.7 Altitude-climate production potential index

表8 年均降雨量-氣候生產潛力指數Tab.8 Annual average precipitation-climate production potential index

3.5 作物生產潛力指數修正

根據影響因素的最優分層方式及在該分層方式下由地理探測器得的解釋力，歸一化計算得海拔因素對光溫生產潛力指數和氣候生產潛力指數的影響權重分別為0.61和0.40，年均降雨量對光溫生產潛力指數和氣候生產潛力指數的影響權重分別為0.39和0.60。

由海拔和年均降雨量的權重和層間生產潛力指數，計算得到研究區所有縣級分等單元的光溫生產潛力指數分布區間為[997,1 217]，氣候生產潛力指數分布區間為[732,955]，分布結果如圖7所示。

圖7 修正后平安區生產潛力指數分布Fig.7 Distribution maps of crop production potential index of Ping’an County after amending

由圖7可以看出，平安區作物光溫生產潛力指數整體由東北向西南逐漸降低，光溫生產潛力指數較高的耕地主要分布在海拔較低的東北部低洼河谷地帶；氣候生產潛力指數與年均降雨量成正相關分布，隨著年均降雨量的增大，氣候生產潛力指數升高。經過修正后的生產潛力指數相對于年度更新評價中使用的生產潛力指數(光溫生產潛力指數1 364，氣候生產潛力指數663)有較大變化。

根據修正后生產潛力指數計算的國家自然等別與平安區耕地質量年度更新評價成果中的國家自然等別進行差值計算，結果如圖8所示。修正后平安區水澆地的國家自然等別整體降低0.4～1.2等，其中降低0.8～1.2等的耕地面積僅占全縣耕地總面積的4.0%，降低0.4～0.8等的耕地面積占22.1%。旱地的國家自然等別提升0～1.2等，等別提升小于0.4等的耕地面積占比最高，為52.8%，提升0.4～1.2等的耕地面積占全縣耕地總面積的21.1%。

圖8 修正前后平安區耕地國家自然等別變化值分布Fig.8 Distribution map of difference of Ping’an County’s national agriculture land physical quality classification before and after amending

3.6 修正模型合理性檢驗

作物的實際產量可以在一定程度上反映作物生產潛力指數修正結果的合理性，借助于平安區測土配方施肥項目布設的產量監測點數據，監測點分布情況如圖9所示。

圖9 產量監測點分布Fig.9 Distribution map of yield monitoring point

利用Pearson模型對實際產量與修正的作物生產潛力指數進行相關性分析，驗證修正結果的合理性。分析結果表明，83個產量監測點的作物實際產量與生產潛力指數相關系數為0.91，呈極顯著相關，修正結果與原作物生產潛力指數相比更具客觀性。

4 結論

(1)基于青海省海東市平安區地形和氣候特征，以平安區及其周邊共70個區縣的氣象站點作為研究樣本，利用地理探測器的解釋力分析功能，建立平安區生產潛力指數修正模型，修正了平安區全縣域內所有縣級分等單元的春小麥生產潛力指數。研究結果表明，影響平安區光溫生產潛力指數和氣候潛力指數的主要因素均為海拔和年均降雨量。

(2)在耕地質量等別調查評價中，平安區縣域內水澆地春小麥修正后的光溫生產潛力指數較原光溫生產潛力指數有所下降，使用修正后光溫生產潛力指數計算得到的國家自然等別指數較最新一期耕地質量等別調查評價成果降低0.4～1.2等；縣域內旱地春小麥修正后氣候生產潛力指數均高于原氣候生產潛力指數，使用修正后的氣候生產潛力指數計算得到的國家自然等別指數較最新一期耕地質量等別調查評價成果提升0～1.2等。

(3)將修正的作物生產潛力指數與作物實際產量進行相關性分析，可以得出，修正后的作物生產潛力指數與實際產量具有極顯著相關性，表明該修正指數比原唯一指數更加合理，采用這種修正方法能夠提高耕地質量評價結果的精度。