基于地理探测器的山区作物生产潜力指数修正研究

杨永侠 张 函 郭雅萍 张丽红 桑 婧

(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)将修正的作物生产潜力指数与作物实际产量进行相关性分析，可以得出，修正后的作物生产潜力指数与实际产量具有极显著相关性，表明该修正指数比原唯一指数更加合理，采用这种修正方法能够提高耕地质量评价结果的精度。