基于遥感技术的南水北调水源区(河南段)石漠化遥感调查与评价

许军强，张斌，袁晶，卢意恺

1.河南省航空物探遥感中心，郑州450053；2.河南省地质科学研究所，郑州450001；3.河南省地矿局第五地质勘查院，郑州450001

0 引言

目前，遥感技术在岩溶石漠化信息提取、空间分析等领域得到了广泛应用，并取得了一系列显著成果。如陈希等[1]利用SAR数据结合光学影像融合对喀斯特地区石漠化等级划分，提取石漠化信息，结果表明，多源遥感数据支持下的石漠化信息提取方法比传统提取方法更加准确和有效。周迪[2]等基于影像纹理和地形数据辅助面向对象方法进行喀斯特地区石漠化信息的提取，可有效地减少因复杂地形导致石漠化信息提取结果“椒盐化”现象，提取精度明显优于基于像素的监督分类和非监督分类方法。夏学齐等[3]提出“植被线”概念，并定义“石漠化几何指数”(GRI)，客观稳定地提取石漠化程度。涂杰楠等[4-5]通过归一化植被指数与植被覆盖度进行岩溶石漠化遥感信息提取，得到两期石漠化空间分布特征，并对两期植被覆盖度数据进行叠加，分析生态治理区石漠化的演变趋势。刘芳等[6]基于通过对喀斯特地区岩石热特性探究，提出基于热红外遥感的石漠化监测方法，进行石漠化等级划分并分析石漠化变化。甘敬民等[7]运用多因子综合分析、专家打分和层次分析等方法，建立多因素石漠化强度等级评价及石漠化演化敏感性分析地理模型，得出石漠化等级及石漠化演化敏感性空间分布结果。熊康宁等[8]基于GIS/RS调查了贵州省岩溶石漠化的分布现状和特征。陈飞等[9]利用CA-Markov预测方法，探寻研究区1990—2016年并预测了2016—2021年石漠化时空演变规律。郭丽琴等[10]选取TM和高分一号、资源三号和遥感二号2期卫星遥感数据为数据源，采用RS和GIS技术，获取了石漠化面积变化情况，探索石漠化时空演变特征及成因分析。李松等[11]采用遥感变化检测的方法对TM/ETM+数据进行石漠化信息提取和精度检验，取得了较好的效果。李建存等[12]利用相隔10年的3期影像数据，进行碳酸盐岩遥感解译和信息提取及野外核查，生成碳酸盐岩分布图、石漠化分级分布图以及石漠化趋势图等，为石漠化的防治提供依据。

1 生态环境背景

综合丹江口水库(河南辖区)汇水区域和碳酸盐岩的分布选取南水北调水源区(河南段)的范围。南水北调水源区(河南段)主要位于河南省南阳市西南部，南阳盆地西缘，包括淅川南部、内乡和邓州西南部一隅。东邻内乡县、邓州市，北依西峡县，西与陕西省交界，南与湖北省相接，地处32°35′～33°08′N， 111°08′～111°44′E之间，总面积约2 584 km2(图1)。研究区属北亚热带向暖温带过渡的季风性气候区，年平均气温15.8℃。年降雨量在391.3～1 423.7 mm之间，多年平均降雨量为804.3 mm，多年月平均降雨量68.2 mm，区内雨量充沛。

图1 研究区交通位置图Fig.1 Traffic location map of research area

研究区地形为西北高、东南低，呈马蹄状。西北部地貌为侵蚀中山与低山区，中部为丘陵区，东南部为岗地与河流阶地缓冲区。区内出露地层从老到新依次是震旦系、寒武系、奥陶系及石炭系，碳酸盐岩分布广泛，有灰岩、白云质灰岩、白云岩及大理岩等。广泛分布的各类碳酸盐岩为岩溶作用提供了物质条件，成为石漠化的主要物质基础。

2 研究方法

2.1 波谱特征分析

