农田是相对封闭的人工生态系统，作物生长的必要条件基本都可以通过农田管理进行控制，土壤养分是其中的关键要素之一，养分含量的高低能在很大程度上影响作物长势及作物最终的产量，合理控制土壤养分含量对作物增产、农民增收、环境保护及土壤肥力保持等都有重要的意义。变量施肥是精准农业的主要内容之一(蒙继华 等，2011)，变量施肥处方图是当前变量施肥的主要依据，而处方图的生成离不开土壤速效养分信息的准确获取。

当前土壤养分信息获取主要有传统地面调查、土壤光谱和植被指数3种方式，传统地面调查法通过实地采样化验的方式获取采样点土壤养分含量(Theng和Tate，1989)，该方法可以在点上保持较高的精度，但在进行面尺度外推时存在代表性差、成本高、时效性差等缺点(Lamsal 等，2009)。土壤的反射光谱是土壤表层颗粒水分、矿物质等组分的综合体现(Molan 等，2014；吴昀昭 等，2003)，对土壤表层养分含量具有一定的指示作用，可以利用地面高光谱数据对土壤中的有机质、全氮、全磷等含量进行反演(Ehsani 等，2001；刘焕军 等，2010；徐丽华和谢德体，2012)。随着高光谱遥感特别是机载高光谱技术的发展，高光谱影像在土壤养分反演中得到应用，使基于该数据的土壤养分快速大面积估算成为可能(Selige 等，2006；Hbirkou 等，2012)。土壤光谱法需要裸土期的高光谱数据建立反演模型，对模型运行时间和所需遥感数据获取时间要求高，且反演主要针对土壤表层，应用价值受限。植被指数法是通过连接土壤养分和作物长势参数从而利用遥感数据间接估算土壤养分(Meng 等，2015；潘瑜春 等，2007)，该方法的理论基础是养分缺失引起的作物长势特征的变化，而随着遥感技术的发展，利用遥感手段可以精确稳定地反演主要农田参数(蒙继华 等，2010)，结合一定的反演模型可以对耕层土壤速效养分含量进行动态模拟，该方法可以部分解决前两种方法中存在的问题，使反演结果更具实际意义，在近几年成为土壤养分反演研究的热点，但目前长势和养分之间的关系一般采用经验统计模型，模拟的精度和稳定性往往无法保证。总体来说，由于土壤速效养分在耕作层的立体分布、在空间上的异质性及时间上的不稳定性，传统方法和土壤光谱法的应用均受到限制，植被指数法摆脱了另外两种方法在应用时间、范围、成本、结果代表性等方面的问题，是实现土壤速效养分反演的有效方法，但是想要使结果更好的指导施肥，需要对模型进行优化升级，使其在理论基础、精度和稳定性方面得到提升。

作物模型因其生理生化基础和对作物生长的详尽数字模拟在作物长势和产量模拟中具有明显的优势，利用作物模型代替植被指数法中的经验统计关系，可以提高反演的精度和稳定性，实现植被指数法的完善和升级。作物模型可以模拟作物的生理生化过程、土壤水分和养分吸收及运移过程，并且可以克服异常气象条件和自然灾害对作物生长的影响(Gerakis和Ritchie，1998；Soler 等，2007；Ma 等，2013)。由于土壤养分是作物生长的必要条件，多数作物模型都将养分参数列为作物生长模拟的核心参数，其变化会引起模型模拟结果的相应变化(Boogaard 等，1998；Akponikpè 等，2010)，这也是利用作物模型进行土壤速效养分反演的前提。作物模型是在气象站点尺度上的设计开发的农学模型，进行区域外推是会遇到参数标定困难的问题(Dorigo 等，2007)，充分利用遥感数据在空间、时间和精度上的优势是解决这个问题的有效途径，其中数据同化技术发挥了重要的作用。遥感数据同化技术起步于20世纪80年代(Maas，1988)，目前同化算法主要分两种：一类是基于统计估计理论的，如最优插值法OI(Optimal Interpolation)、集合卡尔曼滤波EnKF(Ensemble Kalman Filter)、粒子滤波PF(Particle Filter)等；另一类是基于变分方法的，如三维变分(Three-dimensional variational analysis)、四维变分(Four-dimensional variational analysis)等(姜志伟，2012)，这些方法被应用到CERES-Wheat(Crop Environment Resource Synthesis)、WOFOST(World Food Studies)和STICS(Simulateur multidisciplinaire Pour Les Cultures Standard)模型中，在农情监测中发挥了作用(Dong 等，2013；Huang，2016；Chen和Cournède，2014)。利用数据同化方法使作物模型与遥感技术结合实现了大范围作物生长的动态模拟，也为作物模型与遥感反演的作物参数结合实现土壤养分有效反演创造了条件。

