近些年来，随着国内外对资源环境、气候变化、生态和灾害等对国家或地区粮食生产系统影响等一系列问题定量研究的逐步深入，国家农业科学空间基础数据对上述事关国家粮食安全和应对气候变化等问题研究和科学决策具有更加重要的意义(Qiu 等，2003；Ye 等，2013)。其中，大范围、长时间和高精度农作物分布数据准确获取可为气候变化对国家粮食安全、粮食区域平衡、作物种植制度、生产结构和区域布局等方面影响研究提供重要基础数据，也具有重要现实意义和科学意义(Barnwal和Kotani，2013)。

遥感信息由于具有覆盖面积大、探测周期短、资料丰富和费用低等特点，已经在大范围作物空间分布信息获取和作物空间分布制图中发挥了越来越重要作用。冬小麦是中国主要的粮食作物之一，也是重要的商品粮和战略性储备粮。因此，开展准确获取冬小麦分布信息的作物识别和作物分布制图方法研究不仅可为气候变化、资源环境和粮食安全等领域研究提供重要基础科学数据，而且也对农业种植结构调整、农业政策科学制定等具有重要的现实意义。目前，常用的农作物识别遥感提取方法大体分为3类，第1类为基于地物空间光谱信息的差异进行农作物遥感识别的方法，第2类为基于时序遥感数据提取作物物候特征发展出的农作物遥感识别方法，第3类则是基于多源数据融合的农作物遥感识别方法(许文波和田亦陈, 2005)。

基于光谱特征的农作物遥感识别方法较早是以目视解译为主，之后出现了利用计算机基于分类图像统计进行分类的监督分类和非监督分类法，随着多学科交叉领域的发展，以句法结构分类法为主的各种集成分类方法(如神经网络法、模糊数学法、决策树法和基于混合像元分解的方法等)逐渐得到应用(唐华俊 等，2010)。但由于受卫星影像分辨率限制，导致分类过程中“同物异谱”和“同谱异物”现象难以避免，所以单纯基于地物光谱特性进行复杂种植结构作物分类较难获得理想效果(王娜 等，2017)。

对基于时序遥感数据提取作物物候特征的作物遥感识别方法而言，不同作物随着时相的变化具有明显的变化特征，基于时间序列遥感数据提取作物物候时相变化规律可以实现目标作物的识别，该方法通过分析时序数据中作物关键物候期特征进行作物类型识别，可避免利用单一时相影像数据进行作物空间分布提取时因“异物同谱”等原因导致的错分、漏分等现象。其中，基于时序影像单一特征量(如NDVI、EVI等)和作物物候特征的阈值法，不仅利用了作物随时间的生长变化规律，具有操作简单、运算速度快的特点，并且在作物种植类型较为单一的地区取得了较高的识别精度，因此，已经成为农作物特别是冬小麦识别中最常用的方法之一(蔡学良和崔远来，2009；张健康 等，2012；雷光斌 等，2014；刘吉凯 等，2015；Salahddine 等，2017)。近些年来，国内外学者已经开展较多基于阈值法和光学遥感数据的多种作物(如冬小麦、玉米、水稻、大豆等)分类研究，分类总体精度一般高于80%(Foerster 等，2012；黄青 等，2013；许青云 等，2014；杨闫君 等，2015；梁守真 等，2015)。然而，传统阈值法大多根据农作物关键时间点物候历特征及其“光谱—时序”曲线特征确定区分目标作物与其他地物的最佳时间点(胡琼 等，2015；郝震 等，2017)，最佳区分点的植被指数阈值大小主要通过训练样本光谱统计分析方法确定(蔡学良和崔远来，2009；刘佳 等，2015；张淼 等，2015)，或采用试错法，人工不断调试植被指数阈值的方式，根据经验对各关键物候期阈值不断调整尝试并多次分类，选择分类结果精度较高的一组阈值(张霞 等，2008；黄青 等，2013)。尽管上述方法能够取得较高精度结果，但传统阈值检测方法阈值设置缺乏灵巧性与智能性，不仅需耗费大量的工作量，主观性强，阈值检测结果存在较大偏差，从而阻碍了阈值方法作物识别精度进一步提高，阈值优化设置问题是阈值法作物识别中重要关键科学问题。近年国内外学者对阈值法阈值合理设置进行了一系列研究(Foerster 等, 2012；张健康 等，2012；郝鹏宇 等，2012；王利民 等，2016)，如利用地面样点时序植被指数参考曲线相似性阈值确定法、约束条件下的植被指数特征值最小阈值法和基于阈值分割的阈值筛选法等，但将全局优化思想引入阈值优化设置问题中的相关研究未见报道，且阈值设置准确性和自动化程度仍然有待进一步提高。

基于多源数据农作物遥感识别方法一类是遥感信息与非遥感信息融合的作物识别方法，另一类是多种传感器获取的具有不同时间、空间分辨率和不同波段信息的多源遥感数据融合的作物识别方法。其中，由于农作物面积统计数据在气候变化、国家粮食安全等领域是不可替代的重要基础数据，国内外学者已将统计数据作为非遥感信息与遥感信息进行融合，从而实现了大范围作物空间分布信息提取，这为全球变化与粮食安全研究提供了可靠的作物空间分布基础数据(Barnwal和Kotani, 2013；宋茜 等，2015)。如部分学者基于交叉信息熵原理的作物空间分配模型(SPAM)将遥感信息与农业统计数据进行融合，已在全球尺度及区域尺度上获得了较高精度的作物空间分布结果(Liu 等，2008；刘珍环 等；2013)。但是，目前一般遥感信息与非遥感信息融合的作物空间分布制图技术大多将遥感获取的土地利用、农业灌溉和耕地适宜性等作为辅助信息，而对研究对象作物自身生长机理及其变化规律考虑较少，没有充分地对作物自身遥感信息(特别是时序遥感信息)加以直接应用，从而阻碍了作物空间分布制图精度的进一步提高。近年来，一些学者研究表明，将NDVI等遥感特征参量引入遥感信息与非遥感信息(如统计数据)融合的作物空间分布制图中，降低了分类规则对训练样本的依赖，不仅方法简单，而且容易操作，可有效提高制图精度和效率(Pervez和Brown，2010；Brown和Pervez，2014)。