图2 碳酸盐岩地区主要地物波谱曲线Fig.2 Spectrum curves of main objects in area composed of limestone

借鉴碳酸盐岩地区主要地物波谱[13]，碳酸盐岩地区主要地物类型为3种：①植被；②碳酸盐岩；③土壤(图2)。从图中可以看出，植被的反射率在0.48 μm附近有一个18%的反射波峰，在0.55 μm附近有一个反射率为10%～15%的吸收波谷，即在TM3波段中红光波段为植被吸收波段，在0.76～0.82 μm区间反射率陡增，在0.82～1.3 μm区间植被反射率相对稳定在40%～50%的高反射区间，在1.5 μm附近植被反射率陡降到17%的吸收谷，在1.60 μm升到32%的次反射峰，说明TM4波段中的近红外波段(0.76～0.90 μm)是植被强反射区间，可以区分植被类型和覆盖率的差异。碳酸盐岩波谱曲线相对平缓，在各个波段的反射率无明显差异，TM1、TM2、TM3、TM4、TM5波段在0.40～1.75 μm的反射率基本稳定在20%±，其中1.55～1.75 μm是绝大部分造岩矿物波谱响应的强反射波谱，碳酸盐岩在TM5波段具有较强的识别能力，从1.75～2.35 μm反射率缓慢下降至17%±，TM7波段对造岩矿物的识别能力较弱。土壤的波谱曲线在0.40～2.10 μm区间，反射率整体呈上升趋势，但波动较大，在0.40～0.80 μm区间，土壤反射率呈线性逐渐上升到42%，反射率达到一个波峰，0.80～0.90 μm区间逐渐下降到36%后又持续上升至46%，在TM5波段反射率达到最大53%。在TM7波段升降波动明显。

石漠化在遥感图像上的光谱响应特征可以通过波段DN值计算获得，选取良好的石漠化训练样本对表征不同石漠化侵蚀程度的光谱特征以及石漠化分类能够达到良好的效果。为保证样本选取的准确性与可靠性，以野外踏勘选取的石漠化样本点为主，投影到影像上选取典型的潜在石漠化、轻度石漠化、中度石漠化和重度石漠化的训练样区。基于样本区石漠化光谱DN值，生成不同程度的石漠化光谱特征曲线(图3)。

图3 石漠化程度光谱特征Fig.3 Spectral characteristics of rocky desertification

图3中，从均值的可分性看，重度石漠化均值最大，反射率最高，且与中度石漠化、轻度石漠化、潜在石漠化的反射率有显著差异，中度石漠化均值其次，与轻度石漠化、潜在石漠化的反射率有明显差异，轻度石漠化与潜在石漠化的反射率接近，可分性最弱。从波段特征上看，第4波段对重度石漠化和中度石漠化有较好的区分，第1波段对潜在石漠化识别较好，第2、3波段对中度和轻度石漠化有较好的分区，第5波段对重度石漠化区分较好。

2.2 石漠化信息提取方法

以TM/ETM+、OLI数据为主要数据源，对无石漠化和非石漠化做掩膜处理后的石漠化范围进行石漠化程度信息提取，在不同程度石漠化训练样本的基础上，基于改进增强型植被指数[14]和监督分类的方法提取石漠化侵蚀程度信息。本文采用最大似然法、支持向量机和神经网络分类法，通过改进增强型植被指数与数据源的合并，参与分类，实现石漠化信息提取。

(1)最大似然分类

像元数据x，其对于给定的类别k的似然度Lk被定义为x出现时，其归属于类别k的条件概率的计算公式如下：

Lk=P(k|x)=P(k)×P(x|k)/∑P(i)×P(x/i)

(1)

式中：x为像元的观测数据，它是由待分类像元n个特征变量的观测值所组成的列向量。

对于同一分类类别和相同的训练样本，上式是个常量，可以简化为：

Lk=P(k|x)=P(k)×P(x|k)

(2)