本研究针对当前主要养分反演方法中存在的问题，提出利用作物模型代替经验统计模型实现土壤速效养分有效反演的新方法，实现了耕层土壤速效养分的准确、稳定模拟，提高了养分反演结果的应用价值。本文将对算法理论基础、反演思路、算法流程及应用效果进行介绍。

本研究选择双山基地农场作为研究区，农场位于黑龙江北部，所在区域位于中国最大的春玉米种植区，适合以春玉米为对象开展研究。双山基地农场(48°48′N，125°31′E)属于寒温带大陆性季风气候，冬季较长且降雪量较大，夏季降水集中，全年平均降水量555 mm，有效积温(> 10 ℃)2250 ℃，主要种植春玉米、大豆和小麦，其中春玉米的种植面积超过50%，种植模式为一年一熟，播种时间一般为5月中下旬，成熟时间在9月底至10月初。农场地块大(平均地块面积为850亩)，地块作物种植单一，适合农田尺度的遥感监测与应用。农场的土壤类型单一，上层土壤为黑土，土层厚度1.2 m左右，养分含量高，下层为黄土，沙性土壤，厚度在10 m以上，肥力低，本研究主要面向黑土层进行养分反演研究。

本研究在试验区进行了多年连续观测，主要采用样点的方式采集农田作物和土壤要素，作物要素针对玉米整个生长季5月—10月进行连续观测，主要观测参数包括生物量、LAI、作物单产、作物生育期等，用于标定模型作物参数；土壤要素主要包括土壤特性参数和土壤养分参数的观测，土壤特性参数主要用于模型土壤参数标定，养分观测用于验证养分反演精度，养分参数采集时间为2014年5月14日—21日。作物参数和土壤特征参数的观测在之前的文章中有比较详细的介绍(Cheng 等，2016)，本研究重点介绍土壤速效养分观测。采样按照先布设样方后布设样点的顺序，首先在农场范围布设样方，样方的选取及布设原则是：(1)包括不同类型的地块，尽量在高中低产田均匀分布；(2)样方内部的样点分布根据地块形状和种植结构布设，垂直于垄向分布；(3)样方选取同时考虑了本研究使用的遥感数据分辨率，每个样方统一设计为30 m×30 m的方形区域，样方的布设间隔为150—240 m，即5—8个遥感像元。样点在样方内以斜对角线的形式等间隔选取3个采样点，每个采样点以相同的标准采集土壤样品，以3个点的样品化验均值作为样方采样值。基于以上原则，在玉米种植地块中选择12个实验地块，共71个样方，218个采样点，每个样点分别在垄上和垄间各取一个样，采用土钻采取地面以下20 cm的土壤样本，将两份土样充分混合、烘干、碾碎之后，在实验室进行化验，获得每个样品的有速效氮、有效磷、速效钾含量，最终得到71组有效数据(5组对照数据)，图1显示了实验区的具体位置和土壤养分采样点的分布情况。

图 1 研究区域及2014年采样点分布 Figure 1 Location of study area and the field observation sites in 2014

本研究使用时间序列HJ-1 CCD数据，数据采集自环境与灾害监测预报小卫星星座A、B星(HJ-1A/1B星)，空间分辨率为30 m，包含4个波段，波长分别为蓝光波段B1 (430—520 nm)，绿光波段 B2 (520—600 nm)，红光波段B3(630—690 nm)和近红外波段B4(760—900 nm)。所使用的时间序列HJ-1 CCD如表1所示。

表 1(Table 1)

表 1 双山基地时间序列HJ-1 CCD数据 Table 1 Time-series HJ-1 CCD data of Shuangshan Farm 传感器 采集时间 轨道号 HJ-1A CCD1 2014-05-03 451-52 HJ-1A CCD2 2014-06-06 448-56 HJ-1A CCD1 2014-06-23 453-53 HJ-1B CCD2 2014-09-06 448-56 HJ-1A CCD1 2014-09-18 454-53 HJ-1A CCD2 2014-09-24 452-56 HJ-1B CCD2 2014-10-04 452-56

表 1 双山基地时间序列HJ-1 CCD数据 Table 1 Time-series HJ-1 CCD data of Shuangshan Farm

遥感数据预处理主要包括辐射定标、大气校正、几何精校正，其中，辐射定标使用官网给出的定标系数，大气校正利用ENVI的FLAASH模块及影像相关参数进行校正，几何校正使用DEM及TM影像为底图进行几何精校正，误差控制在0.5个像元以内。同时考虑到7、8月份遥感数据缺失的情况，采用数据融合的方式补足缺失数据，融合基于16天合成的MODIS 250 m植被指数产品(MOD13Q1)与HJ-CCD数据计算得到的NDVI产品进行，采用改进型STARFM时空融合算法(Meng 等，2013)，分别得到7月12号、7月28号和8月13号的30 m分辨率NDVI数据，用于后续LAI的反演及数据同化。

