CN117612038A - 基于无人机影像的矿区植被碳汇精细演算方法 - Google Patents
基于无人机影像的矿区植被碳汇精细演算方法 Download PDFInfo
- Publication number
- CN117612038A CN117612038A CN202311566771.2A CN202311566771A CN117612038A CN 117612038 A CN117612038 A CN 117612038A CN 202311566771 A CN202311566771 A CN 202311566771A CN 117612038 A CN117612038 A CN 117612038A
- Authority
- CN
- China
- Prior art keywords
- mining area
- pixel
- vegetation
- pixels
- image
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000005065 mining Methods 0.000 title claims abstract description 153
- 229910052799 carbon Inorganic materials 0.000 title claims abstract description 67
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 title claims abstract description 66
- 238000004364 calculation method Methods 0.000 title claims abstract description 49
- 238000000034 method Methods 0.000 claims abstract description 39
- 238000012545 processing Methods 0.000 claims abstract description 27
- 238000012952 Resampling Methods 0.000 claims abstract description 23
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 20
- 230000005855 radiation Effects 0.000 claims description 26
- 238000001556 precipitation Methods 0.000 claims description 17
- 230000008859 change Effects 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 5
- 238000005192 partition Methods 0.000 claims description 4
- 208000005156 Dehydration Diseases 0.000 claims description 3
- 238000012937 correction Methods 0.000 claims description 3
- 230000000243 photosynthetic effect Effects 0.000 claims description 3
- 230000029058 respiratory gaseous exchange Effects 0.000 claims description 3
- 244000005700 microbiome Species 0.000 claims description 2
- 239000002689 soil Substances 0.000 claims description 2
- 238000000638 solvent extraction Methods 0.000 claims description 2
- 238000012821 model calculation Methods 0.000 abstract 1
- 238000005259 measurement Methods 0.000 description 9
- 210000002569 neuron Anatomy 0.000 description 5
- 238000012544 monitoring process Methods 0.000 description 4
- 239000002023 wood Substances 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 230000000903 blocking effect Effects 0.000 description 3
- 239000003245 coal Substances 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000003702 image correction Methods 0.000 description 3
- 238000010845 search algorithm Methods 0.000 description 3
- 241001270131 Agaricus moelleri Species 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000009919 sequestration Effects 0.000 description 2
- 241000592183 Eidolon Species 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012625 in-situ measurement Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/17—Terrestrial scenes taken from planes or by drones
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/02—Agriculture; Fishing; Forestry; Mining
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/10—Services
- G06Q50/26—Government or public services
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/40—Extraction of image or video features
- G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
- G06V10/462—Salient features, e.g. scale invariant feature transforms [SIFT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/762—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using clustering, e.g. of similar faces in social networks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/188—Vegetation
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Business, Economics & Management (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Multimedia (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Tourism & Hospitality (AREA)
- Economics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Business, Economics & Management (AREA)
- Strategic Management (AREA)
- Primary Health Care (AREA)
- Marketing (AREA)
- Human Resources & Organizations (AREA)
- Life Sciences & Earth Sciences (AREA)
- Software Systems (AREA)
- Marine Sciences & Fisheries (AREA)
- Mining & Mineral Resources (AREA)
- Artificial Intelligence (AREA)
- Agronomy & Crop Science (AREA)
- Computing Systems (AREA)
- Animal Husbandry (AREA)
- Medical Informatics (AREA)
- Evolutionary Computation (AREA)
- Remote Sensing (AREA)
- Databases & Information Systems (AREA)
- Development Economics (AREA)
- Educational Administration (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于无人机影像的矿区植被碳汇精细演算方法，其方法包括：S1、采用多光谱无人机对目标矿区进行重叠正投影拍摄得到影像数据集并得到目标矿区影像；S2、对目标矿区影像逐像素进行灰度化处理与属性划分；S3、基于改进的双线性插值算法对目标矿区影像按照设定像元大小进行重采样处理，统计得到各个像元的植被覆盖度；S4、采集目标矿区的气象参数输入反距离权重插值模型中，反距离权重插值模型采用距离倒数乘方法进行覆盖所有像元的气象参数的插值处理，基于UAV‑CASA模型计算得到各个像元的碳汇量。本发明能够得到各个像元的碳汇量、目标矿区碳汇量分布、总碳汇量、植被净生态系统生产力总量，对科学评价矿区植被碳汇具有重要意义。
Description
技术领域
本发明涉及采矿、遥感及地理信息领域，尤其涉及一种基于无人机影像的矿区植被碳汇精细演算方法。
背景技术
随着全球气候变化和环境问题的不断加剧，减少碳排放已经逐渐成为世界各国的共识，煤炭行业更是碳减排的重点领域。露天煤矿的开采会对其周边生态环境造成破坏，并影响其所在地区生态系统及植被的固碳能力。因此，有必要建立一种面向矿区的精细植被碳汇计算方法，实现对矿区植被碳汇的精确计算及监测。
目前植被碳汇的计算方法主要为实地测量法和遥感影像反演法。对于矿区的植被碳汇计算，实地测量法具有较高的可靠度和精度，但需要大量时间和劳动力成本，且成本高、覆盖范围小；且由于矿区地域的特殊性，在矿区的复杂地形条件下进行设立观测站点、统计植被的叶面积指数等生态数据的作业较为困难，因此并不能满足矿区植被碳汇长期大范围监测的需求。遥感影像反演法可以大范围计算生态系统植被碳汇，注重大尺度数据反馈与宏观统计，基于卫星遥感影像的矿区植被碳汇计算在大气干扰云遮挡等情况的影响下无有效覆盖研究区域且难以实现精细化数据采集。因此通过上述两种方法进行矿区植被碳汇计算均存在一定程度的缺陷；目前，本领域较为关注矿区等小尺度、地形地势复杂区域的植被碳汇精细演算，同时更看重其演算精度，通过拍照影像(以无人机拍照为例，无人机灵活性强，对矿区等小尺度地形地势复杂区域优势明显，可根据需要调整飞行高度、路线和飞行范围，获取矿区范围的遥感影像，进而实现计算范围的灵活变化；且无人机遥感影像数据获取成本低、适用性广，是计算矿区植被碳汇的较为理想的数据源)缺少卫星遥感影像的光谱波段信息，无法基于拍照影像按照遥感影像进行反演，故此，亟需研究小尺度、复杂区域并依赖拍照影像的植被碳汇演算技术。
发明内容
本发明的目的在于解决现有技术存在的技术问题，提供一种基于无人机影像的矿区植被碳汇精细演算方法，基于无人机正投影拍摄获取影像数据，对影像进行逐像素特定方法的灰度化处理，然后依靠设定的灰度级数划分为两个灰度类，根据类间方差最大值得到最优阈值，按照最优阈值对像素植被或非植被属性划分，按照设定像元大小进行重采样处理并得到各个像元的植被覆盖度，进而计算得到各个像元的碳汇量，得到目标矿区碳汇量分布、总碳汇量、植被净生态系统生产力总量；达到对矿区生态系统碳汇精确计算的目的，对科学评价矿区植被碳汇具有重要意义。
本发明的目的通过下述技术方案实现：
一种基于无人机影像的矿区植被碳汇精细演算方法，其方法包括：
S1、采用多光谱无人机对目标矿区进行重叠正投影拍摄得到影像数据集，基于目标矿区构建坐标系，将影像数据集中影像数据配准并拼接、校正成目标矿区影像；
S2、对目标矿区影像逐像素进行如下灰度化处理：将目标矿区影像中按照像素解析R、G、B三个颜色通道，计算得到像素的灰度值EXG，EXG＝2*G1-R1-B1，G1、R1、B1表示三个颜色通道所对应的数值；按照像素的灰度值EXG进行像素的灰度化处理；
将目标矿区影像的像素设定具有L个灰度级数，先预设阈值K，将目标矿区影像的像素按照阈值K划分为C0、C1两个灰度类，然后计算得到目标矿区影像的平均灰度μ、C0类像素的平均灰度μ0、C1类像素的平均灰度μ1，然后按照如下公式计算类间方差δ2(K)：
δ2(k)＝ω0(μ-μ0)2+ω1(μ-μ1)2；其中ω0是C0类的像素占比比，ω1是C1类的像素占比；
更换预设阈值K并获取类间方差δ2(K)的最大值时的预设阈值K作为最优阈值T；将目标矿区影像中灰度值不小于T的像素点划分为植被，其余划分为非植被；
S3、基于改进的双线性插值算法对目标矿区影像按照设定像元大小进行重采样处理，然后统计得到各个像元的植被覆盖度；
S4、基于坐标系的坐标数据构建反距离权重插值模型，采集目标矿区的气象参数输入反距离权重插值模型中，气象参数包括月平均气温、月总降水、月太阳总辐射，反距离权重插值模型采用距离倒数乘方法进行覆盖所有像元的气象参数的插值处理；基于UAV-CASA模型按照如下方法计算得到各个像元的碳汇量NPP(x，t)：
NPP(x，t)＝APAR(x，t)×ε(x，t)
APAR(x，t)＝SOL(x，t)×FPAR(x，t)×0.5×[FVC×(FPARmax-FPARmin)+FPARmin]
ε(x，t)＝T1(x，t)×T2(x，t)×W(x，t)×ε*；其中x表示像元位置，t表示时间，APAR(x，t)表示像元x在t月的植被光合有效辐射，ε(x，t)为像元x在t月真实的植被光能利用率，SOL(x，t)表示像元x在t月的太阳总辐射量，FPARmax与FPARmin分别取0.95和0.001，T1(x，t)和T2(x，t)分别表示低温和高温下对光能利用效率的胁迫系数，W(x，T)为水分胁迫系数，ε*为理想条件下的最大光能利用率。
为了更好地实现本发明，本发明还包括如下方法：
S5、按照如下公式计算得到各个像元的植被净生态系统生产力NEP：
NEP(x，t)＝NPP(x，t)-RH，其中RH表示土壤微生物的呼吸量；
统计汇总目标矿区所有像元的碳汇量得到碳汇总量；统计汇总目标矿区所有像元的植被净生态系统生产力得到植被净生态系统生产力总量。
优选地，步骤S2中还包括如下方法：
基于植被、非植被类型对目标矿区影像中的像素进行聚类分块，在每个分块中基于SIFT算法找出多个SIFT特征点，按照特征匹配算法在分块中扩大领域进行邻近像素搜索及类型重新判断，若特征匹配算法的中心像素为SIFT特征点，则周围邻域像素以SIFT特征点的植被或非植被进行属性更改，直到周围邻域像素不同于SIFT特征点的属性占比超过A1时停止更改作业；若特征匹配算法的中心像素为非SIFT特征点，以周围邻域像素的属性占比超过A2的属性更改中心像素的属性，否则中心像素的属性不更改。
优选地，在步骤S2中，目标矿区影像的像素的L个灰度级数从0～L-1，第i级灰度像素点数为Ni，i位于0～L-1的范围内，目标矿区影像的像素总数为N，则目标矿区影像中第i级像素点的概率为目标矿区影像的平均灰度计算表达式为：/>C0类像素的平均灰度计算表达式为：/>C1类像素的平均灰度计算表达式为：/>
优选地，在步骤S2中，像素点按照坐标(x1，y1)表示像素点的灰度值为f(x1，y1)，g(x1，y1)表示属性，属性为划分植被或非植被的结果。
优选地，在步骤S3中，根据重采样处理的目标分辨率从目标矿区影像计算宽和高的缩放常数，表达式如下：
其中pWidth和tWidth分别表示目标矿区影像重采样的宽度和高度，xt和yt分别表示进行重采样后像元的横坐标和纵坐标；x和y为重采样后像元对应原始图像的横坐标和纵坐标；
遍历统计像元对应原始图像所有像素的属性，并计算植被属性的占比作为像元的植被覆盖度。
优选地，在步骤S3中，根据遍历位置找到其中一个像元j对应原始图像的位置坐标P(x，y)，找出距离像元j最近的四个像元，先在x方向进行两次插值，再在y方向进行一次插值。
优选地，在步骤S4中，反距离权重插值模型的距离倒数乘方法表达式如下：
Z0为像元所在位置的估计值，Zi为像元所在位置附近气象站点的气象参数，di为附近气象站点与像元所在位置的距离，n为气象站点数目，r为指定的幂数。
优选地，在步骤S1中，坐标系中具有目标矿区地图，在目标矿区地图上生成密集点云并生成三角网，利用密集点云及三角网进行影像数据配准，利用影像拼接工具进行影像数据拼接，在目标矿区地图使用ArcMap软件进行校正处理。
优选地，在步骤S1中，多光谱无人机设置的旁向重叠度为50％～70％，航向重叠度为60％～80％。
本发明较现有技术相比，具有以下优点及有益效果：
(1)本发明基于无人机正投影拍摄获取影像数据，对影像进行逐像素特定方法的灰度化处理，然后依靠设定的灰度级数划分为两个灰度类，根据类间方差最大值得到最优阈值，按照最优阈值对像素植被或非植被属性划分，按照设定像元大小进行重采样处理并得到各个像元的植被覆盖度，进而计算得到各个像元的碳汇量，得到目标矿区碳汇量分布、总碳汇量、植被净生态系统生产力总量；达到对矿区生态系统碳汇精确计算的目的，对科学评价矿区植被碳汇具有重要意义。
(2)本发明采用面向无人机影像的变邻域局部搜索算法对处理后像素的植被与非植被结果进行遍历和判断，对可疑结果像素通过对邻近该像元的植被和非植被特征进行分析和判断，进而对该可疑像素划分结果进行验证，实现对像素计算结果的纠正，提高植被覆盖度的计算精度，进而提高矿区植被碳汇的计算精度。
(3)本发明利用改进的双线性插值算法对无人机影像进行重采样，重采样的像元大小可根据实际需求(设置长度与宽度，或设置像元大小)而进行设定，实现对影像边缘范围的精细化处理，同时能够便于实现植被覆盖度FVC计算以及植被覆盖度多尺度范围内动态调整，实现与重采样目标分辨率更好的拟合。
(4)本发明是基于无人机进行正投影拍摄得到影像数据作为数据源，与传统卫星影像相比，具有成本低、操作简单、获取影像速度快、影像空间分辨率高等优点，且不受云层遮挡等条件的影响；本发明可实现从原始影像到植被碳汇的自动化计算处理，准确计算目标区域尤其以小范围矿区等范围碳汇，具有较强的实用价值，能够利用无人机影像精准的计算在矿区复杂场景下的植被碳汇，能够提升计算矿区植被碳汇的准确性；对科学评价矿区植被碳汇具有重要意义。
附图说明
图1为本发明矿区植被碳汇精细演算方法的流程示意图；
图2为实施例中邻近像素搜索的原理示意图；
图3为实施例中从目标矿区截取一部分区域举例计算植被覆盖度的原理示意图；
图4为实施例举例西湾矿区作为目标矿区的目标矿区影像示意图；
图5为西湾矿区截取举例部分像元显示植被或非植被属性前后影像的示意图；
图6为实施例中西湾矿区的气温插值示意图；
图7为实施例中西湾矿区的降水量插值示意图；
图8为实施例中西湾矿区的太阳辐射插值示意图；
图9为西湾矿区碳汇量NPP(x，t)图标化展示的示意图。
具体实施方式
下面结合实施例对本发明作进一步地详细说明：
实施例
如图1所示，一种基于无人机影像的矿区植被碳汇精细演算方法，其方法包括：
S1、采用多光谱无人机对目标矿区进行重叠正投影拍摄得到影像数据集，优选地，多光谱无人机设置的旁向重叠度为50％～70％，航向重叠度为60％～80％，划定目标矿区内的飞行区域并进行航线规划，设定飞行航高、返航高度、航向重叠度、旁向重叠度、飞行方式、航带条数等参数，飞行高度设置不能高于飞行器默认限高；多光谱无人机可以包括传感器系统、地面监控系统等组成部分，其中传感器系统通常包括小型数码相机、定点摄影技术等部分，地面监控系统与数据接收系统主要负责航线规划与数据接收。基于目标矿区构建坐标系，将影像数据集中影像数据配准并拼接、校正成目标矿区影像。优选地，坐标系中具有目标矿区地图，在目标矿区地图上生成密集点云并生成三角网，利用密集点云及三角网进行影像数据配准，利用影像拼接工具(包括PhotoScan和Pix4dmapper等，具体流程为：打开影像拼接工具软件，新建工作空间Workspace，将影像数据集中影像数据添加到工作空间中，将影像数据在坐标系及目标矿区地图上依次排列，依靠密集点云及三角网进行影像数据配准，然后裁剪与拼合处理)进行影像数据拼接，在目标矿区地图使用ArcMap软件进行校正处理(若生成的拼接影像存在偏差，可利用ArcMap等软件进行影像校正)。
本实施例以陕西省榆林市神木市西湾矿区为目标矿区举例说明，多光谱无人机采用大疆精灵4RTK无人机(相机分辨率为0.05m)对西湾矿区进行重叠正投影拍摄得到高分辨率的影像数据构成影像数据集。选择WGS84经纬度坐标系，选择处理模板，生成密集点云和三角网；在Pix4Dmapper软件中新建工作空间，将获取的无人机影像数据添加到工作空间中，Pix4Dmapper软件开始自动处理图像；处理完毕后，将得到拼接完整的DOM正射影像；将处理后的DOM正射影像利用Arcmap进行影像校正，得到的校正后的正射影像如图4所示。
S2、对目标矿区影像逐像素进行如下灰度化处理：将目标矿区影像中按照像素解析R、G、B三个颜色通道(获取目标矿区影像的RGB颜色空间，得到三个颜色通道)，计算得到像素的灰度值EXG，EXG＝2*G1-R1-B1，G1、R1、B1表示三个颜色通道所对应的数值。按照像素的灰度值EXG进行像素的灰度化处理。
将目标矿区影像的像素设定具有L个灰度级数，先预设阈值K(优选地，L为256，L个灰度级数为0～255，阈值K位于0～255范围内)，将目标矿区影像的像素按照阈值K划分为C0、C1两个灰度类，以L个灰度级数从0～L-1为例，则C0灰度类区间为0～K，C1灰度类区间为K～L-1，然后计算得到目标矿区影像的平均灰度μ、C0类像素的平均灰度μ0、C1类像素的平均灰度μ1，然后按照如下公式计算类间方差δ2(k)：
δ2(k)＝ω0(μ-μ0)2+ω1(μ-μ1)2。其中ω0是C0类的像素占比比，ω1是C1类的像素占比。
在本实施例中，目标矿区影像的像素的L个灰度级数从0～L-1(优选地，L为256，L个灰度级数为0～255)，第i级灰度像素点数为Ni，i位于0～L-1的范围内，i表示灰度级的级数，目标矿区影像的像素总数为N，则目标矿区影像中第i级像素点的概率为目标矿区影像的平均灰度计算表达式为：/>C0类像素的平均灰度计算表达式为：类像素的平均灰度计算表达式为：/>
更换预设阈值K并获取类间方差δ2(K)的最大值时的预设阈值K作为最优阈值T，即求取类间方差δ2(K)的最大值时的阈值K，将此时的阈值K作为最优阈值。将目标矿区影像中灰度值不小于T的像素点划分为植被，其余划分为非植被。
在一些实施例中，在步骤S2中，像素点的坐标记为(x1，y1)，f(x1，y1)表示像素点(x1，y1)的灰度值，g(x1，y1)表示像素点(x1，y1)的属性，属性为划分植被或非植被的结果；这样像素点均按照坐标记录有灰度值、属性，可以便于后续的数据处理及邻近像素搜索处理。
在一些实施例中，基于植被、非植被类型对目标矿区影像中的像素进行聚类分块，在每个分块中基于SIFT算法找出多个SIFT特征点，按照特征匹配算法在分块中扩大领域进行邻近像素搜索及类型重新判断，若特征匹配算法的中心像素为SIFT特征点，则周围邻域像素以SIFT特征点的植被或非植被进行属性更改，直到周围邻域像素不同于SIFT特征点的属性占比超过A1(A1为设置的占比阈值)时停止更改作业；若特征匹配算法的中心像素为非SIFT特征点，以周围邻域像素的属性占比超过A2(A2为设置的占比阈值)的属性更改中心像素的属性，否则中心像素的属性不更改。图2示出了分块中四邻域(包括中心像素、周围邻域像素)、六邻域、九邻域的三种情况，按照上述方法进行目标矿区影像中的像素属性更改处理，提高了植被和非植被在像素级别上的划分精度。
本实施例以陕西省榆林市神木市西湾矿区为目标矿区举例说明，按照上述方法进行灰度化处理，然后通过最优阈值T将目标矿区影像中灰度值不小于T的像素点划分为植被，其余划分为非植被；变邻域局部搜索算法对划分结果进行逐像素遍历并结合被判断像素与其周围像素进行判断，提高划分精度；在西湾矿区中截取包含若干区域所有像素的像元并得出植被或非植被属性划分结果，按照白色为植被、黑色为非植被进行图表化展示，如图5所示，图中选择四个包含若干区域所有像素的像元并得到植被或非植被属性图表结果。
将进行灰度化处理后的影像利用基于面向无人机影像的优化OSTU无人机影像阈值确定最优阈值T。通过将灰度化后的像元值与最优阈值T进行对比，将所有像元进行植被与非植被像元的划分，并采用面向无人机影像的变邻域局部搜索算法对划分结果像元进行逐像元遍历并结合被判断像元与其周围像元进行判断，提高划分精度；根据划分的影像区域中植被像元占总像元的百分比可以求出植被的植被覆盖度。西湾矿区部分区域像元划分结果截取部分区域如图6所示；该图选取西湾矿区四个不同区域划分结果及进行展示，其中有图白色部分为植被区域，黑色部分为非植被区域。
S3、基于改进的双线性插值算法对目标矿区影像按照设定像元大小进行重采样处理，然后统计得到各个像元的植被覆盖度。如图3所示(其中黑色为非植被像素或像元，白色为植被像素或像元)，从目标矿区截取一部分区域(该区域可以是包括一片区域的像素，也可以是包括若干个像元的区域，根据计算植被覆盖度的尺度进行灵活应用)举例，以图3的区域是包括一片区域的像素为例，进行植被像素与非植被像素的统计，植被像素的占比作为该像元的植被覆盖度，同理，若一个像元包括若干个小像元(一定区域像素集合)，则进行在小像元级别上的植被占比统计或者进行小像元级别扩展至像素级别的植被占比统计。本实施例以陕西省榆林市神木市西湾矿区为目标矿区举例说明，参见图5的像元举例，根据划分的影像像元中植被像素占总像素的百分比可以求出该像元的植被覆盖度。
在一些实施例中，根据设定像元大小采用重采样处理方式进行上下尺度之间的数据采样处理；重采样是从高分辨率影像中提取出低分辨率影像的过程，本实施例采用改进的双线性插值算法对目标矿区影像进行重采样，其面向特定的小范围目标矿区影像进行重采样，在进行影像重采样计算前遍历图像所有像元(像元的最小尺度可以为像素)，根据重采样目标像元大小自动更改目标矿区影像的范围，算法根据重采样像元大小与原影像的像元大小和数量进行自动计算，根据计算结果对影像边界进行矩形化处理，使数字图像矩阵进行规则排列，根据边界像元与重采样目标像元大小的数量关系进行判断并进行像元的增减，使之与重采样后像元大小及长宽数量吻合；同时该算法在遍历图像时可以在处理速度和图像质量之间取得平衡，既能够较快地进行图像重采样，同时又能最大限度的保留图像高频细节信息。其中，重采样处理的目标分辨率(输入所需要重采样的目标像元大小，根据目标像元与原始像元大小对遍历后的图像矩阵边缘进行处理，使之矩形化，使其与重采样目标像元大小及长宽数量相吻合)从目标矿区影像计算宽和高的缩放常数，表达式如下：
其中pWidth和tWidth分别表示目标矿区影像重采样的宽度和高度，xt和yt分别表示进行重采样后像元的横坐标和纵坐标；x和y为重采样后像元对应原始图像的横坐标和纵坐标，重采样后像元对应原始图像的坐标记为P(x，y)。
在一些实施例中，根据遍历位置找到其中一个像元j对应原始图像的位置坐标P(x，y)，找出距离像元j最近的四个像元，先在x方向进行两次插值，再在y方向进行一次插值。优选地，若距离P(x，y)点最近的四个像元点分别为Q11＝(x1，y1)、Q12＝(x1，y2)、Q21＝(x2，y1)、Q22＝(x2，y2)，则使用上述四个点，先在x方向进行两次插值，并在y方向进行一次插值，最终得到重采样后的目标像素值；插值公式可以简化如下：
f(x，y)＝(1-u)(1-v)f(x1，y1)+u(1-v)f(x2，y1)+(1-u)vf(x1，y2)+uvf(x2，y2)
f(x1，y1)、f(x2，y1)、f(x1，y2)、f(x2，y2)分别为对应的四个点像素值。
然后遍历统计像元对应原始图像所有像素的属性，并计算植被属性的占比作为像元的植被覆盖度。本实施例以陕西省榆林市神木市西湾矿区为目标矿区举例说明，设置重采样分辨率为30m并作为像元大小，西湾矿区按照像元统计出的植被覆盖度分布。
S4、基于坐标系的坐标数据构建反距离权重插值模型，采集目标矿区的气象参数输入反距离权重插值模型中，气象参数包括月平均气温、月总降水、月太阳总辐射，反距离权重插值模型采用距离倒数乘方法进行覆盖所有像元的气象参数的插值处理。基于UAV-CASA模型按照如下方法计算得到各个像元的碳汇量NPP(x，T)：
NPP(x，t)＝APAR(x，t)×ε(x，t)
APAR(x，t)＝SOL(x，t)×FPAR(x，t)×0.5×[FVC×(FPARmax-FPARmin)+FPARmin]
ε(x，t)＝T1(x，t)×T2(x，t)×W(x，t)×ε*。其中x表示像元位置，t表示时间(以月为单位)，APAR(x，t)表示像元x在t月的植被光合有效辐射(MU·m-2)，ε(x，t)为像元x在t月真实的植被光能利用率(gC/MJ)，SOL(x，t)表示像元x在t月的太阳总辐射量(MJ/(m2·month))，FPARmax与FPARmin分别取0.95和0.001，常数0.5表示植被在波长为0.38-0.71μm处所能利用的太阳有效辐射占太阳总辐射的比例。T1(x，t)和T2(x，t)分别表示低温和高温下对光能利用效率的胁迫系数，W(x，T)为水分胁迫系数，ε*为理想条件下的最大光能利用率。本实施例以陕西省榆林市神木市西湾矿区为目标矿区举例说明，基于UAV-CASA模型按照如下方法计算得到各个像元的碳汇量分布如图9所示，进一步地，可以算出西湾矿区整个区域内的总碳汇量(即将西湾矿区内所有像元的碳汇量相加得到)。
在一些实施例中，反距离权重插值模型的距离倒数乘方法(是一种加权平均插值法，进行确切的或圆滑的方式插值处理，在具体实施时，可以利用ERA5和NECP月平均气温、月总降水、月太阳总辐射等气象参数，还可以结合气象站点实测数据对气温降水等数据进行融合，提高气象数据的精度。反距离权重插值模型采用距离倒数乘方法构建插值温度与站点实测温度之间、插值降水与站点实测降水之间、插值太阳总辐射与站点实测太阳总辐射的关系，能够得到目标矿区范围内各个位置的平均气温、降水和日照辐射月度气象数据)表达式如下：
Z0为像元所在位置的估计值，Zi为像元所在位置附近气象站点的气象参数，di为附近气象站点与像元所在位置的距离，n为气象站点数目，r为指定的幂数。通过ERA5和NECP获取得到月平均气温、月总降水、月太阳总辐射的数据，然后按照上述方法依次进行月平均气温、月总降水、月太阳总辐射的插值处理，得到覆盖所有像元的数据。
本实施例以陕西省榆林市神木市西湾矿区为目标矿区举例说明，利用ERA5和NCEP再分析数据基于反距离权重插值法插值出包含西湾矿区范围的月平均气温和月总降水数据，基于反距离权重插值模型将插值后的数据进行细化，并构建基于ERA5和NECP插值后数据与西湾矿区附近站点实测月平均温度及月降水量数据之间的关系，进而提高气温和降水的准确性，反距离权重插值模型在插值时还可以采用如下表达式进行插值处理：
yi＝f(Ii)
上式中Xi(j＝1，2，3…n)为神经元获得的再分析气象数据信息；θi为
阈值；Wij为第j个神经元到第i个神经元间的连接权值；f(·)为传递函数。通过输入气象数据再分析资料和气象站点实测数据，建立反距离权重插值模型，融合各方面因素输出进行学习分析后的气象数据，有效提高气象数据的精度。太阳辐射数据则采用辐射站点实测数据作为西湾矿区真实月太阳辐射值。进行精度提高后的西湾矿区气温、降水、辐射插值示意图如图6(月平均气温)、图7(月总降水)、图9(月太阳总辐射)所示。
S5、按照如下公式计算得到各个像元的植被净生态系统生产力NEP(单位为gC/(m2·a))：
NEP(x，t)＝NPP(x，t)-RH，其中RH表示土壤微生物的呼吸量(gC/(m2·a))。
统计汇总目标矿区所有像元的碳汇量得到碳汇总量。统计汇总目标矿区所有像元的植被净生态系统生产力得到植被净生态系统生产力总量。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种基于无人机影像的矿区植被碳汇精细演算方法，其特征在于：其方法包括：
S1、采用多光谱无人机对目标矿区进行重叠正投影拍摄得到影像数据集，基于目标矿区构建坐标系，将影像数据集中影像数据配准并拼接、校正成目标矿区影像；
S2、对目标矿区影像逐像素进行如下灰度化处理：将目标矿区影像中按照像素解析R、G、B三个颜色通道，计算得到像素的灰度值EXG，EXG＝2*G1-R1-B1，G1、R1、B1表示三个颜色通道所对应的数值；按照像素的灰度值EXG进行像素的灰度化处理；
将目标矿区影像的像素设定具有L个灰度级数，先预设阈值K，将目标矿区影像的像素按照阈值K划分为C0、C1两个灰度类，然后计算得到目标矿区影像的平均灰度μ、C0类像素的平均灰度μ0、C1类像素的平均灰度μ1，然后按照如下公式计算类间方差δ2(k)：
δ2(k)＝ω0(μ-μ0)2+ω1(μ-μ1)2；其中ω0是C0类的像素占比比，ω1是C1类的像素占比；
更换预设阈值K并获取类间方差δ2(k)的最大值时的预设阈值K作为最优阈值T；将目标矿区影像中灰度值不小于T的像素点划分为植被，其余划分为非植被；
S3、基于改进的双线性插值算法对目标矿区影像按照设定像元大小进行重采样处理，然后统计得到各个像元的植被覆盖度；
S4、基于坐标系的坐标数据构建反距离权重插值模型，采集目标矿区的气象参数输入反距离权重插值模型中，气象参数包括月平均气温、月总降水、月太阳总辐射，反距离权重插值模型采用距离倒数乘方法进行覆盖所有像元的气象参数的插值处理；基于UAN-CASA模型按照如下方法计算得到各个像元的碳汇量NPP(x，t)：
NPP(x，t)＝APAR(x，t)×ε(x，t)
APAR(x，t)＝SOL(x，t)×FPAR(x，t)×0.5×[FVC×(FPARmax-FPARmin)+FPARmin]
ε(x，t)＝T1(x，t)×T2(x，t)×W(x，t)×ε*；其中x表示像元位置，t表示时间，APAR(x,t)表示像元x在t月的植被光合有效辐射，ε(x，t)为像元x在t月真实的植被光能利用率，SOL(x，t)表示像元x在t月的太阳总辐射量，FPARmax与FPARmin分别取0.95和0.001，T1(x，t)和T2(x，t)分别表示低温和高温下对光能利用效率的胁迫系数，W(x，t)为水分胁迫系数，ε*为理想条件下的最大光能利用率。
2.按照权利要求1所述的基于无人机影像的矿区植被碳汇精细演算方法，其特征在于：还包括如下方法：
S5、按照如下公式计算得到各个像元的植被净生态系统生产力NEP：
NEP(x，t)＝NPP(x，t)-RH，其中RH表示土壤微生物的呼吸量；
统计汇总目标矿区所有像元的碳汇量得到碳汇总量；统计汇总目标矿区所有像元的植被净生态系统生产力得到植被净生态系统生产力总量。
3.按照权利要求1所述的基于无人机影像的矿区植被碳汇精细演算方法，其特征在于：步骤S2中还包括如下方法：
基于植被、非植被类型对目标矿区影像中的像素进行聚类分块，在每个分块中基于SIFT算法找出多个SIFT特征点，按照特征匹配算法在分块中扩大领域进行邻近像素搜索及类型重新判断，若特征匹配算法的中心像素为SIFT特征点，则周围邻域像素以SIFT特征点的植被或非植被进行属性更改，直到周围邻域像素不同于SIFT特征点的属性占比超过A1时停止更改作业；若特征匹配算法的中心像素为非SIFT特征点，以周围邻域像素的属性占比超过A2的属性更改中心像素的属性，否则中心像素的属性不更改。
4.按照权利要求1所述的基于无人机影像的矿区植被碳汇精细演算方法，其特征在于：在步骤S2中，目标矿区影像的像素的L个灰度级数从0～L-1，第i级灰度像素点数为Ni，i位于0～L-1的范围内，目标矿区影像的像素总数为N，则目标矿区影像中第i级像素点的概率为目标矿区影像的平均灰度计算表达式为：/>C0类像素的平均灰度计算表达式为：/>C1类像素的平均灰度计算表达式为：/>
5.按照权利要求1或4所述的基于无人机影像的矿区植被碳汇精细演算方法，其特征在于：在步骤S2中，像素点按照坐标(X1，y1)表示像素点的灰度值为f(x1，y1)，g(x1，y1)表示属性，属性为划分植被或非植被的结果。
6.按照权利要求1所述的基于无人机影像的矿区植被碳汇精细演算方法，其特征在于：在步骤S3中，根据重采样处理的目标分辨率从目标矿区影像计算宽和高的缩放常数，表达式如下：
x＝xt*(pWidth/tWidth)
y＝yt*(pWidth/tWidth)；其中pWidth和tWidth分别表示目标矿区影像重采样的宽度和高度，xt和yt分别表示进行重采样后像元的横坐标和纵坐标；x和y为重采样后像元对应原始图像的横坐标和纵坐标；
遍历统计像元对应原始图像所有像素的属性，并计算植被属性的占比作为像元的植被覆盖度。
7.按照权利要求6所述的基于无人机影像的矿区植被碳汇精细演算方法，其特征在于：在步骤S3中，根据遍历位置找到其中一个像元j对应原始图像的位置坐标P(x，y)，找出距离像元j最近的四个像元，先在x方向进行两次插值，再在y方向进行一次插值。
8.按照权利要求1所述的基于无人机影像的矿区植被碳汇精细演算方法，其特征在于：在步骤S4中，反距离权重插值模型的距离倒数乘方法表达式如下：
Z0为像元所在位置的估计值，Zi为像元所在位置附近气象站点的气象参数，di为附近气象站点与像元所在位置的距离，n为气象站点数目，r为指定的幂数。
9.按照权利要求1所述的基于无人机影像的矿区植被碳汇精细演算方法，其特征在于：在步骤S1中，坐标系中具有目标矿区地图，在目标矿区地图上生成密集点云并生成三角网，利用密集点云及三角网进行影像数据配准，利用影像拼接工具进行影像数据拼接，在目标矿区地图使用ArcMap软件进行校正处理。
10.按照权利要求1所述的基于无人机影像的矿区植被碳汇精细演算方法，其特征在于：在步骤S1中，多光谱无人机设置的旁向重叠度为50％～70％，航向重叠度为60％～80％。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311566771.2A CN117612038B (zh) | 2023-11-22 | 基于无人机影像的矿区植被碳汇精细演算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311566771.2A CN117612038B (zh) | 2023-11-22 | 基于无人机影像的矿区植被碳汇精细演算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117612038A true CN117612038A (zh) | 2024-02-27 |
| CN117612038B CN117612038B (zh) | 2024-07-02 |
Family
ID=
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106897679A (zh) * | 2017-02-13 | 2017-06-27 | 长江水利委员会长江科学院 | 一种基于改进模糊c均值聚类的语义变化检测方法及系统 |
| CN111259963A (zh) * | 2019-12-27 | 2020-06-09 | 南京林业大学 | 一种区域植被指标的驱动因素分析方法、装置及存储介质 |
| CN115391882A (zh) * | 2022-08-09 | 2022-11-25 | 中信建筑设计研究总院有限公司 | 一种多元数据融合全息空间系统的构建方法 |
| CN116823896A (zh) * | 2023-06-21 | 2023-09-29 | 中国地质大学（武汉） | 一种高植被覆盖下靶矿区范围预测方法、装置及电子设备 |
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106897679A (zh) * | 2017-02-13 | 2017-06-27 | 长江水利委员会长江科学院 | 一种基于改进模糊c均值聚类的语义变化检测方法及系统 |
| CN111259963A (zh) * | 2019-12-27 | 2020-06-09 | 南京林业大学 | 一种区域植被指标的驱动因素分析方法、装置及存储介质 |
| CN115391882A (zh) * | 2022-08-09 | 2022-11-25 | 中信建筑设计研究总院有限公司 | 一种多元数据融合全息空间系统的构建方法 |
| CN116823896A (zh) * | 2023-06-21 | 2023-09-29 | 中国地质大学（武汉） | 一种高植被覆盖下靶矿区范围预测方法、装置及电子设备 |
Non-Patent Citations (1)
| Title |
|---|
| 吴晓华 等: "煤炭企业绿色矿山建设标准体系研究", 《中国煤炭》, 30 June 2022 (2022-06-30), pages 1 - 6 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114037911B (zh) | 一种顾及生态分区的大尺度森林高度遥感反演方法 | |
| CN102074047B (zh) | 一种高精细城市三维建模方法 | |
| CN102506824B (zh) | 一种城市低空无人机系统生成数字正射影像图的方法 | |
| CN110263717B (zh) | 一种融入街景影像的土地利用类别确定方法 | |
| CN102147250B (zh) | 一种数字线划图测图方法 | |
| CN111626947B (zh) | 基于生成对抗网络的地图矢量化样本增强方法及系统 | |
| CN113920262B (zh) | 一种增强边缘取样与改进Unet模型的矿区FVC计算方法及系统 | |
| CN107688818A (zh) | 一种基于卫星遥感影像特征分析的路径智能选择方法及系统 | |
| CN110046563B (zh) | 一种基于无人机点云的输电线路断面高程修正方法 | |
| CN112270320B (zh) | 一种基于卫星影像校正的输电线路杆塔坐标校准方法 | |
| CN113744249B (zh) | 一种海洋生态环境损害调查方法 | |
| Ramli et al. | Homogeneous tree height derivation from tree crown delineation using Seeded Region Growing (SRG) segmentation | |
| CN111428582B (zh) | 一种利用互联网街景照片计算城市天空开阔度的方法 | |
| CN115424135A (zh) | 一种用于植被提取深度学习的四通道影像处理方法 | |
| CN117387580B (zh) | 一种基于倾斜摄影大比例尺地形图的测绘方法及系统 | |
| CN112254713B (zh) | 一种面向高大密集建筑群的无人机倾斜摄影参数确定方法 | |
| Niederheiser et al. | Dense image matching of terrestrial imagery for deriving high-resolution topographic properties of vegetation locations in alpine terrain | |
| CN117612038B (zh) | 基于无人机影像的矿区植被碳汇精细演算方法 | |
| CN116310901A (zh) | 一种基于低空遥感的泥石流物源动态运移识别方法 | |
| CN117612038A (zh) | 基于无人机影像的矿区植被碳汇精细演算方法 | |
| CN115700777A (zh) | 一种基于无人机和数字地表模型的公路施工阶段预测方法 | |
| Oppong Adu et al. | Surface condition assessment of unpaved roads through the use of unmanned aerial vehicle | |
| Trevoho et al. | Aerial data application for construction of large-scale plans | |
| CN115423975A (zh) | 一种基于可见光影像及深度学习算法的地面高程提取方法 | |
| CN113192194A (zh) | 一种植被区地形图的绘制方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |