本研究针对中国近年来重点发展的创新型区域收入保险缺少第三方实时客观产量数据的问题，引 入了卫星遥感估产技术，探讨其应用模式和适用性。以山东省嘉祥县大豆区域收入保险为例，基于哨兵2号 卫星遥感数据提取大豆种植地块，计算归一化植被指数NDVI和作物生理参数，结合气象卫星遥感数据与实 地抽样测产数据，建立了多参数线性回归模型估算大豆产量。研究结果显示，卫星遥感获取的研究区大豆 种植面积为 124 km2 ，与当地农业局上报的 127 km2 相差 3 km2 ；采用实测地块验证，种植分布地块遥感识别 精度达90%；产量估算结果显示，2018年8月23日大豆结荚期的NDVI和9月7日大豆鼓粒期的NDVI对大豆 单产的解释度最佳，多参数回归模型计算全区平均产量为244,500 kg/km2 ，与常年299,800 kg/km2 相比，体现 了受灾严重的农情；产量估算数据与实测数据之间的回归系数达0.92，可满足应用需求。结果表明，基于哨 兵2号卫星遥感数据能够准确识别研究区大豆种植分布，并能在大豆收获后最快一周完成产量估算，指导保 险公司的理赔工作。

1 引言

农业保险是农业风险管理的重要工具，能够 减少农民的直接经济损失、减轻政府和社会各方 压力、提高农业综合实力，保障农业和农村社会 持续稳定发展，受到各国政府的高度重视。中国自 2007年开始实施保费补贴型政策性农业保险， 主要险种是政策性直接物化成本保险，主保个体 产量。目前中国农业保险保费已居世界第二，仅 次于美国。然而，随着农业保险的发展和土地政 策的改革，政策性直接物化成本保险因存在保额 低、验标验险耗财耗力、道德风险高等缺陷，与保险公司和农业经营户，尤其是新型经营主体抵 抗生产风险的需求存在较大差距。随后中国陆续 进行试点试验和增加了指数保险、价格保险、收 入保险、区域产量、区域收入保险等不同险 种［1-3］ 。其中，区域收入保险是中国近年重点发 展的创新型险种，它综合了区域平均产量与农产 品价格，形成了一个对区域平均收入进行保护的 险种，相对降低了验标验险的成本，还可规避基 差风险和道德风险，在中国山东、山西、辽宁等 地已经开始广泛试点，并得到政府部门和保险公 司的大力推广。

开展区域收入保险需要两个重要的数据支 持，一是收获期的区域平均价格，二是收获期的 区域平均产量。区域的单位一般为县、乡镇、甚 至村级行政单元。其中，区域价格可以参考物价 局的实时价格数据，或者定点价格监测数据。区 域产量需要第三方公布的当年产量数据，但中国 统计局官方发布的产量数据一般会滞后 1~3年。 若按照传统的保险运营模式，依然存在信息不对 称和协议理赔等问题，因而需要更客观及时的估 产方法，这就需要遥感技术的介入。

近年来，随着农业保险对信息技术的需求越 来越高，遥感信息技术以其客观特性逐渐被应用 到农业保险的验标和定损等环节 。遥感技术 在验标环节回答田间为何种农作物、种植在哪 里、是否与承保面积相符等问题其实是遥感作物 识别的问题，已经比较成熟 。遥感技术在理赔 环节也有不少辅助理赔的成功案例 。遥感还 能为创新型指数保险险种提供客观指数 ，在农 险中的应用益处已得到保险各方的认可 。然 而将遥感产量估算结果直接用到农业保险理赔的 案例却不多见。