研究选择WOFOST模型模拟实验区春玉米生长并设计相应的算法实现养分反演，该模型是荷兰瓦赫宁根大学和世界粮食研究中心共同开发研制，是欧盟作物长势监测系统(CGMS)的核心模型(van Ittersum 等，2003；Boogaard 等，2013)。模型主要模拟气候条件、土壤条件、管理条件对作物的影响，通过作物、气候、土壤、养分4个模块模拟同化作用、呼吸作用、蒸腾作用、干物质的分配等生理过程，对作物的生物量、产量、叶面积指数等进行模拟，可以在潜在、水分限制、养分限制3个水平上模拟作物生长与产量形成，其中养分限制产量在潜在产量或者水分限制产量基础上通过考虑土壤速效养分对作物生长的胁迫作用模拟得到，土壤养分模块是计算养分限制产量的核心，可以模拟氮磷钾的吸收、胁迫及其相互影响，得到3种养分共同限制条件下的作物产量，同时养分模块也是本研究中进行养分反演的关键模块，模块优化改造和算法补充主要针对该模块进行。

WOFOST模型在中国的应用已经超过10年(Wu 等，2003)，模型主要的参数已经通过地面试验和遥感数据同化等方式在不同气象、土壤和管理条件下进行了比较充分的标定(Ma 等，2013；Huang 等，2016)，然而很多因素会使作物模型参数发生变化，甚至同一作物不同品种之间的作物参数会有较大的差别，因此有必要对模型的核心参数进行标定。本研究中选择的标定参数包括：作物参数、土壤水分参数、土壤养分参数和气象参数，研究团队前期针对WOFOST模型参数标定进行系统研究，总结了一套完整的标定方法和流程(Cheng 等，2016)，本文沿用这些标定方法。

研究采用集合卡尔曼滤波(EnKF)作为遥感数据同化算法，该算法是由Evensen于1994年根据Epstein的随机动态预测理论提出的一种顺序数据同化算法(Evensen，1994)，它基于集合预报的思想，利用状态集合来表示概率密度函数，使用蒙特卡洛方法预报背景场的误差协方差，再计算分析场的集合(Evensen，2009)。通过引入集合，使得非线性的数据也能够进行同化计算，相比变分算法，又提高了计算效率。EnKF标准分析方程为

式中，集合卡尔曼增益系数：

式中，A_a为分析矩阵，A_f为预报矩阵，A_c为状态变量的协方差矩阵，D_t为观测变量矩阵，D_c为观测变量的协方差矩阵，HA为观测变量的误差协方差矩阵，H为观测算子。从上述方程可以看出，在计算的过程中，同化算法综合考虑了WOFOST模型模拟结果与遥感观测数据的误差。由于对于作物生长内在机理的理解并不是非常完善，因此随着模型的驱动，其受到参数标定不完善影响，误差积累逐渐增大，因此引入同化算法对模型模拟过程进行校正；而遥感数据的误差来源主要有两方面，其一是遥感数据反演LAI等作物生理生化参数所带来的误差，其二是在成像过程中受到云层、水汽等因素影响而带来的误差，针对有云层覆盖的影像，参照前人的研究(Lin 等，2008)，通过设置膨胀系数来避免在同化中后期出现排斥遥感数据的现象，为了不影响同化算法的效率，利用影像绿波段与近红外波段的亮度值提取被云层覆盖区域；此外，与成像质量较好的影像相比，受水汽影响严重的影像整体方差较小，因此引入调节因子增大受水汽影响的影像方差，以避免误差较大的外部数据大量进入模型。

时间序列叶面积指数有效反演是进行遥感数据同化的前提，也是实现养分反演中实际长势参数模拟的关键环节。本文首先根据地面观测数据与遥感反演的NDVI产品构建线性回归模型，模型构建过程中同时考虑到生育期的影响，对模型进行了分段，具体模型构建在论文中已有详细论述(Cheng 等，2016)，模型形式为

式中，DVS是WOFOST模型中生育期的数字化表述，DVS等于0时表示作物生长季的起点，1时代表春玉米处于抽穗初期，2时代表作物生长季的结束，NDVI达到峰值。

