农业遥感综合实习 —— 土地生产潜力研究

1. 农业遥感综合实习——土地生产潜力研究 南京信息工程大学遥感学院

2. 遥感是地球系统研究中应用非常广泛的一门科学与技术，遥感信息因为具有周期性、现实性，宏观性和系统性方面的优势，已广泛应用于地球系统的气象、农业、林业、地质、海洋、气象、水文、军事、环保等领域，正在向产业化、业务化的发展方向不断推进。各个部门对遥感专业人才的需求也日益凸显，对于培养遥感人才的高校面临着如何高质量的培养适合社会需求的遥感人才考验。遥感作为一门技术性很强的专业，加强实习实践环节教学是非常必要的。本研究设计了一个遥感和GIS软件综合应用实例。遥感是地球系统研究中应用非常广泛的一门科学与技术，遥感信息因为具有周期性、现实性，宏观性和系统性方面的优势，已广泛应用于地球系统的气象、农业、林业、地质、海洋、气象、水文、军事、环保等领域，正在向产业化、业务化的发展方向不断推进。各个部门对遥感专业人才的需求也日益凸显，对于培养遥感人才的高校面临着如何高质量的培养适合社会需求的遥感人才考验。遥感作为一门技术性很强的专业，加强实习实践环节教学是非常必要的。本研究设计了一个遥感和GIS软件综合应用实例。

3. 一、实习目的 遥感软件和GIS软件的有机结合，针对农业遥感的重要研究内容土地遥感，通过系统学习一种遥感软件，并结合地理信息系统软件对遥感数据进行处理、分析、应用实现遥感专业实习内容的系统化，通过一个综合实习来实习土地生产潜力评价，并建立土地生产潜力遥感反演模型。

4. 二、原理与方法 选择县级区域基于数字高程模型(DEM)数据和气象站观测资料，在遥感、GIS技术支持下，首先对起伏地形条件下的温度、降水等生态环境因子进行空间插值生成空间网格数据。综合考虑光照、温度、水分、土壤等因子，依据光能利用率模型，通过对光合生产潜力—光温生产潜力—气候生产潜力—土地生产潜力几个阶段的逐步订正来建立适宜县级土地生产潜力评价模型；对典型县级区域进行土地生产潜力评价，最后通过遥感作物长势分析对模型结果进行验证，并建立遥感反演土地生产潜力定量模型。

5. 本实习具体内容： 1、遥感图像处理中的图像校正、图像增强、图像分类、植被指数计算等。 2、定量遥感中涉及到的某个参数的定量反演。 3、数字高程模型中所涉及到的DEM的建立及相关分析。 4、地图学及实习所涉及到的地图数字化、地图制图等。 5、遥感影像投影定义与转换。 6、地理信息系统及其实习所涉及到的GIS空间分析。 7、点—面的各种插值方法。 8、地统计分析与GIS相结合，进行空间统计分析。

6. 二、原理与方法 方法：首先将气象站点的温度，降水等数据进行空间插值，形成栅格数据；依次计算出光合、光温、气候生产潜力；考虑土壤因素计算出土壤有效系数（由于时间的限制，土壤有效系数在本实习中给出结果数据，但在实习中仍然给出计算过程）；根据气候生产力跟土壤有效系数计算耕地区域的土地生产潜力，然后用NDVI指数验证结果，并求出遥感反演模型，通过模型对土地生产潜力进行定量反演（事实上是不能这样反演，但因为时间的限制我们只能通过这种办法让学生进行操作）。

7. 波段合成 图像增强 空间化 太阳辐射 投影转换 几何纠正 气象站点数据 温度 等高线数字化 遥感影像处理 坡度 降水 DEM 数据 高程 解译 耕地 坡向 影像波段运算 TM 数据 植被指数 土壤属性参数 土壤普查数据 土壤养分 • 土壤有效系数获取 地形参 数获取 气象参 数获取 土地生产潜 力模型建立 模型验证 典型区域气候生 产潜力空间分布图 遥感定量反演 图1 综合实习技术路线设计

8. 三、实习仪器与数据 • 计算机，TM遥感影像、土地利用、NDVI数据、DEM、气象站点数据（日照、气温、降水）

