Welcome to Smart Agriculture 中文
Information Processing and Decision Making

Application of Satellite Remote Sensing Yield Estimation Technology in Regional Revenue Protection Crop Insurance: A Case of Soybean

  • CHEN Ailian , 1, 2 ,
  • LI Jiayu 3 ,
  • ZHANG Shengjun 3 ,
  • ZHU Yuxia , 1, 2 ,
  • ZHAO Sijian 1, 2 ,
  • SUN Wei 1, 2 ,
  • ZHANG Qiao 1, 2
Expand
  • 1. Agricultural Information Institute, Chinese Academy of Agricultural Sciences, Beijing 100081, China
  • 2. Key Laboratory of Agricultural Information Service Technology, Ministry of Agriculture and Rural Agriculture, Beijing 100081, China
  • 3. China Pacific Property Insurance Company Limited Shandong Branch, Jinan 250001, China

Received date: 2020-06-05

  Revised date: 2020-09-18

  Online published: 2020-11-18

Highlights

In recent years, revenue protection crop insurance is an innovative insurance that has been prioritized in China. But it still lacks the support of the third-party yield data around crop harvest time. Aiming to provide objective yield data for revenue protection crop insurance, satellite remote sensing production estimation technology was employed to discuss its application mode and applicability. Taking the soybean revenue protection insurance in Jiaxiang county, Shandong province as an example, we first extracted soybean planting plots, calculated vegetation index and crop physiological parameters based on Sentinel-2 satellite images in 2018 . Combining to TRMM precipitation data from TRMM precipitation-monitoring radar satellite and MODIS land surface temperature data from Terra/Aqua satellite and site yield data, we established a multi-parameter linear regression model, and estimated soybean yield per unit area. The crop extraction results showed that the soybean planting area in the study area was 1.24 km2, which was in good agreement with the 1.27 km2 reported by the local agricultural bureau; and with using the actual measurement plots, the remote sensing identification accuracy of the planting distribution plots reached 90%. The yield estimation results showed that the NDVI of the soybean pod stage on August 23 and the leaf area index of the soybean seedling stage on September 7 explained the soybean yield per hectare the best, and the average estimated yield of the whole area was 244,500 kg/m2, which reflects the severely affected agricultural conditions, comparing to 299,800 kg/km2 in previous years.The regression coefficient between the estimated yield data and the measured data reached 0.92, which meet the application needs.With this results, the estimated yield of different towns can be summarized, and the regional yield was present, and was used as the real yield in 2018, multiplying with the average soybean price around October 11 to December 10 from the local price bureau, the real revenue was obtained. Compared the real revenue to the expected revenue in the contract of insurance, the claims work was decided. The results indicated that the Sentinel-2 satellite data could be used to identify the soybean planting distribution in the study area accurately, and to complete the yield estimation as soon as one week after the soybean harvest, which could guide the insurance company's claims work. The whole methodology is capable of aiding the claims work in revenue protection crop insurance.

Cite this article

CHEN Ailian , LI Jiayu , ZHANG Shengjun , ZHU Yuxia , ZHAO Sijian , SUN Wei , ZHANG Qiao . Application of Satellite Remote Sensing Yield Estimation Technology in Regional Revenue Protection Crop Insurance: A Case of Soybean[J]. Smart Agriculture, 2020 , 2(3) : 139 -152 . DOI: 10.12133/j.smartag.2020.2.3.202006-SA002

1 引 言

