基于地级市的区域水稻遥感估产与空间化研究

汤斌　王福民　周柳萍等

摘要：利用遥感技术，可以准确、及时地掌握中国主要水稻生长区的水稻生长状况与产量信息，能促进中国农业政策合理制定、粮食价格的理性宏观调控以及国际国内粮食良性贸易。本试验以江苏省为研究区域，利用MOD09A1数据提取水稻面积，针对行政区域为单位的水稻进行产量估算，全省估产结果精确度均值达到99.46%，各市的精确度基本维持在95%左右。在估产结果的基础上，依据最大植被指数的空间分布，将以行政区域为单位的单产值转化到以单个像元为单元的单产值。这不仅细化了区域水稻估产结果，为更直观地了解水稻田产量信息提供了可能，而且可以用于指导水稻品种改良和栽培技术改进等工作。

关键词：MOD09A1；水稻；面积提取；估产；空间化；江苏省

中图分类号： S127文献标志码： A文章编号：1002-1302（2015）11-0525-04

收稿日期：2014-11-29

项目基金：国家自然科学基金（编号：41371393、51109183）；博士点基金（编号：20110101120036）。

作者简介：汤斌（1989—），男，安徽铜陵人，硕士研究生，主要从事农业遥感方面的研究。E-mail：tltb09125006@163.com。民以食为天，我国是人口大国，粮食供给的稳定性与社会安全息息相关，历史上的社会动乱也都与食物短缺有着密切联系，所以农业是具有重要战略性意义的基础行业。而水稻是我国主要粮食作物之一，全国水稻产量约占粮食作物总产量的40%，稻田面积约占全国粮食作物总面积的28%[1]。这也就有了水稻产量的统计与估算的需求。在当今复杂多变的经济和气候条件下准确、及时地掌握中国及全球重点产粮区的水稻生长状况与产量信息，对于中国农业政策的制定、粮食价格的宏观调控以及国际国内粮食贸易具有重要意义[2]。

传统的水稻产量统计方式，主要有人工监测，通过人的实地调查、逐级汇总来了解水稻种植面积和产量的变化，显然，由于数据的搜集、传输和加工都是手工方式，调差结果大大滞后于农业资源的变化，还容易掺入人为因素而扭曲资源的本来面目，因而不客观[3]。这也就凸显出遥感技术的宏观优势、客观优势、重复优势以及廉价优势[4]。

目前遥感估产方法主要包括经验统计方法、光能利用率方法和机理模型方法，这些方法都有自己的优点以及局限性。光能利用率方法主要是基于资源平衡的观点，假定所有资源对植物生长有平等的限制作用，以植物作用过程和Monteith提出的光能利用率为基础建立的[5-6]。遥感数据较容易获得，也更容易掌握净初级生产力（net primary productivity，NPP）的季节、年际变化；但是缺乏可靠的生理生态基础，光能传递及转换过程中还是存在不少不确定性。机理模型方法是以相似性原理为依据，以分析作物生长发育的物理过程、物理机制和环境条件为手段，设法将作物生长发育、产量形成的规律表述为有关的物理学定律，并用数学语言将这些有关的物理学定律写成数学模型，在一定假设条件下，确定边界条件、简化模型，寻求合适的数学解法，通过模拟试验调控模型的参数，最后建立作物估产的数值模拟模型[7]。机理模型方法虽然适用范围较广，但是参数校正工作却什么繁琐，而且在某些情况下某些参数没办法进行校正。常规的经验统计方法是利用前人的研究，采用一些经验的水稻面积提取公式提取水稻面积并估算作物产量。这种方法操作简单，估产精确度比较高，但所建立的估产模型一般只适用于建模的区域，不适宜推广。