针对基于作物长势特征的植被指数法存在的问题，本研究提出利用作物模型和遥感数据同化来实现土壤速效养分的稳定准确反演。反演过程包括数据准备、可行性分析、数据同化算法设计和养分反演4个部分，具体可以分为地面观实验、遥感数据处理、反演理论可行性分析、数据通话算法构建、模型优化改造、养分反演算法设计和反演精度评价7个步骤，流程如图2所示。其中，地面观测实验基于模型参数标定和养分反演精度评价设计相关农田参量的采样，可以生成参数标定数据集和精度评价数据集，遥感数据处理包括数据下载、预处理和叶面积指数反演，生成数据同化数据集，可行性分析主要从模型结构、参数敏感性等角度分析利用作物模型反演土壤养分是否可行，数据同化算法构建主要是实现EnKF同化算法使时间序列遥感数据参与到模型模拟中，将模拟由点外推至区域，模型优化改造包括模型结构调整、算法优化和算法补充，养分反演和精度评价是在上面工作的基础上设计相应的算法实现土壤速效养分的有效反演并利用地面观测数据进行精度评价。反演理论基础、模型优化改造、反演算法设计和精度评价是本文研究的主要内容。

图 2 养分反演基本流程 Figure 2 The process of soil nutrients estimation

WOFOST模型是基于作物生长模拟和产量预测设计的，模型正向运行时通过输入作物参数、土壤参数和气象参数结合遥感数据同化算法以特定步长模拟作物长势，并输出相应长势参数，在参数标定准确的情况下，输出的应该是实际长势，实际上在田块或者像元尺度上标定所有参数非常困难，而其中的土壤养分参数在这个尺度上具有明显的空间异质性和时间不稳定性，是无法标定的主要参数，可以认为是作物模型模拟长势与实际长势存在差异的主要来源。同时，我们可以通过一定的遥感手段获取像元尺度的作物实际长势，如果能够利用一定技术手段将作物长势作为输入参数，土壤速效养分含量作为输出参数反向运行模型，便可以在遥感数据的参与下实现土壤速效养分的有效反演，具体的反演流程如图3所示。

图 3 土壤养分反演方法流程 Figure 3 The technical process of soil nutrients simulation

作物长势参数遥感估算结果和模型模拟值的差异和模型反向运行是养分反演是否可行的关键，模型的逆向运行可以通过引入相应算法并对模型调用顺序进行调整实现。长势差异的理论基础是研究区当前养分状况对作物生长具有胁迫作用，可以通过参数敏感性分析得到验证，两者都将在后面的不同小节中进行介绍。

本研究首先对模型进行重新编码实现养分模块的动态调用，同时通过引入作物养分吸收方程(宋海星和李生秀，2003)，实现养分吸收量的分生育期模拟，解决养分模块实时调用时产生的养分吸收问题。其次，考虑到施肥量在本研究中更具实际指导意义而且在原始模型中施肥量并没有参与作物养分限制产量的计算，本研究通过单位变换并结合模型中对DVS的模拟将土壤基础养分含量和施肥量结合作为统一参数参与模型模拟。最后，养分参数在不同部位中的运移和作物生长后期(LAI下降阶段)养分胁迫作物模拟也会影响养分模拟结果，其中养分参数的运移会改变养分元素在不同作物器官中的含量，特别在作物生长后期，养分会从营养器官运往存储器官，实时观测元素含量的变化需要付出比较多的人力和物力，本研究采用简单的系数法将作物作为一个封闭的整体来解决运移问题。部分作物生长参数在达到生长顶峰之后会在后面的生育阶段内出现负增长，养分限制作物在该阶段内的影响方式与前期不同，对养分的限制机理进行调整可以解决这个问题。

通过算法设计和模型优化，可以实现作物养分限制生长的实时模拟和模型反向运行，进而通过一定的算法实现土壤速效养分的模拟。

考虑到农田化肥使用量是农田养分管理过程中人为可控的主要参数，本研究选择施肥量分析其变化对作物长势的影响，分析针对整个农场2014年所有玉米地块进行，在整体施肥量在–10%—10%的范围内变化时分别模拟所有玉米地块的单产并统计平均单产的变化情况(统计结果如表2所示)。结果表明农场尺度上单独和整体改变氮肥、磷肥和钾肥施肥量都会引起产量模拟结果的变化，考虑到施肥量在改进之后的模型中是养分参数的一部分，可以看出在土壤基础养分含量发生变化时，产量模拟结果同样会发生相应的变化，因此参数敏感性分析的结果支持了养分胁迫的理论基础，通过整合作物模型和时间序列遥感数据实现土壤速效养分反演在算法、模型结构、理论基础上都是可行的。

表 2(Table 2)