农业保险是农业风险管理的重要工具,能够减少农民的直接经济损失、减轻政府和社会各方压力、提高农业综合实力,保障农业和农村社会持续稳定发展,受到各国政府的高度重视。中国自2007年开始实施保费补贴型政策性农业保险,主要险种是政策性直接物化成本保险,主保个体产量。目前中国农业保险保费已居世界第二,仅次于美国。然而,随着农业保险的发展和土地政策的改革,政策性直接物化成本保险因存在保额低、验标验险耗财耗力、道德风险高等缺陷,与保险公司和农业经营户,尤其是新型经营主体抵抗生产风险的需求存在较大差距。随后中国陆续进行试点试验和增加了指数保险、价格保险、收入保险、区域产量、区域收入保险等不同险种1-3。其中,区域收入保险是中国近年重点发展的创新型险种,它综合了区域平均产量与农产品价格,形成了一个对区域平均收入进行保护的险种,相对降低了验标验险的成本,还可规避基差风险和道德风险,在中国山东、山西、辽宁等地已经开始广泛试点,并得到政府部门和保险公司的大力推广4-7
开展区域收入保险需要两个重要的数据支持,一是收获期的区域平均价格,二是收获期的区域平均产量。区域的单位一般为县、乡镇、甚至村级行政单元。其中,区域价格可以参考物价局的实时价格数据,或者定点价格监测数据。区域产量需要第三方公布的当年产量数据,但中国统计局官方发布的产量数据一般会滞后1~3年。若按照传统的保险运营模式,依然存在信息不对称和协议理赔等问题,因而需要更客观及时的估产方法,这就需要遥感技术的介入。
近年来,随着农业保险对信息技术的需求越来越高,遥感信息技术以其客观特性逐渐被应用到农业保险的验标和定损等环节8,9。遥感技术在验标环节回答田间为何种农作物、种植在哪里、是否与承保面积相符等问题其实是遥感作物识别的问题,已经比较成熟8。遥感技术在理赔环节也有不少辅助理赔的成功案例10,11。遥感还能为创新型指数保险险种提供客观指数3,在农险中的应用益处已得到保险各方的认可8-11。然而将遥感产量估算结果直接用到农业保险理赔的案例却不多见。
从遥感产量估算的技术可行性角度看,作物产量估算的遥感估算在小尺度(县域尺度)的研究文献已经有很多,包括小麦12-15、玉米16和水稻卫星遥感估产17等。相比粮食作物而言,大豆等油料作物的产量估算相对复杂,因此研究也更精细。Yu等18及Maimaitijiang等19基于无人机高光谱数据,Gusso等20基于MODIS的大尺度数据,以及Kross等21基于5 m分辨率的Rapideye数据开展了不同尺度的大豆地上生物量或产量的估算研究。国内大豆的产量估算研究也比较成熟。南京农业大学盖钧镒院士及其团队依托大豆育种研究,开展了许多与产量估算相关的研究,涉及了低空、地面和高空遥感。其中,开展基于地面光谱或无人机遥感的作物生理参数,如叶面积指数、光谱参数及其与生物量的关系研究的有齐波等22、陆国政等23;开展基于地面光谱或无人机遥感的大豆籽粒产量研究的有吴琼等24、张宁25和赵晓庆等26。国内其他作者也开展部分基于卫星遥感的大豆叶面积指数或产量的大豆估产研究,如张淮栋等27基于高分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)所示。
预期收入=(预期价格+现货基差补偿价 格)×预期单产
预期价格使用大连商品交易所保险往年11月和1月到期的黄豆1号合约在6月保险销售前五个交易日的日收盘价均值。现货基差补偿价格即约定为0.5 CNY/kg。预期单产为当地前五年单产数据中(去掉最大与最小值后)的三年平均值,结合实地调研情况,约定预期单产。保障水平是对预期收入的保障,可以是0%~100%,考虑到可能存在的道德风险和逆选择,区域收入保险的保障水平理论不得超过90%。
当实际区域收入低于保障收入时触发理赔(公式(2))。
实际区域收入=实际区域价格×遥感区域产量
理赔时,实际区域价格采用当地物价局的一段时间的大豆平均销售价格。研究区当地大豆销售期始于每年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,Cab28)、光合有效辐射分量(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 相关关系分析,再根据相关系数分析结果选择系数高的因子参与逐步回归,并人工添加气象因子加入回归模型,根据决定系数R 2的变化,确定估产模型进行估产。
在产量的99个点中,除去有云的样点,余下88个样点,进行相关性分析。由于样点不多,用其中4/5作为建模样点,余下1/5作为验证样点。建模样点与验证样的选取步骤如下:用求余函数将每个样点在ArcGIS的数据表要素编号FID,除以5,将余数为0的点提取出来,即得到1/5的验证样点17个,其余4/5的71个样点参与回归建模。

3.6 精度评价

(1)大豆种植情况识别的精度评价。由于大豆种植没有用监督分类,而是用非监督以及决策树分类,最后用种实测地块是否落入识别范围作为精度评价,评价公式如下。
P = × 100 %
其中,P为精度。
(2)产量估算结果的精度评价。
用实测产量的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

精度评价结果如下。
P长势正常大豆 = 85 88 × 100 %  =96.59%
P受旱受害大豆 = 21 29 × 100 %  =72.4%
其中,85为长势正常大豆正确识别地块数量;88为长势正常大豆实测地块数量;21为受旱受害大豆正确识别地块数量;29为受旱受害大豆实测地块数量。总体上,长势正常的大豆,识别正确率很高,为96.59%,受旱受害的大豆由于表现不一,识别率约为72.4%。
由于总的识别面积为123.74 km2,而受旱受害的仅有31.47 km2,占总面积的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之间。
表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)。因此在选用因子时,只选了0823NDVI0823Cw
表 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个气象参数,进行逐步回归分析,得到两个与产量显著相关且相关系数都较高的输入因子,分别为0823NDVI0907LAI,所构建的模型(公式6)对产量的解释度达到46.4%。其模型标准误差约为45 kg。再增加4个变量或者8个变量并未增加解释度,标准误也相差不到1 kg,模型系数摘要见表4
表4 大豆多元线性估产模型摘要