综合以上分析，本研究针对阈值法农作物识别过程中阈值设置存在灵巧性差和自动化程度低等弱点，提出引入全局优化算法实现作物识别模型阈值参数选取的方法，通过选择作物识别总量约束以及利用全局优化算法优化阈值检测方法两个方面进行改进，开展利用全局优化算法改进基于阈值检测的农作物分布制图方法研究，以期实现作物面积提取模型阈值参数快速准确的自动率定，从而进一步提高阈值法农作物识别精度。其中，本文以中国粮食主产区黄淮海地区河北省衡水市景县为典型实验区，以冬小麦为研究对象，以国产高分一号(GF-1)数据为主要遥感数据源，在阈值模型和全局优化算法支持下，开展以下两方面研究：首先，选择合适的总量约束措施，确定全局优化的标准参考；其次，确定合适的全局优化算法对阈值模型中待优化阈值参数进行设置，以提高参数计算效率和制图精度。通过研究，实现阈值法阈值参数自动优化选取，提高阈值选取客观性和准确性，从而进一步提高阈值法冬小麦空间分布制图的工作效率和分类精度，为获取大范围、长时间和高精度农作物空间分布基础农业科学数据提供一定思路借鉴。

本研究实验区域(115°54′E—116°28′E，37°28′N—37°52′N)位于中国北方粮食生产基地黄淮海平原内河北省衡水市景县，该县地势较为平坦，海拔在14.1—25 m。县域面积共1183 km²(图1)，其中，耕地面积约830 km²，以冬小麦和夏玉米为主要农作物，二者为一年两熟轮作制，种植总面积占总耕地面积的80%以上。实验区所在区域属温带半湿润季风气候，大于或等于0 ℃年积温4200 ℃—5500 ℃，年累计辐射量5.0×10⁶—5.2×10⁶ KJ. m^–2，无霜期170—220 d，年均降水量500—900 mm。实验区冬小麦集中种植时间为上年10月上中旬，11月下旬—12月上旬达到冬前分蘖峰值，12月中旬开始停止生长，第二年下年3月上中旬进入返青期，拔节期为3月下旬至4月上中旬，抽穗期为4月下旬到5月上旬，乳熟期为5月中下旬，成熟期为6月上旬。

图 1 实验区地理位置图 Figure 1 Location of experimentation area

