CN114021348A - 一种精细土地利用类型的矿区植被碳汇遥感反演方法 - Google Patents
一种精细土地利用类型的矿区植被碳汇遥感反演方法 Download PDFInfo
- Publication number
- CN114021348A CN114021348A CN202111303598.8A CN202111303598A CN114021348A CN 114021348 A CN114021348 A CN 114021348A CN 202111303598 A CN202111303598 A CN 202111303598A CN 114021348 A CN114021348 A CN 114021348A
- Authority
- CN
- China
- Prior art keywords
- vegetation
- carbon
- land
- ndvi
- area
- 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
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 title claims abstract description 195
- 229910052799 carbon Inorganic materials 0.000 title claims abstract description 195
- 238000005065 mining Methods 0.000 title claims abstract description 128
- 238000000034 method Methods 0.000 title claims abstract description 68
- 230000009919 sequestration Effects 0.000 claims abstract description 79
- 238000004364 calculation method Methods 0.000 claims abstract description 51
- 230000003595 spectral effect Effects 0.000 claims abstract description 16
- 230000005855 radiation Effects 0.000 claims description 45
- 239000002689 soil Substances 0.000 claims description 24
- 239000002028 Biomass Substances 0.000 claims description 22
- 238000011160 research Methods 0.000 claims description 15
- 238000009826 distribution Methods 0.000 claims description 12
- 230000007306 turnover Effects 0.000 claims description 12
- 238000010521 absorption reaction Methods 0.000 claims description 11
- 230000000243 photosynthetic effect Effects 0.000 claims description 11
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 9
- 238000005259 measurement Methods 0.000 claims description 8
- 239000011159 matrix material Substances 0.000 claims description 7
- 238000001556 precipitation Methods 0.000 claims description 7
- 238000002310 reflectometry Methods 0.000 claims description 7
- 238000005070 sampling Methods 0.000 claims description 7
- 238000006243 chemical reaction Methods 0.000 claims description 6
- 244000005700 microbiome Species 0.000 claims description 6
- 238000012544 monitoring process Methods 0.000 claims description 6
- 230000029058 respiratory gaseous exchange Effects 0.000 claims description 5
- 230000003287 optical effect Effects 0.000 claims description 4
- 238000012937 correction Methods 0.000 claims description 3
- 238000005520 cutting process Methods 0.000 claims description 3
- 230000029553 photosynthesis Effects 0.000 claims description 3
- 238000010672 photosynthesis Methods 0.000 claims description 3
- 239000003245 coal Substances 0.000 abstract description 6
- 238000000638 solvent extraction Methods 0.000 abstract description 6
- 230000009466 transformation Effects 0.000 abstract description 2
- 235000019646 color tone Nutrition 0.000 description 5
- 238000010586 diagram Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 239000003086 colorant Substances 0.000 description 2
- 239000011152 fibreglass Substances 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 229910052500 inorganic mineral Inorganic materials 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 239000011707 mineral Substances 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 239000002699 waste material Substances 0.000 description 2
- ATJFFYVFTNAWJD-UHFFFAOYSA-N Tin Chemical compound [Sn] ATJFFYVFTNAWJD-UHFFFAOYSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007850 degeneration Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 235000008216 herbs Nutrition 0.000 description 1
- 238000011835 investigation 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
- 239000008239 natural water Substances 0.000 description 1
- 238000006386 neutralization reaction Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 239000002910 solid waste Substances 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 238000013316 zoning Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Artificial Intelligence (AREA)
- Life Sciences & Earth Sciences (AREA)
- Computer Hardware Design (AREA)
- Geometry (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种精细土地利用类型的矿区植被碳汇遥感反演方法，其方法如下：A、获取高分辨率影像数据，结合高光谱分辨率传感器获取的地物光谱曲线特征进行土地利用类型分类；B、采用波段运算方式对矿区影像数据进行植被参数反演，植被参数包括归一化植被指数、比值植被指数、植被覆盖度；C、结合土地利用类型、临界阈值NDVIlim划分为固碳区类型；D、构建植被固碳参数反演模型对植被固碳区进行反演，分别得到植被净初级生产力。本发明通过对矿区土地利用类型精细分类和植被固碳能力分区，结合遥感影像和地面实测数据构建光能利用率模型，提高了矿区植被碳汇计算的精确程度，能够为煤炭清洁开采转型和土地复垦决策提供科学的数据支撑。
Description
技术领域
本发明涉及采矿领域、生态学领域、遥感及地理信息领域，尤其涉及一种精细土地利用类型的矿区植被碳汇遥感反演方法。
背景技术
在全球气候变化的机遇与挑战下，中国力争实现“碳达峰”和“碳中和”的目标，煤炭行业作为碳减排的重要领域，建立碳监测评估方法体系势在必行。植被碳汇作为矿区生态系统碳汇计算的重要组成部分，当前的计算方法主要为地面实地测量法和遥感影像反演法。通过在矿区设立观测站点并进行实地考察，统计植被的叶面积指数、总生物量和净初级生产力等生态数据，可以得到矿区植被碳汇(参见文献：张成业,李军,雷少刚,杨金中,杨楠.矿区生态环境定量遥感监测研究进展与展望[J/OL].金属矿)。但是矿区的实地数据采集工作效率低下，在矿区复杂地形条件下不易作业，难以满足矿区植被碳汇长期大范围监测的需求，一般主要通过气温、总降水量、地表净辐射等数据，接着实现区域长时序、高频次、大范围的植被碳汇计算。然而矿区生产造成了区域特殊的土地利用类型，不同土地利用类型的最大光能利用率和植被有效光合辐射吸收比例调节系数等参数在遥感植被碳汇反演中是不同的，而当前研究往往没有对其精细化处理。综上所述，现有技术方法不能有效且精确的计算出矿区生态系统的植被碳汇，存在以下3个问题：(1)矿区植被碳汇计算过程缺乏对矿区土地利用类型的精细划分，不能评估不同分区植被碳汇的贡献；(2)没有顾及遥感植被碳汇计算模型参数空间上的差异性，缺乏矿区遥感植被碳汇关键参数的反演模型；(3)在矿区采用统一的遥感植被碳汇计算方法，不能反映矿区不同土地利用类型的精确植被碳汇量。
发明内容
本发明的目的在于克服背景技术中所指出的技术问题，提供一种精细土地利用类型的矿区植被碳汇遥感反演方法，通过对矿区土地利用类型精细分类和植被固碳能力分区，结合遥感影像和地面实测数据构建改进的光能利用率模型，提高了矿区植被碳汇计算的精确程度，实现对矿区减排和生态修复状况的精准监测，能够为煤炭清洁开采方式的转型和土地复垦方案的决策提供科学的数据支撑。
本发明的目的通过下述技术方案实现：
一种精细土地利用类型的矿区植被碳汇遥感反演方法，其方法如下：
A、获取高分辨率影像数据，对高分辨率影像数据按照矿区范围进行裁剪得到矿区影像数据，结合高光谱分辨率传感器获取的地物光谱曲线特征进行土地利用类型分类，土地利用类型的大类包括自然地物大类、农业用地大类、矿区用地大类、复垦地大类、城镇用地大类；自然地物大类包括林地、草地、裸地、水体小类，农业用地大类包括耕地小类，矿区用地大类包括露天采场、中转场地、矿区建筑与设施、排土场小类，复垦地大类包括复垦林地、复垦草地、退化土地小类，城镇用地大类包括城镇用地和道路小类；
B、生态系统植被参数反演：采用波段运算方式对矿区影像数据进行植被参数反演，植被参数包括归一化植被指数、比值植被指数、植被覆盖度；
归一化植被指数反演计算公式如下：
其中ρnir为近红外波段地表反射率，ρred为红光波段地表反射率；
比值植被指数反演计算公式如下：
其中ρnir为近红外波段地表反射率，ρred为红光波段地表反射率；
植被覆盖度反演计算公式如下：
其中，NDVI表示遥感影像中该像元的归一化植被指数，即NDVI值，NDVImin为研究区内完全为裸土的像元NDVI值，NDVImax为研究区纯植被像元的NDVI值；
C、结合土地利用类型，将包括草地小类、林地小类、耕地小类、复垦草地小类、复垦林地小类在内的小类划分为植被固碳区；
根据矿区影像数据中露天采场小类的归一化植被指数按照像元统计量95％区间位置确定存在植被固碳能力的临界阈值NDVIlim；将矿区影像数据除开植被固碳区各个小类所对应的NDVI值分别与临界阈值NDVIlim进行逐像元比较并统计大于临界阈值NDVIlim的像元数量所占比例，然后划分为弱植被固碳区、无植被固碳区，
弱植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例>M1，则该土地利用类型小类划分为弱植被固碳区；
无植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例≤M1，则该土地利用类型小类划分为无植被固碳区；
D、构建植被固碳参数反演模型对植被固碳区进行反演，植被净初级生产力NPP计算公式如下：
NPP＝k1·k2·α·F1·εm+k1·k2·F2·εm-k1·k2·α·F2·εm，其中k1为0.5倍的月太阳总辐射，即植被光合有效辐射；k2为光能利用率受到影响的总胁迫系数，F1为通过NDVI计算得到的FPAR,F2为通过RVI计算得到的FPAR，APAR为植被吸收的光合有效辐射；εm为理想条件下植被的最大光能利用率，α为调整系数；
构建CASA模型并实现弱植被固碳区的植被净初级生产力NPP计算；
对于无植被固碳区不计算或剔除。
本发明矿区植被碳汇遥感反演方法步骤C还可采用如下方法(即步骤C提供第二种优选方法)：
根据归一化植被指数对矿区影像数据结合土地利用类型中的小类通过直方图统计得到NDVI值的区间范围和分布情况，并分割矢量边界区域；根据矿区影像数据中露天采场小类的归一化植被指数按照像元统计量95％区间位置确定存在植被固碳能力的临界阈值NDVIlim；将矿区影像数据中矢量边界区域所对应的NDVI值分别与临界阈值NDVIlim进行逐像元比较并统计大于临界阈值NDVIlim的像元数量所占比例，然后划分为植被固碳区、弱植被固碳区、无植被固碳区；
植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例大于M2，则该矢量边界区域或/和土地利用类型小类为植被固碳区；
弱植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例为大于等于M3且小于等于M2，则该矢量边界区域或/和土地利用类型小类为弱植被固碳区；
无植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例为小于M3，则该矢量边界区域或/和土地利用类型小类为无植被固碳区。
本发明矿区植被碳汇遥感反演方法在步骤D中第一种优选的方法如下：在步骤D中，植被固碳参数反演模型方法如下：
D1、按照如下公式得到植被固碳区的植被净初级生产力真实值NPPtrue：
NPP(true，i)＝(AGBi+BGBi)*turnover；其中，AGB表示地上生物量，BGB表示地下生物量，turnover表示单位面积生物量与实测NPP值的转换系数；i表示植被固碳区不同土地利用类型或矢量边界区域；
D2、结合植被净初级生产力真实值NPPtrue及实测气象数据构建植被固碳参数反演模型，利用最小二乘方法进行迭代反演得到参数εm、α；，然后采用改进的光能利用率模型进行植被净初级生产力NPP的计算，得到植被固碳区的植被净初级生产力NPP。
本发明矿区植被碳汇遥感反演方法在步骤D中第二种优选的方法如下：在步骤D中，植被固碳参数反演模型方法如下：
D3、对矿区影像数据的植被固碳区进行植被生物量采样，按照如下公式得到植被固碳区的植被净初级生产力真实值NPPtrue：
NPP(true，i)＝(AGBi+BGBi)*turnover；其中，AGB表示地上生物量，BGB表示地下生物量，turnover表示单位面积生物量与实测NPP值的转换系数；i表示植被固碳区不同土地利用类型；
D4、通过移动气象站对步骤D3所对应植被固碳区的气象要素进行监测，气象要素包括降水量、气温、太阳辐射、土壤温度、土壤湿度；
D5、基于光能利用率模型按照如下公式进行计算：
植被净初级生产力NPP可以由植被吸收的光合有效辐射和实际光能利用率计算得到，公式如下：
NPP＝APAR×ε；式中，APAR为植被吸收的光合有效辐射，由植被光合有效辐射PAR和植被光合有效辐射的吸收比例FPAR求取，按照陆地植被净初级生产力全球估算模型可用以下公式表示：
APAR＝PAR×FPAR；式中，PAR表示植被接收的光合有效辐射，单位为[MJ/(m2·月)]；FPAR为植被光合有效辐射吸收比例，由归一化植被指数NDVI和比值植被指数RVI经过调整系数α计算得到；PAR和FPAR的计算方法如下：
PAR＝0.5×SOL
FPAR＝αFPARNDVI+(1-α)FPARRVI
FPARNDVI＝FVC×(FPARmax-FPARmin)+FPARmin
式中，SOL为太阳总辐射[MJ/(m2·月)]；NDVImax表示完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值；NDVImin表示完全是裸土或纯土壤像元的NDVI值；NDVI表示影像中各个像元对应的NDVI值；RVImin为比值植被指数的最小值，RVImax为比值植被指数的最大值，RVI为影像中各个像元对应的RVI值；FPARmax和FPARmin取值分别为0.95和0.001；FPARNDVI为通过NDVI计算得到的FPAR,FPARRVI为通过RVI计算得到的FPAR；
ε为植被的光能利用率，其计算公式如下：ε＝Ts1×Ts2×Ws×εm，其中Ts1表示低温对光能利用率的胁迫系数，Ts2表示高温对光能利用率的胁迫系数，Ws为水分对光能利用率的胁迫系数，εm为植被的最大光能利用率；
D6、构建不同土地利用类型植被固碳参数εm和α的植被固碳参数反演模型，植被净初级生产力NPP可通过以下公式进行表示：
NPP＝0.5×SOL×[α·FPARNDVI+(1-α)·FPARRVI]×Ts1×Ts2×Ws×εm
将上式中可通过计算得到系数进行合并，公式可简记为：
NPP＝k1×[α·F1+(1-α)·F2]×k2×εm
即，NPP＝k1·k2·α·F1·εm+k1·k2·F2·εm-k1·k2·α·F2·εm
其中k1为0.5倍的月太阳总辐射，即为植被光合有效辐射PAR，k2为光能利用率受到影响的总胁迫系数，F1为通过NDVI计算得到的FPAR,F2为通过RVI计算得到的FPAR；
对上式进行泰勒公式一阶展开，
计算得到
NPP＝NPP(εm0，α0)+(k1·k2·α0·F1+k1·k2·F2-k1·k2·α·F2)·(εm-εm0)
+(k1·k2·F1×εm0-k1·k2·F2×εm0)·(α-α0)
式中，左侧NPP表示地面采样实测NPP，即植被净初级生产力真实值NPPtrue，NPP(εm0，α0)表示基于光能利用率模型的NPP，其最大光能利用率和调节系数设为初值；采用矩阵形式完成植被固碳参数εm和α的植被固碳参数反演模型构建；
D5、在不同土地利用类型进行植被固碳参数反演模型的求解；其中，εm0的初值使用植被最大光能利用率经验值0.389gC/MJ，α0的初值使用经验值0.5，利用光能利用率模型计算相应NPP(εm0，α0)；系数矩阵中各个参数依据气象数据和遥感反演参数获取，然后将不同土地利用类型获取的实测NPP数据代入模型作为NPPtrue进行反演计算；
利用最小二乘方法，依据VTPV＝min的原理，对不同土地利用类型的植被固碳参数模型进行迭代计算；每次迭代计算得到Δε和Δα后，与其相应初值之和作为下一次迭代的εm0，α0来求取模型NPP；每次迭代结束后需对本次迭代结果和上次迭代结果进行比较，若其差值小于0.0001量级即停止迭代计算，将此时Δε，Δα和其对应的εm0，α0求和，得到植被固碳区中草地、林地、耕地、复垦草地、复垦林地小类所对应的最大光能利用率εm和调整系数α。
进一步，本发明矿区植被碳汇遥感反演方法在步骤D之后还包括如下方法：
E、按照步骤D的方法得到植被固碳区、弱植被固碳区的植被净初级生产力NPP；并按照如下公式得到植被固碳区、弱植被固碳区所对应的植被净生态系统生产力NEP：
NEP＝NPP-RH，其中NEP表示植被净生态系统生产力，NPP表示植被净初级生产力，Rh表示土壤微生物的呼吸量。
优选地，本发明矿区植被碳汇遥感反演方法在步骤E中，若植被净生态系统生产力NEP大于0，表明植被固碳量大于土壤碳排放量，则表现为碳汇；若植被净生态系统生产力NEP大于0，表明植被固碳量小于土壤碳排放量，则表现为碳源。
优选地，本发明矿区植被碳汇遥感反演方法在步骤E中，土壤微生物的呼吸量Rh采用如下公式计算：
RH＝0.22×(exp(0.0912T)+ln(0.3145R+1))×30×46.5％，其中，T表示气温(℃)，R表示降水(mm)。
优选地，本发明矿区植被碳汇遥感反演方法在步骤D中，M1选取65％。
优选地，本发明矿区植被碳汇遥感反演方法在进行光谱曲线特征或/和生态系统植被参数反演之前，矿区影像数据需进行包括图像配准、大气校正在内的预处理。
优选地，本发明矿区植被碳汇遥感反演方法步骤A中高分辨率影像数据来源于Landsat系列卫星遥感影像或高分卫星遥感影像或无人机高分辨率遥感影像或其他类似影像(可由上述各类影像构成多源观测数据)；所述高光谱分辨率传感器包括机载高光谱分辨率光谱传感器、高光谱分辨率光谱仪。
本发明较现有技术相比，具有以下优点及有益效果：
(1)本发明通过对矿区土地利用类型精细分类和植被固碳能力分区，结合遥感影像和地面实测数据构建改进的光能利用率模型，旨在提高矿区植被碳汇计算的精确程度，实现对矿区减排和生态修复状况的精准监测，能够为煤炭清洁开采方式的转型和土地复垦方案的决策提供科学的数据支撑。
(2)本发明首先利用多源观测数据对研究矿区进行土地利用类型精细化分类与识别，基于遥感影像进行矿区植被参数的定量反演，采用NDVI阈值法对矿区进行植被固碳能力区的划分，包括植被固碳区、弱植被固碳区和无植被固碳区；然后基于地面采集数据构建植被固碳区不同土地利用类型的植被固碳参数反演模型，利用最小二乘法求解不同土地利用类型的最大光能利用率和植被有效光合辐射吸收比例调整系数；最后基于遥感影像及气象数据，利用CASA模型计算弱植被固碳区的土地利用分类植被碳汇量，利用反演模型参数在植被固碳区构建不同土地利用类型的光能利用率改进模型并精确计算植被固碳区碳汇量，达到矿区植被碳汇遥感精确计算的目的。
(3)本发明针对矿区建立了基于多源观测的土地利用类型精细化分类及识别方法，构建了基于NDVI阈值的矿区植被固碳能力分区方法，明晰了矿区中植被固碳区、弱植被固碳区和无植被固碳区对于矿区植被碳汇的贡献，为矿区植被固碳精确计算提供了最基础的精细化土地利用分类结果及矿区固碳能力分区信息。
(4)本发明通过实测得到的植被净初级生产力NPPtrue，构建了针对矿区不同土地利用类型的遥感植被固碳参数反演模型，获取了更加合理的针对矿区生态系统的不同土地利用类型植被有效光合辐射吸收比例M-FPAR的调整系数和最大光能利用率，为在矿区不同土地利用类型下进行植被碳汇精确计算提供最关键的参数信息。
(5)本发明基于遥感影像和相关气象数据，在植被固碳区构建了不同土地利用类型的光能利用率改进模型，分别对矿区中植被固碳区和弱植被固碳区的碳汇量进行计算，实现了矿区植被碳汇遥感优化计算的目的。
附图说明
图1为实施例中矿区植被碳汇遥感反演方法的原理流程图；
图2为实施例中改进光能利用率模型的NPP计算流程图；
图3为实施例中矩阵形式的示意图；
图4为实施例中2020年矿区精细土地利用类型分类图；
图5为实施例中2020年矿区NDVI空间分布示意图；
图6为实施例中2020年矿区RVI空间分布示意图；
图7为实施例中2020年矿区植被覆盖度空间分布示意图；
图8为实施例中矿区露天采场NDVI值统计与阈值划分结果图；
图9为实施例中基于土地利用类型的矿区植被固碳能力分区图；
图10为实施例中2020年矿区植被碳汇空间分布示意图。
具体实施方式
下面结合实施例对本发明作进一步地详细说明：
实施例
如图1～图10所示，一种精细土地利用类型的矿区植被碳汇遥感反演方法，其方法如下：
A、获取高分辨率影像数据，高分辨率影像数据可以来源于Landsat系列卫星遥感影像或高分卫星遥感影像或无人机高分辨率遥感影像或其他影像(既可来源于单一影像，也可为多源影像数据，为了提高精确度，优选采用多源观测数据，比如Landsat系列陆地资源卫星遥感影像、国产高分卫星遥感影像等高分辨率影像数据；如果由于天气、传感器等影响因素造成某些区域的影像缺失，则需结合无人机高空间分辨率传感器对卫星影像不能精确解译的区域进行确认和再分析)，对高分辨率影像数据按照矿区范围进行裁剪得到矿区影像数据，结合高光谱分辨率传感器(高光谱分辨率传感器包括机载高光谱分辨率光谱传感器、高光谱分辨率光谱仪或其他设备)获取的地物光谱曲线特征进行土地利用类型分类，土地利用类型的大类包括自然地物大类、农业用地大类、矿区用地大类、复垦地大类、城镇用地大类。
自然地物大类包括林地、草地、裸地、水体小类；林地主要指用于林业生产的地区或天然林区；影像中色调大多呈现为深绿色，多呈现为团状或絮状。草地是指生长草本和灌木植物为主并适宜发展畜牧业的土地；影像中色调大多呈现为淡绿色或黄绿色，标准假彩色影像中大多呈现为红色、淡红色或灰褐色，影像结构均一且纹理较细。裸地是指裸露的地表、植被较为稀疏的土地或指在一年内超过特定时间处于无植被生长，但在植被生长旺季仍存在植被的区域；其表面纹理不均一，呈较浅颜色。水体主要指自然水体以及开采废弃地的积水矿坑，标准假彩色中多呈现为蓝黑色或黑色，有明显的区域边界。
农业用地大类包括耕地小类，耕地是指种植农作物获得粮食或产品的土地。影像颜色随着季节的变化有所不同，一般为淡绿色，几何形状明显，大多呈现为块状或圆弧形，边界清晰，结构细腻，有纹理结构。
矿区用地大类包括露天采场、中转场地、矿区建筑与设施、排土场小类，露天采场是指进行剥离、采矿获取岩石矿物的工作场所；呈现为黑色、深暗色调，采坑内有明显道路，边界清晰且形状规则，有明显均一的条纹和层次感，与周围有明显的凹陷形态。中转场地主要指用于堆放煤矿、煤矸石等开采获得的矿物的场所；影像中大多为黑色，分布在露天采坑内或洗煤厂附近，堆积呈现为圆形团状。排土场是指用于堆放矿物开采过程中堆放固体废物的场所；影像中一般色调为灰色或白色，呈现近似梯形或圆形的图斑和凸起的形态。矿区建筑与设施是指矿区范围内用于矿业加工的建筑和停放大型采矿设备的场所；建筑物轮廓清晰，形状规则，相较于城镇建筑分布更为稀疏，大部分矿山建筑物采用蓝色玻璃钢屋顶或者部分橙色玻璃钢屋顶。
复垦地大类包括复垦林地、复垦草地、退化土地小类，复垦林地主要指复垦为林地的排土场；色调呈现为浅绿色或深绿色，颜色随季节变化，形状规则，纹理明显，明显区别于自然林地。复垦草地主要指复垦为林地的排土场；色调多为浅绿色或灰绿色，夏季植被茂盛期颜色较深，有明显的区域边界，在排土场平坦区域呈现片状分布。退化土地是指复垦的排土场二次退化的土地；几何特征明显，形状多为圆形且较规则，色调多为浅黄色或浅褐色，由于土地出现退化，其地物特征与裸地相似。
城镇用地大类包括城镇用地和道路小类，城镇用地的色调由多种颜色混合，边界清晰，形状规则且几何特征明显，影像结构粗糙；道路的颜色相较于周边更加明亮，呈现较为平行且宽度基本不变的双边缘状态，特征较为一致。
本实施例以选取的研究矿区为锡林浩特市胜利一号矿、西二矿、西三矿、东二矿及其周围10公里为例，作为矿区生态系统植被碳汇计算的研究范围。选取2020年的高分二号1m分辨率的卫星影像进行目视解译，并结合矿区实地情况，采用无人机和现场勘测的方式，对目视解译得到的土地利用类型分类结果进行校正。根据该矿区的实际土地利用情况将该区域划分为包括裸地、城镇用地、耕地、林地、草地、露天采场、水体、复垦草地的8种土地利用类型(参见图4)
B、生态系统植被参数反演：采用波段运算方式对矿区影像数据进行植被参数反演(优选地，在进行光谱曲线特征或/和生态系统植被参数反演之前，矿区影像数据需进行包括图像配准、大气校正在内的预处理)，植被参数包括归一化植被指数(NormalizedDifference Vegetation Index,简称NDVI)、比值植被指数(Ratio vegetation index,简称RVI)、植被覆盖度(Fraction of Vegetation Coverage,简称FVC)；本实施例使用Landsat-5、Landsat-7、Landsat-8卫星影像产品和Sentinel-2A影像产品，在Google EarthEngine(GEE)平进行加载。
归一化植被指数反演计算公式如下：
其中ρnir为近红外波段地表反射率，在Landsat-5和Landsat-7中为第4波段，在Landsat-8中为第5波段；ρred为红光波段地表反射率，在Landsat-5和Landsat-7中为第3波段，在Landsat-8中为第4波段；本实施例举例的研究矿区范围进行2020年归一化植被指数反演，其2020年矿区NDVI空间分布参见图5。
比值植被指数反演计算公式如下：
其中ρnir为近红外波段地表反射率，ρred为红光波段地表反射率；本实施例举例的研究矿区范围进行2020年比值植被指数反演，其2020年矿区RVI空间分布参见图6。
植被覆盖度反演计算公式如下：
其中，NDVI表示遥感影像中该像元的归一化植被指数，即NDVI值，NDVImin为研究区内完全为裸土的像元NDVI值，NDVImax为研究区纯植被像元的NDVI值；本实施例举例的研究矿区范围进行2020年植被覆盖度反演，其2020年矿区植被覆盖度空间分布参见图7。
C、结合土地利用类型，将包括草地小类、林地小类、耕地小类、复垦草地小类、复垦林地小类在内的小类划分为植被固碳区；
根据矿区影像数据中露天采场小类的归一化植被指数按照像元统计量95％区间位置(参见图8，在本实施例举例的研究矿区范围，计算得到NDVIlim＝0.1928)确定存在植被固碳能力的临界阈值NDVIlim；将矿区影像数据除开植被固碳区各个小类所对应的NDVI值分别与临界阈值NDVIlim进行逐像元比较并统计大于临界阈值NDVIlim的像元数量所占比例，然后划分为弱植被固碳区、无植被固碳区，
弱植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例>M1(本实施例的M1选取65％)，则该土地利用类型小类划分为弱植被固碳区；
无植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例≤M1，则该土地利用类型小类划分为无植被固碳区。
根据本实施例的一个实施例，本实施例步骤C还可采用如下方法来划分植被固碳区、弱植被固碳区、无植被固碳区(这样，本发明可以采用两种方法进行固碳区划分)：
根据归一化植被指数对矿区影像数据结合土地利用类型中的小类通过直方图统计得到NDVI值的区间范围和分布情况，并分割矢量边界区域；根据矿区影像数据中露天采场小类的归一化植被指数按照像元统计量95％区间位置(参见图8)确定存在植被固碳能力的临界阈值NDVIlim；将矿区影像数据中矢量边界区域所对应的NDVI值分别与临界阈值NDVIlim进行逐像元比较并统计大于临界阈值NDVIlim的像元数量所占比例，然后划分为植被固碳区、弱植被固碳区、无植被固碳区；
植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例大于M2，则该矢量边界区域或/和土地利用类型小类为植被固碳区；
弱植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例为大于等于M3且小于等于M2，M2、M3一般根据当前矿区进行抽样实测选取，则该矢量边界区域或/和土地利用类型小类为弱植被固碳区；
无植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例为小于M3，则该矢量边界区域或/和土地利用类型小类为无植被固碳区。
D、构建植被固碳参数反演模型对植被固碳区进行反演，植被净初级生产力NPP计算公式如下：
NPP＝k1·k2·α·F1·εm+k1·k2·F2·εm-k1·k2·α·F2·εm，其中k1为0.5倍的月太阳总辐射，即植被光合有效辐射；k2为光能利用率受到影响的总胁迫系数，F1为通过NDVI计算得到的FPAR,F2为通过RVI计算得到的FPAR，APAR为植被吸收的光合有效辐射；εm为理想条件下植被的最大光能利用率，α为调整系数；
构建CASA模型并实现弱植被固碳区的植被净初级生产力NPP计算；
对于无植被固碳区不计算或剔除。
根据本实施例的第一个优选实施例，在步骤D中，植被固碳参数反演模型方法如下：
D1、按照如下公式得到植被固碳区的植被净初级生产力真实值NPPtrue：
NPP(true，i)＝(AGBi+BGBi)*turnover；其中，AGB表示地上生物量，BGB表示地下生物量，turnover表示单位面积生物量与实测NPP值的转换系数；i表示植被固碳区不同土地利用类型或矢量边界区域；
D2、结合植被净初级生产力真实值NPPtrue及实测气象数据构建植被固碳参数反演模型，利用最小二乘方法进行迭代反演得到参数εm、α；，然后采用改进的光能利用率模型进行植被净初级生产力NPP的计算，得到植被固碳区的植被净初级生产力NPP。
根据本实施例的第二个优选实施例，在步骤D中，植被固碳参数反演模型方法如下：
D3、对矿区影像数据的植被固碳区进行植被生物量采样，按照如下公式得到植被固碳区的植被净初级生产力真实值NPPtrue：
NPP(true，i)＝(AGBi+BGBi)*turnover；其中，AGB表示地上生物量，BGB表示地下生物量，turnover表示单位面积生物量与实测NPP值的转换系数；i表示植被固碳区不同土地利用类型；
D4、通过移动气象站对步骤D3所对应植被固碳区的气象要素进行监测，气象要素包括降水量、气温、太阳辐射、土壤温度、土壤湿度；
D5、基于光能利用率模型按照如下公式进行计算：
植被净初级生产力NPP可以由植被吸收的光合有效辐射和实际光能利用率计算得到，公式如下：
NPP＝APAR×ε；式中，APAR为植被吸收的光合有效辐射，由植被光合有效辐射PAR和植被光合有效辐射的吸收比例FPAR求取，按照陆地植被净初级生产力全球估算模型可用以下公式表示：
APAR＝PAR×FPAR；式中，PAR表示植被接收的光合有效辐射，单位为[MJ/(m2·月)]；FPAR为植被光合有效辐射吸收比例，由归一化植被指数NDVI和比值植被指数RVI经过调整系数α计算得到；PAR和FPAR的计算方法如下：
PAR＝0.5×SOL
FPAR＝αFPARNDVI+(1-α)FPARRVI
FPARNDVI＝FVC×(FPARmax-FPARmin)+FPARmin
式中，SOL为太阳总辐射[MJ/(m2·月)]；NDVImax表示完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值；NDVImin表示完全是裸土或纯土壤像元的NDVI值；NDVI表示影像中各个像元对应的NDVI值；RVImin为比值植被指数的最小值，RVImax为比值植被指数的最大值，RVI为影像中各个像元对应的RVI值；FPARmax和FPARmin取值分别为0.95和0.001；FPARNDVI为通过NDVI计算得到的FPAR,FPARRVI为通过RVI计算得到的FPAR；
ε为植被的光能利用率，其计算公式如下：ε＝Ts1×Ts2×Ws×εm，其中Ts1表示低温对光能利用率的胁迫系数，Ts2表示高温对光能利用率的胁迫系数，Ws为水分对光能利用率的胁迫系数，εm为植被的最大光能利用率；
D6、构建不同土地利用类型植被固碳参数εm和α的植被固碳参数反演模型，植被净初级生产力NPP可通过以下公式进行表示：
NPP＝0.5×SOL×[α·FPARNDVI+(1-α)·FPARRVI]×Ts1×Ts2×Ws×εm
将上式中可通过计算得到系数进行合并，公式可简记为：
NPP＝k1×[α·F1+(1-α)·F2]×k2×εm
即，NPP＝k1·k2·α·F1·εm+k1·k2·F2·εm-k1·k2·α·F2·εm
其中k1为0.5倍的月太阳总辐射，即为植被光合有效辐射PAR，k2为光能利用率受到影响的总胁迫系数，F1为通过NDVI计算得到的FPAR,F2为通过RVI计算得到的FPAR；
对上式进行泰勒公式一阶展开(展开图参见图3)，
计算得到
NPP＝NPP(εm0，α0)+(k1·k2·α0·F1+k1·k2·F2-k1·k2·α·F2)·(εm-εm0)
+(k1·k2·F1×εm0-k1·k2·F2×εm0)·(α-α0)
式中，左侧NPP表示地面采样实测NPP，即植被净初级生产力真实值NPPtrue，NPP(εm0，α0)表示基于光能利用率模型的NPP，其最大光能利用率和调节系数设为初值；采用矩阵形式(矩阵形式泰勒公式一阶展开参见图3)完成植被固碳参数εm和α的植被固碳参数反演模型构建；
D5、在不同土地利用类型进行植被固碳参数反演模型的求解；其中，εm0的初值使用植被最大光能利用率经验值0.389gC/MJ，α0的初值使用经验值0.5，利用光能利用率模型计算相应NPP(εm0，α0)；系数矩阵中各个参数依据气象数据和遥感反演参数获取，然后将不同土地利用类型获取的实测NPP数据代入模型作为NPPtrue进行反演计算；
利用最小二乘方法，依据VTPV＝min的原理，对不同土地利用类型的植被固碳参数模型进行迭代计算；每次迭代计算得到Δε和Δα后，与其相应初值之和作为下一次迭代的εm0，α0来求取模型NPP；每次迭代结束后需对本次迭代结果和上次迭代结果进行比较，若其差值小于0.0001量级即停止迭代计算，将此时Δε，Δα和其对应的εm0，α0求和，得到植被固碳区中草地、林地、耕地、复垦草地、复垦林地小类所对应的最大光能利用率εm和调整系数α。本实施例举例研究矿区范围调整系数等反演结果如下表1：
表1研究区精细土地利用分类最大光能利用率和植被有效光合辐射吸收比例的调整系数反演结果
本发明对于弱植被固碳区的植被净初级生产力NPP计算还可以采用方法如下：
(1)气象站点数据反演
获取国家气象科学数据中心的气象站点数据，选择某个矿区生态系统范围所在的省区划内的站点数据进行下载，所需的数据种类包括：月平均气温数据集、月平均降水数据集、月平均辐射数据集，站点数据应当覆盖矿区所在的省界范围，并采用反距离权重法进行空间插值，裁剪得到研究区矿区生态系统的气象数据栅格影像。
(2)改进光能利用率模型的NPP计算
基于矿区范围内的卫星遥感影像，按照土地利用类型精细分类标准生成对应的矢量图层，并进行矿区固碳能力分区划分；利用对植被固碳区中不同土地利用类型进行固碳参数反演，最小二乘得到不同土地利用类型的最大光能利用率εm和植被有效光合辐射吸收比例M-FPAR的调整系数Mα。基于遥感影像及相关气象数据，利用上述反演参数改进的光能利用率模型，对植被固碳区各个土地利用类型分别进行NPP精确计算。对于弱植被固碳区，利用遥感影像和气象数据，采用CASA模型进行计算整个弱植被固碳区的NPP计算；对无植被固碳区进行剔除，不参与计算；参照图2所示，即可求解得到矿区生态系统植被大范围研究区的精确净初级生产力。
根据本发明的一个优选实施例：在步骤D之后还包括如下方法：
E、按照步骤D的方法得到植被固碳区、弱植被固碳区的植被净初级生产力NPP；并按照如下公式得到植被固碳区、弱植被固碳区所对应的植被净生态系统生产力NEP：
NEP＝NPP-RH，其中NEP表示植被净生态系统生产力，NPP表示植被净初级生产力，Rh表示土壤微生物的呼吸量。
根据本发明的一个优选实施例：在步骤E中，若植被净生态系统生产力NEP大于0，表明植被固碳量大于土壤碳排放量，则表现为碳汇；若植被净生态系统生产力NEP大于0，表明植被固碳量小于土壤碳排放量，则表现为碳源。本实施例举例研究矿区范围2020年矿区植被碳汇空间分布参见图10.
根据本发明的一个优选实施例：在步骤E中，土壤微生物的呼吸量Rh采用如下公式计算：RH＝0.22×(exp(0.0912T)+ln(0.3145R+1))×30×46.5％，其中，T表示气温(℃)，R表示降水(mm)。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种精细土地利用类型的矿区植被碳汇遥感反演方法，其特征在于：其方法如下：
A、获取高分辨率影像数据，对高分辨率影像数据按照矿区范围进行裁剪得到矿区影像数据，结合高光谱分辨率传感器获取的地物光谱曲线特征进行土地利用类型分类，土地利用类型的大类包括自然地物大类、农业用地大类、矿区用地大类、复垦地大类、城镇用地大类；自然地物大类包括林地、草地、裸地、水体小类，农业用地大类包括耕地小类，矿区用地大类包括露天采场、中转场地、矿区建筑与设施、排土场小类，复垦地大类包括复垦林地、复垦草地、退化土地小类，城镇用地大类包括城镇用地和道路小类；
B、生态系统植被参数反演：采用波段运算方式对矿区影像数据进行植被参数反演，植被参数包括归一化植被指数、比值植被指数、植被覆盖度；
归一化植被指数反演计算公式如下：
其中ρnir为近红外波段地表反射率，ρred为红光波段地表反射率；
比值植被指数反演计算公式如下：
其中ρnir为近红外波段地表反射率，ρred为红光波段地表反射率；
植被覆盖度反演计算公式如下：
其中，NDVI表示遥感影像中该像元的归一化植被指数，即NDVI值，NDVImin为研究区内完全为裸土的像元NDVI值，NDVImax为研究区纯植被像元的NDVI值；
C、结合土地利用类型，将包括草地小类、林地小类、耕地小类、复垦草地小类、复垦林地小类在内的小类划分为植被固碳区；
根据矿区影像数据中露天采场小类的归一化植被指数按照像元统计量95％区间位置确定存在植被固碳能力的临界阈值NDVIlim；将矿区影像数据除开植被固碳区各个小类所对应的NDVI值分别与临界阈值NDVIlim进行逐像元比较并统计大于临界阈值NDVIlim的像元数量所占比例，然后划分为弱植被固碳区、无植被固碳区，
弱植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例>M1，则该土地利用类型小类划分为弱植被固碳区；
无植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例≤M1，则该土地利用类型小类划分为无植被固碳区；
D、构建植被固碳参数反演模型对植被固碳区进行反演，植被净初级生产力NPP计算公式如下：
NPP＝k1·k2·α·F1·εm+k1·k2·F2·εm-k1·k2·α·F2·εm，其中k1为0.5倍的月太阳总辐射，即植被光合有效辐射；k2为光能利用率受到影响的总胁迫系数，F1为通过NDVI计算得到的FPAR，F2为通过RVI计算得到的FPAR，APAR为植被吸收的光合有效辐射；εm为理想条件下植被的最大光能利用率，α为调整系数；
构建CASA模型并实现弱植被固碳区的植被净初级生产力NPP计算；
对于无植被固碳区不计算或剔除。
2.按照权利要求1所述的一种精细土地利用类型的矿区植被碳汇遥感反演方法，其特征在于：步骤C还可采用如下方法：
根据归一化植被指数对矿区影像数据结合土地利用类型中的小类通过直方图统计得到NDVI值的区间范围和分布情况，并分割矢量边界区域；根据矿区影像数据中露天采场小类的归一化植被指数按照像元统计量95％区间位置确定存在植被固碳能力的临界阈值NDVIlim；将矿区影像数据中矢量边界区域所对应的NDVI值分别与临界阈值NDVIlim进行逐像元比较并统计大于临界阈值NDVIlim的像元数量所占比例，然后划分为植被固碳区、弱植被固碳区、无植被固碳区；
植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例大于M2，则该矢量边界区域或/和土地利用类型小类为植被固碳区；
弱植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例为大于等于M3且小于等于M2，则该矢量边界区域或/和土地利用类型小类为弱植被固碳区；
无植被固碳区为：大于临界阈值NDVIlim像元数量占比的比例为小于M3，则该矢量边界区域或/和土地利用类型小类为无植被固碳区。
3.按照权利要求1或2所述的一种精细土地利用类型的矿区植被碳汇遥感反演方法，其特征在于：在步骤D中，植被固碳参数反演模型方法如下：
D1、按照如下公式得到植被固碳区的植被净初级生产力真实值NPPtrue：NPP(true，i)＝(AGBi+BGBi)*turnover；其中，AGB表示地上生物量，BGB表示地下生物量，turnover表示单位面积生物量与实测NPP值的转换系数；i表示植被固碳区不同土地利用类型或矢量边界区域；
D2、结合植被净初级生产力真实值NPPtrue及实测气象数据构建植被固碳参数反演模型，利用最小二乘方法进行迭代反演得到参数εm、α；，然后采用改进的光能利用率模型进行植被净初级生产力NPP的计算，得到植被固碳区的植被净初级生产力NPP。
4.按照权利要求1或2所述的一种精细土地利用类型的矿区植被碳汇遥感反演方法，其特征在于：在步骤D中，植被固碳参数反演模型方法如下：
D3、对矿区影像数据的植被固碳区进行植被生物量采样，按照如下公式得到植被固碳区的植被净初级生产力真实值NPPtrue：NPP(true，i)＝(AGBi+BGBi)*turnover；其中，AGB表示地上生物量，BGB表示地下生物量，turnover表示单位面积生物量与实测NPP值的转换系数；i表示植被固碳区不同土地利用类型；
D4、通过移动气象站对步骤D3所对应植被固碳区的气象要素进行监测，气象要素包括降水量、气温、太阳辐射、土壤温度、土壤湿度；
D5、基于光能利用率模型按照如下公式进行计算：
植被净初级生产力NPP可以由植被吸收的光合有效辐射和实际光能利用率计算得到，公式如下：
NPP＝APAR×ε；式中，APAR为植被吸收的光合有效辐射，由植被光合有效辐射PAR和植被光合有效辐射的吸收比例FPAR求取，按照陆地植被净初级生产力全球估算模型可用以下公式表示：
APAR＝PAR×FPAR；式中，PAR表示植被接收的光合有效辐射，单位为[MJ/(m2·月)]；FPAR为植被光合有效辐射吸收比例，由归一化植被指数NDVI和比值植被指数RVI经过调整系数α计算得到；PAR和FPAR的计算方法如下：
PAR＝0.5×SOL
FPAR＝αFPARNDVI+(1-α)FPARRVI
FPARNDVI＝FVC×(FPARmax-FPARmin)+FPARmin
式中，SOL为太阳总辐射[MJ/(m2·月)]；NDVImax表示完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值；NDVImin表示完全是裸土或纯土壤像元的NDVI值；NDVI表示影像中各个像元对应的NDVI值；RVImin为比值植被指数的最小值，RVImax为比值植被指数的最大值，RVI为影像中各个像元对应的RVI值；FPARmax和FPARmin取值分别为0.95和0.001；FPARNDVI为通过NDVI计算得到的FPAR，FPARRVI为通过RVI计算得到的FPAR；
ε为植被的光能利用率，其计算公式如下：ε＝Ts1×Ts2×Ws×εm，其中Ts1表示低温对光能利用率的胁迫系数，Ts2表示高温对光能利用率的胁迫系数，Ws为水分对光能利用率的胁迫系数，εm为植被的最大光能利用率；
D6、构建不同土地利用类型植被固碳参数εm和α的植被固碳参数反演模型，植被净初级生产力NPP可通过以下公式进行表示：
NPP＝0.5×SOL×[α·FPARNDVI+(1-α)·FPARRVI]×Td1×Ts2×Ws×εm
将上式中可通过计算得到系数进行合并，公式可简记为：
NPP＝k1×[α·F1+(1-α)·F2]×k2×εm
即，NPP＝k1·k2·α·F1·εm+k1·k2·F2·εm-k1·k2·α·F2·εm
其中k1为0.5倍的月太阳总辐射，即为植被光合有效辐射PAR，k2为光能利用率受到影响的总胁迫系数，F1为通过NDVI计算得到的FPAR，F2为通过RVI计算得到的FPAR；
对上式进行泰勒公式一阶展开，
计算得到
NPP＝NPP(εm0，α0)+(k1·k2·α0·F1+k1·k2·F2-k1·k2·α·F2)·(εm-εm0)+(k1·k2·F1×εm0-k1·k2·F2×εm0)·(α-α0)
式中，左侧NPP表示地面采样实测NPP，即植被净初级生产力真实值NPPtrue，NPP(εm0，α0)表示基于光能利用率模型的NPP，其最大光能利用率和调节系数设为初值；采用矩阵形式完成植被固碳参数εm和α的植被固碳参数反演模型构建；
D5、在不同土地利用类型进行植被固碳参数反演模型的求解；其中，εm0的初值使用植被最大光能利用率经验值0.389gC/MJ，α0的初值使用经验值0.5，利用光能利用率模型计算相应NPP(εm0，α0)；系数矩阵中各个参数依据气象数据和遥感反演参数获取，然后将不同土地利用类型获取的实测NPP数据代入模型作为NPPtrue进行反演计算；
利用最小二乘方法，依据VTPV＝min的原理，对不同土地利用类型的植被固碳参数模型进行迭代计算；每次迭代计算得到Δε和Δα后，与其相应初值之和作为下一次迭代的εm0，α0来求取模型NPP；每次迭代结束后需对本次迭代结果和上次迭代结果进行比较，若其差值小于0.0001量级即停止迭代计算，将此时Δε，Δα和其对应的εm0，α0求和，得到植被固碳区中草地、林地、耕地、复垦草地、复垦林地小类所对应的最大光能利用率εm和调整系数α。
5.按照权利要求1或2所述的一种精细土地利用类型的矿区植被碳汇遥感反演方法，其特征在于：在步骤D之后还包括如下方法：
E、按照步骤D的方法得到植被固碳区、弱植被固碳区的植被净初级生产力NPP；并按照如下公式得到植被固碳区、弱植被固碳区所对应的植被净生态系统生产力NEP：
NEP＝NPP-RH，其中NEP表示植被净生态系统生产力，NPP表示植被净初级生产力，Rh表示土壤微生物的呼吸量。
6.按照权利要求5所述的一种精细土地利用类型的矿区植被碳汇遥感反演方法，其特征在于：在步骤E中，若植被净生态系统生产力NEP大于0，表明植被固碳量大于土壤碳排放量，则表现为碳汇；若植被净生态系统生产力NEP大于0，表明植被固碳量小于土壤碳排放量，则表现为碳源。
7.按照权利要求5所述的一种精细土地利用类型的矿区植被碳汇遥感反演方法，其特征在于：在步骤E中，土壤微生物的呼吸量Rh采用如下公式计算：RH＝0.22×(exp(0.0912T)+ln(0.3145R+1))×30×46.5％，其中，T表示气温(℃)，R表示降水(mm)。
8.按照权利要求1所述的一种精细土地利用类型的矿区植被碳汇遥感反演方法，其特征在于：在步骤D中，M1选取65％。
9.按照权利要求1或2所述的一种精细土地利用类型的矿区植被碳汇遥感反演方法，其特征在于：在进行光谱曲线特征或/和生态系统植被参数反演之前，矿区影像数据需进行包括图像配准、大气校正在内的预处理。
10.按照权利要求1或2所述的一种精细土地利用类型的矿区植被碳汇遥感反演方法，其特征在于：高分辨率影像数据来源于Landsat系列卫星遥感影像或高分卫星遥感影像或无人机高分辨率遥感影像；所述高光谱分辨率传感器包括机载高光谱分辨率光谱传感器、高光谱分辨率光谱仪。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111303598.8A CN114021348B (zh) | 2021-11-05 | 2021-11-05 | 一种精细土地利用类型的矿区植被碳汇遥感反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111303598.8A CN114021348B (zh) | 2021-11-05 | 2021-11-05 | 一种精细土地利用类型的矿区植被碳汇遥感反演方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114021348A true CN114021348A (zh) | 2022-02-08 |
| CN114021348B CN114021348B (zh) | 2022-05-20 |
Family
ID=80061535
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111303598.8A Expired - Fee Related CN114021348B (zh) | 2021-11-05 | 2021-11-05 | 一种精细土地利用类型的矿区植被碳汇遥感反演方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114021348B (zh) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115082274A (zh) * | 2022-06-09 | 2022-09-20 | 贵州师范学院 | 基于卫星遥感的地表植被碳汇估算与交易方法及系统 |
| CN115115224A (zh) * | 2022-06-28 | 2022-09-27 | 中国矿业大学(北京) | 一种矿区植被碳汇分级及变化趋势评价分析方法 |
| CN115131616A (zh) * | 2022-07-20 | 2022-09-30 | 中科三清科技有限公司 | 土地利用类型的分类方法、装置、设备及存储介质 |
| CN115272853A (zh) * | 2022-07-27 | 2022-11-01 | 清华大学 | 基于人工智能技术与大数据的工业废弃地识别方法及产品 |
| CN115310757A (zh) * | 2022-07-01 | 2022-11-08 | 南京信息工程大学 | 一种高精度动态可变分辨率的生态碳汇评估方法及系统 |
| CN115526098A (zh) * | 2022-09-14 | 2022-12-27 | 国家能源投资集团有限责任公司 | 一种矿区地表植被叶面积指数遥感计算方法及电子设备 |
| CN115565063A (zh) * | 2022-03-24 | 2023-01-03 | 中国矿业大学(北京) | 基于气候潜在补偿的矿区植被碳汇贡献量计算与分析方法 |
| CN115602259A (zh) * | 2022-09-27 | 2023-01-13 | 常州大学（Cn） | 生物质液体燃料土地利用变化对环境影响研究方法及系统 |
| CN116341986A (zh) * | 2023-04-07 | 2023-06-27 | 北京一点五度科技有限公司 | 一种多源多时相传感器数据融合的碳汇项目审计与评级方法 |
| CN116429938A (zh) * | 2023-04-12 | 2023-07-14 | 中国科学院地理科学与资源研究所 | 一种基于遥感影像的土地n2o排放量的测算方法及系统 |
| CN117689923A (zh) * | 2023-09-20 | 2024-03-12 | 中国测绘科学研究院 | 一种高精度的土地利用与覆被碳收支变化分类方法 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106446564A (zh) * | 2016-09-29 | 2017-02-22 | 华南农业大学 | 一种植被净初级生产力遥感估算方法 |
| CN109916839A (zh) * | 2019-04-03 | 2019-06-21 | 中国矿业大学 | 一种基于高分遥感影像的露头残煤燃烧中心精准探测方法 |
| CN112257225A (zh) * | 2020-09-16 | 2021-01-22 | 中国科学院地理科学与资源研究所 | 一种适用于高寒草地生态系统的npp计算方法 |
| CN112580982A (zh) * | 2020-12-21 | 2021-03-30 | 北京市测绘设计研究院 | 一种基于多时相遥感和casa模型的生态保护红线实施评估 |
| CN113139901A (zh) * | 2021-04-15 | 2021-07-20 | 青岛地质工程勘察院(青岛地质勘查开发局) | 一种流域尺度植被净初级生产力遥感精细反演方法 |
| WO2021184550A1 (zh) * | 2020-03-19 | 2021-09-23 | 中国科学院南京土壤研究所 | 一种基于土壤传递函数的目标土壤性质含量预测方法 |
| CN113591759A (zh) * | 2021-08-09 | 2021-11-02 | 中国矿业大学(北京) | 一种矿区长时序地表土壤含水量遥感数据生产方法及系统 |
-
2021
- 2021-11-05 CN CN202111303598.8A patent/CN114021348B/zh not_active Expired - Fee Related
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106446564A (zh) * | 2016-09-29 | 2017-02-22 | 华南农业大学 | 一种植被净初级生产力遥感估算方法 |
| CN109916839A (zh) * | 2019-04-03 | 2019-06-21 | 中国矿业大学 | 一种基于高分遥感影像的露头残煤燃烧中心精准探测方法 |
| WO2021184550A1 (zh) * | 2020-03-19 | 2021-09-23 | 中国科学院南京土壤研究所 | 一种基于土壤传递函数的目标土壤性质含量预测方法 |
| CN112257225A (zh) * | 2020-09-16 | 2021-01-22 | 中国科学院地理科学与资源研究所 | 一种适用于高寒草地生态系统的npp计算方法 |
| CN112580982A (zh) * | 2020-12-21 | 2021-03-30 | 北京市测绘设计研究院 | 一种基于多时相遥感和casa模型的生态保护红线实施评估 |
| CN113139901A (zh) * | 2021-04-15 | 2021-07-20 | 青岛地质工程勘察院(青岛地质勘查开发局) | 一种流域尺度植被净初级生产力遥感精细反演方法 |
| CN113591759A (zh) * | 2021-08-09 | 2021-11-02 | 中国矿业大学(北京) | 一种矿区长时序地表土壤含水量遥感数据生产方法及系统 |
Non-Patent Citations (4)
| Title |
|---|
| SANTOS, ETF等: "TRIGONAL MESHES IN DIFFRACTION TOMOGRAPHY WITH OPTIMUM REGULARIZATION: AN APPLICATION FOR CARBON SEQUESTRATION MONITORING", 《JOURNAL OF SEISMIC EXPLORATION》 * |
| 包妮沙等: "国产高分遥感数据估算草地NPP", 《东北大学学报(自然科学版)》 * |
| 尹小君等: "基于Landsat和MODIS数据融合的农牧区NPP模拟", 《农业机械学报》 * |
| 王雪等: "蔚县矿区植被净初级生产力时空变化特征及影响因素", 《生态与农村环境学报》 * |
Cited By (16)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115565063B (zh) * | 2022-03-24 | 2023-04-11 | 中国矿业大学(北京) | 基于气候潜在补偿的矿区植被碳汇贡献量计算与分析方法 |
| CN115565063A (zh) * | 2022-03-24 | 2023-01-03 | 中国矿业大学(北京) | 基于气候潜在补偿的矿区植被碳汇贡献量计算与分析方法 |
| CN115082274A (zh) * | 2022-06-09 | 2022-09-20 | 贵州师范学院 | 基于卫星遥感的地表植被碳汇估算与交易方法及系统 |
| CN115115224A (zh) * | 2022-06-28 | 2022-09-27 | 中国矿业大学(北京) | 一种矿区植被碳汇分级及变化趋势评价分析方法 |
| CN115310757A (zh) * | 2022-07-01 | 2022-11-08 | 南京信息工程大学 | 一种高精度动态可变分辨率的生态碳汇评估方法及系统 |
| CN115310757B (zh) * | 2022-07-01 | 2023-06-06 | 南京信息工程大学 | 一种高精度动态可变分辨率的生态碳汇评估方法及系统 |
| CN115131616A (zh) * | 2022-07-20 | 2022-09-30 | 中科三清科技有限公司 | 土地利用类型的分类方法、装置、设备及存储介质 |
| CN115272853A (zh) * | 2022-07-27 | 2022-11-01 | 清华大学 | 基于人工智能技术与大数据的工业废弃地识别方法及产品 |
| CN115272853B (zh) * | 2022-07-27 | 2023-06-27 | 清华大学 | 基于人工智能技术与大数据的工业废弃地识别方法及产品 |
| CN115526098A (zh) * | 2022-09-14 | 2022-12-27 | 国家能源投资集团有限责任公司 | 一种矿区地表植被叶面积指数遥感计算方法及电子设备 |
| CN115602259A (zh) * | 2022-09-27 | 2023-01-13 | 常州大学（Cn） | 生物质液体燃料土地利用变化对环境影响研究方法及系统 |
| CN116341986A (zh) * | 2023-04-07 | 2023-06-27 | 北京一点五度科技有限公司 | 一种多源多时相传感器数据融合的碳汇项目审计与评级方法 |
| CN116341986B (zh) * | 2023-04-07 | 2024-01-12 | 北京一点五度科技有限公司 | 一种多源多时相传感器数据融合的碳汇项目审计与评级方法 |
| CN116429938A (zh) * | 2023-04-12 | 2023-07-14 | 中国科学院地理科学与资源研究所 | 一种基于遥感影像的土地n2o排放量的测算方法及系统 |
| CN116429938B (zh) * | 2023-04-12 | 2024-01-02 | 中国科学院地理科学与资源研究所 | 一种基于遥感影像的土地n2o排放量的测算方法及系统 |
| CN117689923A (zh) * | 2023-09-20 | 2024-03-12 | 中国测绘科学研究院 | 一种高精度的土地利用与覆被碳收支变化分类方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114021348B (zh) | 2022-05-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114021348B (zh) | 一种精细土地利用类型的矿区植被碳汇遥感反演方法 | |
| Chen et al. | Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements | |
| Milesi et al. | Assessing the impact of urban land development on net primary productivity in the southeastern United States | |
| Mathew et al. | Prediction of surface temperatures for the assessment of urban heat island effect over Ahmedabad city using linear time series model | |
| Deyong et al. | How does the conversion of land cover to urban use affect net primary productivity? A case study in Shenzhen city, China | |
| Hively et al. | Using satellite remote sensing to estimate winter cover crop nutrient uptake efficiency | |
| CN101963664B (zh) | 基于水陆地物分类信息的微波遥感混合像元分解方法 | |
| CN107067334B (zh) | 基于多源遥感数据的农田火灾和秸秆焚烧监测方法及系统 | |
| Liu et al. | Application of a new leaf area index algorithm to China's landmass using MODIS data for carbon cycle research | |
| CN103942420A (zh) | 一种建筑物尺度的太阳能快速估算方法 | |
| CN104915674A (zh) | Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法 | |
| CN116757357B (zh) | 一种耦合多源遥感信息的土地生态状况评估方法 | |
| CN112465332A (zh) | 一种城市人工湿地公园生态地质环境稳定性的评价方法 | |
| CN115565063B (zh) | 基于气候潜在补偿的矿区植被碳汇贡献量计算与分析方法 | |
| CN115204691B (zh) | 基于机器学习和遥感技术的城市人为热排放量估算方法 | |
| CN111160799A (zh) | 一种自然资源数据库构建方法 | |
| CN112365158A (zh) | 一种基于遥感数据的矿山复绿监测与评价方法 | |
| CN116504014A (zh) | 一种基于多维感知的森林火灾监测方法及系统 | |
| Tian et al. | Impacts of anthropogenic and biophysical factors on ecological land using logistic regression and random forest: A case study in Mentougou District, Beijing, China | |
| CN113919227B (zh) | 一种矿区生态时间累积效应点与空间累积范围识别方法 | |
| CN116228041A (zh) | 一种废弃矿山生态修复后碳指标计算方法 | |
| CN114429455A (zh) | 一种基于机载激光点云辅助影像植被三维变化检测方法 | |
| Zhu et al. | UAV flight height impacts on wheat biomass estimation via machine and deep learning | |
| Som-ard | Rice security assessment using geo-spatial analysis | |
| CN116046770A (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 | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20220520 |