从遥感产量估算的技术可行性角度看，作物 产量估算的遥感估算在小尺度 （县域尺度） 的研 究文献已经有很多，包括小麦、玉米 和 水稻卫星遥感估产 等。相比粮食作物而言， 大豆等油料作物的产量估算相对复杂，因此研究 也更精细。Yu 等及 Maimaitijiang 等基于无人机高光谱数据，Gusso 等基于 MODIS 的 大尺度数据，以及Kross等基于5 m分辨率的 Rapideye数据开展了不同尺度的大豆地上生物量 或产量的估算研究。国内大豆的产量估算研究也 比较成熟。南京农业大学盖钧镒院士及其团队依 托大豆育种研究，开展了许多与产量估算相关的 研究，涉及了低空、地面和高空遥感。其中，开 展基于地面光谱或无人机遥感的作物生理参数， 如叶面积指数、光谱参数及其与生物量的关系研 究的有齐波等 、陆国政等；开展基于地面 光谱或无人机遥感的大豆籽粒产量研究的有吴琼 等、张宁等和赵晓庆等。国内其他作者 也开展部分基于卫星遥感的大豆叶面积指数或产 量的大豆估产研究，如张淮栋等基于高分 2 号 NDVI 数据进行大豆估产，发现 NDVI 与大豆 产量的相关性不高，存在一些高产而 NDVI极低 的地块。

综合以上研究，当前大豆产量遥感估算的技 术参考研究充分，其中无人机遥感是应用的前 沿，高光谱遥感是研究的前沿。针对区域收入保 险需要估算的区域产量可能是乡镇甚至是村等更 小的行政单元，承保的业务范围又可能是一个县 或者多个县，且对时效的要求极高，而对投入成 本又有精打细算等问题，本研究选择以卫星遥感 作为数据源，以哨兵 2号卫星数据为主，结合气 象卫星数据和地面测产数据，建立多参数线线性 回归模型估算县级以下尺度的大豆收获产量，并 提供一个遥感估产的行业应用案例，以供农业保 险领域各方参考。

2 研究区概况

研究区位于山东省济宁市西部嘉祥县，东经 116°06′~116°27 ′ ，北纬 35°11′~35°38 ′ ，属黄河 冲 积 平 原 ， 总 面 积 约 975 km2 ， 地 理 位 置 见 下图1。 嘉祥县大豆主要种植于该县北部乡镇。2017 年，该县被列为国家首批五大良种大豆繁育基地 之一。由于气候、地理位置的关系，加上多年制种企业的发展，当地种业发达，育种、繁种都形 成了规模，并辐射和带动了周边地区。 嘉祥县大豆的生育期从 6月中旬 （一般端午 之后，收了麦子等待土壤墒情合适） 陆续开始种 植，至9月下旬或者10月上旬大豆收获结束，整 个生育期100~120天。

3 研究方法

3.1 研究内容和技术路线

区域收入保险的定义如下：在保险期间，按 照划定的保险区域，由于各类自然灾害和非人为 原因的病虫害，导致投保作物所在区域内单产产 量降低，或由于市场供应变化导致投保作物的价 格下降，抑或上述两种情况同时发生造成区域保 险作物的实际收入低于保单约定的保障收入时， 保险人按照合同约定对被保险人进行赔付。 保单约定的单位保障收入亦即投保人按照一 定的保障水平投保的预期收入，公式如下。

预期收入= （预期价格+现货基差补偿价 格） ×预期单产

预期价格使用大连商品交易所保险往年 11 月和 1 月到期的黄豆 1 号合约在 6 月保险销售前 五个交易日的日收盘价均值。现货基差补偿价格 即约定为0.5 CNY/kg。预期单产为当地前五年单产数据中 （去掉最大与最小值后） 的三年平均 值，结合实地调研情况，约定预期单产。保障水 平是对预期收入的保障，可以是 0%~100%，考 虑到可能存在的道德风险和逆选择，区域收入保 险的保障水平理论不得超过90%。

当实际区域收入低于保障收入时触发理赔 （公式 （2））。

理赔时，实际区域价格采用当地物价局的一 段时间的大豆平均销售价格。研究区当地大豆销 售期始于每年 10 月初，虽然大豆易存储，可长 期销售。但考虑到前两个月是销售旺季，故价格 采集时间为每年10月11日到12月10日。实际区 域产量则采用遥感估算所得的估产结果。当遥感 估算的当年产量与当年平均价格形成的实际收入 低于预期收入时，触发理赔。