本研究主要利用国产GF-1高分辨率遥感卫星数据进行冬小麦分布提取。GF-1卫星搭载4台16 m分辨率的多光谱相机(WFV1-WFV4)，光谱范围为0.45—0.89 μm，包含蓝(0.45—0.52 μm)、绿(0.52—0.59 μm)、红(0.63—0.69 μm)和近红外(0.77—0.89 μm)4个波段，多光谱相机组合幅宽达800公里，重访周期为4d。本文数据获取时间为2014年10月至2015年6月，影像来源于中国资源卫星应用中心(http://www.cresda.com/CN/ [2017-11-10])推送至农业部遥感应用中心的GF-1数据。选取影像时，以影像云量较低质量较好，并且影像时间间隔尽量一致为原则。最终，本研究共获得覆盖2015年冬小麦整个生长期的15景较高质量影像，所获取GF-1遥感数据详细信息见列表1。

表 1(Table 1)

表 1 实验区高分一号影像信息列表 Table 1 GF-1/WFV image list in experimentation area 时相序号 传感器类型 轨道号 获取时间 01 WFV2 L1A0000393184 2014-10-16 02 WFV4 L1A0000427927 2014-11-02 03 WFV3 L1A0000454185 2014-11-14 04 WFV2 L1A0000580482 2014-12-04 05 WFV1 L1A0000538637 2014-12-24 06 WFV2 L1A0000580482 2015-01-10 07 WFV3 L1A0000639856 2015-02-08 08 WFV4 L1A0000673553 2015-03-01 09 WFV2 L1A0000725498 2015-03-29 10 WFV1 L1A0000750015 2015-04-14 11 WFV1 L1A0000765287 2015-04-22 12 WFV3 L1A0000781757 2015-05-01 13 WFV4 L1A0000812059 2015-05-18 14 WFV1 L1A0000826550 2015-05-25 15 WFV3 L1A0000850828 2015-06-07

表 1 实验区高分一号影像信息列表 Table 1 GF-1/WFV image list in experimentation area

研究中，GF-1遥感影像预处理过程主要包括正射校正、辐射定标、大气校正、几何精校正等，数据预处理过程主要在ENVI5.3中完成。其中，借助高分一号影像自带的RPC(Rational Polynomial Coefficient)文件对原始影像进行无控制点的有理多项式模型区域网平差正射校正。辐射定标利用中国资源卫星应用中心(http://www.cresda.com/CN/ [2017-11-10])提供的绝对定标系数对经过正射校正后的GF-1影像进行辐射定标获得了影像的辐亮度文件。其中辐射定标公式如下：

式中， $L\left(\lambda \right)$ 为转换后辐亮度(W/(m²∙sr∙μm))， ${\rm{Gain}}$ 为定标斜率(W/(m²∙sr∙μm))， ${\rm{DN}}$ 为观测灰度值， ${\rm{Bias}}$ 为定标截距。其中， ${\rm{Gain}}$ 和 ${\rm{Bias}}$ 均来自中国资源卫星应用中心。

采用ENVI5.3软件中自带的FLAASH模型对影像进行大气校正以获得真实地表反射率数据。最后以15 m分辨率的Landsat 8/OLI影像(https://glovis.usgs.gov/ [2017-11-10])及地面实际采集样点作为参考进行几何精校正，误差控制在一个像元内。

归一化差值植被指数NDVI(Normalized Difference Vegetation Index)是反映植被生长状态及植被覆盖度的最佳指示因子之一，目前已被广泛应用于农作物识别和面积估算等方面(Zhang 等，2017；Saeed 等，2017)。其中，基于时序植被指数数据提取冬小麦面积与类型的研究方法已经比较成熟(Goodwin和Collett，2014；Massey 等，2017)，因此，本研究也将采用时序高分一号NDVI数据开展冬小麦识别研究。其中，通过大气校正、几何精校正等预处理后，GF-1 NDVI具体计算公式如下：

式中， ${\rho _{{\rm_{NIR}}}}$ 为GF-1卫星WFV传感器近红外波段的反射率； ${\rho _{\rm_{R}}}$ 为GF-1卫星WFV传感器红光波段的反射率。NDVI值域范围为[–1, 1]。本研究中部分GF-1遥感卫星NDVI结果如图2所示。

图 2 实验区景县GF-1 NDVI计算结果 Figure 2 GF-1 NDVI data of Jingxian county in 2015

本研究于2015年4月在开展了实验区冬小麦分布提取训练样本、验证样本的信息采集。根据实地调查，由于该地区在冬小麦种植时期主要影响冬小麦识别的地物为林地，因此，实验区共分为冬小麦、林地、其他(裸地、建筑用地等)3个类别。其中，林地及其他在最终分类中划分为非冬小麦类别。实地调查中，地物类型样点分布主要遵循均匀分布原则，并利用差分GPS进行样本点定位，最终，共获得2653个地面实测各种地物类型样点。其中，冬小麦1860个，林地385个，其他类型408个。研究中，随机选择600个样本点作为训练样本，用于提取地物时序特征，剩下的2053个样本点作为验证样本，用于后续提取结果的精度验证，具体样本点数量和分布，如表2和图3所示。

表 2(Table 2)

表 2 实验区地面样本点采集数量 Table 2 Number of samples in Jingxian county 样本类别 验证点采集数量/个 冬小麦 非冬小麦 总计 林地 其他 训练样本 400 100 100 600 验证样本 1460 285 308 2053

表 2 实验区地面样本点采集数量 Table 2 Number of samples in Jingxian county

图 3 实验区野外采集样本点分布 Figure 3 Distribution of sample points in experimentation area

本研究涉及的辅助数据主要包括研究区行政边界矢量图(县、市、省级)、县级农业统计数据以及冬小麦物候信息等。其中，景县2015年冬小麦种植面积县级统计数据主要来源于衡水市统计年鉴(衡水市统计局和国家统计局衡水调查队，2016)，该数据主要用于实验区冬小麦分布提取方法研究中的行政单元内冬小麦种植面积总量控制参数。

本研究在冬小麦时序NDVI参考曲线基础上，利用冬小麦关键物候期特征建立基于阈值法的冬小麦分布提取阈值模型。其中，模型中的阈值参数作为待优化参数。在此基础上，以区域农作物种植面积统计数据为作物种植面积总量全局优化的标准参考，用于和冬小麦面积分布提取模型获取的冬小麦面积建立优化目标函数。然后，通过SCE-UA全局优化算法寻找最小目标函数值。当遥感提取的区域冬小麦面积与区域冬小麦面积统计数据间差异最小时，便可得到区域内冬小麦分类模型最优参数，并在最优参数下输出冬小麦空间分布信息，从而获得高精度的冬小麦空间分布结果。最后，完成对冬小麦空间分布结果的精度评价和对比分析，具体如图4所示。

图 4 研究方法流程图 Figure 4 Flowchart of research

其中，利用SCE-UA算法对阈值模型中待优化参数进行全局优化的主要过程是：首先利用训练样本，确定各待优化阈值参数的上下界；然后，SCE-UA全局优化算法按照其自身搜索法则在上下界中寻找一组阈值参数作为初始值，代入冬小麦分布提取阈值模型中，获得研究区冬小麦空间分布结果，并进行像元统计获得区域冬小麦种植面积。将该利用遥感提取的区域冬小麦种植面积数据代入目标函数中，计算目标函数值(即利用遥感提取的通过像元统计获得的区域冬小麦面积与区域冬小麦面积统计数据间差值)，判断目标函数值是否满足收敛条件：若否，则重新初始化阈值，返回至确定阈值初始值的步骤，重新赋予一组阈值参数初值，并进行后续步骤；若满足收敛条件，则输出该组阈值，即为所求的最优化阈值。

在野外采集的实验区样本数据库支持下，利用训练样本在时序遥感影像中分别提取冬小麦及其他各类地物类型的NDVI时间序列数据。最终，绘制不同地物NDVI均值形成的NDVI时序特征曲线，如图5所示。由于研究中采用高分辨率遥感影像，很难获得时间间隔一致的影像，因此，图5中显示的特征曲线为非等时间间隔特征曲线。由图5中可以看出，结合实验区内物候信息，冬小麦各个生育期在时序影像上的表征为10月上中旬为冬小麦播种期，此时影像中主要表现为裸地光谱信息，整体NDVI值偏低；10月下旬至11月上旬为冬小麦出苗期，NDVI值逐渐升高，直到冬前分蘖期，NDVI值达到比较高的水平，此时出现了冬小麦整个生长季的第一个波峰；之后冬小麦进入越冬期，叶绿素含量降低，光合作用减弱，植株基本停止生长，NDVI值逐渐减小，并达到波谷；3月上旬开始，冬小麦进入返青—拔节期，NDVI值迅速增长，是NDVI曲线斜率最大的时期，并于4月下旬至5月上旬进入冬小麦抽穗期，NDVI达到最大值，此时出现冬小麦整个生长季的第2个波峰；之后进入灌浆乳熟期，NDVI值不断降低，直至冬小麦成熟。

图 5 不同地物GF-1/WFV NDVI时序曲线 Figure 5 GF-1/WFV NDVI time series curve of different features

从图5分析可知，与非冬小麦(如林地和其他地物等)对比，整个生长季内只有冬小麦出现了明显的两次波峰，并且两次波峰的NDVI值均明显大于非冬小麦。在出苗期–分蘖期和抽穗期–成熟期，冬小麦与非冬小麦的时序植被指数特征差异较大。在出苗期到分蘖期，冬小麦NDVI值迅速增大，而林地NDVI值却迅速减小，其他类型NDVI值波动较小。在抽穗期到成熟期，冬小麦NDVI值迅速减小，而林地地物类型NDVI值反而逐渐升高，其他类型NDVI变化幅度较小。因此，利用冬小麦生长期内的关键物候期作为区分冬小麦的识别特征，从而建立冬小麦识别模型。

根据以上对冬小麦及其他各类地物类型的物候特征分析，本研究利用冬小麦出苗期、冬前分蘖期、拔节期、抽穗期和成熟期等关键物候期以及冬小麦双峰波形特征NDVI时序变化规律，确定不同物候期下冬小麦像元需要满足的判别函数和特征参数，从而建立基于时序NDVI和物候特征的冬小麦空间分布提取模型。

(1) 根据冬小麦冬前分蘖期和抽穗期平均NDVI值均达到峰值的特征，建立判别函数：

式中，T₁表示冬小麦第4时相NDVI应满足的阈值，T₂表示冬小麦第11时相NDVI应满足的阈值。

(2) 根据只有冬小麦出苗—分蘖期和拔节期平均NDVI值呈快速增长的特征，建立判别函数：

式中，T₃表示冬小麦第4时相与第2时相NDVI差值应满足的阈值，T₄表示冬小麦第11时相与第9时相NDVI差值应满足的阈值。

(3) 根据只有冬小麦进入越冬期和成熟期平均NDVI值呈快速下降的特征，建立判别函数：

式中，T₅表示冬小麦第4时相与第6时相NDVI差值应满足的阈值，T₆表示冬小麦第14时相与第15时相NDVI差值应满足的阈值；

式(3)到式(5)中，t_i(i=1, 2, $ \cdots $ , c)代表冬小麦整个生育期开始对应的c个时相； ${\rm{NDV}}{{\rm{I}}_{{t_{{i}}}}}$ 表示某一像元在t_i时相对应的NDVI平均值；T_i(i=1, 2, $ \cdots $ , d)表示冬小麦各判别函数对应的阈值，为待优化的参数。本研究中时相数c=15，待优化阈值参数个数d=6。以上3个判别函数组成了衡水市景县冬小麦空间分布提取阈值模型，如图6所示。

图 6 阈值模型冬小麦识别流程 Figure 6 The flowchart of winter wheat identification by threshold method

针对传统阈值法中阈值设定存在的弱点，提出利用全局优化算法对阈值参数进行优化选取以提高参数计算效率和制图精度的方法。目前国内外常用传统的全局优化算法主要有遗传算法(葛继科 等，2008)、模拟退火算法(林慧君和彭宏，2006)和蚁群算法等(段海滨 等，2007)。然而，遗传算法全局优化效果较好，但该算法局部搜索能力较弱，运算量大，运行时间长，且对于参数敏感从而优化结果易受影响。模拟退火算法虽然运行效率高，操作简单，局部搜索能力强，但该算法是一种随机搜索算法，全局搜索能力差，容易陷入局部最优解因而难以获得全局最优解(王雪梅和王义和，1997；卢宇婷 等，2015)。(Duan 等，1994, 1993)针对水文模型率定提出的SCE-UA(Shuffled Complex Evolution-University of Arizona)复合形混合演化算法，综合了遗传算法和控制随机搜索算法的优点，是一种对于非线性复杂模型最优化问题寻找全局最优解最为成功的方法之一(董洁平 等，2012)。并且具有可以同时优化模型中的若干个参数并迅速搜索到全局最优解的特点，收敛速度快，稳定性高，且对初始值不敏感(邢会敏 等，2017)，是参数优选问题中非常有效且应用广泛的方法(Liu 等，2005)。

目前SCE-UA算法大多应用于水文模型参数优化中，由于其综合了多种算法的优点，全局搜索能力强、效率高，且针对非线性问题优化效果良好的特点，在遥感同化作物生长模型的参数优化(闫岩 等，2006；黄健熙 等，2015；刘峻明 等，2018)、陆面过程模式参数同化优化(张生雷 等，2015)等领域中也有了广泛的应用。针对本研究中模型没有具体函数表达式，且具有参数多、非线性强、待优化参数有区间约束的特征，本研究首次尝试引入SCE-UA全局优化算法作为冬小麦分布提取阈值模型参数自动优化的全局优化算法。本研究利用SCE-UA全局优化算法的目的是以行政单元冬小麦种植面积统计数据作为外部优化数据比较对象，当提取县级行政单元冬小麦面积与县级行政单元冬小麦面积统计数据间差异最小时，得到县域内冬小麦分类模型最优参数。最后，在最优参数下输出冬小麦分类算法提取的冬小麦分布信息。SCE-UA 算法进行目标函数最小化问题求解的具体步骤(Duan 等，1993；董洁平 等，2012；张景 等，2014)如下：

(1) 初始化。假定待优化问题是n维问题，选取参与进化的复合型个数p(p≥1)和每个复合型所包含的顶点数目 $m\left({m \geqslant n + 1} \right)$ ，计算样本点数目 $s = p \times m$ 。

(2) 产生样本点。在可行域内随机产生 $s$ 个样本点 ${x_1}, \cdots,\;{x_s},\;$ 分别计算每一点 ${x_i}$ 的函数值 ${f_i} = f\left({{x_i}} \right),\;$ $ i = 1, \cdots, s$ 。

(3) 样本点排序。把 $s$ 个样本点 $\left({{x_i}, {f_i}} \right)$ 按函数值升序排列，排序后仍记为 $\left({{x_i}, {f_i}} \right), i = 1, \cdots, s$ ，其中 ${f_1} \leqslant {f_2} \leqslant \cdots \leqslant {f_s}$ ，记D $ = \left\{ {\left({{x_i}, {f_i}} \right), i = 1, \cdots, s} \right\}$ 。

(4) 划分复合型群体。将D划分为p个复合型 ${{{A}}_{{1}}}{{, }}\cdots\!{{, }}{{{A}}_{\rm{p}}}$ ，每个复合型含有m个点，其中 ${A^k} = \left\{ {\left({x_j^k, f_j^k} \right)}\right.\big|$ $ \left.{\left. {x_j^k = {x_{k + m\left({k - 1} \right)}}, f_j^k = {f_{j + m\left({k - 1} \right)}}, j = 1, \cdots\!, m} \right.} \right\}, k = 1, 2, \cdots, p$ 。

(5) 复合型进化。按CCE(复合型进化算法)分别进化各个复合型。

(6) 复合型混合。把进化后的每个复合型的所有顶点组合成新的点集，再次按函数值f_i升序排列，排序后不妨仍记为D，对D按目标函数的升序进行排列。

(7) 收敛性判断。如果满足收敛条件则停止，否则返回步骤(4)。

模型参数的优化过程就是寻找一组阈值参数，利用该组阈值参数使利用模型提取的值与实测值之间按照目标函数设定的规则达到最吻合的状态。本文中待优化阈值参数的过程如下：

(1) 冬小麦分布提取模型阈值参数上下界确定。统计采集的冬小麦地面训练样本NDVI值，计算冬小麦空间分布提取模型中各个判别函数对应的阈值(T₁, T₂,···, T₆)分布范围(图7)。由于传统阈值法中各关键物候期对应阈值一般采用统计训练样本后得出的最小值或者在最小值附近多次尝试，寻找精度最高的一组，因此本研究阈值参数上下界选择设置在最小值附近。又由于SCE-UA优化算法的最大特点是对待优化参数初始值不敏感，可避免优化过程对先验知识的过分依赖，因此，本研究只需保证最优阈值在确定阈值大致范围内即可，所以最终本研究将统计所得最小值设为下界，其下四分位数设为上界。将各关键物候期对应的下四分位数NDVI值到NDVI最小值分别作为参数优化的上下界。

图 7 冬小麦训练样本关键物候期阈值范围分布 Figure 7 Range of threshold value for training samples in critical growth stage of winter wheat

(2) 全局优化算法SCE-UA主要参数确定。在SCE-UA 算法的具体使用中根据待优化模型的复杂程度和参数个数设置相关参数，其中大部分参数取值主要采用已有研究的默认值(宋星原 等，2009)，但复合型的个数p是唯一需要单独考虑的参数。本研究中，据已发表文献(张景 等，2014)，取p=7。根据算法创始人的推荐(Duan 等，1993)，其余各参数缺省取值为 $\beta = 2n + 1$ ， $m = 2n + 1$ ， $\alpha = 1$ ， $q = n + 1$ 。本研究利用的SCE-UA算法中各参数具体含义及其取值结果见表3。

表 3(Table 3)

表 3 复合型混合演化算法中有关参数的取值 Table 3 Parameters’ values of SCE-UA method in this study 待优化参数个数n 复合型个数p 复合型中顶点数m 子复合型顶点数q 复合型演化次数β 子复合型演化次数α 6 7 13 7 13 1

表 3 复合型混合演化算法中有关参数的取值 Table 3 Parameters’ values of SCE-UA method in this study

(3) 目标函数的建立。本文为了评价利用模型提取值与实测值之间的吻合程度，将区域冬小麦种植面积统计数据作为总量控制参考标准，与根据建立的冬小麦空间分布提取模型，在SCE-UA全局优化算法支持下，提取的冬小麦种植面积的差值的绝对值建立优化目标函数：

式中，目标函数值y是利用遥感提取的区域冬小麦空间分布总量与区域冬小麦面积总量控制参考标准间差异绝对值的和，X_simi是区域内基于遥感利用模型提取的冬小麦空间分布进行像素统计后得到的面积总量，X_obsi是区域内冬小麦种植面积的总量控制参考标准(本研究将区域冬小麦种植面积统计数据作为总量控制参考标准)，n为研究区域个数(本研究中，n=1)。

(4) 迭代停止规则确定。当优化过程出现以下3种情况之一时停止优化过程：1)目标函数在临近5个最优目标函数值之差的绝对值小于0.001，即无法显著改善目标函数值；2)目标函数的最大迭代次数不超过10000次；3)待优化参数连续迭代5次后无显著改变。此时认为目标函数达到最优，迭代停止，输出最优阈值。

