农业保险是农业风险管理的重要工具，能够减少农民的直接经济损失、减轻政府和社会各方压力、提高农业综合实力，保障农业和农村社会持续稳定发展，受到各国政府的高度重视。中国自2007年开始实施保费补贴型政策性农业保险，主要险种是政策性直接物化成本保险，主保个体产量。目前中国农业保险保费已居世界第二，仅次于美国。然而，随着农业保险的发展和土地政策的改革，政策性直接物化成本保险因存在保额低、验标验险耗财耗力、道德风险高等缺陷，与保险公司和农业经营户，尤其是新型经营主体抵抗生产风险的需求存在较大差距。随后中国陆续进行试点试验和增加了指数保险、价格保险、收入保险、区域产量、区域收入保险等不同险种。其中，区域收入保险是中国近年重点发展的创新型险种，它综合了区域平均产量与农产品价格，形成了一个对区域平均收入进行保护的险种，相对降低了验标验险的成本，还可规避基差风险和道德风险，在中国山东、山西、辽宁等地已经开始广泛试点，并得到政府部门和保险公司的大力推广。
开展区域收入保险需要两个重要的数据支持，一是收获期的区域平均价格，二是收获期的区域平均产量。区域的单位一般为县、乡镇、甚至村级行政单元。其中，区域价格可以参考物价局的实时价格数据，或者定点价格监测数据。区域产量需要第三方公布的当年产量数据，但中国统计局官方发布的产量数据一般会滞后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。
图1 研究区地理位置示意图
Fig. 1 Location of the study area
嘉祥县大豆主要种植于该县北部乡镇。2017年，该县被列为国家首批五大良种大豆繁育基地之一。由于气候、地理位置的关系，加上多年制种企业的发展，当地种业发达，育种、繁种都形成了规模，并辐射和带动了周边地区。
嘉祥县大豆的生育期从6月中旬（一般端午之后，收了麦子等待土壤墒情合适）陆续开始种植，至9月下旬或者10月上旬大豆收获结束，整个生育期100~120天。
3 研究方法
3.1 研究内容和技术路线
区域收入保险的定义如下：在保险期间，按照划定的保险区域，由于各类自然灾害和非人为原因的病虫害，导致投保作物所在区域内单产产量降低，或由于市场供应变化导致投保作物的价格下降，抑或上述两种情况同时发生造成区域保险作物的实际收入低于保单约定的保障收入时，保险人按照合同约定对被保险人进行赔付。
保单约定的单位保障收入亦即投保人按照一定的保障水平投保的预期收入，如公式（1）所示。
预期收入=（预期价格+现货基差补偿价格）×预期单产
（2）预期价格使用大连商品交易所保险往年11月和1月到期的黄豆1号合约在6月保险销售前五个交易日的日收盘价均值。现货基差补偿价格即约定为0.5 CNY/kg。预期单产为当地前五年单产数据中（去掉最大与最小值后）的三年平均值，结合实地调研情况，约定预期单产。保障水平是对预期收入的保障，可以是0%~100%，考虑到可能存在的道德风险和逆选择，区域收入保险的保障水平理论不得超过90%。
当实际区域收入低于保障收入时触发理赔（公式（2））。
实际区域收入=实际区域价格×遥感区域产量
（3）理赔时，实际区域价格采用当地物价局的一段时间的大豆平均销售价格。研究区当地大豆销售期始于每年10月初，虽然大豆易存储，可长期销售。但考虑到前两个月是销售旺季，故价格采集时间为每年10月11日到12月10日。实际区域产量则采用遥感估算所得的估产结果。当遥感估算的当年产量与当年平均价格形成的实际收入低于预期收入时，触发理赔。
大豆产量估算主要包括野外实地数据采集、大豆种植地块的遥感识别、大豆产量的遥感估算三项。研究技术路线如图2所示。综合使用卫星导航定位技术和互联网技术采集地块样本，使用基于面向对象的图像处理技术，基于单时相多光谱卫星遥感数据提取了耕地地块，在地块边界限制下，再结合多时相多光谱卫星遥感数据的NDVI，采用决策树分类区分目标作物和其他作物，同时区分目标作物类间差距，亦即受灾情况，得到边界整齐的作物地块。然后，结合专家实地测产数据，在作物收获后最快一周内建立经验统计模型，估算地块产量。
图2 农业保险实施技术路线图
Fig. 2 Technology roadmap of agricultural insurance implementation
3.2 研究数据
研究数据包括实测数据、卫星遥感数据和其他辅助地理信息等。其中，辅助地理信息主要包括土地利用数据和行政区划数据，土地利用数据采用2010年全球30 m的地表覆盖数据（下载地址为：http：//www.globallandcover.com/）。
3.2.1 实测数据
实测数据包括实测地块样本和实地测产数据。在大豆生长的关键期共采集了三次产量，时间分别为2018年8月24日，9月19日和10月8日。
（1）地块样本采集。
地块采集采用中国农业科学院农业信息研究所风险分析与管理中心自主研发的“农业遥感移动采集平台”app（图3），勾画地块四至图形，辅以照片信息和实地考察印象，记录至同一个地块下作为作物地块的属性信息。采集大豆样本的同时，也采集部分其他作物样本用于辅助分类。
图3 “农业遥感移动采集平台”app工作页面
Fig. 3 Application interface of “Mobile Collection Platform for Agriculture Remote Sensing”
本研究共采集地块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所示。
图4 产量样点分布图
Fig. 4 Distribution map of yield samples
3.2.2 遥感数据
遥感数据包括2018年大豆生长期（6~10月）的哨兵2号卫星影像数据、表征气象的测雨雷达（Tropical Rainfall Measuring Mission，TRMM）数据和MODIS地表温度数据。由于需要估算乡镇级产量，县级气象站数据不能满足需求，在本研究中没有被采用。这三个数据均为免费数据。其中，哨兵2号卫星数据和MODIS数据均下载于美国地质调查局（United States Geological Survey，USGS，https：//glovis.usgs.gov/）。
哨兵2号分A星和B星，两颗卫星组成重访周期5天的观测网络。本研究一共下载了11景哨兵2号数据，时相如表1所示。
表1 2018年卫星数据分辨率、时相和数量参数表
Table 1 Resolution, acquisition date and numbers of the satellite data used in 2018
|
卫星名称 |
分辨率/m |
时相/MMDD |
主要用途 |
|
哨兵2号 |
10，20，60 |
0614，0624，0714，0724，0803，0808，0823，0907，0922，0927，1002 |
大豆分类 大豆估产建模 |
|
Terra 卫星 MODIS LST |
1000 |
0601—0928每八天合成 |
大豆估产建模 |
|
TRMM测雨雷达 |
约5000~6000 （0.05°） |
0601—0930每十天合成 |
大豆估产建模 |
选择所有11个时相的哨兵2号卫星数据第8、4、3波段作为红绿蓝波段进行组合，得到RGB假彩色合成影像，示意如下（图5）。
图5 哨兵2号卫星数据8、4、2波段红绿蓝组合图像
Fig. 5 RGB composite images of sentinel-2 satellite data 8，4，2 band
MODIS地表温度使用8天最大值合成产品，分辨率为1 km，下载网址为https://lpdaac.usgs.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，下载于https://gpm.nasa.gov/missions/trmm。对TRMM数据进行逐旬总降雨量合成，得到6~9月4个月，每月三旬，共12个降雨变量，加上整个生长季的最大降雨量、最小降雨量和总降雨量共15个降雨变量。
3.3 大豆地块遥感提取
经过图像对比，发现时相为8月23日的哨兵2号卫星影像上，大豆与其他作物区别最大，采用简译软件多尺度分割算法对其进行面向对象分割，对所得的对象进行人工选取样本类别，结合最大似然法进行监督分类，得到耕地与非耕地。在所得耕地地块的基础上，结合实测地块样本，分析不同作物和不同受灾的大豆在多时相哨兵2号卫星数据的NDVI特征进行大豆识别（如图6和图7所示）
图6 不同作物NDVI曲线
Fig. 6 NDVI curves of different crops
由图6可知，8月8日和8月23日两个时相是区分大豆和其他主要作物的关键时相，基于这两个时相采用ISOData算法进行非监督分类，得到大豆和其他作物的分类图。
不同受灾受害大豆的区分有助于产量估算精度的提高，因而根据样本的标注，如大豆长势如何、是否受灾受害、是否早熟或者晚熟等情况提取出大豆类间的长势差异。
不同灾害长势大豆的NDVI如图7所示。据此，可以看出受旱可收获的大豆，主要体现在8月3日之前NDVI很低，受病虫害大豆在8月8日时相前后开始NDVI都很低，受涝大豆则是8月23日之后NDVI值变低。根据这几个时相，调整阈值，建立决策树，区分不同长势大豆。
图7 不同长势大豆NDVI曲线
Fig. 7 NDVI curves of soybeans with different growing stages
3.4 植被指数和作物生理参数
基于哨兵2号数据计算常用的植被指数和作物生理参数，包括NDVI，叶面积指数（Leaf Area Index，LAI）、叶片中叶绿素含量（Chlorophyll Content，Cab［28］）、光合有效辐射分量（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）大豆种植情况识别的精度评价。由于大豆种植没有用监督分类，而是用非监督以及决策树分类，最后用种实测地块是否落入识别范围作为精度评价，评价公式如下。
P = 正确识别地块实测地块×100%正确识别地块实测地块×100%
（2）其中，P为精度。
（3）产量估算结果的精度评价。
用实测产量的1/5作为精度评价样本，与实测样本进行散点图分析，平均估算精度。
3.7 区域收入保险应用
得到地块单产之后，根据保险公司划定的区域，比如乡镇、村庄、特殊经验区域等进行平均产量统计。保险公司依据平均产量，再结合区域价格形成一个实际的区域收入，与保单约定的单产保障收入进行比较来决定赔付与否。
4 研究结果及分析
4.1 大豆种植情况及精度
大豆分布和受害情况遥感识别结果如图8所示。地块面积统计结果显示，嘉祥县2018年共种植大豆种植面积124 km2，与当地政府部门上报的127 km2统计数据吻合较好。
图8 2018年嘉祥县大豆受害情况分布图
Fig. 8 Soybeans damage distribution map of Jiaxiang county in 2018
精度评价结果如下。
其中，85为长势正常大豆正确识别地块数量；88为长势正常大豆实测地块数量；21为受旱受害大豆正确识别地块数量；29为受旱受害大豆实测地块数量。总体上，长势正常的大豆，识别正确率很高，为96.59%，受旱受害的大豆由于表现不一，识别率约为72.4%。
由于总的识别面积为123.74 km2，而受旱受害的仅有31.47 km2，占总面积的25.4%，因此总体精度可以计算为：
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之间。
表2 NDVI、生理参数、气象参数与产量之间的Pearson相关系数表
Table 2 Pearson correlation coefficients between NDVI, biophisical parameters, meteorological parameter and yield
|
变量 |
与产量相关系数 |
变量 |
与产量相关系数 |
|
0823NDVI |
0.644** |
0922Cw |
0.356** |
|
0907LAI |
0.641** |
0803Cab |
0.351** |
|
0907Cab |
0.638** |
0922NDVI |
0.348** |
|
0907FaPar |
0.633** |
0808Cw |
0.333** |
|
0823FaPar |
0.628** |
0803FCover |
0.320** |
|
0907FCover |
0.618** |
0724FCover |
0.290* |
|
0823FCover |
0.615** |
0724FaPar |
0.289* |
|
0823Cab |
0.611** |
Max07FCover |
0.289* |
|
0907NDVI |
0.606** |
Max07FaPar |
0.288* |
|
0823LAI |
0.603** |
0808Cab |
0.281* |
|
0907Cw |
0.528** |
0808FCover |
0.274* |
|
0922FaPar |
0.486** |
0724LAI |
0.273* |
|
1002_0823 |
-0.482** |
0808FaPar |
0.268* |
|
0823_0624 |
0.480** |
0803FaPar |
0.268* |
|
0922FCover |
0.475** |
0803NDVI |
0.261* |
|
0922LAI |
0.449** |
TRMM2018091 |
0.257* |
|
0823Cw |
0.437** |
20180727_Day_LST |
0.255* |
|
0922Cab |
0.437** |
20180711_Day_LST |
0.237* |
|
0803LAI |
0.402** |
20180828_Night_LST |
0.235* |
|
0927NDVI |
0.393** |
Max07NDVI |
0.232* |
|
0808LAI |
0.376** |
注：** 显著水平0.01，*显著水平0.05
分析2018年8月23日时相与9月7日时相的NDVI与作物生理参数之间的自相关性，得到二者显著正相关，相关系数除了8月23日时相的水分含量（0823Cw）之外，都超过0.59（如表3）。因此在选用因子时，只选了0823NDVI和0823Cw。
表 3 8月23日与9月7日时相NDVI与作物生理参数的相关系数表
Table 3 Pearson correlation coefficients between NDVI and biophical parameters on Aug. 23 and Sep.7
|
因子 |
0823NDVI |
0907NDVI |
0907LAI |
0907FaPar |
0907FCover |
0907Cab |
0907Cw |
|
0823LAI |
0.816** |
0.735** |
0.812** |
0.802** |
0.780** |
0.802** |
0.561** |
|
0823FCover |
0.920** |
0.846** |
0.833** |
0.857** |
0.847** |
0.815** |
0.603** |
|
0823FaPar |
0.908** |
0.837** |
0.849** |
0.867** |
0.851** |
0.836** |
0.633** |
|
0823Cab |
0.817** |
0.734** |
0.834** |
0.819** |
0.793** |
0.833** |
0.623** |
|
0823Cw |
0.474** |
0.467** |
0.556** |
0.516** |
0.479** |
0.562** |
0.593** |
|
0823NDVI |
1.000 |
0.879** |
0.779** |
0.813** |
0.813** |
0.744** |
0.597** |
|
0907NDVI |
0.879** |
1.000 |
0.825** |
0.889** |
0.899** |
0.791** |
0.591** |
|
0907LAI |
0.779** |
0.825** |
1.000 |
0.982** |
0.968** |
0.989** |
0.837** |
|
0907FaPar |
0.813** |
0.889** |
0.982** |
1.000 |
0.995** |
0.960** |
0.797** |
|
0907FCover |
0.813** |
0.899** |
0.968** |
0.995** |
1.000 |
0.940** |
0.772** |
|
0907Cab |
0.744** |
0.791** |
0.989** |
0.960** |
0.940** |
1 |
0.828** |
|
0907Cw |
0.597** |
0.591** |
0.837** |
0.797** |
0.772** |
0.828** |
1 |
注：** 显著水平0.01
4.3 大豆产量模型
在相关系数最高的植被指数或者作物生理参数里选择3~5个参数，同时选择3~5个气象参数，进行逐步回归分析，得到两个与产量显著相关且相关系数都较高的输入因子，分别为0823NDVI和0907LAI，所构建的模型（公式6）对产量的解释度达到46.4%。其模型标准误差约为45 kg。再增加4个变量或者8个变量并未增加解释度，标准误也相差不到1 kg，模型系数摘要见表4。
表4 大豆多元线性估产模型摘要
Table 4 Summary of soybeans yield linear models
|
模型序号 |
模型 |
R2 |
标准误 |
F |
输入变量 |
标准化系数 |
变量显著性 |
|
1 |
Y=-282.356+547.583×0823NDVI （7） |
0.415** |
47.554 |
50.379 |
（常量） |
0.000 |
|
|
0823NDVI |
0.644 |
0.000 |
|||||
|
2 |
Y=-171.365+313.127×0823NDVI+35.154×0907LAI （8） |
0.464** |
45.824 |
30.359 |
（常量） |
0.025 |
|
|
0823NDVI |
0.386 |
0.010 |
|||||
|
0907LAI |
0.354 |
0.013 |
|||||
|
3 |
Y=-187.199+289.032×0823NDVI+30.816×0907LAI+957.102×0823Cw-12.239×NDVI（1002-0823） （9） |
0.471** |
46.210 |
15.136 |
（常量） |
0.051 |
|
|
0823NDVI |
0.340 |
0.049 |
|||||
|
0907LAI |
0.311 |
0.042 |
|||||
|
0823Cw |
0.093 |
0.388 |
|||||
|
NDVI(1002-0823） |
-0.025 |
0.845 |
|||||
|
4 |
Y=2025.482+323.644×0823NDVI+33.527×0907LAI+22.704×K_MOD11_10+8.307×K_MOD11_14-16.884×TRMM091+0.923×K_MOD11_23 （10） |
0.482** |
46.407 |
10.242 |
（常量） |
0.640 |
|
|
0823NDVI |
0.381 |
0.018 |
|||||
|
0907LAI |
0.338 |
0.022 |
|||||
|
K_MOD11_10 |
0.106 |
0.348 |
|||||
|
K_MOD11_14 |
-0.196 |
0.208 |
|||||
|
TRMM091 |
0.163 |
0.207 |
|||||
|
K_MOD11_23 |
0.013 |
0.902 |
|||||
|
5 |
Y=-285.649+509.074×0823NDVI+4.999×0907LAI+1349.340×0823Cw-95.854×NDVI（1002-0823）-309.911×NDVI（0823-0624）+231.163×0922FAPAR （11） |
0.520** |
44.684 |
11.912 |
（常量） |
0.005 |
|
|
0823NDVI |
0.599 |
0.007 |
|||||
|
0907LAI |
0.050 |
0.789 |
|||||
|
0823Cw |
0.131 |
0.219 |
|||||
|
NDVI(1002-0823) |
-0.196 |
0.176 |
|||||
|
NDVI(0823-0624） |
-0.428 |
0.025 |
|||||
|
0922FAPAR |
0.320 |
0.036 |
因此本研究最后采用公式8进行估产。即：
Y=-171.365+313.127×0823NDVI+35.154× 0907LAI （8）
其中，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/亩）相比，体现了受灾严重的农情。
图9 2018年嘉祥县大豆地块产量分级图
Fig.9 Graded map of soybean plot yield in Jiaxiang county in 2018
在大豆产量样本中，用剩余的1/5样本作为验证，得到如图10散点图，二者回归系数达0.92，表明产量估算模型的精度较高。
图10 实测产量与模拟产量验证散点图
Fig. 10 Scatter plot of measured yield and simulated yield
5 结论与讨论
5.1 结论
本研究主要基于哨兵2号卫星识别大豆种植分布及受害情况，计算植被指数和作物生理参数，结合实际测产、气象卫星遥感数据进行多参数回归的大豆产量建模。结果表明哨兵2号卫星数据对大豆地块的识别精度达90%，对产量的估算精度达92%。
但本研究对于农业保险迫切的行业应用需求而言仍存在一些缺陷。首先，花荚期、鼓粒期的影像数据是研究的基本保障，如果这个时期的遥感数据存在云及阴雨天气影响，对于估产结果精度和整个应用方案可能产生影响。其次，研究需要大量的实地样本，主要是产量样本的支持，且需要较多专家知识参与，并且采样方案需要综合考虑成本及效率。另外，保险理赔过程涉及当地农户对遥感的认知和接受程度，还涉及当地的民风和经济状况，还需要从多个角度，多种方法结合、多部门数据一起来定约定理赔。
5.2 讨论
根据以上研究可知，从技术角度看，遥感估产在农业保险中的应用整个研究过程有三个重要步骤，缺一不可。
一是野外采样。一般至少需要两次，第一次在作物生长季中期，主要采集作物的识别样本和长势、受灾受害情况。第二次是在作物收获期，与传统的采样一样，尽量根据作物地块识别结果，全局均匀采样，不同长势均匀覆盖。
二是目标作物遥感识别。对于目标作物识别而言，目标作物地块的形态、位置和边界很重要，与其他作物的时间区分很重要，因此目标作物识别之前的影像时间和波段选择、耕地掩膜和面向对象的图像分割很重要，可以较好地规避传统遥感影像分类的像素孤岛等错误。在目标作物识别过程中，对于类间差距，以及何种长势可能导致绝产的标注也很重要，对于杂粮作物以及其他油料、棉花等作物而言，在机理不明确的经验统计模型中，去除一些异常的长势，如叶子很多却不结荚的大豆，可以提高估算精度和解释度。
三是产量估算。在作物地块分好类之后，选取具有一定生理意义的作物光谱参数，加入气象参数，建立模型，并根据前期长势的标注进行产量的加权调整，使得建模过程更完整更可信。
但是这三个步骤也是区域收入保险运营需要考虑的成本投入。首先，采样步骤一般需要邀请与目标作物相关的农学专家参与，进行长势和产量的测定。其次，基于多光谱卫星遥感的作物识别，对作物关节生长期的卫星数据极其依赖，若赶上天气不好，就缺乏可用的多光谱卫星数据，需要启动其他数据预案，包括付费编程数据和无人机数据来完成目标作物的识别，此时可能会因无多时相遥感数据参与产量建模而废弃遥感估产方法，而只根据作物识别结果增加实地抽样测产来确定区域产量，这都极大增加保险运营成本。
农业保险承保范围广而细，全国所有省份几乎都参与了农业保险，保单细到每家每户，这给农业保险定损核损带来了极大困难，这就是农业保险由传统的个体保险发展到区域保险的主要原因。但即便是区域类型保险，理赔的区域也需要精细到乡镇，甚至村。由此可见，农业区域产量或者收入保险急需借助遥感这种大范围全天候客观可靠的技术来完成产量实时估算的任务。遥感技术作为农业保险急需的技术，是农业保险行业发展的新动力。尽管很难，但是遥感用于农业保险是趋势，农业保险对遥感的需求有增不减。遥感应用于农业保险还需要更多的投入应用、经验总结与凝华。
作者：陈爱莲、李家裕、张圣军、朱玉霞、 赵思健、孙伟、张峭
来源：智慧农业
（内容版权归原作者所有，如有侵权请立即与我们联系，我们将及时处理）