常规的经验统计方法是以行政单元为单位，无法体现水稻单产空间变化的信息。因为水稻田一般是连片的、跨区域的，这样单纯的以行政区划分界线作为边界进行水稻估产结果缺乏水稻单产的空间性，不能直观地了解到不同区域的水稻生产力。为此本研究在传统经验统计方法的基础上，采用一种最大植被指数的方法，将以行政单元为单位的产量转换为以像元为单位的信息，得到每个水稻像元级别的水稻产量图，从而不仅提供了更加丰富的水稻产量信息，而且可以指导水稻品种改良和栽培技术改进等工作。

1研究区域概况

研究区域为江苏省。江苏省地处我国东部沿海，介于东经116°18′～121°57′、北纬30°45′～35°20′之间，地居长江、淮河下游，东临黄海，西连安徽，北接山东，南与上海、浙江接壤，全省面积10.26万 km2，包括13个地级市，主要分成苏北徐州、盐城、连云港、淮安、宿迁5市，苏中扬州、泰州、南通3市，苏南南京、苏州、无锡、常州、镇江5市。江苏是全国地势最低的一个省，绝大部分地区海拔都在50 m以下，全省最高峰云台山玉女峰海拔只有625 m，低山丘陵主要集中在北部和西南部，仅占全省总面积的14.3%。全省覆盖h27v05和h28v05两块。

江苏属亚热带和暖温带地区，气候温和，雨量适中，具有寒暑变化显著、四季分明的特征。水稻是江苏省主要种植作物。近30年来，我国沿海的浙江、福建、广东等地水稻总产都出现了不同幅度的下降，而江苏省却出现了大幅上升。江苏省三大区域的水稻种植面积变化差异明显，苏南地区水稻种植面积大幅度减少，苏中地区水稻种植面积稳中略减，苏北地区水稻种植面积大幅度增加，在水稻种植面积上呈现出明显的“南减北增”的趋势[8]。水稻生长期为5—10月，中稻在5月初播种，6月初移栽，7月初拔节，8月初抽穗，9月末黄熟，10月初收割。单季晚稻在5月中旬播种，6月中旬移栽，7月中下旬拔节，8月末抽穗，10月中旬黄熟，10月末收割[9]。

2数据获取与研究方法

2.1数据获取

2.1.1遥感数据江苏省2004—2007年MOD09A1数据，即MOD09A1.h27v05和MOD09A1.h28v05。下载地址为：http：//ladsweb.nascom.nasa.gov/data/search.html。空间分辨率为500 m，时间分辨率为8 d。MOD09数据的用户说明MOD09_UsersGuide，下载地址为：http：//www.gscloud.cn/userfiles/file/MOD09_UserGuide.pdf。

2.1.2农业数据江苏省2004—2007年各地级市稻谷种植面积和单产，在2005—2008年《江苏省农村统计年鉴》上查询得到。

2.1.3江苏省行政区划图江苏省1 ∶25万的市级行政区划数据，用于提取各行政单元遥感数据、面积提取及估产研究等。

2.2研究方法

2.2.1用于估产的水稻面积提取方法本研究根据肖向明等对东南亚和中国南部的研究得到的关于地表水分指数（land surface water index，LSWI）和归一化插值植被指数（normalized difference vegetation index，NDVI）、增强型植被指数（enhanced vegetation index，EVI）的关系[10]，再以孙华生等对于全国水稻种植面积提取的研究结果为基础[11]，提取江苏省水稻种植区域。因为EVI耦合了抗大气植被指数ARVI和土壤调节植被指数SAVI，对气溶胶和土壤背景影响做了进一步校正[12]，所以本研究选择EVI作为水稻面积提取的植被指数。研究采用时间分辨率为8 d的MOD09A1数据进行植被指数计算和水稻面积提取，用MOD09A1自带的QA云检测波段进行去云处理，得到LSWI和去云过后的EVI，进而提取出水稻区域。

（1）EVI和LSWI计算：水稻面积提取以及估产和最终单产的像元化的植被指数和水指数分别是采用EVI和LSWI。MOD09A1包含7个波段，红光波段、近红外波段、蓝光波段、短波红外波段分别为MOD09A1数据的波段1、2、3、6。具体的计算公式[13-15]如下：