(1) 遥感分类精度验证指标。本研究中，主要采用野外采集样本点对冬小麦空间分布提取结果进行精度验证，该方法是遥感分类和冬小麦面积提取精度评价中常用的检验方式(Conese和Maselli，1992；Janssen和van der Wel，1994；Knick 等，1997)。研究中，利用地面点样本数据，主要采用总量精度、总体精度、Kappa系数、制图精度和用户精度等5个指标对本研究中获取的冬小麦空间分布结果进行精度客观评价。其中，总量精度计算公式为：

式中，K为总量精度(%)，X_simi为基于遥感提取冬小麦空间分布进行像素统计获得的面积总量(hm²)，X_obj为研究区内冬小麦种植面积的统计数据(hm²)。其他指标具体定义及计算方式可参考相关参考文献(Congalton，1991)，此处不再赘述。

(2) 其他分类方法对比。为进一步验证本研究所提阈值模型参数优化提取冬小麦空间分布的方法效果，除进行所提方法自身精度验证外，本研究利用GF-1/WFV遥感数据、地面训练样本数据和地面验证数据，同时开展与其他方法冬小麦空间分布提取精度对比研究。其他方法主要包含传统阈值方法，还包括较为常用且精度较高的冬小麦分布提取方法，如最大似然法(Solberg 等，1994；杨鹏 等，2000；黄健熙 等，2017)、支持向量机等方法(Mazzoni 等，2007；Mathur和Foody，2008；贾坤 等，2011)。有关两种方法的介绍在此不多作介绍，详细可见相关参考文献。