大豆产量估算主要包括野外实地数据采集、 大豆种植地块的遥感识别、大豆产量的遥感估算 三项。研究技术路线如图 2所示。综合使用卫星 导航定位技术和互联网技术采集地块样本，使用 基于面向对象的图像处理技术，基于单时相多光 谱卫星遥感数据提取了耕地地块，在地块边界限 制下，再结合多时相多光谱卫星遥感数据的ND‐ VI，采用决策树分类区分目标作物和其他作物， 同时区分目标作物类间差距，亦即受灾情况，得 到边界整齐的作物地块。然后，结合专家实地测 产数据，在作物收获后最快一周内建立经验统计 模型，估算地块产量。

3.2 研究数据

研究数据包括实测数据、卫星遥感数据和其 他辅助地理信息等。其中，辅助地理信息主要包 括土地利用数据和行政区划数据，土地利用数据 采用 2010 年全球 30 m 的地表覆盖数据 （下载地 址为：http：//）。

3.2.1 实测数据

实测数据包括实测地块样本和实地测产数 据。在大豆生长的关键期共采集了三次产量，时间 分 别 为 2018 年 8 月 24 日 ， 9 月 19 日 和 10 月8日。

（1） 地块样本采集。 地块采集采用中国农业科学院农业信息研究 所风险分析与管理中心自主研发的“农业遥感移 动采集平台”app （图 3），勾画地块四至图形， 辅以照片信息和实地考察印象，记录至同一个地块下作为作物地块的属性信息。采集大豆样本的 同时，也采集部分其他作物样本用于辅助分类。 本研究共采集地块 192 个，其中大豆地块 128 个，涉及总面积 2,983,075 m2 ，相当于总种 植面积 （127 km2 ） 的 2.5%。平均地块大小为 23,305 m2 。其他作物地块 64个。

（2） 产量采样。主要以理论产量采样法进行采集。在农业保险的应用中，一般需要邀请与 作物相关的农学专家以及当地农业技术推广站相 关专家来参与和督导工作。本研究邀请了山东省 济宁市农业科学院大豆研究所的专家 2人和嘉祥 县农技站站长等2人参与。

本研究实际进行了两次采样，分别为：大豆 鼓粒期，采样时间为 2018 年 8 月 24 日至 26 日； 收获期，采样时间为 2018 年 10 月 8 日至 10 日。 其中鼓粒期的产量采样主要是计算豆粒数，即在 地块采样的时候，把已经鼓粒的粒数予以清数和 记录。

收获期产量采样根据以遥感识别的大豆地块 作为样点选取的依据，选取拟采样地块，在实际 采样时以寻访和收割方式并行。其中寻访的主要 原因是：研究区域为大田作物而非定点控制实 验，可能赶上已经收割的田块，若是田间有农 户，则予以寻访。对没有收割的大豆，采用收割 法。在选择的地块上，进行三次重复取样操作。 具体操作步骤为：①量测行距；②数出5 m双行 的株数；③随机连续取下 10 株，带回实验室脱 粒烘干至水分剩余15%，并称百粒重，从而得到 单位面积株数、单位面积总粒数和百粒重三要 素。研究共采得产量样点 99 个，其中鼓粒期数 粒数的样点 19 个，收获期产量测量样点 80 个。 样点分布如下图4所示。

3.2.2 遥感数据

遥感数据包括 2018 年大豆生长期 （6~10 月） 的哨兵 2号卫星影像数据、表征气象的测雨 雷 达 （Tropical Rainfall Measuring Mission，TRMM） 数据和 MODIS 地表温度数据。由于需 要估算乡镇级产量，县级气象站数据不能满足需 求，在本研究中没有被采用。这三个数据均为免 费数据。其中，哨兵 2 号卫星数据和 MODIS 数 据均下载于美国地质调查局 （United States Geo‐ logical Survey，USGS，https：//glovis.usgs.gov/）。 哨兵 2号分 A 星和 B 星，两颗卫星组成重访 周期5天的观测网络。本研究一共下载了11景哨 兵2号数据，时相如表1所示。 选择所有11个时相的哨兵2号卫星数据第8、 4、3 波段作为红绿蓝波段进行组合，得到 RGB 假彩色合成影像，示意如下 （图5）。