式中：RED、NIR、BLUE、SWIR分别是红光波段、近红光波段、蓝光波段和短波红外波段的地表反射率。通过计算可以得到生育期内的EVI值和LSWI值。

（2）QA波段云检测与条件时间序列插值去云处理：MOD09A1数据产品包含像元尺度的质量评估信息（quality assessment，QA），由于便于节省归档空间，加快网络传输效率，QA信息被压缩成二进制[16]。所以通过MOD09_UsersGuide中提供的sur_refl_state_500 m关于500 m分辨率的云与阴影的信息描述建立云掩膜。

虽然2 G数据MOD09已经经过严格的图像处理，在一定程度上去除了云、云阴影、气溶胶对图像质量的影响，但是在江苏省水稻生育期的很多图景中还是能很明显地看到有云层覆盖；所以在一定程度上做进一步的去云处理是十分有必要的。方法是在建立每个时期的云掩膜之后，用条件时间序列插值算法（Conditional Temporal Interpolation Filtering，CTIF）[17]对前期计算所得的EVI值进行去云处理。具体操作步骤为：通过云掩膜确定每个时期研究区域内是否受到云影响，总共分为4种情况：（1）如果未受到影响就保留原本的EVI值；（2）若有云层覆盖，且前1个和后1个时期图像相应像元未受到云影响，则取前、后2个时期像元EVI的平均值代替受到原始受到云影响的EVI。（3）若有云层覆盖，且前、后2个时期有且仅有1个图像相应像元受到云影响，则取前或后未受影响的像元EVI值代替原始像元EVI。（4）若有云层覆盖，且前、后2个时期相应像元都受到云影响，则原始像元视为无效。

（3）水稻种植区域的识别算法：由于地表有众多的建筑物、水体、植被、山体、土壤、水稻田以及其他经济作物，但是水稻田有其鲜明的特点，其在生育期前期会被水体覆盖，而在水稻移栽期之后，水稻开始长出水面，并在抽穗期生长达到最旺盛，水稻几乎覆盖了整个地表。这就为提取水稻种植区域提供了可能，在生育期前期可以与植被区分开来，在生育期后期又可以和永久水体加以区分。

由于研究区域江苏省只种植单季稻，所以本研究只考虑单季稻情况。具体的水稻种植区域识别条件有：（1）根据肖向明等对东南亚水稻种植区域提取的研究以及孙华生等对于全国水稻种植区域提取的研究，并结合江苏省的状况最终确定单季稻的提取算法为LSWI>0.12，EVI<0.26，（LSWI+0.065）>EVI。（2）由于在水稻移栽期之后大约第6～11个8 d 水稻生育期将进入抽穗期，也就是水稻生长最旺盛的时期，水稻几乎覆盖了地表。此时这6个8 d的同一像元EVI平均值需要大于0.35，这进一步去除了永久水体的干扰[18]。（3）由于江苏省地形地势起伏不大，故本研究不考虑数字高程模型（Digital Elevation Model，DEM）；（4）经过CTIF处理后，仍受到云干扰的，被视为无效的像元点，不参与以上操作。

根据《江苏省农村统计年鉴》，江苏省全省水稻的移栽期基本介于6月初到6月中旬，所以定为全年的第145、153、161、169天的4个8 d，这在一定程度上完整覆盖了全省的水稻移栽期。最终通过以上识别条件，并结合这4个8 d的图景得到最终的水稻种植区域。