表 2 农场尺度养分参数敏感性分析结果 Table 2 Soil nutrients parameter sensitivity analysis results at farm scale N P K Yield N P K Yield N P K Yield –10% –10% –10% –1.80% 0% –10% –10% –1.23% 10% –10% –10% –0.67% 0% –1.37% 0% –0.80% 0% –0.25% 10% –0.96% 10% –0.39% 10% 0.15% 0% –10% –0.97% 0% –10% –0.41% 0% –10% 0.12% 0% –0.55% 0% 0.00% 0% 0.53% 10% –0.15% 10% 0.40% 10% 0.92% 10% –10% –0.19% 10% –10% 0.36% 10% –10% 0.87% 0% 0.22% 0% 0.76% 0% 1.27% 10% 0.60% 10% 1.14% 10% 1.65%

表 2 农场尺度养分参数敏感性分析结果 Table 2 Soil nutrients parameter sensitivity analysis results at farm scale

养分模拟算法在模型参数标定、EnKF同化算法设计、模型改造及算法补充的基础上进行设计，首先将整个生长模拟过程分为3个阶段：起始生长模拟阶段(阶段1)、水分限制条件下的生长模拟及土壤速效养分反演阶段(阶段2)、养分限制条件下的生长模拟阶段(阶段3)，阶段1和阶段2的分界点(A点)是养分反演的起点，阶段2和阶段3的分界点(B点)是养分反演的终点，养分反演主要在两个点之间实现。阶段1中，模型首先由作物生长的起点开始模拟作物的播种出苗等生育过程，然后在EnKF同化算法的参与下模拟作物生长到A点，这期间的作物生长是水分胁迫下的生长，由于A点之前的作物生长阶段靠近土壤基肥和种肥的施用时间，同时作物生长早期的作物生长养分需求量较少，该阶段的作物生长受气象和土壤水分的影响更加明显，因此可以用该阶段的水分胁迫生长代表实际作物生长，在A点处的生物量和LAI等接近于水分—养分限值条件下的作物长势参数，这一阶段的主要作用是利用遥感数据同化方法在养分模拟阶段开始前实现作物实际生长的有效模拟，确保在A的模型各参数模拟值得准确性；阶段2中，模型在A点开始在由模型自行模拟作物生长，由于土壤水分模块的参与，该阶段是针对作物水分限制长势模拟，养分模块在B点参与计算，利用反演输入参数模拟作物的水分—养分限制长势，由于所输入的养分参数无法保证在像元尺度上的准确性，此时的水分—养分限制长势与遥感反演模型得到实际作物长势必然存在一定的差别，利用该差别在B点设计一定的养分算法可以反演相应的养分值；阶段3中，利用阶段2反演得到的各个像元的养分值代入到A点，模型在养分模块的参与下开始模拟作物实际的水分—养分限制长势至生长季结束，得到作物实际的产量及生物量参数。养分反演的3个阶段基本流程如图4所示。

图 4 养分反演基本流程 Figure 4 Process of proposed soil nutrients estimation method

在第2阶段中，设计了两种不同的养分反演算法：查找表和构建方程组，其中构建方程组法是通过构建方程组的反向求解获取土壤速效养分含量，LAI和土壤养分分别作为自变量和因变量参与方程组计算，该方法具有反演速度快的优点，可以在很大程度减少反演时间，但是不同养分限制条件下所需要构建的方程组是不一样的，并且部分方程组很难直接求解，虽然可以采用趋势线拟合的方式求解，但是该算法对遥感数据的数量和质量提出了更高的要求，同时考虑到研究区域在7月—8月以多雨天气为主，往往很难满足趋势线拟合对影像质量和数量的要求，一定程度上限制了该方法的应用。查找表法是在不改变养分模块模拟顺序的情况下遍历可能的养分输入参数，寻找与实际情况最接近的养分组合，比对的参数是由模型模拟的养分限值条件下的LAI值和利用前文中建立的LAI遥感反演模型计算得到的LAI值，在正常养分范围内两者差值的最小值所对应的养分含量就是所要反演的养分结果。该方法的优点是算法稳定，对遥感影像的数量和质量要求相对较低，可以在实际应用中进行多年多区域推广，也是在本研究中采用的养分反演核心算法。

WOFOST模型参数标定和数据同化算法构建是本研究中实现作物长势阶段性动态模拟的主要方式，其精度会对后续的养分反演产生直接影响，本研究通过文献资料调研、农场数据搜集、田间观测和遥感反演4种方式对模型的核心参数进行了标定，标定之后的参数如表3所示，参数标定之后的模型可以在农场尺度上准确模拟作物生长，遥感数据同化算法的引入实现了像元尺度作物生长实时模拟。

表 3(Table 3)