9. 四、软件实现步骤 4.1 TM遥感影像处理 TM遥感影像处理，用于目视解译耕地。涉及的技术有：遥感影像校正、增强处理、图像目视解译以及地图数字化、制图等。 4.4.1 遥感图像的几何校正 1）ERDAS图标面板工具条：点击DataPrep图标，→Image Geometric Correction →打开Set Geo-Correction Input File对话框（图2）。 图2 Set Geo-Correction Input File对话框

10. 在Set Geo-Correction Input File对话框（图2）中，需要确定校正图像。 2）图像几何校正的计算模型（Geometric Correction Model） ERDAS提供的图像几何校正模型有7种，具体功能如下： 表1 几何校正计算模型与功能

11. 3）图像校正的具体过程 第一步：显示图像文件（Display Image Files） 首先，在ERDAS图标面板中打开两个视窗 （Viewer1/Viewer2），并将两个视窗平铺放置，操作过程如下： ERDAS图表面板菜单条：Session→Title Viewers 然后，在Viewer1中打开需要校正的图像； 在Viewer2中打开作为地理参考的校正过的图像。 第二步：启动几何校正模块（Geometric Correction Tool） Viewer1菜单条：Raster→ Geometric Correction →打开Set Geometric Model对话框（2） →选择多项式几何校正模型：Polynomial→OK →同时打开Geo Correction Tools对话框（3）和Polynomial Model Properties对话框（4）。

12. 在Polynomial Model Properties对话框中，定义多项式模型参数以及投影参数： →定义多项式次方（Polynomial Order）:2 →定义投影参数：（PROJECTION）:略 →Apply→Close →打开GCP Tool Referense Setup 对话框（5） 图3 Set Geometric Model对话框 图4 Geo Correction Tools对话框

13. 图5 Polynomial Properties对话框 图6 GCP Tool Referense Setup 对话框

14. 第三步：启动控制点工具（Start GCP Tools） 图7 Viewer Selection Instructions 首先，在GCP Tool Referense Setup对话框（图5）中选择采点模式： →选择视窗采点模式：Existing Viewer→OK →打开Viewer Selection Instructions指示器（图2-6） →在显示作为地理参考图像panAtlanta，img的Viewer2中点击左键 →打开reference Map Information 提示框（图2-7）；→OK →此时，整个屏幕将自动变化为如图7所示的状态，表明控制点工具被启动，进入控制点采点状态。

15. 图8 reference Map Information 提示框 第四步：采集地面控制点（Ground Control Point） GCP的具体采集过程： 在图像几何校正过程中，采集控制点是一项非常重要和繁重的工作，具体过程如下：

16. 1、 在GCP工具对话框中，点击Select GCP图表，进入GCP选择状态； 2、 在GCP数据表中，将输入GCP的颜色设置为比较明显的黄色。 3、 在Viewer1中移动关联方框位置，寻找明显的地物特征点，作为输入GCP。 4、 在GCP工具对话框中，点击Create GCP图标，并在Viewer3中点击左键定点，GCP数据表将记录一个输入GCP，包括其编号、标识码、X坐标和Y坐标。 5、 在GCP对话框中，点击Select GCP图标，重新进入GCP选择状态。 6、 在GCP数据表中，将参考GCP的颜色设置为比较明显的红色， 7、 在Viewer2中，移动关联方框位置，寻找对应的地物特征点，作为参考GCP。 8、 在GCP工具对话框中，点击Create GCP图标，并在Viewer4中点击左肩顶巅，系统将自动将参考点的坐标（X、Y）显示在GCP数据表中。 9、在GCP对话框中，点击SelectGCP图标，重新进入GCP选择状态，并将光标移回到Viewer1中，准备采集另一个输入控制点。 10、不断重复1-9，采集若干控制点GCP，直到满足所选定的几何模型为止，尔后，没采集一个InputGCP，系统就自动产生一个Ref. GCP，通过移动Ref. GCP可以优化校正模型。

17. 第五步：采集地面检查点（Ground Check Point） 以上采集的 GCP的类型均为控制点，用于控制计算，建立转换模型及多项式方程。下面所要采集的GCP类型是检查点。（略） 第六步：计算转换模型（Compute Transformation） 在控制点采集过程中，一般是设置为自动转换计算模型。所以随着控制点采集过程的完成，转换模型就自动计算生成。 在Geo-Correction Tools对话框中，点击Display Model Properties 图表，可以查阅模型。 第七步：图像重采样（Resample the Image） 重采样过程就是依据未校正图像的像元值，计算生成一幅校正图像的过程。原图像中所有删格数据层都要进行重采样。 ERDAS IMAGE 提供了三种最常用的重采样方法。略

18. 图像重采样的过程： 首先，在Geo-Correction Tools对话框中选择Image Resample 图标。 然后，在Image Resample对话框中，定义重采样参数； →输出图像文件明（OutputFile）:rectify.img →选择重采样方法（Resample Method）:Nearest Neighbor →定义输出图像范围： →定义输出像元的大小： →设置输出统计中忽略零值: →定义重新计算输出缺省值。 第八步：保存几何校正模式（Save rectification Model） 在Geo-Correction Tools对话框中点击Exit按钮，推出几何校正过程，按照系统提示，选择保存图像几何校正模式，并定义模式文件，以便下一次直接利用。 第九步：检验校正结果（Verify rectification Result） 基本方法：同时在两个视窗中打开两幅图像，一幅是矫正以后的图像，一幅是当时的参考图像，通过视窗地理连接功能，及查询光标功能进行目视定性检验。

19. 4.4.2 遥感图像的增强处理 通过上机操作，了解空间增强、辐射增强几种遥感图象增强处理的过程和方法，加深对图象增强处理的理解。 ERDAS IMAGE图像解译模块主要包括了图像的空间增强、辐射增强、光谱增强、高光谱工具、傅立叶变换、地形分析以及其他实用功能。 空间增强技术是利用像元自身及其周围像元的灰度值进行运算，达到增强整个图像之目的。卷积增强（Convolution）是空间增强的一种方法。 卷积增强（Convolution）时将整个像元分块进行平均处理，用于改变图像的空间频率特征 。卷积增强（Convolution）处理的关键是卷计算子----系数矩阵的选择。该系数矩阵又称卷积核（Kernal）。ERDAS IMAGINE将常用的卷计算子放在一个名为default.klb的文件中，分为3*3，5*5、7*7三组，每组又包括“EdgeDetect/Low Pass/Horizontal/Vertical/Summary”等七种不同的处理方式。具体执行过程如下：

20. ERDAS图标面板菜单条：Main→Image Interpreter→Spatial enhancement→convolution→convolution对话框。 图9 Convolution对话框 几个重要参数的设置： 边缘处理方法：（Handle Edges by）：Reflection 卷积归一化处理：Normalize the Kernel

21. 4.2 遥感图像目视解译 涉及的技术有：遥感图像目视解译，数字化与地图制图等。 主要通过ArcView GIS 3.2数字化提取耕地类型。 4.2.1 数字化 1、建面状图层 ①打开遥感影像 ArcView GIS 3.2—Cancel—New—Add Theme—路径—Image Data Source—遥感影像—Ok ②建立图层 View—New Theme（Polygon）—路径—名称 图10 建立图层

22. 2、数字化提取影像图中的耕地要素，保存要素与属性值。2、数字化提取影像图中的耕地要素，保存要素与属性值。 ①数字化:用面数字画工具提取出耕地类型； 图11 数字化工具 ②打开属性表Open Theme Table，输入属性值； ③保存数字化结果Save edits—路径—名称。

23. 4.2.2 进行图例设计、地图注记 1、对数字化面要素图例进行按图中特征进行设计，并保存图例 ①在图层上双击，打开Edit Legend对话框，设计图例。 ②Edit Legend对话框，用Save保存图例文件。 图12图例设计

24. 2、按着所给图注记进行地图注记。 ①用theme —Stop Editing退出编辑状态。 ②用Text tool图标进行注记。 4.2.3 进行图例设计、地图注记 输出JPEG格式图，要求有图例、指北针、比例尺（单位为米）。 ①设置单位View Properties—meters。 ②view—Layout进入制图窗口，进行图名、图例、指北针、比例尺编辑。 ③File—Export—输出JPEG格式图。 4.3 DEM模型的建立 涉及的技术有：等高线数字化，建立DEM模型及其地形参数分析等；数字化等高线同目视解译耕地数字化类似，只是建立的是线状图层。 地形因子提取 常用的地形因子可以划分为微观地形因子与宏观地形因子两种基本类型。按照提取地形因子差分计算的阶数，又可将地形因子分为一阶地形因子、二阶地形因子和高阶地形因子（图13和14）。