Table 4 Summary of soybeans yield linear models

模型序号 模型 R 2 标准误 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
其中,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 讨论

根据以上研究可知,从技术角度看,遥感估产在农业保险中的应用整个研究过程有三个重要步骤,缺一不可。
一是野外采样。一般至少需要两次,第一次在作物生长季中期,主要采集作物的识别样本和长势、受灾受害情况。第二次是在作物收获期,与传统的采样一样,尽量根据作物地块识别结果,全局均匀采样,不同长势均匀覆盖。
二是目标作物遥感识别。对于目标作物识别而言,目标作物地块的形态、位置和边界很重要,与其他作物的时间区分很重要,因此目标作物识别之前的影像时间和波段选择、耕地掩膜和面向对象的图像分割很重要,可以较好地规避传统遥感影像分类的像素孤岛等错误。在目标作物识别过程中,对于类间差距,以及何种长势可能导致绝产的标注也很重要,对于杂粮作物以及其他油料、棉花等作物而言,在机理不明确的经验统计模型中,去除一些异常的长势,如叶子很多却不结荚的大豆,可以提高估算精度和解释度。
三是产量估算。在作物地块分好类之后,选取具有一定生理意义的作物光谱参数,加入气象参数,建立模型,并根据前期长势的标注进行产量的加权调整,使得建模过程更完整更可信。
但是这三个步骤也是区域收入保险运营需要考虑的成本投入。首先,采样步骤一般需要邀请与目标作物相关的农学专家参与,进行长势和产量的测定。其次,基于多光谱卫星遥感的作物识别,对作物关节生长期的卫星数据极其依赖,若赶上天气不好,就缺乏可用的多光谱卫星数据,需要启动其他数据预案,包括付费编程数据和无人机数据来完成目标作物的识别,此时可能会因无多时相遥感数据参与产量建模而废弃遥感估产方法,而只根据作物识别结果增加实地抽样测产来确定区域产量,这都极大增加保险运营成本。
农业保险承保范围广而细,全国所有省份几乎都参与了农业保险,保单细到每家每户,这给农业保险定损核损带来了极大困难,这就是农业保险由传统的个体保险发展到区域保险的主要原因。但即便是区域类型保险,理赔的区域也需要精细到乡镇,甚至村。由此可见,农业区域产量或者收入保险急需借助遥感这种大范围全天候客观可靠的技术来完成产量实时估算的任务。遥感技术作为农业保险急需的技术,是农业保险行业发展的新动力29,30。尽管很难,但是遥感用于农业保险是趋势,农业保险对遥感的需求有增不减。遥感应用于农业保险还需要更多的投入应用、经验总结与凝华。
1
艾思言. 定西市马铃薯区域产量保险的适用性研究[J]. 甘肃农业, 2018(15): 21-26.

AI S. Study on the applicability of potato yield insurance in Dingxi city[J]. Gansu Agriculture, 2018(15): 21-26.

2
李艳, 陈盛伟. 农作物区域产量保险研究进展[J]. 山东农业科学, 2018, 50(7): 161-166.

LI Y, CHEN S. Research review of crop area yield insurance[J]. Shandong Agricultural Sciences, 2018, 50(7): 161-166.

