基于地质统计学影像纹理的海南矿区荒漠化监测

D01:10.13374j.isml00103x2006.08.001 第28卷第8期 北京科技大学学报 Vol.28 Na 8 2006年8月 Journal of University of Science and Technology Beijing Aug.2006 基于地质统计学影像纹理的海南矿区荒漠化监测 廖楚江)王长耀) 式江)黄志强到 1)中国科学院遥感应用研究所遥感科学国家重点实验室，北京1001012)海南省遥感中心，海口570226 3)广西地勋总院.南宁530007 摘要多年来，由于对钛矿的无序开采，使得海南岛东部出现大面积的土地荒漠化.采用遥感的 手段进行跟踪监测，合理地授子采矿权，组织适当的复垦，是解决当地荒漠化的有效途径.基于不 同沙地类型在地表空间结构上的差异，提出将基于地质统计学的影像纹理应用到荒漠化监测中， 通过变异函数纹理来加大各种不同类别沙地间的区别，提高样本选择的分离度.结果表明.运用变 异函数纹理结合光谱波段的最大似然分类方法能够很好地界定海滩沙地和内陆荒漠地的等级，最 高分类精度达到92.4%，证明了基于地质统计学的影像纹理在实现该地区遥感荒漠化监测方面的 有效性. 关键词荒漠化：沙地：地质统计学；纹理：变异函数 分类号TP751.1:TP79 作为全国钛铁矿最为丰富的海南省，近年来 区不曾有过的问题，因此必须寻求新的解决途径 由于矿产的无序开发，岛东部的文昌、琼海、万宁 来实现该地区的荒漠化监测.本文将以海南文昌 等地.出现了大面积的土地荒漠化，而且这种趋势 市为例，探讨基于地质统计学的影像纹理在实现 在逐步加剧.海南岛隶属热带湿润地区，在荒漠 该地区荒漠化评价中的应用. 化的成因上，与我国西北干旱半干旱地区的沙质 荒漠化存在明显的差别，该地区的荒漠化主要是 1研究区域和数据描述 由于人工过度开采钛矿所造成，并受海洋季风的 文昌市位于海南岛的东北端，东面和北面临 影响逐步形成与西北地区相似的流动沙丘和沙 海，中心经纬度坐标为东经110.717°，北纬 地.为防止土地荒漠化而完全不开采钛矿绝非可 19.617°，图1为研究区域相对于全岛的位置.为 能，而合理有序的开采和复垦的结合完全可以防 便于分析和对比，选取的三幅影像分别为1986- 止土地的大面积荒漠化.因此.必须做好实时监 12-04的TM数据、1999-12-24和2001-12-24的 测，合理地授予采矿权，防止过度无序的开采. ETM十数据（皆为1048像素×805像素），其中 目前，我国的荒漠化研究主要集中在西北地 前两幅影像时间跨度较大，便于分析荒漠化的总 区，其中基于遥感的荒漠化评价方法主要有两种： 体趋势，后两幅影像时间跨度较小，便于分析不同 一是运用图像处理软件，通过监督和非监督分类， 程度荒漠化的细微变化. 直接划分类型和程度：另一种是选择几个基于 对两幅影像初步分析得出，与西北地区荒漠 RS,GIS的指标，给定不同的权重，通过综合来得 化监测不同的是，该地区临海，海滩沙地和内陆荒 出结果.海南是我国东南沿海地区惟一受荒漠化 漠地在图像上呈现几乎相同的光谱特征（图】中 威胁较为严重的省份，研究积累较少，既缺乏足够 影像解译经验指导该地区荒漠化程度的正确分 类，也缺乏充分实地调查建立的适合该地区的荒 海南省 漠化监测指标体系，而且监测中还将遇到其他地 收稿日期：2005-07-21修回日期：2005-1009 基金项目：中国科学院知识创筋工程重大项目(No.KZCX一SW- 0 10 km 20 km 0一02)，国家”863计划资助项目(Na2003AA131130)和r华南 热带作物区生态地质环境变化动态遥感监测”课题 图1研究区域图示 作者简介：廖楚江(1979一)，男，博士研究生 Fig.1 Region for research