25. 图13 依据空间区域范围的坡面因子分类体系

26. 图14 基于提取算法的坡面因子分类体系

27. 1. 坡度的提取 地表面任一点的坡度（Slope）是指过该点的切平面与水平地面的夹角。坡度表示了地表面在该点的倾斜程度。 实际应用中，坡度有两种表示方式方法： （1） 坡度(degree of slope)：既水平面与地形面之间夹角。 （2） 坡度百分比（percent slope）：既高程增量与水平增量之比的百分数。 （3） ArcGIS中坡度的提取过程为： （4） 在Spatial Analyst下拉菜单中选择Surface analysis, 在弹出的下一级菜单中点击Slope，出现Slope对话框，如图15，以下所有设置如图中所示。 （5） 在Input Surface的下拉菜单中选择用来生成坡度的表面； （6） 选择一种坡度表示方法，在此分别用两种方法做了坡度。 （7） 在Z factor栏中设定高程变换系数； （8） 在Output cell size栏中设定栅格大小； （9） 在Output raster栏中指定输出坡度的存放路径与文件名。 （10） 点击OK按钮。

28. 图15 坡度计算

29. 2. 坡向的提取 坡向指地表面上一点的切平面的法线矢量在水平面的投影与过该点的正北方向的夹角。对于地面任何一点来说，坡向表征了该点高程值改变量的最大变化方向。在输出的坡向数据中，坡向值有如下规定：正北方向为0度，按顺时针方向计算，取值范围为0°～360°。 ArcGIS中坡向的提取过程为： （1） 在Spatial Analyst下拉菜单中选择Surface analysis, 在弹出的下一级菜单中点击Aspect，出现Aspect对话框，如图8.40，以下所有设置如图中所示。 （2） 在Input Surface的下拉菜单中选择用来生成坡向的表面； （3） 在Output cell size栏中设定栅格大小； （4） 在Output raster栏中指定输出坡向的存放路径与文件名。

30. 图16 坡向计算

31. 4.4 气象站点数据空间插值 涉及的技术有：点­­—面的插值方法。 在ArcGIS中，表面分析的主要功能有：查询表面值、从表面获取坡度和坡向信息、创建等值线、分析表面的可视性、从表面计算山体的阴影、确定坡面线的高度、寻找最陡路径、计算面积和体积、数据重分类、将表面转化为矢量数据等。在本实习中主要介绍ArcGIS表面分析中的栅格插值，坡度、坡向等基本地形因子的提取等基本分析功能。 一般情况下采集到的数据都是以离散点的形式存在的，只有在这些采样点上才有较为准确的数值，而其它未采样点上都没有数值。然而，在实际应用中却很可能需要用到某些未采样点的值，这个时候就需要通过已采样点的数值来推算未采样点值。这样的一个过程也就是栅格插值过程。插值结果将生成一个连续的表面，在这个连续表面上可以得到每一点的值。栅格插值包括简单栅格表面的生成和栅格数据重采样。

32. 在ArcGIS 9 栅格分析模块中，通过栅格插值运算生成表面主要有三种实现方式：反距离权重插值，样条函数插值和克里格插值，如图17所示。 图17 栅格插值方法

33. 下面以一组土壤元素PH值的插值来逐一说明在Spatial Analyst中三种表面生成插值的实现过程。 1. 反距离权重插值（IDW） （1） 在Spatial Analyst下拉菜单中选择Interpolate to Raster，在弹出的下一级菜单中点击Inverse Distance Weighted命令, 弹出IDW对话框，如图8.25。 （2） 在Input points的下拉菜单中选择被用来进行插值的离散点数据； （3） 在Z value field的下拉菜单中选择要加入的字段； （4） 在Power栏中填入进行插值计算的幂值；幂值就是距离的指数。如幂指数为2时则进行反向距离平方插值。幂指数是一个正实数，其缺省值为2。 （5） 在Search radius type 栏中选择一种搜索半径设置类型； 1) Variable：当选择此项时,搜索半径由下面两Maximum distance。首先在Number of points中输入搜索的最近点的个数（缺省值为12），然后在Maximum distance中输入一个控制距离。如果最近点的个数超出控制距离，则将会以控制距离为限制来选取较少的点；

34. 2) Fixed ，Distance和Minimum number of points。首先在Distance中输入搜索半径距离（缺省值是输出栅格大小的五倍），然后在Minimum number of points中输入控制插值点个数的最小整数值。如果搜索半径距离内的点个数小于插值点个数的最小整数值，则搜索半径自动增大。 （6） Use barriers polyline为可选项，输入一个中断线文件。barriers是在插值中，如有某些地方出现异常，（如某些断裂带），而要求插值时考虑到这样的因素，所设置的选项。它是一个打断表面的线特征。这一线特征没有Z值。悬崖，峭壁，堤岸或某些障碍都是典型的barriers。barriers限制了插值计算，它使得计算只在线的两侧各自进行。而落在线上的点则会同时参与线两侧的计算。 （7） Output cell size：指定输出结果的栅格大小； （8） Output raster：为输出结果指定目录及名称； （9） 点击OK按钮。

35. 2、样条函数插值（SPLINE） 样条函数插值采用两种不同的计算方法，选择Regularized生成一个平滑、渐变的表面，得出的插值结果很可能会超出样本点的取值范围。如选择Tension，它会根据要生成的现象的特征生成一个比较坚硬的表面，得出结果的插值更接近限制在样本点的取值范围内。 同时，计算过程中，需要选权重（weight）。选择Regularized时，他决定了表面最小曲率三次导的权重。权重越高表面越光滑。可能用到的典型值有：0、0.001、0.01、1和5。选择Tension时，他决定了Tension的权重。权重越高，表面越粗糙。可能用到的典型值有：0、1、5和10。样条函数插值过程如下： （1）在Spatial Analyst 下拉菜单中选择Interpolate to Raster, 在弹出的下一级菜单中点击Spline，出现Spline对话框，如图8.27，以下所有设置如图中所示。（2） 在Input points的下拉菜单中选择被用来进行插值的离散点数据；在Z value field的下拉菜单中选择要加入的字段；

36. 图18 Spline对话框

37. （3） 在Spline type中填入样条函数插值的类型； （4） 在Weight栏中填入一个影响插值的特征参数； （5）在Number of points 栏中输入每一个区域内用来估值点的个数。它缺省是12. （6）指定输出结果的栅格大小； （7） 为输出结果指定目录及名称； （8） 点击OK按钮。 3. 克里格插值（Kriging） （1） 在Spatial Analyst下拉菜单中选择Interpolate to Raster, 在弹出的下一级菜单中点击 Kriging，出现Kriging对话框，如图8.29，以下所有设置如图中所示。 （2） 在Input points的下拉菜单中选择被用来进行插值的离散点数据； （3） 在Z value field的下拉菜单中选择要加入的字段； （4） 选择你所需要的克里格方法； （5） 在Semivariogram model的下拉菜单中选择你所需要使用的模型；

38. （6） 点击Advanced Parameters按钮，如果知道可以指定这些参数，另外空间分析模块也将为你估算这些参数； （7） 在Search radius type 栏中选择一种搜索半径设置类型； （8） 指定输出结果的栅格大小； （9） 可选项，是否需要生成预测的标准误差； （10） 为输出结果指定目录及名称； （11） 点击OK按钮。 4. 数据重采样 重采样是栅格数据空间分析中处理栅格分辨率匹配问题的常用数据处理方法。进行空间分析时，用来分析的数据资料由于来源不同，经常会出现栅格不同大小问题，这时为了便于分析，就需要将不同的栅格大小转化为同样栅格大小，这样的一个过程就是栅格数据的重采样过程。栅格数据的重采样主要基于三种方法：最邻近采样（NEAREST），双线性ILINEAR）和三次卷积采样（CUBIC）。