本研究在冬小麦种植面积统计数据总量参数控制下，利用章节3.4所述阈值法参数优化过程，在全局优化算法SCE-UA支持下，对阈值法冬小麦空间分布提取模型的参数进行优化，最终得到实验区景县2015年冬小麦空间分布提取模型最优阈值参数T_i，具体结果如表4所示。其中，当像元满足第4时相NDVI值大于0.46且第11时相的NDVI值大于0.73时，判断该像元为冬小麦；当像元满足第4时相比第2时相NDVI值大0.30以上并且第11时相比第9时相NDVI值大0.44以上时，判断该像元为冬小麦；当像元满足第4时相比第6时相NDVI值大0.39以上并且第14时相比第15时相NDVI值大0.37以上时，判断该像元为冬小麦；剩余其他像元归为非冬小麦。利用上述阈值法冬小麦空间分布提取模型最优参数，可实现实验区景县冬小麦空间分布高精度提取。通过对本研究利用上述冬小麦空间分布提取模型优化参数获取的2015年冬小麦空间分布结果分析可知，提取面积与县域冬小麦面积总量控制参考标准间差异y仅为160 m²。可见，通过全局优化算法和总量控制获得的阈值法冬小麦空间分布提取模型阈值参数确定的准确性和自动化程度均得到较大程度提高。