基于地质统计学影像纹理的海南矿区荒漠化监测 廖楚江1) 王长耀1) 丁式江2) 黄志强3) 1) 中国科学院遥感应用研究所遥感科学国家重点实验室, 北京 100101 2) 海南省遥感中心, 海口 570226 3) 广西地勘总院, 南宁 530007 摘 要 多年来, 由于对钛矿的无序开采, 使得海南岛东部出现大面积的土地荒漠化.采用遥感的 手段进行跟踪监测, 合理地授予采矿权, 组织适当的复垦, 是解决当地荒漠化的有效途径.基于不 同沙地类型在地表空间结构上的差异, 提出将基于地质统计学的影像纹理应用到荒漠化监测中, 通过变异函数纹理来加大各种不同类别沙地间的区别, 提高样本选择的分离度.结果表明, 运用变 异函数纹理结合光谱波段的最大似然分类方法能够很好地界定海滩沙地和内陆荒漠地的等级, 最 高分类精度达到 92.4 %, 证明了基于地质统计学的影像纹理在实现该地区遥感荒漠化监测方面的 有效性. 关键词 荒漠化;沙地;地质统计学;纹理;变异函数 分类号 TP751.1 ;TP 79 收稿日期:2005 07 21 修回日期:2005 10 09 基金项目:中国科学院知识创新工程重大项目( No .KZCX1-SW- 01-02) , 国家“ 863 计划” 资助项目( No.2003AA131130) 和“华南 热带作物区生态地质环境变化动态遥感监测”课题 作者简介:廖楚江( 1979—) , 男, 博士研究生 作为全国钛铁矿最为丰富的海南省, 近年来 由于矿产的无序开发, 岛东部的文昌 、琼海 、万宁 等地, 出现了大面积的土地荒漠化, 而且这种趋势 在逐步加剧.海南岛隶属热带湿润地区, 在荒漠 化的成因上, 与我国西北干旱半干旱地区的沙质 荒漠化存在明显的差别, 该地区的荒漠化主要是 由于人工过度开采钛矿所造成, 并受海洋季风的 影响逐步形成与西北地区相似的流动沙丘和沙 地.为防止土地荒漠化而完全不开采钛矿绝非可 能, 而合理有序的开采和复垦的结合完全可以防 止土地的大面积荒漠化.因此, 必须做好实时监 测, 合理地授予采矿权, 防止过度无序的开采 . 目前, 我国的荒漠化研究主要集中在西北地 区, 其中基于遥感的荒漠化评价方法主要有两种: 一是运用图像处理软件, 通过监督和非监督分类, 直接划分类型和程度 ;另一种是选择几个基于 RS, GIS 的指标, 给定不同的权重, 通过综合来得 出结果.海南是我国东南沿海地区惟一受荒漠化 威胁较为严重的省份, 研究积累较少, 既缺乏足够 影像解译经验指导该地区荒漠化程度的正确分 类, 也缺乏充分实地调查建立的适合该地区的荒 漠化监测指标体系, 而且监测中还将遇到其他地 区不曾有过的问题, 因此必须寻求新的解决途径 来实现该地区的荒漠化监测.本文将以海南文昌 市为例, 探讨基于地质统计学的影像纹理在实现 该地区荒漠化评价中的应用. 1 研究区域和数据描述 文昌市位于海南岛的东北端, 东面和北面临 海, 中 心 经 纬度 坐 标为 东 经 110.717°, 北 纬 19.617°, 图 1 为研究区域相对于全岛的位置.为 便于分析和对比, 选取的三幅影像分别为 1986 12 04的 TM 数据、1999 12 24和2001 12 24的 ETM +数据( 皆为 1 048 像素 ×805 像素) , 其中 前两幅影像时间跨度较大, 便于分析荒漠化的总 体趋势, 后两幅影像时间跨度较小, 便于分析不同 程度荒漠化的细微变化. 图 1 研究区域图示 Fig.1 Region for research 对两幅影像初步分析得出, 与西北地区荒漠 化监测不同的是, 该地区临海, 海滩沙地和内陆荒 漠地在图像上呈现几乎相同的光谱特征( 图 1 中 第 28 卷 第 8 期 2006 年 8 月 北 京 科 技 大 学 学 报 Journal of University of Science and Technology Beijing Vol .28 No.8 Aug.2006 DOI :10.13374/j .issn1001 -053x.2006.08.001

。710 北京科技大学学报 2006年第8期 的亮白色)，直接基于TM辐射波段进行分类无 Y()-E[DN(x+k)-DN(x)]2 (1) 法将二者区分开来，因而也就无法得出荒漠化土 地的实际面积：另外，目前在荒漠化监测中成功应 式中()可以阐述为相距为h的一对像素DN 用的荒漠化四分法，即极重度荒漠化、重度荒漠 值平方差的数学期望的/2，即半变异值，Y(h) 化、中度荒漠化和轻度荒漠化在图像上很难界定 是一个依赖位于x十h和x处的两个像素之间步 出来，而无论是直接从图像上解译，或是求取一定 长矢量h的模和角度的矢量函数.另外需要指出 的荒漠化监测指标如MSA VI或NDVI)后的对 的是变异函数有几个别称，包括半变异函数、样本 比，发现都与其他地区荒漠化监测所描述的特征 变异函数和实验变异函数，为了区分变异函数和 相去甚远？，不能加以借鉴 某个特定步长的变异函数值，在此约定，凡是某个 针对三幅影像的特点，本文制定了不同的研 特定步长的变异函数值，一律称为半变异值，而凡 究侧重点，1986年影像显示荒漠土地面积较小， 以步长为自变量，半变异值为因变量的函数成为 而海滩沙地的范围较大，因此对该幅影像的研究 变异函数 侧重于如何更好地分离海滩沙地和荒漠化土地. 2.2空间变异性的纹理计算 而1999年和2001年影像中荒漠化土地虽较 要将地质统计学纹理信息引入到影像分类过 1986年增加了很多，但海滩由于种植沿海防护林 程中，必须计算一系列半变异值.目前，基于变异 的缘故已经变得不明显，因此对这两幅影像侧重 函数的遥感分类比较成功的应用主要有Miranda 于如何更好地进行荒漠化土地类型分级，基于此， 等基于JES-1数据对巴西热带雨林的分类)和 我们既可以利用前者更为精确地得到十几年来该 Chica-Olmo等基于TM数据对西班牙南部的岩 研究区域荒漠化的趋势，又可以利用后者辅助决 石分类.本文使用的变异函数依据计算波段 策者制定荒漠化治理的方向. 的不同，可分为单波段变异函数和多波段变异函 数其中方向变异函数与绝对变异函数属于前者， 2方法 而交叉变异函数和虚拟变异函数可以用来计算两 波段之间的纹理特征 2.1纹理和变异函数 所谓影像中的纹理特征，指的是相同物体在 (1)方向变异函数.方向半变异函数可以参 影像中群集时，是以一种规律的方式排列出现，而 考式(1)得到： 不会呈现完全均匀的灰阶值，这种排列的方式，可 ()2[DN(x)-DN:( 认为是代表地物的某一种特性，这种特性就是组 (2) 织特征到，因此，纹理分析提供了在影像解译和 其中P()为在某一方向上相距步长h的像素对 分类中与土地覆盖类别相关的重要的分离特征. 的数目，DN()为像素和x:十h处的DN值，k 从统计学的角度看，我们可以用两个与DN 为传感器的波段号. 值相关联的概念表达纹理：局部/全局的变化性和 (2)绝对变异函数.绝对变异函数与方向变 空间自相关，前者通常是通过计算一个局部层次 异函数相似，惟一的区别是用绝对差取代了平方 上的方差来达成，即先计算一个移动窗口的均值， 差 后计算DN值与该均值的偏离程度4：而空间自 P(h) 1 相关是假定在影像中DN值不是完全随机分布 Mh1=P和 |DNk(xa)-DNk(xa+h)川 的，每种土地覆盖类别有其独特的空间变异性（空 (3) 间依赖)结构，Lark指出像素之间的变异函数的 (3)交叉变异函数.交叉变异函数对两个波 值依赖于它们的空间自相关，该值可以用作每种 段之间的联合变异性（交叉相关）进行定量化.j, 土地覆盖类别的纹理描述器？.地质统计学可以 k分别代表两个不同的波段. 将这两方面联合起来建模. 从地质统计学的角度而言，我们可以把遥感 ()=p)>I DN,(x)-DN,( 1 影像的DN值看作为场变量，它既是随机分布的， DN&(xG)-DNx(xa+h) (4) 同时也存在一定的空间相关性，在固有假设条件 (4)虚拟变异函数.虚拟变异函数表达的是 下，这两方面可以通过变异函数（也可称为实验变 交叉增量的变异函数，它取代了上文交叉变异函 异函数连接起来用方程表示为9： 数的协方差

的亮白色) , 直接基于 TM 辐射波段进行分类无 法将二者区分开来, 因而也就无法得出荒漠化土 地的实际面积;另外, 目前在荒漠化监测中成功应 用的荒漠化四分法, 即极重度荒漠化、重度荒漠 化、中度荒漠化和轻度荒漠化在图像上很难界定 出来, 而无论是直接从图像上解译, 或是求取一定 的荒漠化监测指标( 如 MSAVI 或 NDV I) 后的对 比, 发现都与其他地区荒漠化监测所描述的特征 相去甚远[ 1 2] , 不能加以借鉴. 针对三幅影像的特点, 本文制定了不同的研 究侧重点, 1986 年影像显示荒漠土地面积较小, 而海滩沙地的范围较大, 因此对该幅影像的研究 侧重于如何更好地分离海滩沙地和荒漠化土地, 而 1999 年和 2001 年影像中荒漠化土地虽较 1986 年增加了很多, 但海滩由于种植沿海防护林 的缘故已经变得不明显, 因此对这两幅影像侧重 于如何更好地进行荒漠化土地类型分级, 基于此, 我们既可以利用前者更为精确地得到十几年来该 研究区域荒漠化的趋势, 又可以利用后者辅助决 策者制定荒漠化治理的方向. 2 方法 2.1 纹理和变异函数 所谓影像中的纹理特征, 指的是相同物体在 影像中群集时, 是以一种规律的方式排列出现, 而 不会呈现完全均匀的灰阶值, 这种排列的方式, 可 认为是代表地物的某一种特性, 这种特性就是组 织特征[ 3] , 因此, 纹理分析提供了在影像解译和 分类中与土地覆盖类别相关的重要的分离特征. 从统计学的角度看, 我们可以用两个与 DN 值相关联的概念表达纹理 :局部/全局的变化性和 空间自相关, 前者通常是通过计算一个局部层次 上的方差来达成, 即先计算一个移动窗口的均值, 后计算 DN 值与该均值的偏离程度 [ 4] ;而空间自 相关是假定在影像中 DN 值不是完全随机分布 的, 每种土地覆盖类别有其独特的空间变异性( 空 间依赖) 结构, Lark 指出像素之间的变异函数的 值依赖于它们的空间自相关, 该值可以用作每种 土地覆盖类别的纹理描述器 [ 5] .地质统计学可以 将这两方面联合起来建模 . 从地质统计学的角度而言, 我们可以把遥感 影像的 DN 值看作为场变量, 它既是随机分布的, 同时也存在一定的空间相关性, 在固有假设条件 下, 这两方面可以通过变异函数( 也可称为实验变 异函数) 连接起来, 用方程表示为[ 6] : γ( h) = 1 2 E[ DN( x +h) -DN( x)] 2 ( 1) 式中 γ( h)可以阐述为相距为 h 的一对像素 DN 值平方差的数学期望的 1/2, 即半变异值, γ( h) 是一个依赖位于 x +h 和x 处的两个像素之间步 长矢量h 的模和角度的矢量函数 .另外需要指出 的是变异函数有几个别称, 包括半变异函数 、样本 变异函数和实验变异函数, 为了区分变异函数和 某个特定步长的变异函数值, 在此约定, 凡是某个 特定步长的变异函数值, 一律称为半变异值, 而凡 以步长为自变量, 半变异值为因变量的函数成为 变异函数 . 2.2 空间变异性的纹理计算 要将地质统计学纹理信息引入到影像分类过 程中, 必须计算一系列半变异值.目前, 基于变异 函数的遥感分类比较成功的应用主要有 Miranda 等基于 JERS-1 数据对巴西热带雨林的分类[ 7] 和 Chica-Olmo 等基于 TM 数据对西班牙南部的岩 石分类 [ 8] .本文使用的变异函数, 依据计算波段 的不同, 可分为单波段变异函数和多波段变异函 数, 其中方向变异函数与绝对变异函数属于前者, 而交叉变异函数和虚拟变异函数可以用来计算两 波段之间的纹理特征 . (1)方向变异函数.方向半变异函数可以参 考式( 1)得到: γk ( h) = 1 2P( h) ∑ P( h) α=1 [ DNk ( xα) -DNk ( xα+h)] 2 ( 2) 其中 P( h)为在某一方向上相距步长 h 的像素对 的数目, DN( ·)为像素 xi 和xi +h 处的 DN 值, k 为传感器的波段号. ( 2) 绝对变异函数.绝对变异函数与方向变 异函数相似, 惟一的区别是用绝对差取代了平方 差. γk ( h) = 1 2P( h) ∑ P( h) α=1 DNk ( xα) -DNk ( xα+h) ( 3) ( 3) 交叉变异函数.交叉变异函数对两个波 段之间的联合变异性(交叉相关)进行定量化.j, k 分别代表两个不同的波段. γjk ( h) = 1 2P( h) ∑ P( h) α=1 DNj( xα) -DNj( xα+h) × DNk ( xα) -DNk ( xα+h) ( 4) ( 4) 虚拟变异函数.虚拟变异函数表达的是 交叉增量的变异函数, 它取代了上文交叉变异函 数的协方差. · 710 · 北 京 科 技 大 学 学 报 2006 年第 8 期