39. （1) 最邻近采样：此重采样法用输入栅格样后的输出栅格中的每个栅格值, 都是输入栅格数据中真实存在而未加任何改变的值。这种方法简单易用，计算量小，重采样的速度最快。 （2) 双线性采样：此重采样法取待采样点（x，y）点周围四个邻点，在y方向（或x方向）内插两次，再在x方向（或y方向）内插一次，得到（x，y）点的栅格值。 （3) 三次卷积采样：这是进一步提高内插精度的一种方法。它的基本思想是增加邻点来获得最佳插值函数。取待计算点周围相邻的16个点，与双线性采样类似，可先在某一方向上内插，如先在x方向上，每四个值依次内插四次，再根据四次的计算结果在y方向上内插，最终得到内插结果。 在ArcGIS中重采样功能是在ArcToolbox下来实现的。在ArcToolbox中，Data Management TOOls下的Raster子菜单中，选择选择Resample命令（图8.31）。参数说明：

40. （1) Input raster：选择输入栅格数据； （2) Output raster：设置输出栅格数据路径及名称； （3) Output cell size：设置输出栅格单元大小； （4) Resampling technique：选择重采样方法。 1） NEAREST:最邻近距离法； 2） BILINEAR：双线性内插法； 3） CUBIC：立方体法。 图18 重采样参数设置框

41. 4.5土地生产潜力计算 涉及的技术有：地统计分析与GIS相结合，进行空间分析。 “机制法”是应用最为广泛的土地生产潜力计算方法，它根据作物生产力形成的机理，考虑光、温、水、土等自然生态因子，从作物截光特征和光合作用入手，依据作物能量转化及产量形成过程，进行逐步“衰减”估算土地生产潜力，其函数式为： YL = Q·f(Q)·f(T)·f(W)·f(S) = YQ·f(T)·f(W)·f(S) = YT·f(W)·f(S) = YW·f(S) (1) 式中，YL为土地生产潜力；Q为太阳总辐射；f(Q)为光合有效系数；YQ为光合生产潜力；f(T)为温度有效系数；YT为光温生产潜力；f(W)为水分有效系数；YW为气候生产潜力；f(S)为土壤有效系数。

42. 4.5.1 计算光合、光温、气候生产潜力 1.光合生产潜力 YQ = 0.92·Q 式中，YQ指年光合生产潜力（k g/h m2·年）；Q指到达地面的太阳年总辐射（J/cm2·年）；0.92是光合转换系数。具体推算方法见文献[16] 。 我们考虑地形因子影响，则太阳倾斜照射（太阳年总辐射Q1）与太阳正射（太阳年总辐射Q2）的关系为：Q2 = Q1·cos(a) 式中a为坡度。 (1) 在ArcMap下进行投影转化。 <1> 首先在ArcMap下定义南阳市行政区划图层（Arctoolbox->Data Management tools->Projection and Transformation->Define Projection）。

43. 图20 投影定义

44. 注：将南阳市行政区划图层定义为Lambert conditional conic投影，参数为： False_Easting 0 False_Northing 0 Central_Meridian 110 Standard_Parallel_1 25 Standard_Parallel_2 47 Scale_Factor 1 Latitude_Of_Origin 10

45. 地理坐标系统选择1984WGS坐标系统。 Select -> geographic corditional system -> world -> 1984WGS。

46. <2> 然后转换图层投影（矢量图层Arctoolbox-> Data Management tools->Projection and Transformation->feature->projec；栅格图层Data management->Projection ->raster->project）。 注：将南阳市行政区划图层转换为DEM投影方式。

47. (2) 在ArcMap下将南阳市DEM剪切下来 Arctoolbox->spatial analyst tools->extraction->Extract by mask。

48. (3) 在ArcView下将南阳市DEM转化成Grid格式 点击File->extention->spatial analyst->ok；点击标题栏上的Theme->convert to Grid->保存。

49. (4) 在ArcView 下计算光合生产潜力 <1> 计算南阳市地形坡向、坡度 点击File->extention->spatial analyst->ok；点击标题栏上的surface->derive aspect->坡向图->surface->derive slope->坡度图。

50. <2> 计算南阳市光合生产潜力 点击标题栏上的analysis->Map Calculator->双击坡度图->点击cos函数->点击evaluate->坡度的余弦图->点击标题栏上的analysis->Map Calculator->双击坡度的余弦图->点击Abs函数->点击evaluate->坡度的余弦正值->点击标题栏上的analysis->Map Calculator->双击坡度的余弦正值图->乘以年太阳总辐射->点击evaluate->南阳市光合生产潜力分布图。