MODIS 地表温度使用 8 天最大值合成产品， 分辨率为 1 km，下载网址为 . gov/dataset_discovery/modis/。地表温度是地表综 合上行热辐射的一个反应，可以综合反映土壤墒 情、大气显热和潜热等情况，本研究用于表征气 象因子中的温度。采用MODIS的昼夜温度产品， 分辨率为 1 km，时间为 8 天；每 8 天两个变量， 分别为 8 天白昼温度最大值和夜间温度最大值， 从2018年6月1日开始至2018年9月29止，共30 个变量，加上整个生长季白昼极端高温和低温， 夜间极端高温和低温4个变量，共34个变量。 TRMM 数据分辨率为 0.05°，时间分辨率为 3 h，单位为 mm，下载于 . nasa. gov/ missions/trmm。对 TRMM 数据进行逐旬总降雨 量合成，得到6~9月4个月，每月三旬，共12个 降雨变量，加上整个生长季的最大降雨量、最小 降雨量和总降雨量共15个降雨变量。

3.3 大豆地块遥感提取

经过图像对比，发现时相为8月23日的哨兵 2号卫星影像上，大豆与其他作物区别最大，采 用简译软件多尺度分割算法对其进行面向对象分 割，对所得的对象进行人工选取样本类别，结合 最大似然法进行监督分类，得到耕地与非耕地。 在所得耕地地块的基础上，结合实测地块样本，分析不同作物和不同受灾的大豆在多时相哨兵 2 号卫星数据的 NDVI 特征进行大豆识别 （如图 6 和图7所示）。

由图6可知，8月8日和8月23日两个时相是 区分大豆和其他主要作物的关键时相，基于这两 个时相采用 ISOData算法进行非监督分类，得到 大豆和其他作物的分类图。 不同受灾受害大豆的区分有助于产量估算精 度的提高，因而根据样本的标注，如大豆长势如 何、是否受灾受害、是否早熟或者晚熟等情况提 取出大豆类间的长势差异。 不同灾害长势大豆的 NDVI 如图 7 所示。据 此，可以看出受旱可收获的大豆，主要体现在 8 月3日之前NDVI很低，受病虫害大豆在8月8日 时相前后开始 NDVI 都很低，受涝大豆则是 8 月 23日之后NDVI值变低。根据这几个时相，调整 阈值，建立决策树，区分不同长势大豆。

3.4 植被指数和作物生理参数

基于哨兵 2号数据计算常用的植被指数和作 物生理参数，包括NDVI，叶面积指数 （Leaf Ar‐ ea Index，LAI）、叶绿素含量 （Chlorophyll Con‐ tent in the Leaf， CAB）、 光 合 有 效 辐 射 分 量 （Fraction of Absorbed Photosynthetically Active Radiation，FAPAR）、植被盖度分量 （Fraction of Vegetation Cover， FCover） 和 叶 面 水 分 含 量 （Canopy Water Content，CW） 等 5 个生理参数。 其中NDVI用第8波段和第4波段计算。5个生理 参数具体方法采用哨兵 2 号卫星数据分发中心 ——欧空局 （European Space Agency） 提供的软 件 SNAP中的作物生理参数提取方法得到，方法 参考见文献 ［28］。由于7月与8月数据部分有云 覆盖，因而 7月提取的作物生理参数和反射光谱 参数采用 7 月份最大值合成。8 月、9 月和 10 月 是大豆生长到关键期，每一期大豆生长情况相差 较大，因此 8月、9月和 10月数据全使用，并在 样点提取时剔除有云的点。

3.5 统计分析

采用SPSS软件进行作物生理参数、NDVI与 产量的 Pearson 相关关系分析，再根据相关系数 分析结果选择系数高的因子参与逐步回归，并人工添加气象因子加入回归模型，根据决定系数R2 的变化，确定估产模型进行估产。 在产量的 99 个点中，除去有云的样点，余 下 88个样点，进行相关性分析。由于样点不多， 用其中 4/5 作为建模样点，余下 1/5 作为验证样 点。建模样点与验证样的选取步骤如下：用求余 函数将每个样点在 ArcGIS 的数据表要素编号 FID，除以5，将余数为0的点提取出来，即得到 1/5 的验证样点 17 个，其余 4/5 的 71 个样点参与 回归建模。