表 4(Table 4)

表 4 冬小麦识别模型阈值参数优化结果(2015) Table 4 Parameter optimization results of threshold method for winter wheat identification (2015) 阈值参数 T1 T2 T3 T4 T5 T6 下界 0.46 0.69 0.20 0.39 0.34 0.26 上界 0.55 0.78 0.31 0.51 0.43 0.39 优化阈值 0.46 0.73 0.30 0.44 0.39 0.37

表 4 冬小麦识别模型阈值参数优化结果(2015) Table 4 Parameter optimization results of threshold method for winter wheat identification (2015)

根据4.1节中冬小麦空间分布提取阈值法最优模型对实验区景县冬小麦空间分布进行提取，具体空间分布结果如图8(a)所示。从图8(a)中的冬小麦空间分布上看，2015年景县冬小麦主要集中分布于东部和东南部，西北部冬小麦分布较为破碎。通过统计可知，利用全局优化算法和总量控制下阈值法提取景县全县冬小麦种植面积为40917.98 hm²，县级冬小麦种植面积统计数据为40918.00 hm²，两者间相对误差仅为0.01%，即总量精度达到了99.99%，这说明采用作物面积统计数据总量控制和全局优化算法支持下的冬小麦空间分布提取模型实际提取冬小麦总量精度达到很高水平。最后，本研究将分类后的结果与采集地面验证样本计算混淆矩阵，对方法进行精度验证，具体如表5所示。验证结果表明，本文方法在保证较高总量精度的同时，总体精度为97.03%，Kappa系数为0.94，均达到了较高水平。其中，冬小麦和其他类型的制图精度分别为97.48%和96.00%，用户精度分别为98.25%和94.30%，实现了高精度冬小麦空间分布提取的目标。以上结果证明了阈值参数优化法中利用作物面积统计数据总量控制和全局优化算法的有效性。

图 8 不同冬小麦识别方法结果 Figure 8 Classification results of different winter wheat identification methods

表 5(Table 5)

表 5 不同冬小麦识别方法混淆矩阵 Table 5 Confusion matrix of different winter wheat identification methods 分类方法 地物类型 总体 精度/% Kappa 制图 精度/% 用户 精度/% 阈值模型参 数优化方法 冬小麦 97.03 0.94 97.48 98.25 非冬小麦 96.00 94.30 传统阈值法 冬小麦 92.48 0.82 95.30 94.00 非冬小麦 86.00 88.84 支持向量机 冬小麦 96.88 0.93 95.41 100.00 非冬小麦 100.00 91.12 最大似然法 冬小麦 94.60 0.88 94.00 98.18 非冬小麦 96.00 87.43

表 5 不同冬小麦识别方法混淆矩阵 Table 5 Confusion matrix of different winter wheat identification methods

为进一步验证本研究所提阈值模型参数优化提取冬小麦空间分布的方法效果，并从多个方面进一步证明本研究方法的有效性和可靠性，本研究开展所提阈值模型参数优化方法与传统阈值方法、支持向量机和最大似然法等方法冬小麦提取精度对比研究，各种主要方法获得的冬小麦空间分布结果如图8所示，具体各分类方法分类结果精度评价如表5所示。其中，在传统阈值法中，本研究阈值确定主要通过对采集的冬小麦所有训练样本NDVI值进行统计，确定冬小麦各关键物候期对应的最小NDVI值作为该时相阈值。支持向量机法以及最大似然法均利用ENVI5.1软件中监督分类模块对冬小麦抽穗期GF-1/WFV单幅影像(2015-04-22)进行分类。其中，最大似然法分类采用ENVI5.1软件系统默认值(王增林和朱大明，2010)；支持向量机核函数选择最为常用的高斯径向基核函数RBF(Radial Based Function)，内核函数(Gamma in Kernel Function)值为输入影像波段数的倒数，本研究输入影像为高分一号，其波段数为4，因此，内核函数取值为0.25；惩罚参数(Penalty Parameter)反映对误差的宽容度，这个值设置越高，说明越不能容忍出现误差，但是将造成负面影响产生严重的“椒盐”现象，根据经验设置为100(骆剑承 等，2002)。

