CN116452023B - 基于低频微波雷达vod数据的公里级碳储量评估方法 - Google Patents

基于低频微波雷达vod数据的公里级碳储量评估方法 Download PDF

Info

Publication number
CN116452023B
CN116452023B CN202211664301.5A CN202211664301A CN116452023B CN 116452023 B CN116452023 B CN 116452023B CN 202211664301 A CN202211664301 A CN 202211664301A CN 116452023 B CN116452023 B CN 116452023B
Authority
CN
China
Prior art keywords
vegetation
data
water content
model
resolution
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.)
Active
Application number
CN202211664301.5A
Other languages
English (en)
Other versions
CN116452023A (zh
Inventor
樊磊
李月臣
常中兵
李萍
冉林鑫
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southwest University
Original Assignee
Southwest University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Southwest University filed Critical Southwest University
Priority to CN202211664301.5A priority Critical patent/CN116452023B/zh
Publication of CN116452023A publication Critical patent/CN116452023A/zh
Application granted granted Critical
Publication of CN116452023B publication Critical patent/CN116452023B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N33/00—Investigating or analysing materials by specific methods not covered by groups G01N1/00 - G01N31/00
    • G01N33/0098—Plants or trees
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N33/00—Investigating or analysing materials by specific methods not covered by groups G01N1/00 - G01N31/00
    • G01N33/18—Water
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/02—Systems using the reflection of electromagnetic waves other than radio waves
    • G01S17/06—Systems determining position data of a target
    • G01S17/46—Indirect determination of position data
    • 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
    • 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/764—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using classification, e.g. of video objects
    • G06V10/765—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using classification, e.g. of video objects using rules for classification or partitioning the feature space
    • 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/77—Processing image or video features in feature spaces; using data integration or data reduction, e.g. principal component analysis [PCA] or independent component analysis [ICA] or self-organising maps [SOM]; Blind source separation
    • G06V10/774—Generating sets of training patterns; Bootstrap methods, e.g. bagging or boosting
    • 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)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • Theoretical Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Human Resources & Organizations (AREA)
  • Chemical & Material Sciences (AREA)
  • Strategic Management (AREA)
  • Economics (AREA)
  • Multimedia (AREA)
  • Medicinal Chemistry (AREA)
  • Educational Administration (AREA)
  • Electromagnetism (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Marketing (AREA)
  • Immunology (AREA)
  • Biochemistry (AREA)
  • Analytical Chemistry (AREA)
  • Pathology (AREA)
  • General Business, Economics & Management (AREA)
  • Development Economics (AREA)
  • Food Science & Technology (AREA)
  • Tourism & Hospitality (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Computing Systems (AREA)
  • Databases & Information Systems (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Quality & Reliability (AREA)
  • Wood Science & Technology (AREA)
  • Botany (AREA)
  • Primary Health Care (AREA)
  • Mining & Mineral Resources (AREA)
  • Marine Sciences & Fisheries (AREA)
  • Animal Husbandry (AREA)

Abstract

本发明提出了一种基于低频微波雷达VOD数据的公里级碳储量评估方法，包括以下步骤：S1，采集数据：基于低频微波卫星获取的植被含水量数据、基于激光雷达卫星获取的树高数据、基于哨兵1号获取的后向散射数据、基于MODIS卫星获取的植被指数和森林覆盖度数据、基于多源遥感卫星数据合成的植被碳储量参考底图和基于地面观测获取的样地生物量数据；S2，构建高分辨率植被含水量降尺度模型；S3，采用经验函数拟合植被含水量与碳储量底图，构建基于公里级植被含水量的碳储量动态估算模型，实现植被碳储量动态监测。本发明能对植被地上碳储量的时序动态变化进行监测；且能辨识景观破碎度高、生态修复区离散、地形复杂等特点区域的植被地上碳储量的空间异质性，从而能科学评估空间生态修复碳汇能力。

Description

基于低频微波雷达VOD数据的公里级碳储量评估方法
技术领域
本发明涉及自动化遥感信息提取相关技术领域，特别是涉及一种基于低频微波雷达VOD数据的公里级碳储量评估方法。
背景技术
当前，森林清查手段被认为最可靠的碳储量观测手段，但其代表性低且需要耗费大量的人力、物力，并不适用于全区域的碳储量动态估算。遥感技术具有监测范围广、成本低的特点，在大尺度地上碳储量估算中具有不可替代的优势。现有的遥感碳储量估算手段主要包括：(1)激光雷达(LiDAR)手段被认为最有效的估算手段。但是过去十年，无任何面向碳储量的在轨激光雷达卫星，导致该手段无法应用于碳储量的动态监测；(2)光学遥感手段具有较高的时空分辨率，但是多云雨的天气导致光学遥感数据缺失严重，并且光学遥感对中、高植被区域的碳储量变化不敏感。(3)主动微波遥感(SAR)具备了全天侯工作的优势，但是其对高植被区域的碳储量变化存在明显的饱和现象。
近年来，随着低频被动微波卫星的发展，由于其与碳储量的物理联系明确、不受云雨天气干扰可全天候工作、具有较强的穿透能力而对中、高碳储量变化敏感、时间分辨率高等的优势，已广泛应用于全球碳储量的动态监测中，并已服务于全球碳循环实时监测系统。但是现有的低频被动微波数据存在空间分辨率低(仅25公里)的问题，导致其无法应用于具有碳汇分布具有景观破碎度高、生态修复区离散、地形复杂等特点的区域。
发明内容
本发明旨在至少解决现有技术中存在的技术问题，特别创新地提出了一种基于低频微波雷达VOD数据的公里级碳储量评估方法。
为了实现本发明的上述目的，本发明提供了一种基于低频微波雷达VOD数据的公里级碳储量评估方法，包括以下步骤：
S1，采集数据：基于低频微波卫星获取的植被含水量数据、基于激光雷达卫星获取的树高数据、基于哨兵1号获取的后向散射数据、基于MODIS卫星获取的植被指数和森林覆盖度数据、基于多源遥感卫星数据合成的植被碳储量参考底图和基于地面观测获取的样地生物量数据；
S2，构建高分辨率植被含水量降尺度模型；
S2-1，构建高分辨率植被含水量先验知识库：
通过哨兵数据反演的植被含水量(VOD)以及树高数据(FH)，得到高分辨率植被含水量空间分布信息；
将所述高分辨率植被含水量空间分布信息结合MODIS叶面积指数和MODIS植被覆盖度，构建高分辨率植被含水量先验知识库；
S2-2，将粗分辨率的植被含水量的月观测数据、高分辨率植被含水量先验知识库输入高分辨率植被含水量降尺度模型，对粗分辨率的植被含水量的月观测数据进行降尺度，通过所述模型获取公里级植被含水量月观测数据；
S3，采用经验函数拟合植被含水量与碳储量底图，构建基于公里级植被含水量的碳储量动态估算模型，实现植被碳储量动态监测；
所述公里级植被含水量通过高分辨率植被含水量降尺度模型获得。
进一步地，所述高分辨率植被含水量降尺度模型表示为：
VODHigh＝f(VOD,FH,LAI,TC)
其中，VODHigh为基于随机森林模型估算植被含水量；
f()为随机森林算法；
VOD为哨兵数据反演的植被含水量；
FH为树高；
LAI为MODIS叶面积指数；
TC为MODIS植被覆盖度。
采用自助抽样法对所述随机森林模型进行训练：
首先，从原始数据集中利用bootstrapping方法有放回的随机抽取k个新样本集，并用新样本构建K颗分类回归树；其次，假设有n个特征，在每棵树的每个节点选择Mtry个特征，其中Mtry<n，计算每个特征的信息熵，并通过概率值大小选择分类能力最强的特征进行节点的分裂；最后，直接将生成的多棵树组成随机森林对新数据进行分类或回归，分类结果采用简单多数投票机制，最终的分类决策方式如下：
其中，H(x)标识组合分类算法，hi表示单个决策树分类算法，Y为目标变量。
随机森林模型性能根据大数定律可得到泛化误差的收敛上界：
其中k表示随机森林中树的数目。随机森林中树增多的同时，模型泛化误差将会逐渐趋于一个上面公式的上界。
进一步地，在训练高分辨率植被含水量降尺度模型时，将模型中的参数Ntree和Mtry的值分别设定为500和2，其中Ntree为随机森林所包含的决策树数目，Mtry为每个决策树包含的节点个数。
一般Ntree默认为500，Mtry默认为logN。当Ntree设定为500时，OOB误差变化趋于稳定，故本Ntree值设定为500；而随着Mtry值增大，OOB误差随之增大，故Mtry值设定为2。通过对Ntree和Mtry的值进行合理设定，有助于提高模型预测精度。
进一步地，所述经验函数为：
当监测区域为热带地区时，采用线性回归模型或arctangent回归模型作为经验函数；
采用线性回归模型来拟合高分辨率植被含水量和碳储量之间的关系为：
AGC＝a*VOD+b
采用arctangent回归模型来拟合高分辨率植被含水量和碳储量之间的关系为：
式中，VOD为高分辨率植被含水量；
a、b、c、d是回归参数；
当监测区域为温带地区时，采用logistic回归模型作为经验函数；
采用logistic回归模型来拟合高分辨率植被含水量和碳储量之间的关系为：
式中，AGC为植被碳储量；
VOD为高分辨率植被含水量；
a、b、c、d是回归参数。
进一步地，在采集数据后对植被含水量数据进行预处理：
剔除植被含水量模型反演误差TB-RMSE大于8K的植被含水量数据，用于剔除受RFI影响的植被含水量数据；
进一步地，对植被含水量数据进行预处理还包括：
利用土地利用数据剔除城镇、裸地和水体等非植被比例大于10％的升、降轨植被含水量数据，所述土地利用数据包括地区的河流、湖泊、湿地中的一种或者任意组合；用于避免非植被地表类型对植被含水量观测造成的影响。
进一步地，还包括：利用模型反演误差TB-RMSE对升、降轨最优植被含水量观测进行融合，解决升、降轨植被含水量数据受地形及射频干扰影响程度不同的问题；
(1)首先分别计算每个月内升、降轨植被含水量数据及其模型反演误差TB-RMSE的均值，如果同月内升轨和降轨植被含水量数据均值差大于0.05，且升轨或降轨植被含水量数据的TB-RMSE均值大于5K，则移除此季度内的所有升轨或降轨植被含水量数据；
(2)如果一天中同时存在升轨和降轨数据，则通过比较RMSE-TB来选择升轨和降轨中质量较高的观测作为该日的植被含水量观测，从而融合同一天的升、降轨数据，生成质量较高的植被含水量日观测数据；
(3)计算每个季度内的融合后植被含水量日观测数据的均值mean和标准差STD，如果植被含水量日观测数据数值超出mean±2STD，则移除该植被含水量日观测数据。
通过以上过滤和融合步骤，有效剔除了地区地形和射频干扰对植被含水量数据的影响，选择了最优的升、降轨数据获取质量较高的植被含水量日观测数据，从而融合生成了高质量的植被含水量日合成数据。
进一步地，所述高分辨率植被含水量降尺度模型中，需要采用VOD-碳储量反演模型将植被含水量转为碳储量，VOD-碳储量反演模型的构建步骤如下：
首先将水云模型和Ulaby模型耦合得到初步土壤-植被散射模型，再结合地区的土壤水分数据、后向散射数据，构建地区的土壤-植被散射模型；
然后基于土壤-植被散射模型对哨兵1号后向散射数据进行土壤-植被含水量信息分离，获取公里级植被含水量；
土壤-植被散射模型的方程如下：
其中，表示雷达接收到的总体后向散射系数；
A是模型的经验系数，与植被的类型以及雷达参数有关；
V1是对植被冠层的描述量；
VOD表示植被含水量；
θ表示入射角；
SM表示土壤水分；
C代表干燥土壤的后向散射系数；
D代表雷达对土壤水分含量的敏感性；
根据所述土壤-植被散射模型反演出基于哨兵数据的高分辨率植被含水量VOD，得到VOD-碳储量反演模型：
所述高分辨率为小于等于1km的空间分辨率。
反演的后向散射数据可通过土壤-植被散射模型生成植被含水量数据，公里级植被含水量辅助信息表征了植被含水量的时空变化特征。
进一步地，在采集数据后对MODIS叶面积指数进行预处理：
(1)借助MODIS数据质量控制文件，剔除受到云雨天气干扰的像元；
(2)通过最大值合成法将叶面积指数数据生产为年尺度数据；
(3)对年数据进行拼接、重投影和裁剪预处理，获取地区MODIS叶面积指数数据；
(4)将MODIS叶面积指数进行重采样，使其与植被碳储量底图具有相同的投影和空间分辨率。
进一步地，在采集数据后对地面样地观测数据的筛选及质量控制：筛选出森林覆盖率的标准差σTC<15％的数据。
综上所述，由于采用了上述技术方案，本发明具有以下优点：
1.对植被地上碳储量的时序动态变化进行监测；
2.能辨识景观破碎度高、生态修复区离散、地形复杂等特点区域的植被地上碳储量的空间异质性，从而能科学评估空间生态修复碳汇能力。
本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
图1是VOD升轨观测中RFI干扰分布，图1中的(a)为该像元检测到的RFI时间序列；图1中的(b)为该像元的VOD受RFI影响的时间序列。
图2是VOD降轨观测中RFI干扰分布，图2中的(a)为该像元检测到的RFI时间序列；图2中的(b)为该像元的VOD受RFI影响的时间序列
图3是树高数据与三套独立数据验证结果，图3中的(a)为超过100万个的GEDI脚点数据；图3中的(b)为33平方千米的无人机激光雷达数据；图3中的(c)为近6万条的森林清查数据。
图4是预处理后的哨兵1号后向散射数据，图4中的(a)为基于哨兵1号卫星的后向散射原始数据；图4中的(b)为预处理后数据。
图5是MODIS叶面积指数的观测质量分布状况示意图。
图6是通过谷歌影像目视解译筛选地面观测样点示意图。
图7是高质量植被含水量数据质量控制与融合流程图。
图8是植被含水量升、降轨融合数据展示图。
图9是重庆地区高质量融合植被含水量示意图。
图10是基于哨兵数据的高分辨率植被含水量辅助数据示意图。
图11是高分辨率植被含水量降尺度流程图。
图12是随机森林模型变量重要性示意图。
图13是随机森林模型的交叉验证精度示意图。
图14是重庆地区2015-2021年公里级植被含水量示意图。
图15是基于植被含水量的公里级植被碳储量估算流程图。
图16是基于不同碳储量底图构建的植被含水量与碳储量拟合关系示意图。
图17是基于不同碳储量底图反演的公里级植被碳储量示意图。
图18是重庆地区2015年公里级植被碳储量示意图。
图19是基于地面样点的公里级植被碳储量数据精度验证示意图。
图20是基于碳储量参考底图的公里级植被碳储量数据交叉验证示意图。
图21是公里级植被碳储量数据与植被指数显著性(P<0.05)相关系数空间分布。
具体实施方式
下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。
本发明基于低频微波遥感的碳储量动态监测原理，以低频被动微波植被含水量为基础，借助多源高分辨率植被含水量空间分布辅助信息，得到公里级植被地上碳储量动态监测方法。
1.数据收集与预处理
1.1数据信息
数据主要包括基于低频微波卫星获取的植被含水量、基于激光雷达卫星获取的树高、基于主动微波(哨兵1号)卫星获取的后向散射、基于MODIS卫星获取的植被指数和森林覆盖度、基于多源遥感卫星数据合成的植被碳储量参考底图和基于地面观测获取的样地生物量数据等，如表1所示。
表1
数据集 产品/卫星传感器 时间 空间分辨率
低频植被含水量 SMOS VOD 2015–2021 25km
哨兵后向散射 哨兵1号 2015–2021 10m
激光雷达树高 激光雷达 2019 30m
MODIS叶面积指数 MOD15A2 2015–2021 1km
MODIS植被覆盖度 MOD44B 2015–2021 1km
植被碳储量底图 文献收集 2015 1km
植被碳储量底图 多源数据融合 2015 1km
地面样地数据 野外调查 2015 85个
基于低频被动微波卫星的植被含水量产品：植被含水量(VOD)的空间分辨率为25km，时间分辨率为1–3天。该产品基于欧空局的SMOS卫星L波段(1.4GHz)的升轨和降轨亮温数据，采用SMOS-IC算法反演生成。由于SMOS-IC算法在反演时不需要过多的辅助数据(如再分析土壤水分数据和光学植被指数数据等)输入，因此独立于其他植被指数和植被碳储量产品，该数据集提供了全球逐日的植被含水量和土壤湿度数据。
该VOD产品具有全天候、全天时、不受云雨干扰的连续观测能力，且在中高植被区对碳储量不存在明显的饱和现象(饱和点上限高于100Mg C/ha)。与光学植被指数相比，微波植被含水量对植被的叶和茎组分都比较敏感，而且不易受到大气气溶胶和云层的影响。光学植被指数仅反映植被冠层顶部的状况，而微波植被含水量由于具有较强的穿透能力，其不仅代表叶片的理化状况，也包含了树干的含水量和结构信息，而且不易出现信号饱和现象。基于以上优势，微波植被含水量已被广泛应用于全球不同区域的碳储量动态估算研究中，包括热带地区、欧洲地区、亚马逊地区以及我国西南森林地区等。
基于激光雷达卫星的树高产品：该树高产品基于全球生态系统动力学调查卫星(GEDI)和冰、云和陆地高程卫星(ICESat-2)的全波形激光雷达数据，结合高分辨率植被、环境和地形数据，通过深度学习框架下的空间内插模型生成，空间分辨率为30m，时间为2019年。
基于主动微波哨兵1号卫星的后向散射产品：主动微波卫星的后向散射观测包括植被和地表的直接散射、经植被衰减的地表散射、植被内部的体散射、植被与地表相互作用的二次散射等。其中经植被衰减的地表散射与植被含水量密切相关，相关研究已经证明了主动微波卫星的后向散射数据反演植被含水量的潜力。
本发明使用了基于C波段的合成孔径雷达哨兵1号卫星(Sentinel-1)提供的后向散射产品(干涉宽幅模式地距多视Level-1级数据)，空间分辨率为10m，时间分辨率为6天。哨兵1号卫星凭借其丰富的观测数据、高重访周期(6/12天)、高分辨率(10m)的优势，可为高分辨率植被含水量反演提供辅助信息。
本发明所选用的哨兵1号数据是干涉宽幅模式地距(GRD)多视Level-1级数据产品，其相对于SLC数据而言，GRD产品是经过多视处理并投影至地距的聚焦数据，其像元灰度值代表幅度信息，具有相同的方位向和距离向分辨率，降低了斑点噪声和几何分辨率，图像质量更高。
基于MODIS卫星的植被指数与植被覆盖度：基于MODIS卫星的叶面积指数是森林生态系统的一个重要结构参数，表征叶片的疏密程度和冠层结构特征，影响着植被冠层内的光合、呼吸和蒸腾作用等生理生化过程，是估算植被生态过程的重要变量。基于MODIS卫星的植被覆盖度产品代表了植被冠层覆盖度，与植被光合作用、固碳过程密切相关。因此叶面积及植被覆盖度均与植被碳储量变化有密切联系，可用于间接验证的植被碳储量变化。
基于MODIS卫星的叶面积指数(Leaf Area Index，LAI)用于本发明的版本为MOD15A2，空间分辨率为500m，提供了全球每8天尺度的叶面积指数数据。
基于MODIS卫星的植被覆盖度(Percent Tree Cover，TC)数据由MOD44B植被覆盖度产品(Vegetation Continuous Fields，VCF)提供，空间分辨率为250m，时间分辨率为年尺度。植被覆盖度代表了植被冠层覆盖率的百分比，定义为高度大于或等于5m的树冠遮挡的比例。
重庆植被碳储量参考底图：6种已在国际顶级期刊发布的植被地上碳储量静态产品(Aboveground Carbon,AGC)作为本发明的植被碳储量底图，用于本发明中VOD-植被碳储量模型的构建。这6种植被碳储量底图分别来自Saatchi et.al.(2011)、Baccini et.al.(2012)、Santoro et.al.(2021)、Su et.al.(2016)和Chang et.al.(2021)，以下分别称为“Saatchi”、“Baccini”、“CCI”、“Su”、“Saatchi-WT”和“Saatchi-RF”数据。上述产品的空间分辨率均为1km，产品单位为生物量密度(Mg/ha)，误差小于15％。
上述产品的算法均是基于ICESat卫星的激光雷达(LiDAR)点云反演的植被冠层结构参数，结合气候、地形和光学(MODIS、Landsat等)遥感产品、地面观测数据，在异速生长方程或其他经验模型的框架下，反演获取地上生物量。
其中，“Saatchi-WT”和“Saatchi-RF”植被碳储量底图由是基于线性加权和随机森林回归两种数据融合算法，通过融合多源数据(森林清查数据、样地调查数据、遥感反演数据)生产。该产品在中国具有较高的精度(误差小于12％)，与地面观测和森林清查数据具有较高的一致性。
基于IPCC地上生物量密度与地上碳储量密度转换关系，将上述生物量底图(单位为Mg/ha)，乘以转换系数(0.5)，反演为碳储量底图(单位为Mg C/ha)。
重庆地面碳储量地面样地数据：本发明采用的重庆地区地面碳储量样地由本团队自主获取，数据采集时间为2015年。重庆地区共布设具有代表性的森林样地85个，每个样地面积为1000m2。
为确保地面样地数据能充分代表重庆地区森林资源状况，在布设样地时结合重庆森林资源清查数据和植被类型分布图，充分考虑重庆地区森林类型特点(包括森林类型、林龄和起源等)、各类森林的面积、蓄积构成及其地域分布权重，设置各类森林的调查样地数，并确定所设样地的空间分布。
在野外观测中，每个调查样点内设置3个重复样地，样地间隔至少100m。在每个样地中，样地面积为1000m2，将样地进一步划分成10个(10m×10m)样方，然后以样方为单元进行林木调查、测定与记录。在每个样方内，测量整个样方内所有胸径≥5cm的树木，记录胸径、树高等信息。结合重庆地区样地优势树种，选取与样地气候条件相似，同种或同科属树种的生物量方程，计算各树种生物量。
W＝0.0548(D2H)0.8545 (2.1)
W＝0.0414(D2H)0.9345 (2.2)
W＝0.0461(D2H)0.6109 (2.3)
式中，W为林木生物量(kg)；D为林木胸径(cm)；H为林木树高(m)。式2.1、2.2和2.3分别代表重庆地区马尾松、石栎和栓皮栎3种树种的生物量方程。
采集优势树种干、枝、叶、根样品，测定植物组织碳含量，计算样方尺度的碳密度。碳密度计算公式如下：
式中，表示植物组织碳库(gCm-2)，OR∈(干，枝，叶，根)，fi,j(DBH,H)表示i株j植被(T表示树种，S表示林下植被)植物器官生物量异速方程(g干物质量/株)；CFi,j表示i株j植被器官碳含量(g kg-1)，AP为样地面积(m2)。
1.2数据处理：
1.2.1低频被动微波植被含水量的质量控制
利用误差TB-RMSE定量分析植被含水量产品受电磁干扰(RFI)影响的程度，TB-RMSE为原数据本身提供的一个表征不确定性的数据。根据已有研究结果和重庆市数据特征，将TB-RMSE大于8K的数据定义为受RFI干扰的VOD观测。利用该条件，分别对VOD升轨、降轨的日观测数据开展质量控制，剔除受RFI影响的VOD数据。
以重庆森林地区的VOD升轨观测为例，图1中的(a)显示了该像元在2016年检测到的射频干扰(RFI)的分布，其中第200-300天期间遭受了RFI干扰(TB-RMSE>8K)，因此对该像元的开展质量控制，对该段时间内受到RFI影响的观测进行剔除，从而获得高质量的VOD观测，如图1中的(b)。
图2中的(a)显示了重庆地区2016年灌木VOD降轨观测的射频干扰(RFI)情况，图2中的(b)显示了该观测对应的VOD观测受RFI影响的时间序列分布。利用本发明建立的质量控制标准，可对该像元受到RFI影响的VOD观测进行有效剔除。
此外，针对非植被地表类型对植被含水量(VOD)观测造成的影响，本发明利用重庆市河流、湖泊和湿地等土地利用数据，剔除了城镇、裸地和水体等非植被比例大于10％的VOD升、降轨数据。
1.2.2基于激光雷达卫星的树高的精度验证。
树高产品来源于Liu et.al.，其估算的全国森林冠层平均高度为15.90m。通过与三套独立的验证数据集(超过100万个的GEDI脚点数据、33km2的无人机激光雷达数据和近6万条的森林清查数据)相比，该产品具有较高的精度可满足本发明的应用需求，具体如图3所示：三套数据验证精度分别为R2＝0.55,RMSE＝5.32m；R2＝0.58,RMSE＝4.93m；R2＝0.60,RMSE＝4.88m；。
1.2.3哨兵1号后向散射的预处理。
基于SNAP数据处理平台，本发明对哨兵1号后向散射数据，开展数据预处理：
(1)轨道校正：轨道偏差纠正；
(2)辐射定标：辐射定标将传感器输出的信号(DN值)通过定量关系式转化为对应的地表辐射亮度值；
(3)多视处理、滤波处理：去除影像斑点噪声；
(4)地形校正：消除图像的几何畸变，引入30m分辨率的数字高程模型(DEM)作为地形参照文件开展地形纠正，从而得到高质量的后向散射数据，如图4所示。
1.2.4MODIS叶面积指数的数据预处理。
(1)针对MODIS叶面积指数产品在重庆地区易受云雨天气而造成数据质量低的问题，借助MODIS数据质量控制(QA)文件，剔除了受到云雨天气干扰的像元，QA＝0表示数据质量优，QA>0表示像元数据质量存在问题，观测质量分布状况如图5所示。
(2)为减少云、气溶胶、视角等因素的影响，通过最大值合成法将叶面积指数数据生产为年尺度数据。
(3)对每年的MODIS叶面积指数进行拼接、重投影和裁剪等几何纠正预处理，获取重庆地区MODIS叶面积指数产品。
(4)将MODIS叶面积指数重采样为1km，使其与植被碳储量底图具有相同的投影和空间分辨率。
1.2.5地面观测数据质量控制。
对地面样地观测数据的筛选及质量控制流程如下：
(1)根据实地探勘，选择样点未受到严重的人为干扰，且在1km像元具有代表性；
(2)为减少地面观测样点(约1ha)和遥感像元(1km)之间的空间尺度不匹配造成的误差，选择位于均质像元(homogeneous pixel)内的样点。均质像元是指遥感像元内30m分辨率的森林覆盖率的标准差小于15％(σTC<15％)，其中15％的阈值通过高分Google影像的目视解释来确定，如图6所示。
1.2.6对上述图像均执行以下操作：
首先获取同一时刻拍摄同一位置的K张图像，K为大于或者等于10的整数，对同一时刻拍摄同一位置的K张图像进行融合，执行以下步骤：
(1)对每张图像均执行以下操作：
其中，Piexmn表示第m行第n列像素点的颜色值；m＝1、2、3、……、M，M表示纵向像素点个数；n＝1、2、3、……、N，N表示横向像素点个数；
Piex(m-1)(n-1)表示第m-1行第n-1列像素点的颜色值；
Piex(m-1)n表示第m-1行第n列像素点的颜色值；
Piex(m-1)(n+1)表示第m-1行第n+1列像素点的颜色值；
Piexm(n-1)表示第m行第n-1列像素点的颜色值；
Piexm(n+1)表示第m行第n+1列像素点的颜色值；
Piex(m+1)(n-1)表示第m+1行第n-1列像素点的颜色值；
Piex(m+1)n表示第m+1行第n列像素点的颜色值；
Piex(m+1)(n+1)表示第m+1行第n+1列像素点的颜色值；
Piex′mn表示处理后的第m行第n列像素点的颜色值；
Piex(m-1)(n-1)、Piex(m-1)n、Piex(m-1)(n+1)、Piexm(n-1)、Piexm(n+1)、Piex(m+1)(n-1)、Piex(m+1)n、Piex(m+1)(n+1)是Piexmn的八个相临近像素点的颜色值；由于Piex00、Piex01、Piex02、Piex03、……、Piex0(N+1)，Piex00、Piex10、Piex20、Piex30、……、Piex(M+1)0，Piex0(N+1)、Piex1(N+1)、Piex2(N+1)、Piex3(N+1)、……、Piex(M+1)(N+1)，Piex(M+1)0、Piex(M+1)1、Piex(M+1)2、Piex(M+1)3、……、Piex(M+1)(N+1)，即是这些像素点不存在，因此不计算Piex11、Piex12、Piex13、……、Piex1N，Piex11、Piex21、Piex31、……、PiexM1，Piex1N、Piex2N、Piex3N、……、PiexMN，PiexM1、PiexM2、PiexM3、……、PiexMN，即是这些像素点。只计算以下像素点：
p1～p12表示像素点颜色值调节系数；优选的，p1＝p4＝p7＝p10，p2＝p5＝p8＝p11，p3＝p6＝p9＝p12，更优的，p1＝p4＝p7＝p10＝p3＝p6＝p9＝p12＝1，p2＝p5＝p8＝p11＝2，此时：
对每张图像执行以上操作后，得到处理后的图像；
(2)对每张处理后的图像执行以下操作：
其中，Piex′m′n′,1表示处理后的第1张图像中的第m′行第n′列像素点的颜色值；m′＝1、2、3、……、M，M表示纵向像素点个数；n′＝1、2、3、……、N，N表示横向像素点个数；
Piex′m′n′,2表示处理后的第2张图像中的第m′行第n′列像素点的颜色值；
Piex′m′n′,3表示处理后的第3张图像中的第m′行第n′列像素点的颜色值；
Piex′m′n′,K表示处理后的第K张图像中的第m′行第n′列像素点的颜色值；
K表示同一时刻拍摄图像的张数；
Piex″m′n′表示第一图像中的第m′行第n′列像素点的颜色值；
对每张处理后的图像执行以上操作后，得到第一图像；
(3)
其中，Piex′m″n″,k表示处理后的第k张图像中的第m″行第n″列像素点的颜色值；m″＝1、2、3、……、M，M表示纵向像素点个数；n″＝1、2、3、……、N，N表示横向像素点个数；k＝1、2、3、……、K；
Piex″m″n″表示第一图像中的第m″行第n″列像素点的颜色值；
Piexk表示第k张图像的颜色值；
执行以上操作后，得到K张图像的颜色值；
(4)将Piex1、Piex2、Piex3、……、PiexK按照从小到大的顺序依次排列，得到：Piex′1、Piex′2、Piex′3、……、Piex′K，Piex1表示第1张图像的颜色值；Piex2表示第2张图像的颜色值；Piex3表示第3张图像的颜色值；PiexK表示第K张图像的颜色值；Piex′1表示Piex1、Piex2、Piex3、……、PiexK按照从小到大的顺序依次排列后位于第1个(首位)，Piex′2表示Piex1、Piex2、Piex3、……、PiexK按照从小到大的顺序依次排列后位于第2个，Piex′3表示Piex1、Piex2、Piex3、……、PiexK按照从小到大的顺序依次排列后位于第3个，Piex′K表示Piex1、Piex2、Piex3、……、PiexK按照从小到大的顺序依次排列后位于第K个(末位)；
(5)取排序后的前k′个所对应的处理后的图像，k′为大于或者等于2且小于或者等于7，优选的k′＝2，执行以下操作：
表示Piex′k″所对应的处理后的图像；m″′＝1、2、3、……、M，M表示纵向像素点个数；n″′＝1、2、3、……、N，N表示横向像素点个数；k′＝2、3、4……、7；
Piex″′m″′n″′表示融合后的第m″′行第n″′列像素点的颜色值；
执行以上操作后，得到图像的最终颜色值；该颜色值为灰度值，若颜色值为RGB时，将RGB转换为灰度值，其转换公式为：
RGB灰度值＝r*R`m`n+g*G`m`n+b*B`m`n，
其中，R`m`n表示第`m行第`n列像素点的红色值；`m＝1、2、3、……、M，M表示纵向像素点个数；`n＝1、2、3、……、N，N表示横向像素点个数；
G`m`n表示第`m行第`n列像素点的绿色值；
B`m`n表示第`m行第`n列像素点的蓝色值；
r表示红色值系数；r+g+b＝1；
g表示绿色值系数；
b表示蓝色值系数。
综上所述，对低频微波植被含水量产品进行了质量控制，剔除了电磁干扰、地形和水体等因素的影响；对哨兵1号后向散射产品进行了辐射定标、多视处理、滤波处理、地形校正等预处理；对MODIS叶面积指数产品进行了质量控制，剔除了受云雨影响的观测；对地面样地观测数据进行了筛选和质量控制，以保障用于本发明的地面样地在公里级具有代表性。通过数据预处理和质量控制过程保证了数据的质量和精度，为植被含水量降尺度模型和碳储量估算模型提供了可靠的数据源。
2.植被地上碳储量遥感监测模型构建
2.1构建低频被动微波植被含水量高质量数据集
低频微波VOD产品在重庆地区主要面临两方面的问题导致其存在较大的不确定性：<1>在重庆地区受到电磁干扰(RFI)、水体、非植被像元的影响，导致植被含水量产品在部分地区存在观测异常，造成高、低估植被含水量的问题。该问题在数据预处理阶段已解决。<2>由于被动微波卫星天线孔径的倾斜，从升、降轨道获取的植被含水量日观测产品受重庆地形及射频电磁干扰的情况不一致，导致植被含水量升、降轨数据的质量不同。针对上述问题，提出植被含水量融合技术，从而获取高质量的植被含水量日观测数据集。具体流程如图7所示：
(1)针对升、降轨植被含水量数据受地形及射频干扰影响程度不同的问题，利用模型反演误差(TB-RMSE)对升、降轨最优VOD观测进行融合。首先分别计算每个月内升、降轨VOD产品及其模型反演误差(TB-RMSE)的均值。如果同月内升轨和降轨植被含水量产品均值差大于0.05(表明其中一轨存在观测误差，因为同月内植被含水量变化不会超过0.05)，且升轨或降轨植被含水量产品的TB-RMSE均值大于5K(说明该像元持续受到RFI的影响)，则移除此季度内的所有升轨或降轨植被含水量产品。
(2)在此基础上，如果一天中同时存在升轨和降轨数据，则通过比较RMSE-TB来选择升轨和降轨中质量较高的观测作为该日的VOD观测(RMSE-TB越低，VOD观测质量越高)，从而融合同一天的升、降轨数据，生成质量较高的植被含水量日观测产品。
(3)计算每个季度(3个月)内的融合后植被含水量日观测产品的均值(mean)和标准差(standard deviation，STD)。如果植被含水量日观测产品数值超出mean±2STD(即为系统观测误差造成的异常值)，则移除该植被含水量日观测产品。
通过以上过滤和融合步骤，有效剔除了重庆市地形和射频干扰对植被含水量产品的影响，选择了最优的升、降轨数据获取质量较高的植被含水量日观测产品，从而融合生成了高质量的植被含水量日合成数据。图8展示了重庆灌木地区的单一像元经过上述流程后的VOD融合时间序列，ASC和DESC分别代表升、降轨VOD观测；图9为经过上述流程获取的2015年重庆地区VOD空间分布)。
2.2基于哨兵1号的高分辨率植被含水量辅助数据获取
基于水云模型(Water Cloud Model，WCM)，根据重庆地区植被覆盖的实际情况，结合重庆地区土壤水分数据、哨兵1号后向散射数据，利用Ulaby模型计算土壤后向散射贡献，构建适宜于重庆地区的土壤-植被散射模型。
2.2.1重庆地区土壤-植被散射模型构建
在植被覆盖地表情况下，后向散射主要来自于土壤和植被两部分。水云模型将植被层假定为一个各向同性散射体，把地表后向散射描述为植被直接散射与经过植被双程衰减的土壤散射之和。水云模型可以比较准确的描述植被覆盖地区的雷达微波散射机制，为得到植被覆盖地区的土壤后相散射系数的真实值提供了切实可行的方案。由于水云模型简洁的表达方式以及合理的微波散射机制假设，被广泛应用于地表上层生物量反演、植被覆盖下土壤水分反演等研究。水云模型可以表示为：
其中，分别代表雷达接收到的总体后向散射系数、植被层的后向散射系数和土壤层的后向散射系数；τ2代表植被层的双向衰减因子(透过率)；θ是信号入射角；A是模型的经验系数，代表全植被覆盖下冠层的后向散射系数，与植被的类型以及雷达参数有关；V1是对植被冠层的描述量，通常设为1。
在植被茂密的区域(very dense vegetation regions，VDV)，植被冠层的透过率可以近似为零(τ2＝0)，这时土壤的散射项可以忽略不计。则公式3.1可以简化为：
其中VDV像元的A(以下称为A0)可以通过来计算，计算公式为：
选取第t天所有VDV的像元，取其前95％处的A0值来获取第t天像元尺度的A(假设第t天所有像元的A0(t)值一样)。
在上述过程中VDV像元的选取采用NDVI阈值法：
NDVIVDV>NDVI75％(t) (3.7)
式中，NDVIVDV为VDV的NDVI的阈值，取其像元NDVI时间序列前75％处的值。
2.2.2重庆地区土壤散射模型构建
我们首先选择可以忽略植被影响的区域构建土壤散射模型。
(1)选择裸土像元
结合MODIS的NDVI数据和ESA CCI的土地利用数据来确定重庆地区的裸土像元和裸土期。定义为该像元的土地利用类型为裸地或者NDVI<0.1。
(2)裸地像元土壤模型参数的计算
本发明选择Ulaby模型来构建土壤散射模型，Ulaby模型不依赖地表粗糙度参数的影响，可以直接建立土壤水分与后向散射系数之间的经验关系：
其中，为裸露地表后向散射系数(dB)，C和D是Ulaby模型待求解的参数，C代表干燥土壤的后向散射系数，D代表雷达对土壤水分含量的敏感性
在对C和D的求解过程中存在两种情况：一是针对后向散射系数与土壤水分SM呈现显著性相关的像元，可通过/>和SM之间的线性关系，来计算像元尺度的C和D。计算过程中须保证正相关性(R>0)，且p-value<0.05；二是针对后向散射系数/>与土壤水分SM呈现不显著相关的像元，首先找出土壤水分极低的像元(SM<0.1m3/m3)来直接计算C，即C为该段时期的后向散射系数/>的均值。
(3)重庆地区土壤模型参数的校准
以裸土地区的模型参数(C和D)作为先验知识，利用随机森林回归算法获取重庆地区C和D的预测模型(预测因子包括土壤性质数据和地形数据)，进而基于该模型来实现重庆地区土壤散射模型构建。
2.2.3高分辨率植被含水量辅助信息获取
基于土壤-植被散射模型对哨兵1号后向散射数据进行土壤-植被含水量信息分离，获取公里级植被含水量辅助信息，该信息表征了植被含水量的时空变化特征。
将土壤散射模型代入公式3.1，构建了土壤-植被散射模型，模型的方程如下：
将模型转换成反演VOD的模型(公式(3.11))，将哨兵1号获取的后向散射系数和入射角θ、土壤水分SM、土壤散射模型参数C和D、植被散射参数A，V1(通常等于1)代入到公式3.11中，反演出基于哨兵数据的高分辨率植被含水量VOD，如图10所示。
2.3高分辨率植被含水量降尺度模型
低频微波植被含水量产品其较粗的空间分辨率(25km)，导致单一像元内代表了不同生态过程的混合(如森林砍伐、森林退化和植树造林等)，无法满足国土空间修复中的精细尺度植被碳变化监测的需求。针对上述问题，本发明拟借助哨兵1号获取的植被含水量空间分布辅助信息及树高数据(表征了植被含水量的空间分布)，构建高分辨率植被含水量空间分布的先验知识库，进而利用机器学习算法对粗分辨率的植被含水量的月观测产品进行降尺度，获取重庆市公里级植被含水量月观测产品。具体流程如图11所示：
2.3.1高分辨率植被含水量先验知识库构建
在降尺度过程中，高分辨率辅助信息可以刻画复杂地表条件下植被含水量的空间异质性，因此合理地引入高分辨率辅助信息是准确获取高分辨率VOD产品的关键。
激光雷达信号能够穿透整个植被冠层且不受信号饱和的影响，可提供高分辨率的植被垂直结构信息，已经成为提高森林地上碳储量估算精度的重要手段。哨兵1号卫星具有高重访周期、高分辨率的优势，其反演的后向散射数据可通过土壤-植被散射模型生成植被含水量产品，成为提供高分辨率森林结构辅助信息的理想数据源。因此，本发明借助哨兵1号和激光雷达树高等数据提供的高分辨率植被含水量空间分布信息，结合MODIS叶面积指数产品和植被覆盖度产品，构建高分辨率辅助信息先验知识库，辅助粗分辨率VOD降尺度。
2.3.2高分辨率植被含水量辅助信息获取
(1)模型构建
现有的降尺度方法通常可以分为经验和物理模型，其中基于经验回归的方法结构简单、适应性强，且只需要借助于遥感数据即可实现降尺度，是当前应用最广泛的降尺度方法之一。
在构建经验回归关系时，现有研究通常采用机器学习的方法来解决遥感中的非线性回归问题。其中随机森林模型是一种集合学习方法，其与传统的回归树不同，随机森林模型通过反复迭代训练决策树来最小化损失函数，通过结合多种决策树生成预测模型，因此具有较好的可解释性和鲁棒性。
随机森林算法的基本思想为：首先，从原始数据集中利用bootstrapping方法有放回的随机抽取k个新样本集，并用新样本构建K颗分类回归树；其次，假设有n个特征，在每棵树的每个节点选择Mtry个特征，其中Mtry<n，计算每个特征的信息熵，并通过概率值大小选择分类能力最强的特征进行节点的分裂；最后，直接将生成的多棵树组成随机森林对新数据进行分类或回归，分类结果采用简单多数投票机制，最终的分类决策方式如下：
其中，H(x)标识组合分类算法，hi表示单个决策树分类算法，Y为目标变量。
随机森林模型性能根据大数定律可得到泛化误差的收敛上界：
其中k表示随机森林中树的数目。随机森林中树增多的同时，模型泛化误差将会逐渐趋于一个上面公式的上界。
随机森林由一系列的决策树构成，每棵树都由可放回的随机重复采样(bootstrapping，称为自助抽样法)训练样本生成，使得一些样本可能被多次使用。同时，每颗树在分裂节点时随机选择样本变量，使用所有特征的一个子集来拆分每个决策树中的每个节点。最终，通过加权平均每个决策树的预测作为最终的预测结果。以上所述的样本随机性和特征随机性，使得随机森林不容易出现过拟合问题，具有很好的抗噪声能力，对比其他算法具有一定的优势。
在随机森林模型支持下，在低分辨率(25km)尺度下建立植被含水量月产品与植被含水量辅助信息先验知识库提供的植被含水量的非线性关系，构建植被含水量降尺度模型。
VODHigh＝f(VOD,FH,LAI,TC) (3.14)
式中，VODHigh为基于随机森林模型估算植被含水量，VOD为哨兵数据反演的植被含水量，FH为激光雷达树高产品，LAI为MODIS叶面积指数产品，TC为MODIS植被覆盖度产品。
为提高模型预测精度，在进行随机森林模型训练时，需要对模型中的两个重要参数Ntree和Mtry进行优化。其中Ntree为随机森林所包含的决策树数目，默认为500；Mtry为每个决策树包含的节点个数，默认为logN。
在本研究中，以质量控制后的低频被动微波植被含水量产品为因变量，以高分辨率植被含水量先验知识库中哨兵数据反演的植被含水量(VOD)、激光雷达树高产品(FH)，MODIS叶面积指数产品(LAI)和MODIS植被覆盖度产品(TC)为自变量，采用bootstrapping(自助抽样法)的方法对随机森林模型进行训练。根据袋外数据(out of bag,OOB)误差的平均精度下降程度，来判断该参数的重要性。通过实验发现，当Ntree设定为500时，OOB误差变化趋于稳定，故本研究中Ntree值设定为500；而随着Mtry值增大，OOB误差随之增大，故Mtry值设定为2。
随机森林算法可以对特征变量的重要程度进行分析评价，在提高模型精度的同时减少数据的冗余和处理工作量。图12显示了各个高分辨率辅助数据在构建随机森林模型时的变量重要性，其中重要性得分越高，说明该变量对分类结果的影响和贡献就越大。图中，TC为植被覆盖度，SAR为哨兵植被含水量，FH为激光雷达树高产品，LAI为叶面积指数；从图中可以看出，MODIS植被覆盖度和哨兵植被含水量产品占主要贡献。
(2)模型验证
本发明采用样本外检验(out-of-sample test)的方法评估了植被含水量降尺度模型的精度。首先将所有数据随机分为训练数据(75％的观测数据，即in-sample)和验证数据(25％的观测数据，即out-of-sample)，其中验证数据用来评估高分辨率植被含水量降尺度模型的精度。
高分辨率植被含水量降尺度模型的精度使用以下两个指标进行评估：调整后的决定系数(R2)以及粗分辨率微波植被含水量产品与降尺度模型估算的植被含水量数据之间的均方根误差(RMSE)。精度评价显示如图13所示，本发明构建的高分辨率植被含水量降尺度模型的R2和RMSE分别为0.97和0.07，表明该模型具有满意的精度，可用于生产公里级VOD产品。
2.3.3重庆地区公里级植被含水量产品
基于微波植被含水量观测数据和高分辨率辅助信息先验知识库，借助于高分辨率植被含水量降尺度模型，反演获得重庆市2015–2021年公里级植被含水量产品。从图14可以看出，反演获得的公里级植被含水量产品准确刻画了重庆地区植被生长情况。总体上，植被含水量较高的区域主要分布在重庆东北、东南和四山区域，重庆西北地区和主城区植被含水量较低。
2.4基于高分辨率植被含水量的公里级碳储量动态估算模型
2.4.1高分辨率植被碳储量动态估算模型
以重庆市植被碳储量参考底图数据为响应变量，以高分辨率植被含水量年际产品为预测变量，建立植被含水量与碳储量底图之间的反演关系，构建基于公里级植被含水量的碳储量动态估算模型，并通过VOD-碳储量反演模型将植被含水量VOD转为碳储量，参见公式3.11。高分辨率植被碳储量动态估算模型如图15所示。
已有研究采用不同的经验函数拟合VOD与地上碳储量关系，来实现植被碳储量动态监测，例如线性回归模型、arctangent回归模型和logistic回归模型。研究表明，线性模型和arctangent回归模型均适用于热带地区，logistic回归模型多用于温带地区。因此，在本发明中，采用logistic回归模型来拟合高分辨率植被含水量和碳储量之间的关系：
式中，AGC为植被碳储量，单位为Mg C/ha；VOD为高分辨率植被含水量；a、b、c、d是回归参数，其中a和d的单位是Mg C/ha，b和c是无量纲的量，具体的模型参数参见表2。
本发明利用空间相关性(R2)和均方根误差(RMSE)以评价基于植被含水量的碳储量模型精度，分别表示为：
式中AGCref和AGCestimate分别为参考底图和基于VOD反演的植被碳储量，和分别是参考底图和基于VOD反演的植被碳储量平均值，N是观测数据的总数。
根据图16的结果表明：基于植被含水量的公里级碳储量估算模型均具有较高的精度(R2＝0.62–0.81，RMSE＝9.19–14.82Mg C/ha)。基于R2和RMSE两个评价指标，本团队生产两个参考底图(“Saatchi-WT”和“Saatchi-RF”数据)与VOD具有更高的相关性(R2分别为0.81和0.75，RMSE分别为9.19Mg C/ha和10.57Mg C/ha)。
值得注意的是，本发明使用的6个碳储量参考底图(Saatchi、Baccini、CCI、Saatchi-WT、Saatchi-RF和Su)由于其本身存在不确定性，不同参考碳储量数据之间碳储量的估算也不一致。考虑到不同碳储量参考底图的不确定性。为最大程度减小上述碳储量底图的不确定性，使用所有6个碳储量底图数据来训练碳储量动态估算模型，首先得到6个碳储量估算模型。然后使用这6个碳储量估算模型生成6个碳储量产品，如图17所示；并使用6个碳储量的中值作为最终的2015–2021年植被碳储量年际产品，如图18所示。
表2基于不同碳储量底图的碳储量估算模型拟合参数(式3.15)
AGC a b c d R2 RMSE
Saatchi 111.2 22.6 0.3 3.0 0.73 13.48
CCI 102.0 2.6 0.8 -6.3 0.63 14.38
Saatchi-WT 95.5 19.8 0.3 2.7 0.68 14.82
Saatchi-RF 97.2 20.6 0.3 2.9 0.75 10.57
Su 107.2 19.4 0.3 6.1 0.81 9.19
Baccini 147.5 9.7 0.4 -7.3 0.72 12.61
2.4.2公里级碳储量产品精度评价
直接验证、交叉验证、间接验证等三种验证方法用于验证重庆公里级碳储量年际产品。精度评价指标包括：调整后的决定系数(R2)以及均方根误差(RMSE)。
(1)基于地面样地数据的直接验证
根据地面样地获取的碳储量数据，在公里尺度上对本发明生产的碳储量产品开展直接评价。结果如图18所示：植被碳储量年际产品与地面观测数据具有较高的一致性，其中R2达到了0.60且RMSE仅为10.92Mg C/ha，可满足公里级尺度碳储量的动态监测要求。
基于地面样点的公里级植被碳储量产品精度验证，结果如图19所示，该产品在植被碳密度较高时(大于50Mg C/ha)，并没有明显的饱和现象。证明该产品克服了现有碳储量产品在浓密森林存在的信号饱和问题，可准确监测不同碳密度条件下的植被碳储量空间分布。
(2)基于碳储量底图的交叉验证
基于碳储量参考底图，利用bootstrap随机抽样方法，进行1000次重复随机抽样，其中抽样率＝80％，即80％数据用于训练公里级碳储量模型，其余20％数据用于精度评价。结果如图20所示：本发明建立的模型估算结果与碳储量参考底图呈现更高的一致性，R2和RMSE分别为0.83和6.20Mg C/ha。碳储量估算模型有效克服了信号饱和问题，可监测到碳密度大于100Mg C/ha时碳储量的变化。而且，碳储量估算模型在较低碳密度较低时(低于20MgC/ha)对碳储量也具有较强的拟合能力。
(3)基于光学植被指数的间接验证
开展时间尺度上的公里级植被碳储量产品与光学遥感植被指数产品(如MODIS叶面积指数和植被覆盖度)相关分析，对碳储量产品进行间接评估。其中，Pearson相关系数(R)用于间接评价植被碳储量产品：
式中,VI代表了光学植被指数(MODIS叶面积指数和植被覆盖度)，AGC本发明发展的碳储量年际产品。和/>分别代表了光学植被指数和植被碳储量(AGC)的平均值。p值用来定义相关性的显著性水平，本研究中使用p<0.05表示显著性。
间接验证结果如图21所示：碳储量产品与MODIS叶面积指数LAI和植被覆盖度TC的相关性分别达到0.63和0.67，证明本发明的植被碳储量监测方法能够准备的捕捉植被的动态变化。从R值的空间分布来看，85％的LAI和84％的植被覆盖度变化均与植被碳储量呈一致的趋势(R>0)。
尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。

Claims (5)

1.一种基于低频微波雷达VOD数据的公里级碳储量评估方法，其特征在于，包括以下步骤：
S1，采集数据：基于低频微波卫星获取的植被含水量数据、基于激光雷达卫星获取的树高数据、基于哨兵1号获取的后向散射数据、基于MODIS卫星获取的植被指数和森林覆盖度数据、基于多源遥感卫星数据合成的植被碳储量参考底图和基于地面观测获取的样地生物量数据；
在采集数据后对植被含水量数据进行预处理：剔除植被含水量模型反演误差TB-RMSE大于8K的植被含水量数据，用于剔除受RFI影响的植被含水量数据；
还包括：利用模型反演误差TB-RMSE对升、降轨最优植被含水量观测进行融合，解决升、降轨植被含水量数据受地形及射频干扰影响程度不同的问题；
(1)首先分别计算每个月内升、降轨植被含水量数据及其模型反演误差TB-RMSE的均值，如果同月内升轨和降轨植被含水量数据均值差大于0.05，且升轨或降轨植被含水量数据的TB-RMSE均值大于5K，则移除此季度内的所有升轨或降轨植被含水量数据；
(2)如果一天中同时存在升轨和降轨数据，则通过比较RMSE-TB来选择升轨和降轨中质量较高的观测作为该日的植被含水量观测，从而融合同一天的升、降轨数据，生成质量较高的植被含水量日观测数据；
(3)计算每个季度内的融合后植被含水量日观测数据的均值mean和标准差STD，如果植被含水量日观测数据数值超出mean±2STD，则移除该植被含水量日观测数据；
S2，构建高分辨率植被含水量降尺度模型；
S2-1，构建高分辨率植被含水量先验知识库：
通过哨兵数据反演的植被含水量(VOD)以及树高数据(FH)，得到高分辨率植被含水量空间分布信息；
将所述高分辨率植被含水量空间分布信息结合MODIS叶面积指数和MODIS植被覆盖度，构建高分辨率植被含水量先验知识库；
S2-2，将粗分辨率的植被含水量的月观测数据、高分辨率植被含水量先验知识库输入高分辨率植被含水量降尺度模型，对粗分辨率的植被含水量的月观测数据进行降尺度，通过所述模型获取公里级植被含水量月观测数据；
S3，采用经验函数拟合植被含水量与碳储量底图，构建基于公里级植被含水量的碳储量动态估算模型，实现植被碳储量动态监测；
所述经验函数为：
当监测区域为热带地区时，采用线性回归模型或arctangent回归模型作为经验函数；采用线性回归模型拟合高分辨率植被含水量和碳储量之间的关系为：
AGC＝a*VOD+b
采用arctangent回归模型来拟合高分辨率植被含水量和碳储量之间的关系为：
式中，VOD为高分辨率植被含水量；
a、b、c、d是回归参数；
当监测区域为温带地区时，采用logistic回归模型作为经验函数；
采用logistic回归模型来拟合高分辨率植被含水量和碳储量之间的关系为：
式中，AGC为植被碳储量；
VOD为高分辨率植被含水量；
a、b、c、d是回归参数；
所述公里级植被含水量通过高分辨率植被含水量降尺度模型获得；
高分辨率植被含水量降尺度模型表示为：
VODHigh＝f(VOD,FH,LAI,TC)
其中，VODHigh为基于随机森林模型估算植被含水量；
f()为随机森林算法；
VOD为哨兵数据反演的植被含水量；
FH为树高；
LAI为MODIS叶面积指数；
TC为MODIS植被覆盖度；
所述高分辨率植被含水量降尺度模型中，需要采用VOD-碳储量反演模型将植被含水量转为碳储量，VOD-碳储量反演模型的构建步骤如下：
首先将水云模型和Ulaby模型耦合得到初步土壤-植被散射模型，再结合地区的土壤水分数据、后向散射数据，构建地区的土壤-植被散射模型；
然后基于土壤-植被散射模型对哨兵1号后向散射数据进行土壤-植被含水量信息分离，获取公里级植被含水量；
土壤-植被散射模型的方程如下：
其中，表示雷达接收到的总体后向散射系数；
A是模型的经验系数，与植被的类型以及雷达参数有关；
V1是对植被冠层的描述量；
VOD表示植被含水量；
θ表示入射角；
SM表示土壤水分；
C代表干燥土壤的后向散射系数；
D代表雷达对土壤水分含量的敏感性；
根据所述土壤-植被散射模型反演出基于哨兵数据的高分辨率植被含水量VOD，得到VOD-碳储量反演模型：
所述高分辨率为小于等于1km的空间分辨率。
2.根据权利要求1所述的一种基于低频微波雷达VOD数据的公里级碳储量评估方法，其特征在于，在训练高分辨率植被含水量降尺度模型时，将模型中的参数Ntree和Mtry的值分别设定为500和2，其中Ntree为随机森林所包含的决策树数目，Mtry为每个决策树包含的节点个数。
3.根据权利要求1所述的一种基于低频微波雷达VOD数据的公里级碳储量评估方法，其特征在于，对植被含水量数据进行预处理还包括：利用土地利用数据剔除包括城镇、裸地和水体的非植被比例大于10％的升、降轨植被含水量数据，所述土地利用数据包括地区的河流、湖泊、湿地中的一种或者任意组合。
4.根据权利要求1所述的一种基于低频微波雷达VOD数据的公里级碳储量评估方法，其特征在于，在采集数据后对MODIS叶面积指数进行预处理：
(1)借助MODIS数据质量控制文件，剔除受到云雨天气干扰的像元；
(2)通过最大值合成法将叶面积指数数据生产为年尺度数据；
(3)对年数据进行拼接、重投影和裁剪预处理，获取地区MODIS叶面积指数数据；
(4)将MODIS叶面积指数进行重采样，使其与植被碳储量底图具有相同的投影和空间分辨率。
5.根据权利要求1所述的一种基于低频微波雷达VOD数据的公里级碳储量评估方法，其特征在于，在采集数据后对地面样地观测数据的筛选及质量控制：筛选出森林覆盖率的标准差σTC<15％的数据。
CN202211664301.5A 2022-12-23 2022-12-23 基于低频微波雷达vod数据的公里级碳储量评估方法 Active CN116452023B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211664301.5A CN116452023B (zh) 2022-12-23 2022-12-23 基于低频微波雷达vod数据的公里级碳储量评估方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211664301.5A CN116452023B (zh) 2022-12-23 2022-12-23 基于低频微波雷达vod数据的公里级碳储量评估方法

Publications (2)

Publication Number Publication Date
CN116452023A CN116452023A (zh) 2023-07-18
CN116452023B true CN116452023B (zh) 2023-09-26

Family

ID=87134439

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211664301.5A Active CN116452023B (zh) 2022-12-23 2022-12-23 基于低频微波雷达vod数据的公里级碳储量评估方法

Country Status (1)

Country Link
CN (1) CN116452023B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117496363A (zh) * 2023-08-28 2024-02-02 广东省国土资源测绘院 融合主被动微波遥感的作物产量估算方法、系统及设备
CN117456351B (zh) * 2023-10-08 2024-05-17 宁波大学 一种星空地协同的滨海湿地盐沼植被碳储量估算方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110646587A (zh) * 2019-09-29 2020-01-03 武汉大学 结合多源遥感数据的高分辨率农业干旱监测方法及装置
CN113205475A (zh) * 2020-01-16 2021-08-03 吉林大学 基于多源卫星遥感数据的森林高度反演方法
CN114781011A (zh) * 2022-04-07 2022-07-22 武汉大学 一种像素级全球森林碳储量高精度计算方法及系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110646587A (zh) * 2019-09-29 2020-01-03 武汉大学 结合多源遥感数据的高分辨率农业干旱监测方法及装置
CN113205475A (zh) * 2020-01-16 2021-08-03 吉林大学 基于多源卫星遥感数据的森林高度反演方法
CN114781011A (zh) * 2022-04-07 2022-07-22 武汉大学 一种像素级全球森林碳储量高精度计算方法及系统

Also Published As

Publication number Publication date
CN116452023A (zh) 2023-07-18

Similar Documents

Publication Publication Date Title
Silva et al. Comparison of small-and large-footprint lidar characterization of tropical forest aboveground structure and biomass: a case study from Central Gabon
Nelson et al. Estimating Siberian timber volume using MODIS and ICESat/GLAS
Santoro et al. Retrieval of growing stock volume in boreal forest using hyper-temporal series of Envisat ASAR ScanSAR backscatter measurements
CN116452023B (zh) 基于低频微波雷达vod数据的公里级碳储量评估方法
Bwangoy et al. Wetland mapping in the Congo Basin using optical and radar remotely sensed data and derived topographical indices
Quegan et al. Multitemporal ERS SAR analysis applied to forest mapping
CN114387516B (zh) 一种针对复杂地形环境下中小田块的单季稻sar识别方法
Ahmad et al. Satellite remote sensing and GIS-based crops forecasting & estimation system in Pakistan
CN114819737B (zh) 公路路域植被的碳储量估算方法、系统及存储介质
Milenković et al. Assessing Amazon rainforest regrowth with GEDI and ICESat-2 data
CN116484712A (zh) 植被区域地表温度重建方法、装置、电子设备及存储介质
Dashpurev et al. A cost-effective method to monitor vegetation changes in steppes ecosystems: A case study on remote sensing of fire and infrastructure effects in eastern Mongolia
CN109212553A (zh) 无人机LiDAR和随机森林提取银杏生物物理特性的方法
Nasirzadehdizaji et al. Application of sentinel-1 multi-temporal data for crop monitoring and mapping
Traore et al. Assessing the inter-relationship between vegetation productivity, rainfall, population and land cover over the Bani River Basin in Mali (West Africa)
Guo et al. Estimating aboveground biomass using Pléiades satellite image in a karst watershed of Guizhou Province, Southwestern China
Hameid et al. The relationship between vegetation and rainfall in central Sudan
Kumar Retrieval of Forest Parameters from Envisat ASAR Data for Biomass Inventory in Dudhwa National Park, UP Indiana
Millard Development of methods to map and monitor peatland ecosystems and hydrologic conditions using Radarsat-2 Synthetic Aperture Radar
Al-Ali et al. Potionential of Spectral Indices for Halophyte Vegetation Cover Detection in Arid and Salt-Affected Landscape
White Improving capacity for large-area monitoring of forest disturbance and recovery
Novresiandi et al. C-band dual-polarization synthetic aperture radar application for peat depth classification: a case study in Siak Regency, Riau Province, Indonesia
Gudmundsdottir Detection of potential arable land with remote sensing and GIS: a case study for Kjósarhreppur
Horrocks An analysis of multispectral unmanned aerial systems for salt marsh foreshore land cover classification and digital elevation model generation
Di Biase A proposal for designed-based mapping of forest resources

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
GR01 Patent grant