2.2.2水稻估产方法在分别对2004—2007年4年水稻种植区域进行提取之后，就进入了估产的阶段。本研究是以2004—2006年3年水稻移栽期后生育期作为研究时段，即全年第161～305天，以CTIF法去云处理后的EVI作为研究的自变量，采用多元线性回归中的逐步回归作为估产方法，对2007年的江苏省各地级市的水稻产量进行估产。具体的操作步骤：（1）在提取出水稻种植区域以后，建立水稻区域掩膜，与全年第161～305天经过CTIF去云处理后的EVI图景进行提取，得到水稻有效区域内的EVI值；（2）在ArcGIS中用市级行政区划对水稻有效区域进行裁剪，将无数据的像元点设为no data，再计算各个地级市区域内的EVI均值；（3）将计算得到的2004—2006年3年的全年第161～305天江苏省各个地级市的EVI均值与其相对应的实际水稻单产统计值带入预测分析软件（Predictive Analytics Software，PASW），以EVI均值为自变量，实际水稻单产为因变量，建立多元线性回归模型，采用逐步回归，最终建立回归模型，并选择最佳模型；（4）在选取最佳估产模型之后，对2007年进行估产，并统计江苏省各个地级市单产的估计误差。

2.2.3水稻单位面积产量空间化由于原来的水稻估产基本上都是以行政区划分界线作为边界，估产结果是以行政单位作为基础，得到各个行政区域内的水稻单产值。即以一个水稻单产值作为整个行政区的水稻单产，从而忽略了区域内水稻单产的空间变化，难以获得区域内水稻单产的空间差异，为此本研究采用一种最大植被指数方法将行政区域内的单产转换为区域单产，计算公式[19]如下：

Yieldp=Yieldavg×EVIpEVIavg。 （3）

式中：Yieldp、Yieldavg分别为单像元点的水稻单产和地级市水稻单产均值，EVIp、EVIavg分别为单像元点的EVI值和各地级市的EVI均值。

具体的操作步骤：（1）江苏省的水稻抽穗期集中在8月初到8月末，主要是全年的第217、225、233、241天之间的4个8 d，对这4个8 d做最大值合成处理（Maximum Value Composition，MVC），得到EVIp；（2）再导入ArcGIS中，计算各个地级市水稻区域的平均EVI值，即EVIavg；（3）通过公式（3）计算得到最终基于像元级的水稻单产空间分布图。

3结果与分析

3.1水稻面积提取

经过植被指数EVI和水指数LSWI的计算，QA波段云检测，CTIF时间序列去云处理过后，得到2004—2007年的水稻区域面积。提取的江苏省水稻田主要分布在苏中、苏北地区，这与实际情况相符。这4年提取的全省水稻像元点分别为102 994、81 389、96 812、75 236个，相对误差分别为18.5%、-11.8%、2.4%、-15.58%。针对估产的稻田面积提取，更需要的是大块稻田被提取出来，2004—2006年的稻田分布都较均匀，符合实际，但是，2007年苏南地区水稻田提取面积过少，不能达到要求。

改革开放以后，苏南的南京、无锡、苏州、常州和镇江都快速进入城市化，城市建设和经济发展都在江苏省名列前茅，这使得苏南地区水稻田面积大幅减少。不过在改革开放的进程中，江苏省的水稻总产量却有着大幅增加，除了与水稻田生产力增加有关外，还与江苏省大力扶持苏北、苏中的水稻种植有关。“国家粮食丰产科技工程”从2004年开始启动江苏省水稻项目，其中80%的扶持计划都是针对苏北、苏中地区；但是，2007年苏南地区水稻田面积减幅较大的原因，还有可能是CTIF去云处理后，去除了很多像元，导致提取的稻田面积大幅减少。针对这一情况，保留2007年苏北苏中稻田区域，用2006年苏南5市的水稻区域代替2007年苏南水稻区域。改进后的水稻分布更准确合理，稻田面积相对误差变为1.6%。

3.2水稻估产模型的确定及估产结果的评估

在PASW中采用多元线性逐步回归法建立估产模型。由于t为偏回归系数为0的假设检验的t值，具有较好预测效果的变量t值应大于2或者小于-2。Sig为F值大于F临界值的概率，方差分析结果表明，当回归方程包含不同的自变量时，其显著概率值均小于0.001，即拒绝回归系数均为0的原假设。共线性的判断指标主要是容忍度（Tolerance）接近于1，方差膨胀因子VIF值都不大，此时自变量之间可视为不存在共线性问题[20]。所以最终选定估产模型为：