式中：先验概率P(k)已知，求像元x与类别k的似然度，P(x|k)是唯一的未知变量。要计算x与类别k的似然度，必须先给出总体的概率密度函数P(x|k)。对于未知类别X，它属于类别A的概率为Pa，它属于类别B的概率Pb，且Pa>Pb，所以未知类别X属于类别A。

(2)支持向量机分类

支持向量机是以结构风险最小化为原则，在线性分类的基础上提出并应用到线性分类器中。

对非线性可分的训练样本点，用一个非线性函数ψ将训练样本点空间映射到一个高维的特征空间并进行线性分类。若可以找到映射ψ，则内积运算(xi·x)可用(ψ(xi) ·ψ(x))代替。用函数K(xi，x)=(ψ(xi)·ψ(x))表示内积运算而不直接利用函数ψ，称函数为核函数。

最优分类超平面可表示为：

(3)

在ENVI中的SVM经常使用的核函数为径向基核函数：

(4)

在ENVI软件下，选择菜单栏中的Classification模块中的support Vector Machine，核函数类型选择径向基核函数。

(3)神经网络分类

神经网络是模仿人的大脑神经元行为活动，具有分布式并行信息处理特征的一种数学模型。它是一个复杂而又庞大的系统，系统内存在大量的像人脑神经元一样的节点，节点之间相互连接，数据信息在神经网络中不断的输入和输出，从而达到数据处理的目的。

BP人工神经网络算法的分类步骤描述如下：

① 对权值Wij置初值。对每一层的权值Wij的初始值设为较小的随机数，其中Wi，n+1=-θ(θ为阈值)；

② 给定训练样本集，输入训练样本和相应的目标输出，即输入向量X和期望输出Y；

(5)

⑤ 更新权值Wij和阈值θ。其计算表达式为：

(6)

式中：η为学习速率，一般取值为0.2～0.6。

当使用动量时，更新权值Wij数学表达式可采用如下：

(7)

=Wij(t)-Wij(t-1)

(8)

式中：α为动量系数，0<α<1。

⑥ 在得出各层各个权系数之后，判断误差函数E是否收敛所给出的学习精度ε(E≤ε)，如果满足学习精度要求，则算法结束；如果未满足要求，则返回③继续执行。选择ENVI4.8 Classification模块下的Neural Net做神经网络分类，分类结果如图4所示。

a.最大似然分类图;b.支持向量机分类;c.神经网络分类图。图4 石漠化分类图Fig.4 Rocky desertification classification

2.3 分类精度评价

评价石漠化强度分级结果的分类精度时，通过建立的解译标志选取分类训练样本，把选取的样本随机分为两部分，其中70%用于分类，30%用于精度评价，提高精度评价的可靠性，分类精度和Kappa系数如表1所示。

表1 分类精度评价

3 石漠化遥感监测与演变分析

利用ArcGIS10.1空间叠加分析技术对两期石漠化程度进行叠加分析处理，分别以1996、2006和2006年、2016年石漠化程度信息完成两期的石漠化变化趋势信息提取。基本了解1996—2006年和2006—2016年两个动态时期内的石漠化侵蚀程度变化(表2、图5)，掌握各石漠化程度的稳定、改善或恶化情况(图6)。从表2、图5可以看出，1996—2006年的变化趋势是以稳定状态为主，局部石漠化程度处于趋好的状态，轻微改善面积为209.47 km2，明显改善面积22.63 km2；2006—2016年的变化趋势依然以稳定状态为主，生态地质状况总体良好，经野外调查发现，部分严重石漠化区域经治理已有明显改善情况，但局部石漠化依然有恶化现象。

表2 1996—2006年、2006—2016年二期石漠化面积变化趋势统计表

Table 2 Statistic table of trend of rocky desertification area in 1996—2006 and 2006—2016

类型1996—2006年/km22006—2016年/km2严重加剧22.7375.15加剧48.27144.80稳定1 323.621 189.29轻微改善209.47167.40明显改善22.6330.00

