一种地表生境因子量化解析方法

1.本发明属于数字景观场地勘探领域，具体涉及一种基于高分遥感影像与半 监督递归分类技术的地表生境因子解析方法和设备。


背景技术：

2.生境因子，指生态系统中存在着的具有生态影响力的自然要素或人工要 素，其空间形态存在着面状、线状、点状分布，其空间分布位置关系根据其相 互关系及具体环境呈现出交织混杂的特征。对生境因子的高精度解析不仅对维 护和保护生态系统中局部生态格局、推动区域生态格局持续改善、整体向上发 展有着关键性的作用，也能提高生态环境建设的科学性与生态性，助力实现高 质量的景观环境生态文明建设。
3.当前风景环境建设生境因子划分存在着以下这些问题：
4.(1)生境因子数据参数获取与解析的精度、效率与科学性有待提升，呈 现出小数据分析大空间、单一数据分析复杂空间、综合数据难以处理等情况。
5.运用相关生态技术、遥感技术对风景环境建设生境因子解析以定性研究为 主，量化解析程度有待深入；其生境因子等数据多数依靠人工实地调研收集或 依靠中低分辨率遥感影像(分辨率为30m以下)解析，数据收集过程繁琐，其 精度难在中
‑
微观尺度的生态环境空间分析中有待提升；复合生态参数数据统筹 与综合处理较为困难，多数研究仅涉及部分参数，难以整体、科学地把握生境 因子空间整体属性。
6.(2)生境因子研究多基于现有模型与框架，以宏观尺度为研究对象，运 用固定参数指数进行数据解析，聚焦单一生态敏感参数指数，对多源复合生境 因子数据统筹解析、综合分析的整合性与系统性研究较少，整体解析系统运行 效率、精度及全局性有待提升。
7.综上可知，目前生境因子研究层面常以区级、市级、省级宏观尺度居多， 其向中
‑
微观尺度进一步下探分析尺度与解析精度能力有待加强；研究对象多为 某一生境因子或生态参数；应对复合生境因子解析时参数选呈现偏重选取典型 参数分析，在应对各类中
‑
微观尺度不同的生态环境时缺乏对其特殊性的考量； 在研究对象或单位主要以基于遥感影像的选点为主，在像元层面进行生态空间 解析；分析水平与解析精度难以支撑复杂的生态环境状况与变化，深入指导中
‑ꢀ
微观层面的生境因子解析研究。
8.在指导生态建设实践中，现有分析常为运用精度为30m以下的中低分辨率 遥感影像，针对的分析尺度多为区域及其以上空间尺度，在中
‑
微观实际环境开 发建设项目中，这类中低精度分析存在着生境因子混杂、模糊等情况，较难指 导实际建设；生态数据获取需人工调研标注，人工勘测成本高，数据更新速度 难以应对高速建设的变化与需求；基于该尺度下的研究对象数据单元多为网格 点状分布的像元，落实在实际空间中即为抽象的数据点，对生态空间内蕴含的 生境因子及其对象缺少具体落点；解析过程多以现有模型进行运算处理，在应 对动态变化情况的数据补充缺少途径，其“数据采集
‑
处理
‑
分析”环节整体动 态流程尚待进一步整合发展。
9.

技术实现要素：