表 3 WOFOST模型主要作物及土壤参数标定值 Table 3 The value of major input crop and soil parameters for WOFOST 参数名称 参数含义 数值 单位 TDWI 初始干物质总重 120 kg/ha TSUM1 出苗到开花需要的积温 890 °C·d TSUM2 开花到成熟需要的积温 710 °C·d CVL 叶片的干物质转化效率 0.66 kg/kg CVO 存储器官的干物质转化效率 0.80 kg/kg CVR 根的干物质转化效率 0.74 kg/kg CVS 茎的干物质转化效率 0.65 kg/kg FRTB 干物质分配到根的比例 0—0.40 kg/kg FOTB01 地上干物质分配到存储器官的比例 (DVS=0.1) 0 kg/kg FOTB05 地上干物质分配到存储器官的比例(DVS=0.5) 0 kg/kg FOTB11 地上干物质分配到存储器官的比例(DVS=1.1) 0.23 kg/kg FOTB17 地上干物质分配到存储器官的比例(DVS=1.7) 0.76 kg/kg FLTB01 地上干物质分配到叶片的比例(DVS=0.1) 0.75 kg/kg FLTB05 地上干物质分配到叶片的比例(DVS=0.5) 0.70 kg/kg FLTB11 地上干物质分配到叶片的比例(DVS=1.1) 0.31 kg/kg FLTB17 地上干物质分配到叶片的比例(DVS=1.7) 0.20 kg/kg FSTB01 地上干物质分配到茎的比例(DVS=0.1) 0.25 kg/kg FSTB05 地上干物质分配到茎的比例(DVS=0.5) 0.30 kg/kg FSTB11 地上干物质分配到茎的比例(DVS=1.1) 0.46 kg/kg FSTB17 地上干物质分配到茎的比例(DVS=1.7) 0.04 kg/kg TMPFTB 最大CO2同化速率随温度的变化率 0.42 — RML 叶片的维持呼吸消耗系数 0.04 — RMO 存储器官的维持呼吸消耗系数 0.02 — RMR 根的维持呼吸消耗系数 0.01 mm RMS 茎的维持呼吸消耗系数 0.01 day NBASE 土壤基础速效氮含量 308 mg/kg PBASE 土壤基础有效磷含量 37 mg/kg KBASE 土壤基础速效钾含量 174 mg/kg NF 氮肥施肥量 251 kg/ha PF 磷肥施肥量 129 kg/ha KF 磷肥施肥量 141 kg/ha SMTAB-1 土壤体积含水量(pF= –1) 0.410 cm3/cm3 SMTAB01 土壤体积含水量(pF=1) 0.395 cm3/cm3 SMTAB02 土壤体积含水量(pF=2) 0.265 cm3/cm3 SMTAB06 土壤体积含水量(pF=6) 0.084 cm3/cm3 SMFCF 田间持水量 0.295 cm3/cm3 SMW 凋萎点土壤含水量 0.084 cm3/cm3 SM0 土壤饱和含水量 0.410 cm3/cm3 RDMCR 土壤允许的最大根系深度 1.8 cm CONTAB 土壤渗透系数(pF=0—4.2) –7.2—2.0 cm3/cm3 NMEANSO 存储器官平均氮元素的浓度 0.0154 — PMEANSO 存储器官平均磷元素的浓度 0.0047 — KMEANSO 存储器官平均钾元素的浓度 0.0040 — NMEANVE 营养器官平均氮元素的浓度 0.0086 — PMEANVE 营养器官平均磷元素的浓度 0.0023 — KMEANVE 营养器官平均钾元素的浓度 0.0128 —

表 3 WOFOST模型主要作物及土壤参数标定值 Table 3 The value of major input crop and soil parameters for WOFOST

研究在养分反演算法应用之前，分析了作物模型对关键作物参数的模拟精度，选择单产、LAI和关键生育期普遍出现时间为指标，分析结果如表4所示，结果表明在模型关键参数标定和EnKF数据同化算法构建的基础上，作物模型可以对整个生长季内的主要生长参数进行准确模拟。同时，在养分反演过程作物模型主要应用在阶段1和阶段2进行作物长势参数模拟，时间跨度要小于整个生长季，考虑到误差累计作用，在A点处的参数模拟精度通常要高于表4中的结果，可以满足阶段2中养分反演的精度要求。

表 4(Table 4)

表 4 WOFOST模型主要作物参数模拟精度分析 Table 4 Simulation accuracy of major crop growth parameters using WOFOST 参数名称 观测结果 模拟结果 R2 Yield 9718.35 9784.56 0.61 LAI(DVS=0.5) 1.69 1.81 0.79 LAI(DVS=1.1) 3.97 4.12 0.53 LAI(DVS=1.7) 2.95 3.04 0.64 出苗期 5—29 5—25 — 扬花期 7—25 7—27 — 成熟期 9—18 9—22 —

表 4 WOFOST模型主要作物参数模拟精度分析 Table 4 Simulation accuracy of major crop growth parameters using WOFOST

在双山基地农场对养分反演算法进行了应用和验证，首先在2014年5月份对双山基地开展了地面调查，得到了玉米地块的养分观测值，然后利用模型在对应时间进行养分估算，并利用测值对养分模拟结果进行精度分析和算法优选。

从养分算法设计的原理可以看出，改变B点的值会改变遥感数据的使用，同时引起养分反演针对的作物生育期的变化，而养分对作物生长的限制作用会随着作物的生育期发生变化的，例如，速效氮对叶绿素及叶片的生长作用比较明显，这种作用在作物生长的前期比较明显，后期会因为叶片饱和、相互遮挡等因素受到弱化，在遥感数据上体现的不如前期明显。因此改变B点的值会在一定程度上影响养分模拟结果的精度。本研究中在固定A点的基础上，改变B点的对应时间，分别模拟养分反演结果，利用实测数据分析不同时间的模拟精度，寻找不同养分对应的最佳生育期。

实验分别设定土壤基础含量与施肥量合并之后的速效氮输入参数在0—1000 mg/kg 的范围内以1为步长遍历，0—200 mg/kg 的范围内以0.5为步长，遍历有效磷的值，在0—800 mg/kg 的范围内以1为步长，遍历速效钾的值。同时，考虑到野外观测的时间、作物生长特性及遥感数据同化的需求，A点设定在DVS=0.4分别设定B点的对应时间为儒略日(DOY)=174—274范围内以10天为步长变化。遥感数据通过构建归一化植被指数(NDVI)时间序列曲线，提取对应时间的NDVI值并模拟得到实际状态下LAI值。最后以相对偏差(MRD)和相关系数(R²)为指标分析相应的精度。

速效氮反演的精度呈现先上升后下降再上升的趋势，在第一个峰值即DOY=193处，相关系数达到最大(R²=0.68，MRD=7.45%)，在DOY=203处，相对偏差达到最小(R²=0.59，MRD=6.73%)。考虑到指标的指示价值和数值上的差异，可以将DOY=193作为速效氮模拟的最佳时间节点。有效磷反演的精度呈相似的变化，在第一个峰值即DOY=203处，相关系数和相对偏差精度均达到最佳(R²=0.74，MRD=6.17%)，将DOY=203作为速效磷模拟的最佳时间节点。

速效钾的反演精度分析采用相同的方法，结果表明速效钾的精度较低，R²的最大值仅为0.14，没有达到极显著相关的水平，这与速效钾本身的不稳定、养分对作物生长影响类型和程度等有关，钾在作物生长过程中对茎的生长作用比较明显，具有增强作物抗倒伏能力等效果，但是作用的程度要明显低于速效氮和有效磷(李合生，2012)。当3种养分同时产生生长胁迫的时候，速效钾的养分胁迫效果很容易被另外两种养分掩盖，因此要直接刻画这种胁迫比较困难。本研究尝试在消除氮和磷胁迫作用的基础上开展速效钾的反演研究，首先对速效氮和有效磷进行反演，将反演结果带入到A点，在养分模块的参与下模拟氮磷胁迫作用下的作物生长，在B点使用相同的反演算法对速效钾做反演并进行精度分析，分析分别在速效氮单独作用、有效磷单独作用、两者共同作用下的速效钾反演结果，结果表明3种情况下有效钾的反演精度均在DOY=223处达到最大，并且消除速效氮单独作用的精度最高(R²=0.53，MRD=9.97%)。精度分析结果显示，利用该方法对速效氮、有效磷和速效钾(消除速效氮胁迫)模拟均达到了比较高的精度(图5)。

图 5 2014年5月份双山基地速效养分模拟精度分析结果 Figure 5 Analysis results of soil available nutrients simulation of Shuangshan Farm at May, 2014

同时，利用分析得到的最佳监测期对双山基地农场的土壤速效养分进行了监测(图6)，从结果中可以看出，3种养分在农场和田块内部均呈现了一定的空间异质性，该方法具备了一定的养分空间分布模拟能力，使反演结果对变量施肥具有更好的指导作用。

图 6 2014年5月份双山基地土壤速效养分监测结果 Figure 6 Soil available nutrients map of Shuangshan Farm at May, 2014