从本研究阈值模型参数优化法与传统阈值法分类结果对比图(图8)中可以看出，传统阈值法虽操作简单，运算速度快，但由于阈值设定仅根据经验并采用简单统计方法确定，无法准确去除其他背景地物的影响，导致分类结果冬小麦种植面积总量远高于实际冬小麦种植面积总量，且存在大量破碎的小斑块，将非冬小麦错分为冬小麦的现象严重。支持向量机法和最大似然法尽管分类精度也达到一定水平，但上述两种方法主要基于单一时相影像进行分类，辨识光谱相近的两种或多种冬小麦时易受“同谱异物”等原因的影响，一定程度阻碍分类结果精度的进一步提高。此外，与本研究阈值模型参数优化法结果相比，支持向量机法和最大似然法获取的冬小麦种植面积总量明显高于实际种植面积总量的问题，这也体现了本研究所提阈值模型参数优化法采用作物面积统计数据总量控制进行优化参数的优势。

由各分类方法结果精度评价结果看(表5)，从总体精度、Kappa、制图精度和用户精度等指标衡量，传统阈值法精度均低于本研究所提阈值模型参数优化法、支持向量机和最大似然分类方法，而支持向量机和最大似然法除了总体精度和Kappa系数较高外，用户精度也达到了较高的精度水平，证明支持向量机和最大似然法也是分类效果较优的作物分类算法，这与很多学者的结论是一致的(孙坤和鲁铁定，2017；王小明 等，2013)。另外，从表5也可看出，本研究方法分类结果比传统阈值法、最大似然法和支持向量机法的总体精度和Kappa系数均有所提高，且制图精度也达到最高水平。其中，总体精度分别提高4.55%、2.43%和0.15%，Kappa系数分别提高0.12、0.06和0.01，制图精度分别提高2.18%、3.48%和2.07%，尽管本研究所提阈值参数优化法的用户精度比支持向量机略低，也达到98%以上，可以满足用户精度需求。相较于传统阈值法而言，本研究所提阈值参数优化法的用户精度提高了4.25%，也充分表明了本研究对传统阈值法改进的有效性。尽管支持向量机法在本实验中用户精度高于本研究方法，但是由于该算法对训练样本依赖性强，训练样本数目越大，算法计算复杂程度越高。因此对于大规模训练样本，利用该算法分类运算量大，分类耗时长，因此导致该方法在大尺度区域应用时，前期工作量大，计算效率低，而且若在没有充足训练样本的情况下，分类精度容易偏低，使得该方法在大范围、长时间农作物空间分布制图中应用难度增加。最大似然法分类对于地面训练样本同样具有很强的依赖性，该算法要求采集的训练样本光谱特征值服从正态分布，但由于实际应用中遥感信息的统计分布很难满足正态分布特点，因此当区域内选择的样本分布较为离散从而未能满足最大似然法的假设分布条件，或者样本的选取不具有代表性，就会导致分类结果精度下降，这在一定程度上也增加了区域应用中取得较高分类精度的难度。

综合看，所提阈值参数优化法不仅比传统阈值法分类精度有较大提高，而且与支持向量机和最大似然法等相比，所提分类方法在保证了冬小麦种植面积总量精度的基础上，同时分类空间精度也达到较高水平，可以满足大范围、长时间农业科学基础数据中农作物空间分布制图的精度需求；并且本方法继承了传统阈值法提取冬小麦空间分布操作简单、计算速度快、能较好地适应当地条件的优点，对于开展大尺度区域作物空间分布制图有良好的适用性。其中，尽管面积提取方法的模型和输入数据时相略有不同，但是，通过对比足以证明本研究所提作物面积提取方法常用方法间的精度都达到了较高水平，且本研究的大部分精度指标均高于上述对比方法，从而进一步证明本研究作物面积提取方法的有效性和可靠性。

此外，为进一步对本研究所提阈值模型参数优化方法和传统阈值法、最大似然法、支持向量机法的作物识别结果进行更加细致的空间差异分析，本研究将所提阈值模型参数优化方法与其他方法的作物空间分布结果之间进行空间叠加对比(图9)，并对结果进行分析。从冬小麦分类结果差异的空间分布看，本研究所提阈值模型参数优化方法与其他方法间分类结果差异较大的区域主要分布在景县研究区西北部，如图9中的局部放大区域所示，主要原因是由于该县西北部地区冬小麦种植集中连片区域相对较少，地块较为破碎，导致纯净、混合像元共存现象增加(朱爽和张锦水，2016)，从而导致西北部地区各分类方法冬小麦识别结果差异较大。由冬小麦分类结果间差异统计结果可知，传统阈值法分类结果面积总量远高于本研究方法，分类结果与本研究方法空间分布不一致的比例达到17.59%，其原因是传统阈值法将非冬小麦错分为冬小麦的现象较为严重，导致分类结果中冬小麦面积总量远高于真实总量，而本研究方法识别结果与常用较优的支持向量机法和最大似然法作物识别结果空间一致比例均达到90%以上，这说明本研究所提方法不仅能够改进传统阈值法阈值设置主观性强、智能性低等弱点，而且分类结果精度得到较大水平提高。

图 9 本研究方法与其他分类结果空间叠加对比 Figure 9 Spatial overlay results of classification results by the study method with other methods

通过以上分析可知，与传统阈值法相比，本研究方法利用SCE-UA全局优化算法对阈值参数进行优化，一定程度上克服了传统阈值法阈值确定主观性和不确定性都较强的缺点，并且一定程度上避免了以往阈值模型中人工试错确定阈值需耗费大量工作量的缺点，从而实现了阈值法阈值自动优化的冬小麦空间分布提取，提高了阈值法分类中阈值选取的智能性与准确性，较大提高了作物空间分布制图的精度和工作效率。同时，本研究方法与常用的支持向量机法和最大似然法等分类方法相比，不仅继承了传统阈值法提取冬小麦空间分布操作简单、计算速度快、能较好地适应当地条件的优点，作物空间制图精度也达到较高水平，实现了高精度冬小麦空间分布提取的目标，为开展大范围、长时间农业科学基础数据中农作物空间分布制图提供了方法借鉴。