Yield=7 494.581-0.406×EVI217+0.778×EVI273。 （4）

式中：Yield为水稻单产值，EVI217、EVI273分别为全年第217、273天。

最终选择的估产模型以全年第217、273天作为因变量，这2个时期分别是8月上旬中旬和10月初，黄敬峰等利用盆栽试验和小区试验研究证实了孕穗期到抽穗期为建立水稻遥感估产模型的最佳时期[21]。模型选择的全年第217天，江苏全省绝大部分地区正处在孕穗期与抽穗期之间，恰好是水稻估产的最佳时期，而全年第273天全省则是介于黄熟期左右，植被指数也很高；所以从理论角度来看，所选估产模型符合条件。

最终的估产结果评估见表1：各个地级市的估产误差在-7%～10%之间浮动，全省的平均估产误差为0.54%。正常情况下，由于CTIF去云处理，删除了部分像元，再加上气象等因素的影响，整体的估产结果应该偏低，但是全省结果却偏高。通过数据可知，单产估计偏高的分别为无锡、苏州、南通、扬州、镇江和泰州。除了泰州位于苏北外，其他地级市都位于江苏省中南部，而2007年这些地区提取的水稻面积都相对较少，这样也就造成了单位面积产量的增加。但是全省估产结果精确度均值达到99.46%，各地市的精确度基本维持在95%左右，所以本方法对于大范围的水稻估产能达到很高精度，具有一定指导意义。

3.3水稻单位面积产量空间化

在将江苏省第217、225、233、241天之间的4个8 d的EVI做MVC处理之后，得到EVIp，再做各个地级市的平均值得到EVIavg，再通过公式3完成2007年基于每个像元的水稻单位面积产量图。基于地级市的水稻单产图与空间化的产量图如图1。

江苏省的单产在改革开放后总体上呈上升趋势，在全省范围内，稻田单产最高的是苏南地区，其次是苏中和苏北地区，但是这之间的差距也在渐渐缩小，并且有赶超之势。从图1中的江苏省地级市单产图可以看出，2007年，江苏省单产最高的市是南通市，整体上，苏南、苏中地区的单产比苏北高，这个结果符合实际。而基于像元的空间化产量图则应该呈现苏北较低，苏中、苏南较高的情形。而图中显示的正是如此，单产较高的点主要分布在南通、泰州、盐城、连云港、无锡等地，这与江苏省单产统计数据相符。

综上所述，这种水稻单位面积产量空间化的方法是可行的，而且结果也有一定准确度。空间化后的稻田单产图相对于原来以市为单位的单产图更加直观，能迅速了解更小区域范围的稻田生产力。主要作用有：对于具体了解水稻产区单位面积产量具有参考价值；对于分析全省水稻针对城市扩张、人口增长下的增产潜力具有辅助功效；对于全省范围内的水稻品种改良和栽培技术改进具有指导意义；对于更大范围内稻田的增产与城市发展相协调具有借鉴作用。

4小结

本研究采用MOD09A1数据，计算增强型植被指数EVI和地表水分指数LSWI，通过MOD09自带的质量评估信息QA建立云掩膜，用时间序列插值法CTIF对EVI进行去云处理；再结合去云过后的EVI和LSWI依次进行水稻面积提取，估产和水稻单产的空间化操作。各个地级市的水稻单产预测的相对误差在-7%～10%之间，全省估产的精确度达到9946%，各市的精确度在95%左右。这对于农业政策的制定、粮食价格的宏观调控都具有一定的指导意义。而水稻单产具体到像元程度时，则能更直观地了解全省范围各地水稻生产率状况，更方便指导水稻品种改良和栽培技术改进等工作，使水稻产量与城市扩张和人口增长相协调，促进和谐社会的发展。