本文针对传统土壤速效养分反演中存在的问题，提出利用作物模型代替作物长势胁迫法中的植被指数经验统计模型，从而实现土壤速效养分稳定准确获取，该方法首先利用遥感数据同化技术将时间序列遥感数据同化到模型中，实现像元尺度作物参数的分阶段准确模拟，然后针对养分反演过程中的对模型结构和算法的要求对模型进行了优化和算法补充，在此基础上设计了土壤速效养分反演算法实现了土壤中速效氮、有效磷和速效钾的模拟，最后该方法在双山基地农场进行了应用并结合地面观测数据对反演效果进行了评价，结果表明，该方法针对3种速效养分的反演均到达了比较高的精度。

相比于其他土壤养分反演方法，新方法的优势主要体现在反演结果的应用价值、模拟的稳定性、模拟时间和经济成本3个方面。应用价值方面，该方法是针对耕层土壤速效养分模拟设计的，相比表层土壤总养分的估算，整个耕作层使结果更有代表性，速效养分模拟对地块地力评价和养分缺失分析更有价值，可以确保反演结果在指导精确变量施肥中的有效应用；在稳定性方面，新方法利用WOFOST作物模型代替了先前方法中普遍使用的统计模型，作物模型在理论基础、过程机理描述方面具有优势，可以有效克服统计模型在时间、空间外推时稳定性差的缺点；在成本方面，新方法采用基于养分生长胁迫的思路，能够克服积雪、作物覆盖等的影响，方法的应用时间可以扩展到整个生长季，实现了土壤养分的实时快速模拟，能够有效降低时间成本，同时，本研究中使用的是中等空间分辨率的多光谱HJ-CCD数据，数据覆盖范围广，时间分辨率高，获取免费，可以有效控制方法应用的经济成本。

同时，新方法当前还存在一定的局限性，具体来说，主要可以分为3个方面：(1)本研究的理论前题是养分胁迫对作物生长的影响，通过分析胁迫引起的LAI模拟值的变化实现土壤速效养分含量的反演，因此其他可能会对LAI值产生干扰的因素有必要进行分析并通过一定的方式不同胁迫进行分离，根据研究区的实际情况，将会对本研究的LAI模拟值产生影响的要素分为两种类型：类型一，在农场尺度上相对均一的要素，该部分要素包括大部分的气象条件和农田管理措施。其中，气象条件主要包括日照、温度等作物生长关键要素，管理措施主要包括播种、施肥、灌溉等田间管理措施，本研究选择的研究区双山基地的实验地块分布在16 km×18 km的范围之内，气象条件相对均一。同时，田间管理措施是在农场统一调控之下进行的，播种时间相差不大，施肥、灌溉等田间管理措施基本一致，由于这部分的参数在WOFOST中均有所体现，因此不会对养分反演的精度造成太大的影响；类型二，在田块之间甚至田块内部存在差异的要素，包括气象胁迫和土层厚度，本研究中气象胁迫主要是指在降水量、温度、日照、农田坡度、作物长势等条件的共同作用下引发的干旱或者洪涝灾害造成的作物生长抑制的情况，在田块甚至亚田块尺度上的胁迫程度往往会有所差异，在气象灾害发生时，气象胁迫作用甚至会超过养分胁迫成为影响作物生长的主导因子，会对养分反演造成比较大的干扰，如何消除气象胁迫的影响，使模型可以在发生气象胁迫时更好的模拟水分限制下作物生长是解决问题的关键，也是下一阶段研究的重要内容之一。由于水土流失和风蚀作用，研究区域内的黑土层厚度在空间上呈现了比较大的差异，土壤保持较好的地块黑土层厚度一般在1.5 m以上，部分侵蚀严重的地块黑土层厚度不足半米，在相同土壤速效养分含量情况下，两者可供给的养分量会有较大的差别，作物长势也会有相应的变化，这会导致反演精度降低，在后续研究中需要考虑让更多作物与土壤参数进入到养分模拟中；(2)本研究反演的土壤养分结果是A点的实时含量，可以较好的用于指导追肥，但是想要应用于指导基肥施肥，还需要研究相应的养分变换算法；(3)该方法使用的WOFOST模型的土壤模块是静态模块，不涉及养分在土壤中的运移、转换、流失和补给等方面的描述，这在一定程度上行会影响反演结果的精度和稳定性，消除这类影响当前比较有效方法是与其他土壤养分模型进行耦合，需要在后面的研究中加强这方面的研究。

新方法实现了土壤速效养分的快速稳定模拟，并保持了较高的精度，随着后续研究中针对当前问题的解决，新方法在指导变量施肥、提高农田生产效率、确保农田可持续发展方面可以发挥更大的支撑作用。

志 谢 此次研究的野外实验的数据获取得到了原沈阳军区双山农副业基地和中国科学院东北地理与农业生态研究所的支持，同时荷兰瓦赫宁根大学提供WOFOST模型的源代码和部分参数值对本研究有比较大的帮助，在此对以上单位表示衷心感谢。