3
张峭, 李越, 郑茗曦. 农业指数保险的发展、应用与建议[J]. 农村金融研究, 2018(6): 14-20.

ZHANG Q, LI Y, ZHENG M. Development, application and recommendations of agricultural index insurance[J]. Rural Finance Research, 2018(6): 14-20.

4
汪必旺, 张峭. 美国农作物收入保险规模演变及展望[J]. 农业展望, 2018, 14(6): 34-40.

WANG B, ZHANG Q. Scale evolution of the U.S. crop revenue insurance and its prospect[J]. Agricultural Outlook, 2018, 14(6): 34-40.

5
庹国柱, 张峭. 论我国农业保险的政策目标[J]. 保险研究, 2018(7): 7-15.

TUO G, ZHANG Q. On the policy objectives of agricultural insurance in China[J]. Agricultural Research, 2018(7): 7-15.

6
张峭. 农产品价格"保险+期货"的实践探索[N]. 粮油市场报, 2017-07-29 (B 03).

ZHANG Q. Practical exploration of agricultural product prices "insurance + futures" [N]. Grain and Oil Market News, 2017-07-29 (B 03).

7
李越, 张峭, 王克. 中国农业保险保障水平现状、问题与建议(下)[N]. 中国保险报, 2016-11-27(004).

LI Y, ZHANG Q, WANG K. The current situation, problems and suggestions of China's agricultural insurance protection level (Part 2)[N]. China Insurance News, 2016-11-27(004).

8
刘振功. 基于遥感技术的农业保险业务模式创新研究[D]. 济南: 山东大学, 2016.

LIU Z. Research on the innovation of agricultural insurance business model based on remote sensing technology [D]. Jinan: Shandong University, 2016.

9
葛瑞华. 农业保险中信息技术的应用现状与问题研究[D]. 成都: 成都理工大学, 2014.

GE R. Research on the application status and problems of information technology in agricultural insurance[D]. Chengdu: Chengdu University of Technology, 2014

10
蒙继华, 付伟, 徐晋, 等. 遥感在种植业保险估损中的应用[J]. 遥感技术与应用, 2017(2): 238-246.

Meng J, Fu W, Xu J, et al. Application of remote sensing in planting industry insurance loss estimation[J]. Remote Sensing Technology and Application, 2017(2): 238-246.

11
祁鑫. 遥感技术应用于农业保险业务模式[J]. 创新农技服务, 2017, 34(14): 171.

QI X. Remote sensing technology applied to agricultural insurance business model[J]. Innovative Agricultural Technology Service, 2017, 34(14): 171.

12
HUANG J, MA H, SEDANO F, et al. Evaluation of regional estimates of winter wheat yield by assimilating three remotely sensed reflectance datasets into the coupled WOFOST-PROSAIL model[J]. European Journal of Agronomy, 2019, 102: 1-13.

13
FRANCH B, VERMOTE E F, SKAKUN S, et al. Remote sensing based yield monitoring: Application to winter wheat in United States and Ukraine[J]. International Journal of Applied Earth Observation and Geoinformation, 2019, 76: 112-127.

14
JIN X, KUMAR L, LI Z, et al. A review of data assimilation of remote sensing and crop models[J]. European Journal of Agronomy, 2018, 92: 141-152.

15
SCHWALBERT R A, AMADO T J C, NIETO L, et al. Forecasting maize yield at field scale based on high-resolution satellite imagery[J]. Biosystems Engineering, 2018, 171: 179-192.

16
苏伟, 侯宁, 李琪, 等. 基于Sentinel-2遥感影像的玉米冠层叶面积指数反演[J]. 农业机械学报, 2018, 49(1): 151-156.

SU W, HOU N, LI Q, et al. Retrieving leaf area index of corn canopy based on Sentinel-2 remote sensing image[J]. Transactions of the CSAM, 2018, 49(1): 151-156.

17
彭代亮, 黄敬峰, 王福民, 等. 基于分区的水稻总产遥感估算研究[J]. 中国科学:信息科学, 2010, 40(S1): 184-195.

PENG D, HUANG J, WANG F, et al. Research on remote sensing estimation of total rice yield based on subregions[J]. Science China: Information Science, 2010, 40(S1): 184-195.