3.6 精度评价

（1） 大豆种植情况识别的精度评价。由于大 豆种植没有用监督分类，而是用非监督以及决策 树分类，最后用种实测地块是否落入识别范围作 为精度评价，评价公式如下。

（2） 产量估算结果的精度评价。 用实测产量的 1/5 作为精度评价样本，与实 测样本进行散点图分析，平均估算精度。

3.7 区域收入保险应用

得到地块单产之后，根据保险公司划定的区 域，比如乡镇、村庄、特殊经验区域等进行平均 产量统计。保险公司依据平均产量，再结合区域 价格形成一个实际的区域收入，与保单约定的单 产保障收入进行比较来决定赔付与否。

4 研究结果及分析

4.1 大豆种植情况及精度

大豆分布和受害情况遥感识别结果如图 8所 示。地块面积统计结果显示，嘉祥县 2018 年共 种植大豆种植面积 124 km2 ，与当地政府部门上 报的127 km2 统计数据吻合较好。精度评价结果如下。

其中，85 为长势正常大豆正确识别地块数 量；88 为长势正常大豆实测地块数量；21 为受 旱受害大豆正确识别地块数量；29为受旱受害大 豆实测地块数量。总体上，长势正常的大豆，识 别正确率很高，为 96.59%，受旱受害的大豆由 于表现不一，识别率约为72.4%。 由于总的识别面积为 123.74 km2 ，而受旱受 害的仅有 31.47km2 ，占总面积的 25.4%，因此总 体精度可以计算为：

P总体=72.4%×0.254+97.6%×(1-0.254)=90.45%

4.2 产量影响因子

根据 Pearson 相关分析结果 （表 2），与产量 关系最密切的是 8 月 23 日的 NDVI，其次是 9 月 7 日的作物生理参数，相关系数均大于 0.6。与 产量相关性显著的气象因子有四个，分别为： ①9 月上旬降雨 TRMM2018091，表征鼓粒期到 成熟初期的降雨；②7 月 11 日后 8 天白天最大值 合成值 20180727_Day_LST，表征始花期的极高 温；③7 月 27 日后 8 天白天温度最大值合成值 20180711_Day_LST， 表 征 盛 花 期 的 极 高 温 ； ④ 8 月 28 日 后 8 天 夜 间 温 度 最 大 值 合 成 值 20180828_Night_LST 鼓粒期的极高温等四个因 子相关系数在0.23~0.26之间。

分析 2018年 8月 23日时相与 9月 7日时相的 NDVI与作物生理参数之间的自相关性，得到二 者显著正相关，相关系数除了8月23日时相的水 分含量 （0823Cw） 之外，都超过0.59 （如表3）。 因 此 在 选 用 因 子 时 ， 只 选 了 0823NDVI 和 0823Cw。

4.3 大豆产量模型

在相关系数最高的植被指数或者作物生理参数里选择 3~5 个参数，同时选择 3~5 个气象参 数，进行逐步回归分析，得到两个与产量显著相 关 且 相 关 系 数 都 较 高 的 输 入 因 子 ， 分 别 为 0823NDVI 和 0907LAI，所构建的模型 （公式 6） 对产量的解释度达到 46.4%。其模型标准误差约 为45 kg。再增加4个变量或者8个变量并未增加 解释度，标准误也相差不到1 kg，模型系数摘要 见表4。 因此本研究最后采用公式8进行估产。即： Y=-171.365+313.127×0823NDVI+35.154× 0907LAI,其中，Y 为 666.7 m2 的产量，0823NDVI 为 8 月23日NDVI，表示盛花始荚期的归一化植被指 数；0907LAI为9月7日叶面积指数，表示盛荚鼓 粒期的叶面积指数。

4.4 大豆区域产量估算结果及评价

根据大豆受灾长势监测分级结果以及实地采 样依据可知，受旱受害大豆以及受涝大豆为绝产 大豆，在估产中，将这两类地块先分离出去，最 后得到的区域产量再乘以每个乡镇绝产地块所占 比例，得到平均单产等级图，如下图9所示。整个 区域的平均产量为 244,500 kg/km2 （163 kg/亩）， 与常年 299,800 kg/km2 （200 kg/亩） 相比，体现 了受灾严重的农情。