a.1996—2006年；b.2006—2016年。图5 石漠化变化趋势柱状图Fig.5 Histogram of rocky desertification change trend from 1996—2006 and 2006—2016

a.1996—2006年；b. 2006—2016年。图6 石漠化面积变化趋势图Fig.6 Chart of change trend of rocky desertification area

4 讨论

在GIS平台下，首先分别基于地层岩性、地形地貌、工程地质和水文地质等因素不同类别的确定性系数CF函数值来构建影响因素的栅格图层，然后利用GIS平台下空间分析模块中的栅格计算式对各统计单元进行因子叠加分析。通过定义石漠化评价指数(stony desertification index，SDI)来进行评价石漠化的易发性。石漠化评价指数公式如下：

SDI=∑RiCFi

(9)

式中：SDI为评价单元的石漠化易发性指数；CFi为该单元各影响因素的CF值；Ri为根据层次分析法得到的各影响因素的权重。SDI也是一个处于区间[-1，1]的指标，SDI值越高，说明石漠化易发性越高，反之，SDI值越低，石漠化易发性越低。

分级方法采用GIS空间分析的栅格重分类工具，按照自然间距分类(natural break)方法，将石漠化易发性结果栅格图分为4级，即石漠化高易发区(0.07～0.63)、石漠化易发区(-0.28～0.07)、石漠化低易发区(-0.64～-0.28)和石漠化不易发区(-0.86～-0.64)(图7)。通过对研究区石漠化易发性综合分区评价，石漠化易发性共分为高易发区、易发区、低易发区和不易发区4个分区，4个分区的空间分布和地层走向大体一致，说明地层岩性在石漠化易发性综合分区评价中的影响最高，其权重系数在影响因子中最高也得以证明。

图7 石漠化易发性综合分区图Fig.7 Comprehensive regionalization map of rocky desertification susceptibility

高易发区主要分布在研究区的中部和北部，中部的高易发区呈片状特征，面积较大，此区域以奥陶系地层为主，岩性主要为白云质灰岩和灰岩等，北部的高易发区呈条带状特征，面积较小，此易发区域以寒武系、奥陶系地层为主，岩性主要为灰岩、白云质灰岩和含泥质条带灰岩等，高易发区的整体分布与该区域的岩性具有高石漠化发生率相吻合。易发区面积较大，主要集中在研究区北部，以石炭系地层为主，岩性主要为白云质灰岩，同样具有较高的石漠化发生率。低易发区面积较小，主要集中在研究区西北部、中部及西南部，分布地层有震旦系、志留系及寒武系地层，主要岩性为白云质大理岩、白云岩、泥灰岩和砂砾岩等，此区域的石漠化发生率较低。不易发区的分布与第四系地层吻合，处于平原地貌，岩性为砂砾石和黏土等，石漠化发生率最低。

5 结论

(1)利用改进增强型植被指数和监督分类相结合的方法提取石漠化侵蚀程度信息。通过采用最大似然法、支持向量机和神经网络分类均能有效提取研究区石漠化程度信息，最大似然法提取石漠化总体精度最高， 达到84.93， Kappa系数为0.792 3。

(2)利用空间叠加分析技术完成对两期石漠化变化趋势监测。分别对1996、2006、2016年石漠化程度信息进行石漠化空间叠加分析。1996—2006年轻微改善面积为209.47 km2，明显改善面积22.63 km2；2006—2016年的变化趋势依然以稳定状态为主，生态地质状况总体良好。

(3)基于石漠化评价指数进行石漠化易发性分区评价。以地层岩性、地形地貌、工程地质和水文地质等因素不同类别的不确定性系数及权重构建评价指数，得出石漠化高易发区、易发区、低易发区和不易发区。4个综合分区指数分别为0.07～0.63、-0.28～0.07、-0.64～-0.28和-0.86～-0.64， 4个分区的空间分布和地层走向大体一致，表明该指数可定量表达石漠化易发性评价。