18
YU N, LI L, SCHMITZ N, et al. Development of methods to improve soybean yield estimation and predict plant maturity with an unmanned aerial vehicle based platform[J]. Remote Sensing of Environment, 2016, 187: 91-101.

19
MAIMAITIJIANG, MAITINIYAZI, SAGANVASIT, SIDIKEPAHEDING, et al. Vegetation Index Weighted Canopy Volume Model (CVMVI) for soybean biomass estimation from unmanned aerial system-based RGB imagery[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2019, 151: 27-41.

20
GUSSO A, DUCATI J R, VERONEZ M R, et al. Spectral model for soybean yield estimate using MODIS/EVI data[J]. International Journal of Geosciences, 2013, 4: 1233-1241.

21
KROSS A, MCNAIRN H, LAPEN D, et al. Assessment of RapidEye vegetation indices for estimation of leaf area index and biomass in corn and soybean crops[J]. International Journal of Applied Earth Observation and Geoinformation, 2015, 34: 235-248.

22
齐波, 张宁, 赵团结, 等. 利用高光谱技术估测大豆育种材料的叶面积指数[J]. 作物学报, 2015, 41(7): 1073-1085.

QI B, ZHANG N, ZHAO T, et al. Prediction of leaf area index using hyperspectral remote sensing in breeding programs of soybean[J]. Acta Agronomica Sinica, 2015, 41(7): 1073-1085.

23
陆国政, 杨贵军, 赵晓庆, 等. 基于多载荷无人机遥感的大豆地上鲜生物量反演[J]. 大豆科学, 2017, 36(1): 41-50.

LU G, YANG G, ZHAO X, et al. Inversion of soybean fresh biomass based on multi-payload unmanned aerial vehicles (UAVs)[J]. Soybean Science, 2017, 36(1): 41-50.

24
吴琼, 齐波, 赵团结, 等. 高光谱遥感估测大豆冠层生长和籽粒产量的探讨[J]. 作物学报, 2013, 39(2): 309-318.

WU Q, QI B, ZHAO T, et al. A tentative study on utilization of canopy hyperspectral reflectance to estimate canopy growth and seed yield in soybean[J]. Acta Agronomica Sinica, 2013, 39(2): 309-318.

25
张宁. 基于GreenSeeker主动遥感的大豆NDVI/RVI分析与籽粒产量估测的初步研究[D]. 南京: 南京农业大学, 2013.

ZHANG N, A preliminary study on soybean NDVI/RVI analysis and grain yield estimation based on GreenSeeker active remote sensing [D]. Nanjing: Nanjing Agricultural University, 2013.

26
赵晓庆, 杨贵军, 刘建刚, 等. 基于无人机载高光谱空间尺度优化的大豆育种产量估算[J]. 农业工程学报, 2017, 33(1): 110-116.

ZHAO X, YANG G, LIU J, et al. Estimation of soybean breeding yield based on optimization of spatial scale of UAV hyperspectral image[J]. Transactions of the CSAE, 2017, 33(1): 110-116.

27
张淮栋, 陈争光, 张成龙, 基于高分二号-NDVI的大豆遥感估产的时相选择[J]. 湖北农业科学, 2018, 57(6): 103-108.

ZHANG H, CHEN Z, ZHANG C, Time selection of soybean yield estimation based on Gaofen- 2 -NDVI[J]. Hubei Agricultural Sciences, 2018, 57(6): 103-108.

28
JACQUEMOUD S, VERHOEF W, BARET F, et al. PROSPECT+SAIL models: A review of use for vegetation characterization[J]. Remote Sensing of Environment, 2009, 113: 56-66.

29
赵思健, 张峭, 陈敬敏. 大数据视角下的农业保险创新与提升[J]. 保险理论与实践, 2017(12): 22-42.

ZHAO S, ZHANG Q, CHEN J. Innovation and improvement of agricultural insurance from the perspective of big data[J]. Insurance Theory and Practice, 2017(12): 22-42

30
赵思健, 张峭. 农业保险科技的现状、体系与展望[J]. 中国保险, 2018(7): 28-35.

ZHAO S, ZHANG Q. The current situation, system and prospect of agricultural insurance technology[J]. China Insurance, 2018(7): 28-35.

Outlines

/