(1) 研究方法应用前景。与一般农作物面积遥感早期识别和动态监测不同，本文的主要目的为获取大范围、长时间和高精度农作物空间分布基础农业科学数据进行方法探索。由于农作物空间分布基础农业科学数据要求的时效性低且一般针对历史时段进行作物空间分布制图，所以尽管统计数据具有一定的滞后性，本方法在农作物空间分布基础农业科学数据制图中仍具有广泛的应用前景。与一般基于遥感信息和非遥感信息进行多源数据融合农作物空间分布制图方法相比，本方法在保证数量精度较高的情况下，也能获得较高的空间位置精度。因此，本研究方法也可为农作物面积农业统计数据降尺度恢复重建和大范围区域冬小麦空间分布制图研究提供一定技术参考。

(2) 研究方法的不足与改进。目前基于全局优化算法支持下的阈值模型参数优化法提取冬小麦分布信息取得了较好结果，但是研究中仍然存在一些问题或不足。第一，由于受云雨天气的影响，研究中较难获得等时间间隔且云量低的影像，因此，目前研究仅采用了非等时间间隔影像；第二，本研究目前只在县域范围进行研究，且研究对象为种植结构简单的冬小麦，利用的数据为较高空间分辨率的光学遥感数据，且仅采用了单一遥感特征参数NDVI进行阈值法冬小麦面积提取；第三，本研究采用统计数据作为全局优化的总量控制标准参考，由于统计数据具有一定的滞后性，导致本文目前仅用过去年份的数据开展农作物空间分布制图方法应用验证研究，而对本方法适用的时间范围扩展性有待进一步开展研究；此外，本文目前只侧重于全局优化算法支持阈值检测的冬小麦空间分布制图方法研究，而对方法精度的可能影响因素与控制措施研究存在一定不足，如不同物候特征时间组合、统计数据误差等因素对作物空间分布制图精度的影响。同时，本研究中SCE-UA全局优化算法的参数设置目前采用了该算法推荐的默认参数，并未开展优化算法参数变化对本文研究结果的影响研究。

针对方法存在的可能问题或不足，在基于总量控制和全局优化的前提下，今后在以下几个方面进行深入研究：(1)利用NDVI最大值合成方式获取等时间间隔的遥感数据，开展对冬小麦分布提取精度影响研究；(2)利用其他中低、高空间分辨率遥感数据在更大区域进行冬小麦分布提取的有效性、结果精度验证和优化参数结果应用及评价等研究。另外，在现有计算框架下进一步讨论应用其他遥感特征参量(如增强型植被指数EVI、宽范围动态植被指数WDRVI等)进行冬小麦空间分布制图研究；(3)进一步探究利用本研究方法获得研究区历史长时序阈值模型参数规律，通过对算法阈值标准化，开展将作物识别模型历史阈值用于现在和未来作物空间分布监测及预测的方法研究，明确阈值应用的时间扩展性和有效性，从而使得本方法具有更广泛的适用性；(4)开展统计数据自身误差、物候特征组合以及其他总量控制(如基于抽样技术确定区域总量)的全局优化标准参考等对本研究方法的精度影响研究，从而提高本研究方法的精度，降低可能影响因素对方法精度的不确定性影响，从而进一步增强本方法的适用性；(5)进一步开展SCE-UA算法不同参数设置对本研究方法的精度影响研究。

为获取大范围、长时间和高精度农作物空间分布基础农业科学数据，本研究针对阈值法农作物识别过程中阈值设置存在灵巧性差和自动化程度低等弱点，在中国粮食生产基地黄淮海地区典型实验区内，开展了利用全局优化算法改进基于阈值检测的冬小麦分布制图方法探索研究。主要结论如下：

(1) 在改进传统阈值方法基础上，本文引入冬小麦种植面积统计数据作为全局优化的总量控制参考标准，提出了基于时序NDVI遥感数据的阈值模型参数优化冬小麦空间分布制图方法，成功提取了区域冬小麦作物空间分布信息。通过地面验证表明，本研究所提阈值模型参数优化方法冬小麦提取分类结果精度均达到较高水平，可以满足大范围、长时间农业科学基础数据中农作物空间分布制图的精度需求。与传统的阈值法、最大似然和支持向量机等分类方法相比，本研究所提阈值模型参数优化法区域冬小麦作物分类总体精度和Kappa系数分别都有所提高，其中，总体精度分别提高4.55%、2.43%和0.15%，Kappa系数分别提高0.12、0.06和0.01，这证明了本研究所提基于统计数据总量控制的全局优化算法阈值模型参数优化作物分布信息提取方法的有效性和可行性，可获得高精度冬小麦作物空间分布制图结果。

(2) SCE-UA全局优化算法可实现阈值法冬小麦分布提取模型中阈值自动寻优的目的。本文首次尝试将SCE-UA全局优化算法引入冬小麦分布提取阈值模型进行参数自动优化，将遥感提取的实验区冬小麦空间分布总量与实验区冬小麦面积总量控制参考标准间差异绝对值的和作为目标函数，对待优化的阈值参数进行自动优化。结果表明，利用SCE-UA全局优化算法对阈值模型阈值参数进行优化选取，在一定程度上避免了传统阈值模型阈值确定通过人工试错方法特有的主观性和不确定性，不仅大大提高了工作效率，并且提高了冬小麦提取阈值模型参数确定的精确性和自动化程度，对提高中国冬小麦空间分布制图精度和自动化水平具有重要意义。

(3) 从本文开展的作物面积统计数据总量约束下全局优化阈值的冬小麦分布制图方法和精度验证研究可知，本方法由于采用统计数据作为总量控制参考标准，在保证总量精度的基础上，位置精度也达到较高水平，可以获得总量精度和位置精度双优的作物空间分布制图结果。另外，本文所提阈值模型参数优化冬小麦分布制图方法继承了传统阈值法提取作物空间分布操作简单、计算速度快、区域适用性强等优点，不仅作物空间制图精度可达到较高水平，实现了高精度冬小麦空间分布制图的目标，而且在大尺度区域作物空间分布制图中具有良好的适用性，为开展大范围、长时间农业科学基础数据中农作物空间分布制图提供了方法借鉴。