在大豆产量样本中，用剩余的 1/5 样本作为 验证，得到如图 10 散点图，二者回归系数达0.92，表明产量估算模型的精度较高。

5 结论与讨论

5.1 结论

本研究主要基于哨兵 2号卫星识别大豆种植 分布及受害情况，计算植被指数和作物生理参 数，结合实际测产、气象卫星遥感数据进行多参 数回归的大豆产量建模。结果表明哨兵 2号卫星 数据对大豆地块的识别精度达90%，对产量的估 算精度达92%。 但本研究对于农业保险迫切的行业应用需求 而言仍存在一些缺陷。首先，花荚期、鼓粒期的 影像数据是研究的基本保障，如果这个时期的遥 感数据存在云及阴雨天气影响，对于估产结果精 度和整个应用方案可能产生影响。其次，研究需 要大量的实地样本，主要是产量样本的支持，且 需要较多专家知识参与，并且采样方案需要综合 考虑成本及效率。另外，保险理赔过程涉及当地 农户对遥感的认知和接受程度，还涉及当地的民 风和经济状况，还需要从多个角度，多种方法结 合、多部门数据一起来定约定理赔。

一是野外采样。一般至少需要两次，第一次 在作物生长季中期，主要采集作物的识别样本和 长势、受灾受害情况。第二次是在作物收获期， 与传统的采样一样，尽量根据作物地块识别结 果，全局均匀采样，不同长势均匀覆盖。

二是目标作物遥感识别。对于目标作物识别 而言，目标作物地块的形态、位置和边界很重 要，与其他作物的时间区分很重要，因此目标作 物识别之前的影像时间和波段选择、耕地掩膜和 面向对象的图像分割很重要，可以较好地规避传 统遥感影像分类的像素孤岛等错误。在目标作物 识别过程中，对于类间差距，以及何种长势可能 导致绝产的标注也很重要，对于杂粮作物以及其 他油料、棉花等作物而言，在机理不明确的经验 统计模型中，去除一些异常的长势，如叶子很多 却不结荚的大豆，可以提高估算精度和解释度。

三是产量估算。在作物地块分好类之后，选 取具有一定生理意义的作物光谱参数，加入气象参数，建立模型，并根据前期长势的标注进行产 量的加权调整，使得建模过程更完整更可信。

但是这三个步骤也是区域收入保险运营需要 考虑的成本投入。首先，采样步骤一般需要邀请 与目标作物相关的农学专家参与，进行长势和产 量的测定。其次，基于多光谱卫星遥感的作物识 别，对作物关节生长期的卫星数据极其依赖，若 赶上天气不好，就缺乏可用的多光谱卫星数据， 需要启动其他数据预案，包括付费编程数据和无 人机数据来完成目标作物的识别，此时可能会因 无多时相遥感数据参与产量建模而废弃遥感估产 方法，而只根据作物识别结果增加实地抽样测产 来确定区域产量，这都极大增加保险运营成本。

农业保险承保范围广而细，全国所有省份几乎都参与了农业保险，保单细到每家每户，这给 农业保险定损核损带来了极大困难，这就是农业 保险由传统的个体保险发展到区域保险的主要原 因。但即便是区域类型保险，理赔的区域也需要精细到乡镇，甚至村。由此可见，农业区域产量 或者收入保险急需借助遥感这种大范围全天候客 观可靠的技术来完成产量实时估算的任务。遥感 技术作为农业保险急需的技术，是农业保险行业 发展的新动力［29,30］ 。尽管很难，但是遥感用于农 业保险是趋势，农业保险对遥感的需求有增不 减。遥感应用于农业保险还需要更多的投入应 用、经验总结与凝华。

作者：陈爱莲、李家裕、张圣军、朱玉霞、 赵思健、孙 伟、张 峭

来源：智慧农业

迎大家在评论区交流分享

更多精彩资讯，请持续关注数字农业分会官方账号