Vol.28 No.8 廖楚江等：基于地质统计学影像纹理的海南矿区荒漠化监测 711。 Ah)=2P万2[DN(-DN十M 1 了证明变异函数区分几种不同沙地类型的能力， 我们在1999年4,3,2波段合成影像上为第二波 (5) 段选择了3个类别的小样本，分别为海滩沙地（近 2.3变异函数参数 海亮白色区域)、极重度荒漠地（内陆亮白色区域） 地质统计学影像纹理是通过计算移动窗口内 和轻度荒漠地（亮白色周边的红白色区域），利用 的相邻像元的半变异值得到的，计算方法如图2 方向变异函数公式进行计算，得到的变异函数图 所示.在计算的过程中需考虑计算窗口大小、像 如图3所示.在研究变异函数图时，一般关注三 元取样距离（步长）、计算方向和波段组合等几个 个参数9，即变差距离、基台值和块金值（如 因子9.计算窗口（移动视窗）大小，主要是依据 图4)，从变异函数的定义可以看出，当一对样本 在计算每个像元的纹理值时，需要同时考虑到该 距离增加时，所对应的变异函数值一般也增加，当 像元周围多大的范围来决定，一般是在反复实验 变异函数达到一个平台时的样本间距离叫做变差 的基础上得到的，一个理想的窗口不能太大，否则 距离或相关尺度.变差距离表明当样本间距等于 将使得相邻的土地覆盖类型互相干扰，同时也不 或大于此距离时，样本就变得完全独立，在变差距 能太小，否则将减小变异函数分类器的使用范围. 离变异函数所达到的平台叫做基台值.块金值是 像元取样距离（步长）依据窗口大小的改变而改 在极短的样本距离(h≈O)之间变异函数从原点 变，它代表了像元与像元之间的关系，图2中窗口 的跳升值（不连续性），它是由于样本误差和短距 大小为3，步长可为1或者2，一般情况下，计算结 离的变异性引起的. 果取所有步长的平均值.考虑计算方向是因为地 1000r 潜在荒漠化 表景物的排列顺序与方向并不固定，因此在计算 800 ··士人A 纹理特征时，需考虑计算东一西向(0)、南一北向 A (90°)、东北一西南向(45°)，西北一东南向(135°)4 个方向的纹理当空间方向性不明显时，可将4个 警 400 严重荒漠化 ■ 方向的变异函数值加总取平均，即成为所谓等向 200 海滩沙地 变异元.波段组合是针对交叉变异函数和虚拟变 异函数的一个参数，选择不同的波段组合，将对地 5 10 步长，h 表景物的区分产生重要影响. 窗口移动方向 图3不同类别变异函数图 1 Fig 3 Variograms of different classes 10198636154 6457534763 3×3 189 窗口 52524651101 基台值 5054658596 9999839489 ( 光谱影像 纹理影像 图2地质统计学纹理计算示意图 Fig 2 Caleul ation of geostatisti cal texture 块金值 变差距离 2.4荒漠地纹理特征 步长h 由于荒漠化土地的各个等级和海滩沙地在各 图4变异函数图参数意义 波段的光谱值都十分接近，因此在训练样区选择 Fig 4 Parameter meaning of variogram 时极容易造成误判，也因此造成在正式分类时很 难将它们区分开.事实上，它们在空间结构上存 从图3中可以看出，海滩沙地的变异函数曲 在很大的差异，显然海滩的地表结构与内陆的荒 线非常平缓，当步长距离大于1时呈现小至于无 漠化土地的地表结构明显不同，而不同等级的荒 的空间自相关而极重度荒漠地变异函数在步长 漠化土地也有一定的结构差异，如极重度荒漠化 大于2时才变得平缓，其基台值也比海滩沙地更 土地沙丘轮廓明显，峰谷显著，而对于程度稍轻的 大，说明极重度荒漠地比海滩沙地有着更大的变 重度荒漠化土地，则显示为脊线轮廓明显9，为 异性，轻度荒漠地显示最大的变异性，变异函数曲

γjk ( h) = 1 2P ( h) ∑ P( h) α=1 [ DN j( xα) -DN k ( xα+h] ( 5) 2.3 变异函数参数 地质统计学影像纹理是通过计算移动窗口内 的相邻像元的半变异值得到的, 计算方法如图 2 所示 .在计算的过程中需考虑计算窗口大小、像 元取样距离(步长) 、计算方向和波段组合等几个 因子 [ 8] .计算窗口(移动视窗) 大小, 主要是依据 在计算每个像元的纹理值时, 需要同时考虑到该 像元周围多大的范围来决定, 一般是在反复实验 的基础上得到的, 一个理想的窗口不能太大, 否则 将使得相邻的土地覆盖类型互相干扰, 同时也不 能太小, 否则将减小变异函数分类器的使用范围. 像元取样距离( 步长) 依据窗口大小的改变而改 变, 它代表了像元与像元之间的关系, 图 2 中窗口 大小为 3, 步长可为 1 或者 2, 一般情况下, 计算结 果取所有步长的平均值.考虑计算方向是因为地 表景物的排列顺序与方向并不固定, 因此在计算 纹理特征时, 需考虑计算东-西向( 0°) 、南-北向 (90°) 、东北-西南向( 45°), 西北-东南向( 135°) 4 个方向的纹理, 当空间方向性不明显时, 可将 4 个 方向的变异函数值加总取平均, 即成为所谓等向 变异元.波段组合是针对交叉变异函数和虚拟变 异函数的一个参数, 选择不同的波段组合, 将对地 表景物的区分产生重要影响. 图 2 地质统计学纹理计算示意图 Fig.2 Calculation of geostatisti cal texture 2.4 荒漠地纹理特征 由于荒漠化土地的各个等级和海滩沙地在各 波段的光谱值都十分接近, 因此在训练样区选择 时极容易造成误判, 也因此造成在正式分类时很 难将它们区分开 .事实上, 它们在空间结构上存 在很大的差异, 显然海滩的地表结构与内陆的荒 漠化土地的地表结构明显不同, 而不同等级的荒 漠化土地也有一定的结构差异, 如极重度荒漠化 土地沙丘轮廓明显, 峰谷显著, 而对于程度稍轻的 重度荒漠化土地, 则显示为脊线轮廓明显 [ 8] .为 了证明变异函数区分几种不同沙地类型的能力, 我们在 1999 年 4, 3, 2 波段合成影像上为第二波 段选择了 3 个类别的小样本, 分别为海滩沙地(近 海亮白色区域) 、极重度荒漠地(内陆亮白色区域) 和轻度荒漠地( 亮白色周边的红白色区域) , 利用 方向变异函数公式进行计算, 得到的变异函数图 如图 3 所示.在研究变异函数图时, 一般关注三 个参数 [ 9] , 即 变差距离 、基台值 和块金值 ( 如 图 4), 从变异函数的定义可以看出, 当一对样本 距离增加时, 所对应的变异函数值一般也增加, 当 变异函数达到一个平台时的样本间距离叫做变差 距离或相关尺度.变差距离表明当样本间距等于 或大于此距离时, 样本就变得完全独立, 在变差距 离变异函数所达到的平台叫做基台值.块金值是 在极短的样本距离( h ≈0) 之间变异函数从原点 的跳升值( 不连续性), 它是由于样本误差和短距 离的变异性引起的. 图 3 不同类别变异函数图 Fig.3 Variograms of different classes 图 4 变异函数图参数意义 Fig.4 Parameter meaning of variogram 从图 3 中可以看出, 海滩沙地的变异函数曲 线非常平缓, 当步长距离大于 1 时呈现小至于无 的空间自相关, 而极重度荒漠地变异函数在步长 大于 2 时才变得平缓, 其基台值也比海滩沙地更 大, 说明极重度荒漠地比海滩沙地有着更大的变 异性, 轻度荒漠地显示最大的变异性, 变异函数曲 Vol.28 No.8 廖楚江等:基于地质统计学影像纹理的海南矿区荒漠化监测 · 711 ·

。712 北京科技大学学报 2006年第8期 线在步长大于13个像素时才趋于平稳. 为逻辑与运算符.这一过程可以在ENVI软件中 完成，结果将得到包括海滩和内陆荒漠的图像从 3 结果与分析 而将研究的重点集中在如何区分海滩沙地和内陆 为了将注意力集中在沙地（海滩沙地和内陆 荒漠地以及荒漠化土地分级上来 荒漠化土地)研究上，首先我们对非沙地类型进行 3.1研究变量选择 掩膜处理，在非沙地类型中，只有废弃耕地和矿山 为了突出几种不同类型沙地的光谱和形态特 等少数地物的光谱反射特性在可见光波度可能与 征，并减少用于纹理分析的变量，需要首先对TM 沙地具有相似性，但可以利用红外波段(7波段) 6个波段进行主成分分析. 的反射值及植被指数值加以区分7,1g.在处理过 由于主成分的选择直接牵涉到后续纹理计算 程中分三次查询运算和一次逻辑运算来实现。具 的效果，以交叉变异函数17×17窗口为例，取海 体计算过程如下： 滩沙地，内陆荒漠地和轻度荒漠地三种类别的半 非沙地(1)=C<RGB彩色合成图像中的光 变异函数值进行分析，对比了两种主成分分析方 谱反射特性值C, 法的优劣：一种是直接对6个波段求取第一和第 非沙地(2)=R,<红外波段图像中的光谱反 二主成分：另一种是特征主成分分析方法1！.该 射特性值<R, 方法将6个波段分成两组：第一组为可见光波段 非沙地(3)=N<改进型土壤调节植被指数 (TM1,2,3):第二组为红外波段(TM4,5,7),随 (MSAVIK<N 后，求出每一组最有代表性的主成分. 非沙地=非沙地(1)and非沙地(2)and非沙 从图5可以看出，使用特征主成分分析方法 地(3). 两种类型地物特征差异更为明显，因此对两幅图 其中，C,C,R,R,N,和N1分别代表不同图像 像采用特征主成分分析，以可见光波段组的第一 中非沙地光谱反射特性值的最小与最大值；and 主成分（记为P℃1）和红外波段组的第一主成分 (a)特征主成分交叉变异函数 150[(间六波段主成分交叉变异函数 40 严重荒漠地 严重荒漠地 海滩沙地 0 海滩沙地。一·一、 潜在荒漠地 5 9 13 13 9 公 步长 步长 图5两种主成分交叉变异元分析方法对比 Fig.5 Comparison of two principle component analyses (记为PC2)作为纹理计算的基础. 影像，而(b)和(d则分别为两个时期5X5窗口的 3.2纹理波段求取 PC1和PC2的方向变异函数纹理波段和二者的 对三个不同时期的影像.分别采用3×3.5× 交叉变异函数纹理波段合成的影像.从图中可以 5和77窗口来计算PC1和PC2的变异函数纹 看到：()中的内陆中心的亮白色与沿海的亮白色 理波段，计算的单变量和多变量变异函数包括方 无法分离，而对应的(b)中将它们分离成了三种不 向变异函数()、绝对变异函数(M)、虚拟变异函 同的灰阶，很好地实现了海滩沙地与荒漠地的分 数(PV)和交叉变异函数(CV),结果将采用步长 离：（©中不同的荒漠化程度很难从这些几乎完全 为1个像素(30m)在4个方向(NS,E-W,N45E 相同的亮白色中分离出来，而用纹理波段合成的 和N45W)的平均值. (d山则将它们分离成了四种不同的灰阶，分别对应 从每幅影像我们将计算得到18个纹理波段， 不同程度的荒漠化 即三种不同窗口下PC1和PC2的各两个单变量 3.3精度分析 (V和M)变异函数纹理波段，以及它们的两个多 采用最大似然分类法进行分类.由于强调的 变量变异函数(CV和PV)变异函数纹理波段.图 重点不同，1986年影像拟分为三类即极重度荒 6中列出的4幅图中，(a)和(c)分别为截取的 漠化、重度荒漠化和其他（包括海滩沙地）；而 1986年和1999年TM4.3.2波段合成的假彩色 1999年和2001年影像拟分为5类，即极重度荒

线在步长大于 13 个像素时才趋于平稳 . 3 结果与分析 为了将注意力集中在沙地(海滩沙地和内陆 荒漠化土地)研究上, 首先我们对非沙地类型进行 掩膜处理, 在非沙地类型中, 只有废弃耕地和矿山 等少数地物的光谱反射特性在可见光波度可能与 沙地具有相似性, 但可以利用红外波段( 7 波段) 的反射值及植被指数值加以区分[ 7, 10] .在处理过 程中分三次查询运算和一次逻辑运算来实现, 具 体计算过程如下 : 非沙地( 1) =Cs <RGB 彩色合成图像中的光 谱反射特性值<Cl , 非沙地( 2) =R s <红外波段图像中的光谱反 射特性值 <R l, 非沙地( 3) =N s <改进型土壤调节植被指数 (MSAVI) <N l, 非沙地=非沙地( 1) and 非沙地( 2) and 非沙 地( 3) . 其中, Cs, Cl, R s, R l, Ns 和 N l 分别代表不同图像 中非沙地光谱反射特性值的最小与最大值;and 为逻辑与运算符 .这一过程可以在 ENVI 软件中 完成, 结果将得到包括海滩和内陆荒漠的图像, 从 而将研究的重点集中在如何区分海滩沙地和内陆 荒漠地以及荒漠化土地分级上来. 3.1 研究变量选择 为了突出几种不同类型沙地的光谱和形态特 征, 并减少用于纹理分析的变量, 需要首先对 TM 6 个波段进行主成分分析. 由于主成分的选择直接牵涉到后续纹理计算 的效果, 以交叉变异函数 17 ×17 窗口为例, 取海 滩沙地, 内陆荒漠地和轻度荒漠地三种类别的半 变异函数值进行分析, 对比了两种主成分分析方 法的优劣 :一种是直接对 6 个波段求取第一和第 二主成分;另一种是特征主成分分析方法[ 11] .该 方法将 6 个波段分成两组 :第一组为可见光波段 ( TM1, 2, 3) ;第二组为红外波段( TM4, 5, 7), 随 后, 求出每一组最有代表性的主成分. 从图 5 可以看出, 使用特征主成分分析方法, 两种类型地物特征差异更为明显, 因此对两幅图 像采用特征主成分分析, 以可见光波段组的第一 主成分( 记为PC1) 和红外波段组的第一主成分 图 5 两种主成分交叉变异元分析方法对比 Fig.5 Comparison of two principle component analyses (记为 PC2) 作为纹理计算的基础. 3.2 纹理波段求取 对三个不同时期的影像, 分别采用 3 ×3, 5 × 5 和 7 ×7 窗口来计算 PC1 和 PC2 的变异函数纹 理波段, 计算的单变量和多变量变异函数包括方 向变异函数( V ) 、绝对变异函数( M) 、虚拟变异函 数( PV) 和交叉变异函数( CV) , 结果将采用步长 为 1 个像素( 30 m) 在 4 个方向( N-S, E-W, N45E 和 N45W) 的平均值. 从每幅影像我们将计算得到 18 个纹理波段, 即三种不同窗口下 PC1 和 PC2 的各两个单变量 ( V 和 M) 变异函数纹理波段, 以及它们的两个多 变量变异函数( CV 和 PV) 变异函数纹理波段 .图 6 中列出的 4 幅图中, ( a ) 和( c) 分别为截取的 1986 年和 1999 年 TM 4, 3, 2 波段合成的假彩色 影像, 而( b) 和( d) 则分别为两个时期 5 ×5 窗口的 PC1 和 PC2 的方向变异函数纹理波段和二者的 交叉变异函数纹理波段合成的影像 .从图中可以 看到 :( a) 中的内陆中心的亮白色与沿海的亮白色 无法分离, 而对应的( b) 中将它们分离成了三种不 同的灰阶, 很好地实现了海滩沙地与荒漠地的分 离;( c) 中不同的荒漠化程度很难从这些几乎完全 相同的亮白色中分离出来, 而用纹理波段合成的 ( d) 则将它们分离成了四种不同的灰阶, 分别对应 不同程度的荒漠化. 3.3 精度分析 采用最大似然分类法进行分类 .由于强调的 重点不同, 1986 年影像拟分为三类, 即极重度荒 漠化、重度荒漠化和其他( 包括海滩沙地) ;而 1999 年和 2001 年影像拟分为 5 类, 即极重度荒 · 712 · 北 京 科 技 大 学 学 报 2006 年第 8 期

Vol.28 No.8 廖楚江等：基于地质统计学影像纹理的海南矿区荒漠化监测 。713。 (d) 图6单纯光谱波段与纹理波段合成对比(a)1986年第4波段：(b)1986年P2纹理：(91999年第4波段（d199年PC2纹理 Fig.6 Comparison between combination of spectral bands and texture bands:(a)Band 4(1986);(b)texture of PC2(1982):(c) Band 4(1999):(d)texture of PC2 (1999) 漠化、重度荒漠化、中度荒漠化、轻度荒漠化和其 阵，矩阵中精度为生产者精度（即某一类别中被正 他（包括海滩沙地）.以土地利用图为参照，对 确分类的样本数与属于该类别的检验样本总数的 1986年影像共选取398个样本，其中训练样本 比).我们注意到：尽管7×7窗口在极重度荒漠 279个，检验样本119个；1999年和2001年影像 化和重度荒漠化的分类精度比5×5窗口要高，但 各取521个样本，其中训练样本364个，检验样本 是对于海滩沙地的分类精度却偏低.对于1986 157个.训练样本的质量采用样本分离度进行测 年影像，重点是海滩沙地的剥离，因此选择5×5 度.表1和表2显示的是两幅影像最终选取的样 窗口作为该幅影像的纹理计算窗口：对于1999年 表1基于纹理波段的影像训练样本分离度矩阵表 影像，由于强调的是荒漠地的分级，显然其中7X (1986) 7窗口的各项分类精度要高（见表4）.最终两幅 Table 1 Trans formed divergence matrix for training sam 影像的总体分类精度（检验样本中被正确分类的 ples 样本数与总样本数的比)分别为89.25%和 类别 极重度 重度 其他 9240%,分别比单纯采用光谱波段分类提高了 极重度 0 1880 2000 7.1%和8.3%. 重度 1880 0 2000 表3分类精度矩阵(1986) 其他 2000 2000 0 Table 3 Matrix for dassification accuracy (1986) 视窗 极重度 重度 其他 表2 基于纹理波段的影像训练样本分离度矩阵表 3X3 88.17 87.32 91.75 (1999) 5X5 88.86 8824 9473 Table 2 Transformed divergence matrix for training sam- 7X7 89.73 9086 9269 ples(1999) 类别 极重度 重度 中度 轻度 其他 表4分类精度矩阵(1999) 极重度 0 1770 1980 2000 2000 Table 4 Matrix for dassification accuracy (1999) 重度 1770 0 1830 1800 2000 视窗 极重度 重度 中度 轻度其他 中度 1980 1830 0 1887 2000 3X3 85.48 85.06 87.19 8820 85.31 轻度 2000 1800 1887 0 2000 5X5 87.33 8840 89.15 91.94 9448 其他2000 2000 2000 2000 0 7X7 89758873 9023 925187.62 本分离度矩阵（分离度取值范围是0~2000），其 4 中表1是1986年基于5×5窗口的纹理选取的训 讨论 练样本分离度矩阵，表2是1999年基于7×7窗 图7显示的是根据三幅影像得到的分类结果 口的纹理选取的训练样本分离度矩阵.可以看 图.可以看出：1986年影像己经将海滩的干扰基 到，在这两种情况下，两幅影像的样本分离度都很 本消除掉了，而1999年和2001年影像的分类结 高（大于1500为较优）.分类精度矩阵如表3和 果则详细阐明了荒漠化的各个级别.根据分类结 表4所示，其中表3为1986年影像的分类精度矩 果得到1986年总荒漠化面积18461hm2,而

图6 单纯光谱波段与纹理波段合成对比.(a) 1986年第 4波段;(b) 1986 年PC2 纹理;( c) 1999 年第4 波段;( d) 1999 年 PC2纹理 Fig.6 Comparison between combination of spectral bands and texture bands:( a) Band 4 ( 1986) ;( b) texture of PC2 ( 1982) ;( c) Band 4 ( 1999) ;( d) texture of PC2 ( 1999) 漠化 、重度荒漠化 、中度荒漠化 、轻度荒漠化和其 他( 包括海滩沙地) .以土地利用图为参照, 对 1986 年影像共选取 398 个样本, 其中训练样本 279 个, 检验样本 119 个;1999 年和 2001 年影像 各取 521 个样本, 其中训练样本 364 个, 检验样本 157 个.训练样本的质量采用样本分离度进行测 度.表 1 和表 2 显示的是两幅影像最终选取的样 表1 基于纹理 波段的 影像训练 样本分 离度矩阵 表 ( 1986) Table 1 Transformed divergence matrix for training sam￾ples 类别 极重度 重度 其他 极重度 0 1 880 2 000 重度 1 880 0 2 000 其他 2 000 2 000 0 表2 基于纹理 波段的 影像训练 样本分 离度矩阵 表 ( 1999) Table 2 Transformed divergence matrix for training sam￾ples ( 1999) 类别 极重度 重度 中度 轻度 其他 极重度 0 1 770 1 980 2 000 2 000 重度 1770 0 1830 1800 2 000 中度 1 980 1 830 0 1 887 2 000 轻度 2 000 1 800 1 887 0 2 000 其他 2 000 2 000 2 000 2 000 0 本分离度矩阵( 分离度取值范围是 0 ～ 2 000) , 其 中表 1 是 1986 年基于 5 ×5 窗口的纹理选取的训 练样本分离度矩阵, 表 2 是 1999 年基于 7 ×7 窗 口的纹理选取的训练样本分离度矩阵.可以看 到, 在这两种情况下, 两幅影像的样本分离度都很 高( 大于 1 500 为较优) .分类精度矩阵如表 3 和 表 4 所示, 其中表 3 为 1986 年影像的分类精度矩 阵, 矩阵中精度为生产者精度( 即某一类别中被正 确分类的样本数与属于该类别的检验样本总数的 比) .我们注意到 :尽管 7 ×7 窗口在极重度荒漠 化和重度荒漠化的分类精度比 5 ×5 窗口要高, 但 是对于海滩沙地的分类精度却偏低.对于 1986 年影像, 重点是海滩沙地的剥离, 因此选择 5 ×5 窗口作为该幅影像的纹理计算窗口 ;对于 1999 年 影像, 由于强调的是荒漠地的分级, 显然其中 7 × 7 窗口的各项分类精度要高( 见表 4) .最终两幅 影像的总体分类精度( 检验样本中被正确分类的 样本数 与总 样本 数的 比) 分 别为 89.25 %和 92.40 %, 分别比单纯采用光谱波段分类提高了 7.1 %和 8.3 %. 表 3 分类精度矩阵( 1986) Table 3 Matrix for classification accuracy ( 1986) 视窗 极重度 重度 其他 3×3 88.17 87.32 91.75 5×5 88.86 88.24 94.73 7×7 89.73 90.86 92.69 表 4 分类精度矩阵( 1999) Table 4 Matrix for classification accuracy ( 1999) 视窗 极重度 重度 中度 轻度 其他 3×3 85.48 85.06 87.19 88.20 85.31 5×5 87.33 88.40 89.15 91.94 94.48 7×7 89.75 88.73 90.23 92.51 87.62 4 讨论 图 7 显示的是根据三幅影像得到的分类结果 图.可以看出:1986 年影像已经将海滩的干扰基 本消除掉了, 而 1999 年和 2001 年影像的分类结 果则详细阐明了荒漠化的各个级别 .根据分类结 果得到 1986 年总荒漠化面积 1 846.1 hm 2 , 而 Vol.28 No.8 廖楚江等:基于地质统计学影像纹理的海南矿区荒漠化监测 · 713 ·