参考文献：

[1]中国农业年鉴编委会.中国农业年鉴：2000[M]. 北京：中国农业出版社，2001：264.

[2]陈劲松，黄健熙，林珲，等. 基于遥感信息和作物生长模型同化的水稻估产方法研究[J]. 中国科学：信息科学，2010，40（增刊1）：173-183.

[3]陆登槐. 应用高新技术监测农业资源[J]. 农村实用工程技术，1993（12）：3.

[4]陆登槐. 农业遥感的应用效益及在我国的发展战略[J]. 农业工程学报，1998，14（3）：69-75.

[5]崔霞，冯琦胜，梁天刚. 基于遥感技术的植被净初级生产力研究进展[J]. 草业科学，2007，24（10）：36-42.

[6]陈雯，孙成明，刘涛，等. 植被NPP估算的模型方法与机理[J]. 安徽农学通报，2014，20（8）：30-33，117.

[7]张建华. 作物估产的遥感-数值模拟方法[J]. 干旱区资源与环境，2000，14（2）：82-86.

[8]佴军，张洪程，陆建飞. 江苏省水稻生产30年地域格局变化及影响因素分析[J]. 中国农业科学，2012，45（16）：3446-3452.

[9]赵锐，汤君友，何隆华. 江苏省水稻长势遥感监测与估产[J]. 国土资源遥感，2002（3）：9-11.

[10]Xiao X M，Boles S，Liu J Y，et al. Mapping paddy rice agriculture in southern China using multi-temporal MODIS images[J]. Remote Sensing of Environment，2005，95（4）：480-492.

[11]孙华生. 利用多时相MODIS数据提取中国水稻种植面积和长势信息[D]. 杭州：浙江大学，2009.

[12]邓睿，黄敬峰，王福民，等. 基于中分辨率成像光谱仪（MODIS）数据的水稻遥感估产研究——以江苏省为例[J]. 中国水稻科学，2010，24（1）：87-92.

[13]Huete，R A，Liu H Q，Batchiy K，et al. A comparison of vegetation indices global set of TM images for EOS-MODIS[J]. Remote Sensing of Environment，1997，59：440-451.

[14]Xiao X，He L，Salas W，et al. Quantitative relationships between field-measured leaf area index and vegetation index derived from VEGETATION images for paddy rice fields[J]. International Journal of Remote Sensing，2002，23（18）：3595-3604.

[15]Xiao M X，Boles，Frolking S，et al. Observation of flooding and rice transplanting of Paddy rice fields at the site to landscape scales in China using vegetaion sensor data[J]. International Journal of Remote Sensing，2002，23：3009-3022.

[16]王正兴，王亚琴. 利用LDOPE工具解码MODIS陆地产品质量信息[J]. 遥感技术与应用，2013，28（3）：459-466.

[17]Groten S E. NDVI-crop monitoring and early yield assessment of Burkina Faso[J]. Intemational Joumal of Remote Sensing，1993，14（8）：1495-1515.

[18]Sun H S，Huang J F，Huete A R，et al. Mapping paddy rice with multi-date moderate-resolution imaging spectroradiometer （MODIS） data in China[J]. Journal of Zhejiang University-Science A，2009，10（10）：1509-1522.

[19]Cai X L，Sharma B R，sensing I R. Census and weather data for an assessment of rice yield，water consumption and water productivity in the Indo-Gangetic river basin[J]. Agricultural Water Management，2010，97：309-316.

[20]卢纹岱. SPSS统计分析[M]. 北京：电子工业出版社，2010.

[21]黄敬峰，王人潮，蒋亨显，等. 基于GIS的浙江省水稻遥感估产最佳时相选择[J]. 应用生态学报，2002，13（3）：290-294.李静，宋飞虎，浦宏杰，等. 苹果控制排湿压力微波干燥模型研究[J]. 江苏农业科学，2015，43（11：529-532.