10.针对现有技术的不足，本发明提出了一种地表生境因子量化解析方法。
11.本发明的目的可以通过以下技术方案实现：
12.一种地表生境因子量化解析方法，包括以下步骤：
13.步骤1：获取目标地区的地表卫星影像；
14.步骤2：对所述目的地区的地表景观要素进行分类；
15.步骤3：依据基地地域典型生境因子对象类型，测试并选定其分割尺度参数， 以进行生境因子对象分类；
16.步骤4：提取每一个所述生境因子对象的波段数值，
17.步骤5：根据步骤3所得的生境因子参数，构建数据集，并获得高置信度样 本库；
18.步骤6：运用最大似然法模型对高置信度数据样本进行递归最优特征选取。
19.可选地，所述步骤1中获取所述卫星影像后，对所述卫星影像进行预处理：
20.步骤1.1：对原始影像数据进行辐射定标，以消除传感器本身的误差；
21.步骤1.2：进行大气校正，以消除大气折射与散射对地表景观要素反射的影 像；
22.步骤1.3：根据变换系数矩阵参数，进行正射校正；
23.步骤1.4：将正射校正后的遥感影像进行融合处理，根据实验所需选取合适 的遥感影像进行解析。
24.可选地，所述步骤3中基于多尺度分割函数模型分割所述生境因子对象， 且所述分割函数模型的参数配置为分割尺度、分割层数、颜色权重、形状权重、 平滑度、紧致度中的一种或多种。
25.可选地，所述步骤4中的波段数值为red波段数值、green波段数值、blue 波段数值、nir波段数值、bright波段数值、ndvi波段数值、ndbi波段数值、 mndwi波段数值、glcm波段数值中的一种或多种。
26.可选地，所述步骤5包括以下步骤：
27.获取所述步骤4中各个波段数值的最大值、最小值、平均值与标准差；
28.根据已知生境因子类型选取机器学习多维数据样本库；
29.运用主成分分析法对生境因子机器学习多维数据样本库数据库数据降维处 理；在多次训练迭代后，获取每一类生境因子高影响参数，构建高置信度样本库；
30.可选地，所述景观要素分为植被要素、水体要素、土壤要素、建成要素中的 一种或多种。
31.1本发明采用多源高分辨率遥感影像融合进行生境因子的提取，实现了各类 立地条件下，高效率、高维度、高刷新率的生境因子数据采集。
32.2本发明利用面向对象的智能分类技术，将生境因子解析单元由以往的点 (像元)层面提升至具有生态参数属性的生态对象层面，降低了以往解析中存在 的椒盐现象，以生境因子对象为数据集合，实现了生境因子对象的多维度数据整 体分析。
33.3本发明通过主成分分析降维减少了生境因子解析低关联度参数，构建了生 境因子降维参数数据库，提高了生境样地影像解析速度与效率。
34.4本发明提出一种生境因子高置信度样本数据构建方法，解析了多种生境因 子的光谱参数，通过主成分分析降维减少了低关联度参数，构建生境因子高置信 度参数数据库，降低了遥感影像多维参数的干扰，提高了高分辨率遥感影像解析 速度与效率，使得光谱相似的生境因子能够区分。
35.5本发明提出了一种极大似然法递归特征半监督分类学习分析方法，结合高 置信度样本库学习，实现了较少样本标记下生境因子样本自适应泛化解析，实现 了复合参数条件下多种生境因子的高效率、高准确度解析。
附图说明
36.下面结合附图对本发明作进一步的说明。
37.图1为本技术的一些示例的目标地区的多光谱影像；
38.图2为本技术的一个示例的目标地区的全色影像；
39.图3为本技术的一个示例的目标地区的正射影像；
40.图4为本技术的一个示例的目标地区的融合影像；
41.图5为本技术的一个示例的目的地区及其周边遥感影像图；
42.图6为本技术的一个示例的目的地区的地表景观要素分割层；
43.图7为本技术的一个示例的目的地区的地表景观要素分割层；
44.图8为本技术的一个示例的目的地区的地表景观要素分割层；
45.图9为本技术的一个示例的目的地区的地表景观要素ndvi解译图；
46.图10为本技术一个示例的目的地区的地表景观要素ndbi解译图；
47.图11为本技术一个示例的地表景观要素mndwi解译图；
48.图12为本技术一个示例的地表景观要素glcm解译图。
49.图13为本技术一个示例的目标地区的遥感影像地表景观要素的土地覆盖光 谱信息；
50.图14为本技术一个示例的地表生境因子pca降维递归分类流程图；
51.图15为本技术的一个示例中的地表生境因子分类解析结果图。
具体实施方式
52.下面结合说明书附图对本发明的具体技术方案作进一步地描述。
53.采用本发明的方法进行生境因子解析包括七个步骤，分别为多源遥感影像 预处理、生境因子对象分类、生境因子对象多层多尺度参数调节聚类、生境因子 对象库参数解析、生境因子对象库参数建立、生境因子解析参数pca降维、高 置信度数据库样本迭代分类。
54.1)多源遥感影像预处理
55.多源遥感影像预处理步骤，采用本发明的方法进行处理的原始高分辨率分 遥感影像其来源包括卫星遥感影像或无人机遥感影像。本发明采用高分辨率卫 星多光谱遥感影像，其空间分辨率达3.2米，含有四个波段，分别是蓝波段(0.45
‑ꢀ
0.52μm),绿波段(0.52
‑
0.59μm),红波段(0.63
‑
0.69μm),近红外波段(0.77
‑ꢀ
0.89μm)。本发明的一些具体的示例中，处理的遥感影像如图1、图2所示。
56.首先，对原始影像数据进行辐射定标，以消除传感器本身的误差。辐射定标 通过将遥感影像的原始像元数字量化值(digital number，dn)进行大气反射率 转换，以获取传感器所采集的地表景观要素真实反射数据。依照实际需求可分为 绝对定标和相对定标，本例采用绝对定标，其转换公式如式下：
57.lλ＝gain*dn+offset
58.式中：lλ为辐射亮度值，gain为各波段的增益,offset为各波段的偏置,增 益、偏移信息于遥感影像的头文件中获取。
59.其次，进行大气校正，以消除大气折射、散射等现象对地表景观要素反射的 影像。本发明的示例中采用flaash大气校正模型与改进的modtran辐射传 输算法，利用遥感影像像素光谱上的特征取代影像成像时的实测大气参数来对 成像区域大气属性进行反演，消除气溶胶和水蒸气的散射效应，得到精确度较高 的地表反射率，在后续解析步骤中可达到良好的普适性与较高的精度。经过 flaash大气校正模型处理后获得具有地表实际反射率的基地遥感影像，如图3 所示。
60.同时，根据rpc(rational polynomial coefficients，变换系数矩阵)参数， 进行正射校正，投影为wgs84/utm。遥感正射校正通过选取地面控制点，结合 影像范围内的数字高程模型(dem)数据，进行投影差改正和倾斜改正，修正遥 感影像因地形变化引起的像素点偏移以将遥感影像重采样成正射影像，如图4所 示。
61.最后，将正射校正后的遥感影像进行融合处理，根据实验所需选取合适的遥 感影像进行解析。本示例中采用全色波段(2.1m)与多光谱波段(5.8m)遥感 影像，为了提高遥感影像的精度与要素识别能力，増强遥感影像数据的可用性， 对全色波段(2.1m)与多光谱波段(5.8m)遥感影像进行影像融合。运用his变 换融合能够提高影像的空间分辨率，且较好地保留地表景观要素的光谱特性，提 高遥感影像的精度和解译性。经his变换融合后的zy
‑
3多光谱影像空间分辨 率为2.1m，波段分别为蓝波段(blue)、绿波段(green)、红波段(red)和近 红外波段(nir)。为了提高地表景观要素解译速度，减少冗余计算，将融合后的 影像截取基地及其周边部分进行分类处理，如图5所示。
62.2)基地生境因子对象分类
63.根据实地调研情况结合遥感影像解译技术特征与地表景观要素遥感影像图 像表现，首先将基地松岭头目标地区的地表景观要素分为植被要素、水体要素、 土壤要素、建成要素四大类，在此基础上将所有地表生态敏感要素分为乔木、灌 木、地被、滩涂、道路、裸地、水体、茶园、竹林和建筑10类。(表1)
64.表1松岭头地区地表生态敏感要 素分类表
[0065][0066]
3)生境因子对象多层多尺度参数调节聚类
[0067]
依据基地地域典型生境因子对象类型，测试并选定其分割尺度参数，以进行 生境因子对象分类。
[0068]
本案采用基于多尺度分割函数(multiresolution segmentation)模型，运用 分割尺度、分割层数、颜色权重、形状权重、平滑度、紧致度共6个参数完成生 境因子对象分割。分割规则的确定主要考虑根据生境因子对象要素的光谱特点、 几何形状、拓扑特征来对影像进行实用性的分割，并形成一个有a、b、c三级 层次的多尺度遥感影像对象等级网络，各层参数如下。
[0069]
层a的分割参数：尺度为100、颜色权重0.9、形状权重为0.4、平滑度为 0.5、紧致度为0.5；层b的分割参数：尺度为50、颜色权重0.5、形状权重为 0.5、平滑度为0.3、紧致度为0.3；层c的分割参数：尺度为30、颜色权重0.9、 形状权重为0.1、平滑度为0.5、紧致度为0.2。根据多次实验解析，以上参数组 合可有效分割试验场地内生境因子。
[0070]
分割尺度不宜过大，以免造成地物破碎，考虑到松岭头地区水体与非水体交 错程度较高，宜从宏观尺度进行水体与非水体的区分，其分割尺度调整为100； 层b主要区分乔木、灌木、地被、裸地、滩涂，在继承层a分割的基础上，结合 乔木、灌木、地被、裸地、滩涂等具有自然演变特性地的形状与纹理特征，分割 尺度参数为50，将自然要素分开；层c主要区分竹林、茶园、建筑、道路这些 具有人工干预特性地地表景观要素，通过解析其所呈现的光谱、形状、紧致度、 平滑度等特征，并与现有场地基础资料进行比较，有效提取农业植被与建成景观 信息要素，如图6～8所示。表2为各个层次分割的尺度参数及区分度指数。
[0071]
表2：多层分割尺度及区分度参数
[0072][0073]
4)生境因子对象库参数解析
[0074]
对步骤2)获得的生境因子对象参数解析，分别获取其参数。
[0075]
提取每一个生境因子对象的red波段数值。计算公式如下：
[0076][0077]
其中，x为当前计算生境因子对象，i为当前生境因子对象x中一个像素点， n为当前计算生境因子对象x中所有像素点集合，为当前计算子超像素x 中一个像素点i的红波段值。k
x
(red)为当前计算生境因子对象x的red波段数值。
[0078]
提取每一个生境因子对象的green波段数值。计算公式如下：
[0079][0080]
其中，x为当前计算生境因子对象，i为当前生境因子对象x中一个像素点， n为当前计算生境因子对象x中所有像素点集合，为当前计算生境因 子对象x中一个像素点i的红波段值。k
x
(green)为当前计算生境因子对象x的 green波段数值。
[0081]
提取每一个生境因子对象的blue波段数值。计算公式如下：
[0082][0083]
其中，x为当前计算生境因子对象，i为当前生境因子对象x中一个像素点， n为当前计算生境因子对象x中所有像素点集合，为当前计算生境因子 对象x中一个像素点i的红波段值。k
x
(blue)为当前计算生境因子对象x的blue 波段数值。
[0084]
提取每一个生境因子对象的nir波段数值。计算公式如下：
[0085][0086]
其中，x为当前计算生境因子对象，i为当前生境因子对象x中一个像素点， n为当前计算生境因子对象x中所有像素点集合，为当前计算生境因子对 象x中一个像素点i的红波段值。k
x
(nir)为当前计算生境因子对象x的nir波段 数值。
[0087]
提取每一个生境因子对象的bright波段数值。计算公式如下：
[0088][0089]
其中，x为当前计算生境因子对象，i为当前生境因子对象x中一个像素点， n为当前计算生境因子对象x中所有像素点集合，为当前计算生境因 子对象x中一个像素点i的红波段值。k
x
(bright)为当前计算生境因子对象x的 bright波段数值。
[0090]
计算每一个生境因子对象的归一化差异水体指数mndwi，计算公式如下：
[0091][0092]
其中，x为当前计算生境因子对象，i为当前计算生境因子对象x中一个像素 点，n
(x)为当前计算生境因子对象x中所有像素点集合，为当前计算 生境因子对象x中一个像素点i的绿波段值，为当前计算生境因子对象x 中一个像素点i的mir波段值，k
x
(mndwi)为当前计算生境因子对象x的归一化差 异水体指数值，
‑
1<＝k
x
(mndwi)<＝1，正值表示为水体区域，否则为陆地区 域。
[0093]
计算每一个生境因子对象的归一化植被指数ndvi，计算公式如下：
[0094]
所述生境因子对象的归一化植被指数k
x
(ndvi)计算公式如下：
[0095][0096]
其中，x为当前计算生境因子对象，i为当前计算生境因子对象x中一个像素 点，n为当前计算生境因子对象x中所有像素点集合，为当前计算生境因 子对象x中一个像素点i的红波段值，为当前计算生境因子对象x中一个 像素点i的近红外nir波段值，k
x
(ndvi)为当前计算生境因子对象x的归一化植 被指数值，如图9所示。
[0097]
计算每一个生境因子对象的归一化建筑指数ndbi，计算公式如下：
[0098][0099]
其中，x为当前计算生境因子对象，i为当前计算生境因子对象x中一个像素 点，n为当前计算生境因子对象x中所有像素点集合，为当前计算生境因 子对象x中一个像素点i的红波段值，为当前计算生境因子对象x中一个 像素点i的近红外nir波段值，k
x
(ndvi)为当前计算生境因子对象x的归一化植 被指数值，如图10所示。
[0100]
计算每一个生境因子对象的归一化差异水体指数mndwi，计算公式如下：
[0101][0102]
其中，x为当前计算生境因子对象，i为当前计算生境因子对象x中一个像元， n(x)为当前计算生境因子对象x中所有像元集合，为当前计算生境因 子对象x中一个像元i的绿波段值，为当前计算生境因子对象x中一个像 元i的mir波段值，k
x
(mndwi)为当前计算生境因子对象x的归一化差异水体指数 值，
‑
1<＝k
x
(mndwi)<＝1，正值表示为水体区域，否则为陆地区域，如图11 所示。
[0103]
5)生境因子对象库参数建立
[0104]
根据步骤3)所得的生境因子参数，构建数据集。
[0105]
根据场地基地内生境因子类型构建生境因子数据库。对研究区域内的生境 因子参数解析，得每一块生境因子对象的red波段最大值、最小值、均值、标 准差；green波段最大值、最小值、均值、标准差；blue波段最大值、最小值、 均值、标准差；nir波段最大值、最小值、均值、标准差；bright波段均值；归 一化植被指数(ndvi)最大值、最小值、均值、标准差；归一化差异水体指数 (mndwi)最大值、最小值、均值、标准差；归一化建筑指数(ndbi)最大
值、 最小值、均值、标准差；灰度共生矩阵(glcm)最大值、最小值、均值、标准 差；对象形态参数：紧致度值；长宽比值、纹理特征值。解译图如图12所示。
[0106]
根据已知生境因子类型选取机器学习多维数据样本库。
[0107]
基于融合后多光谱遥感影像图像，结合人工现场勘验所得松岭头地区地表 生境因子数据，运用人眼视觉肌理特征及先验知识对场地内生境因子进行进一 步校正，选取试验地松岭头地区乔木、灌木、地被、竹林、茶园、水体、建筑、 道路、裸地、滩涂十类生境因子要素对象，每类对象数量为50
‑
200不等，构建 生境因子机器学习多维数据样本库。样本库数据标准如图13所示。
[0108]
6)生境因子解析参数pca降维
[0109]
运用主成分分析法对生境因子机器学习多维数据样本库数据库数据降维处 理。在多次训练迭代后，获取每一类生境因子高影响参数，构建高置信度样本库。
[0110]
运用最大似然法模型对高置信度数据样本进行递归最优特征选取。
[0111]
首先配置解析波段为：red波段、green波段、blue波段、nir波段、ndvi 波段、ndbi波段、mndwi波段、glcm波段。
[0112]
其中，运用主成分分析法(principal component analysis,pca)基于特征值 分析对生境因子多维数据样本库数据库数据降维处理，其具体的方法可以采用 以下步骤。
[0113]
(1)输入生境因子多维数据样本集，确定降低空间维度数。
[0114]
输入生境因子多维数据样本集d
i
，
[0115]
d
i
＝{x1,x2,x3,
…
,x
m
}
[0116]
确定生境因子多维数据低维空间维度数d
′
,经过参数选取实验，d
′
取t＝95％ 时，使下式成立的最小值。
[0117][0118]
(2)标准化矩阵数据，计算生境因子多维数据样本各变量的均值与标准差。
[0119]
对所有生境因子多维数据样本数据进行中心化：
[0120][0121]
(3)相关系数矩阵函数，计算特征向量和特征值。
[0122]
计算所有生境因子多维数据样本协方差矩阵xx
t
，对生境因子多维数据样本 协方差矩阵xx
t
做特征值分解。
[0123]
取低维空间维度数d
′
特征值所对应的特征向量w1,w2,w3,
…
,w
d
′
。
[0124]
(4)输出生境因子多维数据样本投影矩阵w
i
[0125]
w
i
＝{w1,w2,w3,
…
,w
d
′
}
[0126]
根据生境因子多维数据样本投影矩阵w
i
，构建多类生境因子影响参数，构 建生境因子多维数据高置信度样本库。
[0127]
7)高置信度数据库样本mle迭代分类
[0128]
运用极大似然估计方法(maximum likelihood estimate，mle)对降维后 的生境因子数据样本进行递归最优特征选取。
[0129]
首先对生境因子降维数据集进行预处理，根据专家先验识别与实地采样调 研结果在生境因子集中抽取训练案例，将生境因子降维数据集分为标注数据与 未标注数据，标注数据配置乔木、灌木、地被、竹林、茶园、水体、建筑、道路、 裸地、滩涂十类生境因子，将未标注数据与高置信度样本混合，以备递归分类器 进行机器学习；
[0130]
其次，设置递归分类器迭代条件，输出结果类别为：10，剔除分数为：0.01， 先验概率权重：equal。当输入数据不满足该条件时，将其职于高置信度样本库 继续迭代训练；当输入数据满足该条件时，输出初级分类器。
[0131]
然后，对初级分类器分类解析结果进行精度验证，采用混淆矩阵算法(error matrix based on samples,embs)建立误差矩阵进行测度评估与回归分析， 选取总体精度(oa)和kappa系数2个测度评价指标(生产者精度(pa)、用户 精度(ua)用于参考)，相关算法预计如下：
[0132][0133]
式中：oa为总体分类精度，r为误差矩阵的行或列数，x
ii
为某类别在矩阵主 对角线上的数目，n为参考样点总数。
[0134][0135]
式中k是kappa系数，r是误差矩阵的行数，x
ii
是对角线上i行i列(即某类生 境因子的正确分类数)的值，x
i+
和x
+i
分别为第i行和第i列的和，n为参考样点 总数。
[0136]
迭代升级初级分类器，设置分类精度理想预期oa≥90％，kappa≥85％，若 初级分类器泛化分类结果不满足精度验证条件，则进行分类器优化，即递归分 类器重新开始迭代，直到初级分类器泛化分类精度达到精度验证条件；若初级 分类器泛化分类结果满足精度验证条件，则输出生境因子最终分类器，如图14 所示。
[0137]
最后，运用生境因子最终分类器对实验样地生境因子进行泛化分类解析， 输出分类解析结果，如图15所示。
[0138]
在本说明书的描述中，参考术语“一个实施例”、“示例”、“具体示例”等 的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于 本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述 不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者 特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
[0139]
以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业 的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中 描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明 还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。