CN114970184B - 同步反演高分辨率人为co2排放与自然co2通量的同化方法及系统 - Google Patents
同步反演高分辨率人为co2排放与自然co2通量的同化方法及系统 Download PDFInfo
- Publication number
- CN114970184B CN114970184B CN202210642075.4A CN202210642075A CN114970184B CN 114970184 B CN114970184 B CN 114970184B CN 202210642075 A CN202210642075 A CN 202210642075A CN 114970184 B CN114970184 B CN 114970184B
- Authority
- CN
- China
- Prior art keywords
- flux
- assimilation
- data
- natural
- artificial
- 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
Links
- 230000004907 flux Effects 0.000 title claims abstract description 119
- 238000000034 method Methods 0.000 title claims abstract description 94
- 230000001360 synchronised effect Effects 0.000 title claims abstract description 19
- 238000005457 optimization Methods 0.000 claims abstract description 40
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 16
- 239000002028 Biomass Substances 0.000 claims abstract description 4
- 238000002485 combustion reaction Methods 0.000 claims abstract description 4
- 238000002360 preparation method Methods 0.000 claims abstract description 4
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 claims description 76
- 229910052799 carbon Inorganic materials 0.000 claims description 76
- 238000006243 chemical reaction Methods 0.000 claims description 21
- 238000004458 analytical method Methods 0.000 claims description 19
- 239000013598 vector Substances 0.000 claims description 18
- 238000004088 simulation Methods 0.000 claims description 17
- 238000004364 calculation method Methods 0.000 claims description 16
- 238000005516 engineering process Methods 0.000 claims description 14
- 230000004807 localization Effects 0.000 claims description 13
- 239000011159 matrix material Substances 0.000 claims description 13
- 230000008569 process Effects 0.000 claims description 13
- 230000000241 respiratory effect Effects 0.000 claims description 12
- 230000029058 respiratory gaseous exchange Effects 0.000 claims description 10
- 238000003908 quality control method Methods 0.000 claims description 8
- 238000012795 verification Methods 0.000 claims description 8
- PMGQWSIVQFOFOQ-YKVZVUFRSA-N clemastine fumarate Chemical compound OC(=O)\C=C\C(O)=O.CN1CCC[C@@H]1CCO[C@@](C)(C=1C=CC(Cl)=CC=1)C1=CC=CC=C1 PMGQWSIVQFOFOQ-YKVZVUFRSA-N 0.000 claims description 6
- 238000000354 decomposition reaction Methods 0.000 claims description 6
- 230000004927 fusion Effects 0.000 claims description 6
- 230000029553 photosynthesis Effects 0.000 claims description 6
- 238000010672 photosynthesis Methods 0.000 claims description 6
- 238000000342 Monte Carlo simulation Methods 0.000 claims description 5
- 238000012216 screening Methods 0.000 claims description 4
- 230000005856 abnormality Effects 0.000 claims description 3
- 238000009499 grossing Methods 0.000 claims description 3
- 238000001556 precipitation Methods 0.000 claims description 3
- 238000000638 solvent extraction Methods 0.000 claims description 3
- 238000003860 storage Methods 0.000 claims description 3
- 239000002352 surface water Substances 0.000 claims description 3
- 230000017105 transposition Effects 0.000 claims description 2
- CURLTUGMZLYLDI-UHFFFAOYSA-N Carbon dioxide Chemical compound O=C=O CURLTUGMZLYLDI-UHFFFAOYSA-N 0.000 description 270
- 229910002092 carbon dioxide Inorganic materials 0.000 description 135
- 239000001569 carbon dioxide Substances 0.000 description 126
- 230000006870 function Effects 0.000 description 29
- 230000005540 biological transmission Effects 0.000 description 7
- 230000008859 change Effects 0.000 description 7
- GVVPGTZRZFNKDS-JXMROGBWSA-N geranyl diphosphate Chemical compound CC(C)=CCC\C(C)=C\CO[P@](O)(=O)OP(O)(O)=O GVVPGTZRZFNKDS-JXMROGBWSA-N 0.000 description 7
- 238000006386 neutralization reaction Methods 0.000 description 7
- 230000000243 photosynthetic effect Effects 0.000 description 7
- 238000010168 coupling process Methods 0.000 description 6
- 230000008878 coupling Effects 0.000 description 5
- 238000005859 coupling reaction Methods 0.000 description 5
- 238000003745 diagnosis Methods 0.000 description 5
- 230000005855 radiation Effects 0.000 description 5
- 238000013459 approach Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000011835 investigation Methods 0.000 description 4
- 230000009467 reduction Effects 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 101100136092 Drosophila melanogaster peng gene Proteins 0.000 description 2
- 230000009471 action Effects 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 229930002875 chlorophyll Natural products 0.000 description 2
- 235000019804 chlorophyll Nutrition 0.000 description 2
- ATNHDLDRLWWWCB-AENOIHSZSA-M chlorophyll a Chemical compound C1([C@@H](C(=O)OC)C(=O)C2=C3C)=C2N2C3=CC(C(CC)=C3C)=[N+]4C3=CC3=C(C=C)C(C)=C5N3[Mg-2]42[N+]2=C1[C@@H](CCC(=O)OC\C=C(/C)CCC[C@H](C)CCC[C@H](C)CCCC(C)C)[C@H](C)C2=C5 ATNHDLDRLWWWCB-AENOIHSZSA-M 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 238000012804 iterative process Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- 208000030853 Asthma-Chronic Obstructive Pulmonary Disease Overlap Syndrome Diseases 0.000 description 1
- 208000005156 Dehydration Diseases 0.000 description 1
- 241000196324 Embryophyta Species 0.000 description 1
- 241000283986 Lepus Species 0.000 description 1
- 240000000233 Melia azedarach Species 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 239000000443 aerosol Substances 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000002790 cross-validation Methods 0.000 description 1
- 230000001351 cycling effect Effects 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000007599 discharging Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000007716 flux method Methods 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 230000008014 freezing Effects 0.000 description 1
- 238000007710 freezing Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000012886 linear function Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 230000001932 seasonal effect Effects 0.000 description 1
- 239000000779 smoke Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000004800 variational method Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
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
-
- 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
- Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
- Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
- Y02P90/80—Management or planning
- Y02P90/84—Greenhouse gas [GHG] management systems
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供了一种同步反演高精度、高分辨率人为CO2排放与自然CO2通量的同化方法及系统。该方法包括:制备先验人为CO2排放、生物质燃烧的数据,驱动WRF‑GHG模型前向模拟CO2浓度场和自然CO2通量场;基于所设置的模型模拟区域范围内的CO2通量观测数据,对WRF‑GHG模型模拟的自然CO2通量结果进行优化,并对嵌套于WRF‑GHG模型中的自然生态CO2诊断模型VPRM模块进行参数优化;基于尺度匹配后的地基和多源卫星CO2浓度数据,采用POD4DVar数据同化算法,实现地基‑多源卫星联合同化,进而实现对人为CO2排放的反演核算;生成区域高时空分辨率数据产品集。本发明提供的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法及系统能够解决其无法区分人为CO2排放和自然CO2通量的技术难题。
Description
技术领域
本发明涉及大气科学、地理科学、全球变化科学技术领域,特别是涉及一种同步反演高分辨率人为CO2排放与自然CO2通量的同化方法及系统。
背景技术
中国政府明确提出2030年前碳达峰、2060年前碳中和的目标,从相对减排到绝对减排,进而零排放,成为雄心勃勃的“中国方案”与路线图。我国践行碳中和承诺会导致二氧化碳(CO2)的排放碳源和生态系统碳汇发生明显的变化,这种变化可以指示碳中和各项政策和各项计划实施的有效性,其中涉及与国际科学前沿和面向国家经济社会高质量发展需求的关键科学技术问题。为了科学评价碳中和各项政策和各项计划实施的有效性,急需甄别人为CO2排放和自然CO2通量的方法技术。量化全球和区域CO2通量非常重要。陆地生态系统复杂多样、观测数据时空尺度和精度上差异大等原因导致各种方法估算的陆地生态系统CO2源汇不确定性大(Piao et al.,2012;Schimel et al.,2015;Arneth et al.,2017;Bastos et al.,2020)。通过全球和区域尺度CO2收支估算不确定性的评价(Bastos etal.,2020;Le Quéréet al.,2018)发现,区域尺度上的不确定性相对于全球尺度要大(Martin et al.,2019;Hutchins et al.,2017;Fischer et al.,2017;Gately andHutyra,2017;Oda et al.,2018)。未来减排决策和缓解气候变化策略的制定中离不开高精度、高分辨率CO2源汇信息。因此,降低CO2源汇估算的不确定性,提高CO2通量的反演精度是全球变化研究的热点之一(Pillai et al.,2016;Turnbull et al.,2015;Ye et al.,2018;Wang et al.,2019;Oda et al.,2018;Wang et al.,2020)。
当前CO2通量的估算方法包括基于陆地观测的“自下而上(Bottom-up)”方法(于贵瑞等,2011;Tian et al.,2011;刘纪远等,2003;Asefi-Najafabady et al.,2014;Jianget al.,2016;Kondo et al.,2020)和基于大气观测的“自上而下(Top-down)”方法(Peterset al.,2007;Deng et al.,2016;Zhang et al.,2014;Peng et al.,2015;Wang etal.,2019;Smith et al.,2020;陈报章,张慧芳,2015)。“自上而下”的方法已成为目前CO2循环研究的主流方向(Peters et al.,2005,2007;Zhang et al.,2014;Tian et al.,2014;Peng et al.,2015;Wang etal.,2019)。自上而下的方法受观测数据数量和模型效率等因素的制约,其反演结果分辨率低、不确定性大(Bréon et al.,2015)。此外,大部分CO2通量同化反演方法是在假定自然或人为CO2通量的其中一方不存在误差的前提下去优化另一方,这种假设一定程度上增加所反演的CO2通量的不确定性并降低了其应用价值。传统上采用自下而上的物料平衡方法估算国家和行政单元尺度人为CO2的年总排放量,但该方法存在不确定性大且滞后性明显的问题。因此,迫切需要发展区域大气CO2反演法,以实现自然CO2通量和人为CO2排放的在线同步反演核算。
早期的大气反演模型并不是为估算人为碳排放而设计和开发的,而是在假设人为碳排放无误差(模型中固定不变)的情况下,利用大气CO2浓度观测来优化陆地或海洋自然CO2通量(Peters et al.,2007,2010;Jiang et al.,2016;Saeki et al.,2013;Zhang etal.,2014)。这种“人为碳排放无误差”假设存在着很大的缺陷性,因为大气CO2浓度变化是近地表自然和人为碳通量经过大气混合后的结果,设定人为碳排放无误差(实际上基于调查清单法估算出的人为CO2排放存在着极大的不确定性)必定会增大自然碳通量(陆地或海洋)估算的不确定性(Gurney et al.,2005)。近期有研究者尝试应用大气碳反演系统计算人为CO2排放量,其做法是:因为大气CO2浓度变化主要受“自然CO2通量”和“人为碳排放”影响,因此在假定存在一个具有相同“自然CO2通量”而没有“人为CO2排放”的背景区(background region)存在的前提下,利用背景区剔除自然CO2通量的影响,计算出大气CO2浓度增量△CO2。在理论上,△CO2被认为是由人为CO2排放量导致的是合理的。这样就可以通过降低模拟和观测的CO2浓度差值来纠正人为CO2排放估算的误差(Kort et al.,2012;Schneising et al.,2013;Wong et al.,2015;Bréon et al.et al.2015;Pillai et al.,2016;Martin et al.,2019)。选定一块自然条件相似(有相似的自然生态系统碳通量)而人烟稀少(近似本地人为源排放为0)的背景区来计算CO2浓度增量(△CO2),以区分出人为碳排放。这种改进后的大气CO2反演法被称之为一种人为碳排放反演法(Kort et al.2012;Schneising et al.2013;Wong et al.2015;Bréon et al.et al.2015;Pillai etal.2016;Ye et al.2020)。采用该方法已成功反演出法国(Staufer et al.et al.2016)、洛杉矶(Wong et al.,2015)、柏林(Hase et al.,2015)、韩国(Shim et al.2019)和北京(Feng et al.,2019)等城市人为CO2排放。但这种不受人为CO2排放影响的理想背景区的“假设”条件在地球上很难找到,实际上很难找到绝对的“背景区”来计算人为碳源导致的CO2浓度增量。因此,这种方法的应用价值十分有限。目前基于大气反演模型达到方法估算人为碳排放的结果仍存在着很大的不确定性,而且不能区分自然CO2通量和人为CO2排放,这是评估碳中和行动有效性急需解决的技术问题。
发明内容
本发明要解决的技术问题是提供一种同步反演高分辨率人为CO2排放与自然CO2通量的同化方法及系统,能够解决其无法区分人为CO2排放和自然CO2通量的技术难题。
为解决上述技术问题,本发明提供了一种同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,所述方法包括:解析遥感影像数据MODIS09A1数据,产生增强型植被指数和地表水指数数据集作为VPRM输入数据;制备先验人为CO2排放、生物质燃烧的数据,驱动WRF-GHG模型前向模拟CO2浓度场和自然CO2通量场;基于所设置的模型模拟区域范围内的CO2通量观测数据,对WRF-GHG模型模拟的自然CO2通量结果进行优化,并对嵌套于WRF-GHG模型中的自然生态CO2诊断模型VPRM模块进行参数优化;基于尺度匹配后的地基和多源卫星CO2浓度数据,采用POD4DVar数据同化算法,实现地基-多源卫星联合同化,进而实现对人为CO2排放的反演核算;生成区域高时空分辨率数据产品集,区域高时空分辨率数据产品集包括:CO2浓度场、自然CO2通量场、人为CO2碳排放场;对4维的区域高时空分辨率数据产品集进行准确性验证和模型鲁棒性验证;构建区域碳同化反演系统,实现同步同化目标区域高分辨、高精度人为CO2排放与自然CO2通量。
在一些实施方式中,对VPRM模块进行的参数优化包括:数据质量控制;根据站点通量足迹模型剔除周边人为源对通量观测数据产生较大影响的观测值;λ和PAR0的参数优化。
在一些实施方式中,数据质量控制包括:当夜间摩擦风速低于临界摩擦风速时观测的NEE数据剔除;剔除那些超出所设阈值,存在明显异常的NEE观测数据;剔除那些与降水同时段的NEE观测数据;对NEE观测数据进行统计进行3倍标准差的剔除筛选。
在一些实施方式中,λ和PAR0的参数优化包括:利用夜间的NEE观测值通过蒙特卡洛的方法得到呼吸参数α和β;利用α和β和站点温度观测值计算得到白天的生态系统呼吸Re,由白天的NEE减去Re得到GEE,;再利用VPRM中GEE的计算公式以及对应的Tscale、Wscale、Pscale、EVI、PAR值反演拟合得到λ和PAR0。
在一些实施方式中,根据站点通量足迹模型剔除周边人为源对通量观测数据产生较大影响的观测值,包括:采用涡度相关CO2通量的footprint模式对CO2通量观测数据进行筛选,剔除那些受人为源影响较为严重的CO2通量观测数据。
在一些实施方式中,地基和多源卫星CO2浓度数据的尺度匹配包括:卫星柱浓度尺度转换、多源数据融合及局地化技术。
在一些实施方式中,卫星柱浓度尺度转换包括:卫星柱浓度垂直尺度转换、多源卫星柱浓度的水平空间匹配和尺度转换。
在一些实施方式中,卫星柱浓度垂直尺度转换包括:把碳同化系统模拟的分层浓度插值到GOSAT/GOSAT-2、OCO-2/OCO-3或TanSat卫星的分层高度,以碳卫星先验廓线、气压加权函数、平均核函数为基础,采用Rodgers&Connor的算法,计算与卫星柱浓度相匹配的模型模拟柱浓度。
在一些实施方式中,多源卫星柱浓度的水平空间匹配和尺度转换包括:OCO-2/OCO-3卫星一次打点8个,1秒打3次,对OCO-2/OCO-3在同1秒内的XCO2观测数据求平均。
在一些实施方式中,局地化技术包括:分区;计算回归系数βi;计算权重αmin;计算结合观测空间代表性的权重αmin,obs;引入局地化函数后,执行同化分析。
在一些实施方式中,多源数据融合包括:进行分层浓度与柱浓度的转换;利用集合卡尔曼平滑EnKS方法对观测浓度进行统一同化。
此外,本发明还提供了一种同步反演高分辨率人为CO2排放与自然CO2通量的同化系统,所述系统包括:一个或多个处理器;存储装置,用于存储一个或多个程序;当所述一个或多个程序被所述一个或多个处理器执行,使得所述一个或多个处理器实现根据前文所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法。
采用这样的设计后,本发明至少具有以下优点:
本发明使用WRF-GHG区域模式,选择华北平原为试验区,结合禹城地基站点观测数据和GOSAT与OCO-2卫星柱浓度数据,运用POD4DVar的方法构建了区域高精度地基-多源卫星联合碳同化系统,实现了自然CO2通量和人为CO2排放的同步同化反演。
(1)在线耦合模型对VPRM模型参数优化明显。在WRF-GHG模式框架下实现了生态诊断模型VPRM与大气传输模型WRF的在线耦合,基于禹城通量观测数据优化VPRM模型参数,接着从站点和区域尺度验证、评价模拟的生态系统CO2交换(GEE)和生态系统呼吸(Re)。结果显示优化效果显著,说明参数优化之前模式对生态系统净CO2通量(NEE)存在明显低估问题。在区域尺度上把WRF-GHG碳同化系统参数优化后模拟的NEE(5km分辨率)与CarbonTracker-China的NEE(1°x1°分辨率)反演值进行对比,结果表明两套数据季节变化基本一致,但是4月份WRF-GHG模型模拟的NEE为碳汇,CarbonTracker-China反演的NEE为弱的碳源,从NEE通量大小的贡献量来看CarbonTracker-China的NEE反演值,总体低于WRF-GHG区域碳同化系统的模拟值,主要原因是WRF-GHG同化系统模拟的NEE精度更高,更能代表研究区域的CO2通量特征,能捕捉到农作物更多的生态信息。
(2)WRF-GHG同化前后CO2浓度模拟值结果验证。采用未参与同化的GOSAT和OCO-2卫星柱浓度数据和地基CO2浓度观测数据对同化前后的CO2浓度模拟值进行了对比分析,同化前后CO2浓度模拟值与GOSAT卫星柱浓度数据对比发现,在1月份和4月份相关性大约提高了17.7%,均方根误差(RMSE)平均降了66.7%;与OCO-2卫星柱浓度观测信息对比,其系统模拟值比同化之前效果有很大提高,相关性大约提高了17%,RMSE降低了大约1ppm,平均偏差大约降了37.2%,误差的另一个测量指标IOA大约提高了30%。同化后的站点模拟值与观测值之间的相关性提高了40%,同化后的平均偏差大约降了51.87%。
(3)WRF-GHG同化前后人为CO2通量估算精度明显提高。在城市尺度上,以济宁、邯郸、淄博3个排放量大的城市为例进行了CO2排放量验证统计,并与其他独立数据,即生态环境部规划院的中国城市2012年CO2排放数据集进行对比分析,结果显示排放量基本一致,并发现其他基于传统方法的排放数据存在明显低估问题。在区域尺度上同化系统反演的人为碳通量同清华大学的MEIC通量产品进行对比,结果显示同化系统反演的人为碳通量更能体现季节的变化规律,在基于本发明技术计算反演值排放量上总体高于MEIC排放数据。
(4)基于本发明,实现了自然CO2通量和人为CO2排放的同步同化反演并优化,提高了传统人为CO2排放清单空间分辨率粗放和严重时间滞后性问题,将成为正在实施的碳中和路径的不确定性评估的重要工具,为国家和省、市碳中和行动有效性评估提供科技支撑。
附图说明
上述仅是本发明技术方案的概述,为了能够更清楚了解本发明的技术手段,以下结合附图与具体实施方式对本发明作进一步的详细说明。
图1是同步同化人为排放与自然碳通量总技术路线图;
图2是耦合生态系统模型(VPRM)-大气传输模型优化碳关键参数技术流程;
图3是耦合生态系统模型(VPRM)关键参数优化技术流程图;
图4是CO2柱浓度尺度转换技术示意图;
图5是WRF-GHG模式中新增模块的代码结构;
图6是区域碳同化系统运行流程图;
图7是示例嵌套模拟区域设置;
图8是禹城站周边区域地表覆被;
图9是同化系统运行周期示意图;
图10是不同集合下区域碳同化系统CO2浓度模拟值与观测值对比;
图11是WRF默认MODIS地表覆被类型与LUCC 2010年地表覆被类型分布图;
图12是WRF默认MODIS地表覆被数据与LUCC 2010年地表覆被数据对比分析直方图;
图13是VPRM模式优化输出流程图;
图14是春季GPP优化前后与观测值对比分析图;
图15是春季Re优化前后与观测值对比分析图;
图16是春季NEE优化前后与观测值对比分析图;
图17是夏季GPP优化前后与观测值对比分析图;
图18是夏季Re优化前后与观测值对比分析图;
图19是夏季NEE优化前后与观测值对比分析图;
图20是秋季GPP优化前后与观测值对比分析图;
图21是秋季Re优化前后与观测值对比分析图;
图22是秋季NEE优化前后与观测值对比分析图;
图23是冬季GPP优化前后与观测值对比分析图;
图24是冬季Re优化前后与观测值对比分析图;
图25是冬季NEE优化前后与观测值对比分析图;
图26是2016年四个月份(1、4、7和10月份)中午12:00区域尺度(5km)平均光合作用强度;
图27是2016年四个月份(1、4、7和10月份)中午12:00区域尺度(5km)平均呼吸作用强度;
图28是2016年四个月份(1、4、7和10月份)中午12:00区域尺度(5km)平均地表温度;
图29是2016年四个月份(1、4、7和10月份)GPP、Re和NEE对研究区域贡献量的月均值;
图30是2016年四个月份(1、4、7和10月份)NEE对研究区域的贡献量;
图31是同化系统反演的区域高分辨率(5km)人为源碳通量;
图32是城市尺度排放量对比分析;
图33是区域人为源碳通量统计分析。
具体实施方式
以下结合附图对本发明的优选实施例进行说明,应当理解,此处所描述的优选实施例仅用于说明和解释本发明,并不用于限定本发明。
技术方案的目标是发展同化算法,解决其无法区分人为CO2排放和自然CO2通量的技术难题。核心技术包括如下两个部分:
(1)本发明通过正交4维变分(POD4DVar:POD-based four-dimensionalvariation data assimilation method)方法构造并求解人为CO2排放源和CO2浓度的雅克比矩阵,实现目标函数最小化,采取多尺度交叉验证方法,实现CO2通量和CO2浓度双同化,并基于站点通量观测数据、采用蒙特卡洛参数优化方法对生态诊断模型参数进行优化。
(2)目前碳同化系统只是简单地把陆地生态模型当做一个先验通量的提供者,而没有实现大气传输模型与陆地生态系统模型的耦合,只对净碳通量(NEE)进行同化反演,因而无法对两个碳通量分量生态系统初级生产力(GPP:gross primary productivity)和生态系统呼吸(Re:ecosystem respiration)进行直接反演,更无法对生态系统模型中关键参数(光合参数和呼吸参数)进行优化。该同化系统基于生态诊断模型和大气传输模型的在线耦合模式WRF-VPRM(WRF-GHG模式中针对大气CO2模块)实现自然碳通量(GEE与Re)和人为CO2排放同步同化。
本发明总体技术方案如图1所示。
(1)解析遥感影像数据MODIS09A1数据,产生增强型植被指数(EVI)和地表水指数(LSWI)数据集作为VPRM输入数据;
(2)制备先验人为CO2排放、生物质燃烧等数据,驱动WRF-GHG模型前向模拟CO2浓度场和自然CO2通量场;
(3)基于所设置的模型模拟区域范围内的CO2通量观测数据,对WRF-GHG模型模拟的自然CO2通量结果进行优化,并对嵌套于WRF-GHG模型中的自然生态CO2诊断模型VPRM模块进行参数优化,参数优化后的自然源CO2通量误差很小并为人为CO2排放的同化反演提供条件;
(4)基于尺度匹配后的地基和多源卫星CO2浓度数据,采用POD4DVar数据同化算法,实现地基-多源卫星联合同化,进而实现对人为CO2排放的反演核算;
(5)生成区域高时空分辨率数据产品集(CO2浓度场、自然CO2通量场、人为CO2碳排放场);
(6)对(5)的4维数据集(三维立体空间+时间)进行准确性验证和模型鲁棒性验证;
(7)最后构建区域碳同化反演系统(RCAS),实现同步同化目标区域高分辨、高精度人为CO2排放与自然CO2通量。
所涉及的四个方面的核心方法及关键技术叙述如下。
3.2.1生态系统模型-大气传输模型耦合技术
基于碳卫星太阳诱导叶绿素荧光(SIF)数据、13CO2同位素观测数据以及生态系统碳通量观测数据,采用卡尔曼滤波(Kalman Filter)和马尔可夫链—蒙特卡罗(MarkovChain-Monte Carlo)同化方法,实现陆面生态模型与大气传输模型耦合,对陆面生态模型VPRM中碳关键过程进行参数优化,实现生态系统CO2通量分量(生态系统初级生产力GEE和生态系统呼吸Re)及关键碳过程参数(最大光能利用率LUEmax、呼吸温度敏感参数(Q10)等)的同步优化。
陆面生态系统模型的优化是获取高精度先验自然碳通量的基础。利用站点CO2通量及叶绿素荧光(SIF)观测资料,结合野外调查材料,采用蒙特卡洛算法和资料同化技术(EnKF),实现对陆面生态模型VPRM中碳关键过程进行参数的优化、对比和验证,实现陆面生态模型在站点尺度的优化(图2)。
3.2.2陆面生态系统模型(VPRM)参数优化方法
植被光合呼吸模型(VPRM)是一种基于光能利用效率(LUE)的陆地生态诊断模型,它是Mahadevan在VPM模型的基础上发展而来,与VPM相比,VPRM增加呼吸项(Re)及其相关的呼吸参数α和β(在冰点温度时呼吸基础率),以及用于反映光合有效辐射与光合作用关系的半饱和值参数PAR0。该研究以WRF模型耦合VPRM生态诊断模型为基础,基于通量观测数据结合模拟自然碳通量值,反演光合、光合敏感参数,使用优化后的模型参数预测中国区域不同时空分辨率陆表碳通量分量的时空格局,基于多尺度地面观测与卫星产品,交差验证模拟结果,综合评价国内外现有碳同化系统和碳循环模型及参数,并为进一步改善WRF-VPRM模拟区域大气CO2浓度的时空变化打下坚实的基础。揭示自然与人为因素对大气CO2浓度的相对历史贡献,分析其驱动机制与未来的变化趋势。
如图3所示,在VPRM模型中,NEE的计算主要包括两部分:其中一部分为总生态系统CO2交换量GEE,由光照驱动计算;另一部分为生态系统呼吸项Re,由温度驱动计算。其具体表达式如下:
式中,λ(μmol CO2m-2s-1/(μmol PARm-2s-1))代表最大光能利用率(或最大光量子效率),该参数在无水分压迫和温度适中的情况下对C3植物而言为1/6左(Mahadevan等.,2008),PAR0(μmolm-2s-1)代表光合作用达到半饱和值时对应的光合有效辐射。PAR(μmolm-2s-1)代表光合有效辐射,EVI为增强型植被指数,PAPARPAV代表植被所吸收的光合有效辐射和光合辐射之间的比值(Xiao等.,2004)。Tscale、Wscale和Pscale代表对光合作用的影响函数,分别代表温度、水分胁迫和叶面性状,其计算公式如下:
式中,T代表气温,单位是(℃),Tmin、Tmax和Topt分别代表光合作用所需要的最小、最大和最适温度。当空气温度低于Tmin,Tscale不在采用公式(4)计算,直接取值为0。LSWI代表地表水分指数,LSWImax是每个格网点生长季内最大的LSWI值。植被处在不同的生长期阶段,其Pscale的值不同,模式中对于常绿林Pscale在全年中的值设为1.0。Tscale、Wscale和Pscale3个函数值有一定的取值范围,其范围为0.0~1.0。公式(2.1)最后1项代表呼吸过程对NEE的贡献。多数生态诊断模型把呼吸项看成是温度的指数函数,但是VPRM模式将呼吸项简化为温度的线性函数。其中α(μmol m-2s-1/℃)、β(μmol m-2s-1/℃)可根据通量塔观测数据进行优化调整。
为了使VPRM能够更好的对研究区域进行生态碳通量的模拟,需要对VPRM参数进行优化。参照ChinaFLUX通量数据后处理方法(http://www.chinaflux.org)以及Falge等(2001)和Zhang等(2006)的方法对通量塔检测的NEE数据进行一系列预处理,主要包括运用相关插值算法对缺失数据进行插补、坐标轴旋转、摩擦风速大小的质量控制、储存项计算和密度纠正计算等。参数优化时选择未插补和经过严格摩擦风速筛选的数据以避免数据处理不当引入误差。
具体地,采用如下方案进行参数优化:
(1)数据质量控制:当夜间摩擦风速低于临界摩擦风速时观测的NEE数据剔除;剔除那些超出所设阈值,存在明显异常的NEE观测数据;剔除那些与降水同时段的NEE观测数据;对NEE观测数据进行统计进行3倍标准差的剔除筛选。
(2)根据站点通量足迹模型(Footprrint模式)(Chen等.,2009)剔除周边人为源对通量观测数据产生较大影响的观测值。
(3)参数优化:利用夜间的NEE观测值通过蒙特卡洛的方法得到呼吸参数α和β。利用α和β和站点温度观测值计算得到白天的生态系统呼吸Re,由白天的NEE减去Re得到GEE,本文提到的GEE数值大小和GPP一致,只是符号相反,下同。再利用VPRM中GEE的计算公式以及对应的Tscale、Wscale、Pscale、EVI、PAR值反演拟合得到λ和PAR0
以华北平原应用案例为例,参与参数优化的禹城站通量观测数据需要严格的质量控制,除了上面提到的常规的质量控制以外,本研究采用涡度相关CO2通量的footprint模式对CO2通量观测数据进行筛选,剔除那些受人为源影响较为严重的CO2通量观测数据。
3.2.3卫星柱浓度尺度转换、多源数据融合及局地化技术
(1)CO2浓度观测数据
项目申请依托中国精度CO2浓度全球大气本底站、我国区域背景站、校验站数据,收集整理的全球站点(https://gml.noaa.gov/ccgg/obspack/)和多源卫星观测数据分别来自GOSAT(下载路径:https://disc.gsfc.nasa.gov/datasets?keywords=GOSAT&page=1,ACOS-GOSAT V9r)、GOSAT-2(下载路径:https://prdct.gosat-2.nies.go.jp/aboutdata/howtoaccessdata.html.en#01_howto)、OCO-2(下载路径:https://disc.gsfc.nasa.gov/OCO-2,OCO2_L2_Lite_FP V10r)、OCO-3(下载路径:https://disc.gsfc.nasa.gov/OCO-3,OCO3_L2_lite_FP EarlyR)和TanSat(下载路径:http://data.nsmc.org.cn/portalsite/Data/Satellite.aspx)。利用多源卫星观测开展碳源汇估算以及大气CO2反演法优化研究时,要求碳卫星XCO2测量精度达到1ppm(-0.3%)。因此,我们从GOSAT、OCO-2及TanSat卫星观测数据中过滤出误差<1ppm的数据进入碳同化反演系统。
(2)卫星柱浓度垂直尺度转换
在将卫星数据引入碳同化系统之前,首先要对卫星柱浓度数据与大气传输模型模拟浓度间的尺度进行转化,实现观测浓度与模拟浓度间的尺度匹配。目前已有多种方法用于实现模拟浓度的尺度转换,按照GOSAT/GOSAT-2、OCO-2/OCO-3、TanSat卫星参数及本项目的尺度要求,选择卫星核函数尺度转化法(图4),实现模拟浓度与卫星柱浓度的尺度匹配。尺度转换过程如下:把碳同化系统模拟的分层浓度插值到GOSAT/GOSAT-2、OCO-2/OCO-3或TanSat卫星的分层高度,以碳卫星先验廓线、气压加权函数、平均核函数为基础,采用Rodgers&Connor的算法(Rodgers CD,Connor BJ(2003)Intercomparison of remotesounding instruments.J Geophys Res.doi:10.1029/2002JD002299),计算与卫星柱浓度相匹配的模型模拟柱浓度:
其中,是转化后的模拟柱浓度;hT是卫星(GOSAT/GOSAT-2、OCO-2/OCO-3、TanSat)提供的气压加权函数;Xa是卫星(GOSAT/GOSAT-2、OCO-2/OCO-3、TanSat)的先验廓线;A分别是GOSAT/GOSAT-2、OCO-2/OCO-3、TanSat的平均核函数;Xh是大气传输模型在卫星分层高度插值获取到的模拟浓度。对于卫星本身不能提供核函数、先验廓线等信息的,我们将借助其它方法(如TCCON柱浓度转换法)或其它模型参数(如RCAS同化系统提供的先验廓线信息)按公式(2.1)进行转化。
(3)多源卫星柱浓度的水平空间匹配和尺度转换
多源碳卫星的空间分辨率各异,GOSAT、GOSAT-2、OCO-2/OCO-3和TanSat空间分辨率分别为10.5km×10.5km、9.7km×9.7km、1.29km×2.25km和2km×2km,需要根据模式输出格网的大小进行空间分辨率统一化处理。根据模式分辨率选择与GOSAT和GOSAT-2空间分辨率接近的10km×10km作为多源碳卫星空间匹配标准。OCO-2/OCO-3和TanSat分辨率远大于10km×10km,需进行空间尺度匹配,其中,OCO-2/OCO-3卫星一次打点8个(1.29km×8个≈10km),1秒打3次(24个点,2.25km×3次≈7km),我们对OCO-2/OCO-3在同1秒内(<=24个点,注意部分点会由于气溶胶污染(一般光学厚度AOD>0.3时)被过滤掉)的XCO2观测数据求平均,以尽量使得卫星观测点向空间匹配标准靠近。对于TanSat,我们采取类似于OCO-2/OCO-3的空间处理方式。综合水平空间匹配和垂直尺度转换的模拟浓度计算公式如下:
其中,i代表的是第i个卫星观测,n为1秒内获取的卫星观测总数(对于OCO-2/OCO-3,在假设没有卫星点被过滤掉的情况下,其n=24)。式中其它参数说明请参考公式(1)。
(4)多源观测的自适应局地化方案
由于站点、卫星等多源观测数据的空间分辨率及时空代表性不相同,其观测影响周围状态量的空间范围也不相同。本发明技术结合观测数据类型、格网距离信息采用分级集合滤波技术来研制考虑多源观测数据同化影响范围的局地化方案。局地化技术及过程如下:
1)分区:为减小计算量、提高计算效率,将整个模拟区域均匀划分为m个小区(每个方块中包含有n个成员,所以集合成员总数为m×n),每个小区的中心点选定为“局地化中心”。
2)计算回归系数βi:对模拟区分区后,开始计算各小区的回归系数。任一个观测对小区i的回归系数βi计算公式如下:
式中分子σx,y表示由第i组集合计算得到的x和y先验样本协方差,分母σy,y表示观测的先验方差。分子σx,y和分母σy,y对应为数据同化中的PHT和HPHT(H为观测算子,P为背景误差协方差矩阵,上标T表示转置,x表示状态变量,y表示观测)。
3)计算权重αmin:假设忽略其它的误差源,回归系数βi的真值应从一个回归系数样本的相同分布中随机获得。样本的不确定性反映了计算得到增量的不确定性。因此,我们定义回归置信因子α,使得增量值和利用准确计算得到的增量值之间的均方根误差数学期望最小化:
即最小化公式:
等价于:
对α求导并令其为0:
最终求解得到置信因子αmin,即作为观测对状态量的影响权重系数:
4)计算结合观测空间代表性的权重αmin,obs:根据空间代表性,对站点、卫星观测数据进行空间代表性分级,并给定代表性系数ξ。定义站点观测代表性系数ξ为1,GOSAT/GOSAT-2代表性系数ξ为0.7,OCO-2/OCO-3/TanSat的代表性系数ξ为0.5。更新后的权重αmin,obs如下:
αmin,obs=αmin×ξ (3.8)
5)在同化系统中的技术实现方法:αmin,obs为观测对模式状态量的局地化权重系数。引入局地化函数后,同化分析方程变为:
该方法按照观测与背景场模式状态量间误差的相关性给予不同等级的局地化影响权重,因此,将这种分级集合滤波同化技术是高效的自适应局地化方案。
(5)卫星-地基CO2数据融合技术
卫星监测资料具有覆盖面广、时空分辨率高的特点。它的引进改善现有的大气CO2反演模型观测数据不足的现状。卫星-地基CO2数据联合同化方法的关键技术是如何把模拟的CO2浓度转化为与卫星浓度同一个尺度的柱平均浓度,以方便模型进行对比、分析及数据同化。我们利用公式(3.10)(Rodgers et al.(2003)进行分层浓度与柱浓度的转换。
其中,是我们转化后的模拟平均柱浓度;hT是卫星提供的气压加权函数;Xa是卫星提供的先验廓线;A是卫星提供的平均核函数;Xh是大气传输模型在与卫星相对应的分层高度插值获取到的模型浓度。
将分层模拟浓度转化为柱浓度后,利用集合卡尔曼平滑EnKS方法对观测浓度(包含地基观测和卫星观测)进行统一同化(图4)。
整个技术路线分2步走:
(1)根据观测点(x,y,t,z)的坐标(x,y)、时间(t)、高度(z)信息提取相应的模拟浓度,首先判断当前的观测数据类型:如果是站点观测点,系统直接采样提取相应浓度数据;如果不是站点数据,则要根据公式(3.10)对大气传输模型模拟的分层浓度进行转化,计算出与观测数据相对应的柱平均浓度。
(2)再通过最小化观测与模拟浓度间的差值来同化CO2的浓度与通量。
3.3.4POD4DVar高效CO2同化方法
本研究将本征正交分解法(POD,Proper Orthogonal Decomposition)应用到四维变分同化(4DVar,Four-Dimensional Variational Data assimilation)算法中(Tian等.,2009,2011),结合WRF-GHG模型,构建了基于本征正交分解的四维变分同化反演模型,在很大程度上降低了同化反演所需计算机资源,提高了运算效率和估算精度。该模型考虑了CO2滞后特征,可同化地基、卫星等多源CO2浓度观测数据,反演得到全球和区域等不同尺度的碳通量。
POD4DVar是一种基于本征正交分解的显式四维变分方法,它运用蒙特卡洛的方法生成类似于集合卡尔曼滤波方法(EnKF)的四维样本集合。由于POD正交基具有最小二乘意义上的最优性,将POD方法应用于四维空间的预报集合提取正交基,它可以捕获预报集合的更多能量,更好的表征四维变量的空间结构和时间演变特征。该方法无需积分伴随模式,在运行和维护上都较为简便。
由于传统的4DVar数据同化方法需要求取切线和伴随线性算子,然而对于非线性的模型算子求取切线和伴随算子是非常困难的。POD4DVar能够把隐式的优化问题显式化,用一组基向量捕捉数据的时空演变特征。这种方法不仅简化了数据同化的过程,而且保留了传统4DVar的两大优势:(1)该方法可以提供时间平滑性约束(2)在一个同化窗口中可以把多个时次的观测数据引入。尤其当初始场误差和预报模型不确定性明显的模式运用POD4DVar方法进行数据同化,其同化效果尤为明显。
传统的4Dvar,通过模型的模拟值与观测值之间的差异构造代价函数(公式4.1),对代价函数求导获取后验通量
其中T代表矩阵转置,b为背景场CO2人为碳通量值,k代表CO2浓度观测时间序列,Hk代表观测算子,矩阵B和R分别代表通量背景误差和浓度观测误差协方差。代表同化窗口初始时刻的控制变量(即先验通量)。在代价函数中/>通过公式(4.2)与/>联系在一起。很显然代价函数关于/>的梯度值很难获取,然而POD4DVar方法大大简化了其计算过程。在POD4DVar同化算法中,假设在同化窗口(0,T)中有S个时间序列,用蒙特卡洛方法在t=t0时刻对初始人为碳通量数据添加扰动,生成N个初始场集合/>n=1,2,…,N。通过公式(4.2)获得所用同化窗口内的,所有集合的状态变量:
当扰动样本足够大时,样本集合可以基本表征状态变量的时空特征,所有的状态变量的扰动信息能被扩展到一个有限维的向量空间同样最后的目标分析资料/>(i=0,1,2,…S-1),在整个同化窗口内能被表示为:
当集合数N增加到一定数量,分析状态向量能够被集合空间所表示假设Ω的基向量线性空间用/>(n=1,2,…K,K≤N)表示,分析状态向量/>能够表达为基向量/>的线性组合:
用公式(4.4)和(4.5)重写公式(4.1),由β=(β1...βK)T替换,因此状态变量在代价函数中被显示化,梯度计算在很大程度上被简化。切线模型和伴随模型不再需要。最小代价函数公式(4.1)用变量β显式处理之后,其能用一般的优化算法计算。但是与EnKF不同的是,此算法在每一次的同化循环中,初始状态变量在初始时刻都应该添加扰动,并且最后的优化结果/>只能获得一组。接下来的重点就是怎样获取合适的基向量。
N个集合状态向量的平均值表示如下:
根据构建一个新的偏差向量集合:
其构建的矩阵A(M×N),M=Mg×Mv×S,其中Mg,Mv分别表示模式空间格点的数目和模型变量的数目(这里仅代表人为碳通量)。为了计算本征正交分解(POD)的数目,需要解决特征值问题:
(AAT)M×MV=λV... (4.8)
在实际中M≥N,因此为了解决高效方便的求解问题,将其转化为N×N的特征值问题,转化方式如下:
ATAV=λV
AATAV=AλV
AAT(AV)=λ(AV)... (4.9)
此时转化为求取N×N的特征值问题:
TVk=λkVk,k=1,...,N... (4.10)
其中T=(ATA)N×N,Vk是特征向量V的第k列,λk是特征值λ第k行。由于T为对称半正定的,因此得到一组降序排列的特征值λ1≥λ2≥…λN≥0。非零特征值λk(1≤k≤N)被挑选用来标准正交化,并且POD模式被给出,通过
在四维空间中可以用POD重写,其表示方式为:
其中P为POD模式的数目,其大小的确定按照下面的公式:
其中,0<γ<1表示由降维空间Dp=span{φ1,...φp}所捕获的总能量的百分数。为了捕获POD基的大部分能量,允许γ在1附近取值。根据公式(4.4)和(4.9)重构代价函数(4.1),控制变量β=(β1,β2,...,βp)T代替初始时刻的状态变量如此控制变量在代价函数中被显式处理,在求取代价函数最小值的过程中,切线算子和伴随模式不在需要。求取代价函数的最小值问题,转化成用一般的优化算法求取β=(β1,β2,...,βp)T的问题。然而β的求取仍然需要迭代过程,依旧有很高的计算量。此问题被解决按照如下流程:
构建POD模式矩阵:
Φ=(φ1,φ2,…,φp)... (4.13)
其中φj=(φj(t0),φj(t1),...,φj(tS-1))T,j=1,2,...,P转化公式(4.13)为:
Φ=(φ1,Φ2,...,ΦS-1)T... (4.14)
其中Φk=(φ1(tk),φ2(tk),...,φp(tk)),公式(4.11)被重写如下:
其中β=(β1,β2,...,βp)T。用公式(4.15)重写代价函数:
其中,Hj是切线观测算子,CO2浓度观测向量其中mj是观测向量的个数,可用通过扰动的方式使某一时刻的观测值产生N个集合,其表示为:
式中,N个集合的扰动平均值为0,用矩阵表示为:
Ej=(ε1,ε2,...,εN),观察误差矩阵能够按照公式(4.18)进行评估。
背景误差协方差B的构造方式类似于观测误差协方差R的构造方式。
由于是对称矩阵,代价函数的梯度通过简单的计算就很容易获得,
为了求取代价函数最小值,需要,即:
经过以上的一系列转化之后,公式(4.20)可以不需要迭代过程直接计算获取,这就大大简化了计算量,提高了同化反演效率。
以上所述,仅是本发明的较佳实施例而已,并非对本发明作任何形式上的限制,本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰,均落在本发明的保护范围内。
Claims (12)
1.一种同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,其特征在于,包括:
解析遥感影像数据MODIS09A1数据,产生增强型植被指数和地表水指数数据集作为VPRM输入数据;
制备先验人为CO2排放、生物质燃烧的数据,驱动WRF-GHG模型前向模拟CO2浓度场和自然CO2通量场;
基于所设置的模型模拟区域范围内的CO2通量观测数据,对WRF-GHG模型模拟的自然CO2通量结果进行优化,并对嵌套于WRF-GHG模型中的自然生态CO2诊断模型VPRM模块进行参数优化;
基于尺度匹配后的地基和多源卫星CO2浓度数据,采用POD4DVar数据同化算法,实现地基-多源卫星联合同化,进而实现对人为CO2排放的反演核算;
生成区域高时空分辨率数据集,区域高时空分辨率数据集包括:CO2浓度场、自然CO2通量场、人为CO2碳排放场;
对4维的区域高时空分辨率数据产品集进行准确性验证和模型鲁棒性验证;
构建区域碳同化反演系统,实现同步同化目标区域高分辨、高精度人为CO2排放与自然CO2通量;
采用POD4DVar数据同化算法,实现地基-多源卫星联合同化,进而实现对人为CO2排放的反演核算,包括:
结合WRF-GHG模型构建基于本征正交分解的四维变分同化反演模型;
将POD方法应用于四维空间的预报集合提取正交基;
当扰动样本足够大时,样本集合基本表征状态变量的时空特征,所有的状态变量的扰动信息被扩展到有限维的向量空间
在四维空间中将用POD重写,其表示方式为:
其中,为目标分析状态向量;/>为N个集合状态向量的平均值;p为POD模式的数目;βj为第j个观测向量的控制变量;
构建POD模式矩阵:
Φ=(φ1,φ2,…,φp)
基于φj=(φj(t0),φj(t1),…,φj(ts-1))T,j=1,2,…,P,转化公式为:
Φ=(Φ1,Φ2,..,ΦS-1)T
基于Φk=(φ1(tk),φ2(tk),…,φp(tk)),将用POD重写的公式重写如下:
基于β=(β1,β2,...,βp)T重写代价函数:
其中,上标T代表矩阵转置;下标b为背景场;β为状态变量;Ф为POD模式矩阵;J(β)为重写的代价函数;为同化窗口t0时刻的控制变量;/>为CO2通量的背景场;Φ0为同化窗口初始时刻的POD矩阵;Φj为与第j个观测向量相匹配的POD矩阵;B为通量背景误差;m为模拟区域均匀划分数;Hj为切线观测算子;/>为同化窗口tj时刻的控制变量;Rj为浓度观测误差协方差;/>为CO2浓度观测向量,表达式为/>其中,mj是观测向量/>的个数,可用通过扰动的方式使某一时刻的观测值产生N个集合。
2.根据权利要求1所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,其特征在于,对VPRM模块进行的参数优化包括:数据质量控制;根据站点通量足迹模型剔除周边人为源对通量观测数据产生较大影响的观测值;植物光合与呼吸过程的参数优化。
3.根据权利要求2所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,其特征在于,数据质量控制包括:当夜间摩擦风速低于临界摩擦风速时观测的NEE数据剔除;剔除那些超出所设阈值,存在明显异常的NEE观测数据;剔除那些与降水同时段的NEE观测数据;对NEE观测数据进行统计进行3倍标准差的剔除筛选。
4.根据权利要求2所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,其特征在于,λ和PΛR0的参数优化包括:利用夜间的NEE观测值通过蒙特卡洛的方法得到呼吸参数α和β;利用α和β和站点温度观测值计算得到白天的生态系统呼吸Re,由白天的NEE减去Re得到GEE,;再利用VPRM中GEE的计算公式以及对应的Tscale、Wscale、Pscale、EVI、PAR值反演拟合得到λ和PΛR0。
5.根据权利要求2所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,其特征在于,根据站点通量足迹模型剔除周边人为源对通量观测数据产生较大影响的观测值,包括:
采用涡度相关CO2通量的footprint模式对CO2通量观测数据进行筛选,剔除那些受人为排放CO2源影响较为严重的观测数据。
6.根据权利要求1所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,其特征在于,地基和多源卫星CO2浓度数据的尺度匹配包括:卫星柱浓度尺度转换、多源数据融合及局地化技术。
7.根据权利要求6所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,其特征在于,卫星柱浓度尺度转换包括:卫星柱浓度垂直尺度转换、多源卫星柱浓度的水平空间匹配和尺度转换。
8.根据权利要求7所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,其特征在于,卫星柱浓度垂直尺度转换包括:
把碳同化系统模拟的分层浓度插值到GOSAT/GOSAT-2、OCO-2/OCO-3或TanSat卫星的分层高度,以碳卫星先验廓线、气压加权函数、平均核函数为基础,采用Rodgers&Connor的算法,计算与卫星柱浓度相匹配的模型模拟柱浓度。
9.根据权利要求7所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,其特征在于,多源卫星柱浓度的水平空间匹配和尺度转换包括:
OCO-2/OCO-3卫星一次打点8个,1秒打3次,对OCO-2/OCO-3在同1秒内的XCO2观测数据求平均。
10.根据权利要求6所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,其特征在于,局地化技术包括:分区;计算回归系数βi;计算权重αmin;计算结合观测空间代表性的权重αmin,obs;引入局地化函数后,执行同化分析。
11.根据权利要求6所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法,其特征在于,多源数据融合包括:
进行分层浓度与柱浓度的转换;
利用集合卡尔曼平滑EnKS方法对观测浓度进行统一同化。
12.一种同步反演高分辨率人为CO2排放与自然CO2通量的同化系统,其特征在于,包括:
一个或多个处理器;
存储装置,用于存储一个或多个程序;
当所述一个或多个程序被所述一个或多个处理器执行,使得所述一个或多个处理器实现根据权利要求1至11任意一项所述的同步反演高分辨率人为CO2排放与自然CO2通量的同化方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210642075.4A CN114970184B (zh) | 2022-06-07 | 2022-06-07 | 同步反演高分辨率人为co2排放与自然co2通量的同化方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210642075.4A CN114970184B (zh) | 2022-06-07 | 2022-06-07 | 同步反演高分辨率人为co2排放与自然co2通量的同化方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114970184A CN114970184A (zh) | 2022-08-30 |
CN114970184B true CN114970184B (zh) | 2024-04-02 |
Family
ID=82960875
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210642075.4A Active CN114970184B (zh) | 2022-06-07 | 2022-06-07 | 同步反演高分辨率人为co2排放与自然co2通量的同化方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114970184B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115808501B (zh) * | 2022-11-23 | 2024-06-21 | 武汉大学 | 一种基于oco-2卫星反演发电厂碳排放强度的方法 |
CN115907178B (zh) * | 2022-11-30 | 2023-12-15 | 中国地质大学(武汉) | 一种净生态系统co2交换量的预测方法 |
CN116758999A (zh) * | 2023-08-14 | 2023-09-15 | 中科三清科技有限公司 | 动态植被碳通量模型参数确定方法、装置和电子设备 |
Citations (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
FR2953021A1 (fr) * | 2009-11-26 | 2011-05-27 | Tanguy Griffon | Methode de mesure des emissions hebdomadaires et annuelles d'un gaz a effet de serre sur une surface donnee |
CN104834823A (zh) * | 2015-05-18 | 2015-08-12 | 中国科学院地理科学与资源研究所 | 基于卫星-地基co2数据联合同化的碳源汇估测方法 |
CN111257241A (zh) * | 2020-01-20 | 2020-06-09 | 中国科学院地理科学与资源研究所 | 一种deei的基于卫星观测的大气二氧化碳浓度反演算法 |
CN111487216A (zh) * | 2020-05-06 | 2020-08-04 | 北京中科锐景科技有限公司 | 一种二氧化碳通量反演方法、系统 |
CN111723482A (zh) * | 2020-06-17 | 2020-09-29 | 南京大学 | 一种基于卫星co2柱浓度观测反演地表碳通量的方法 |
CN111881569A (zh) * | 2020-07-24 | 2020-11-03 | 中国科学院大气物理研究所 | 二氧化碳柱浓度的反演方法、装置、存储介质及电子设备 |
CN111914388A (zh) * | 2020-06-06 | 2020-11-10 | 中南大学 | 基于时平均最低浓度的城区自然因素co2浓度的计算方法 |
CN112597651A (zh) * | 2020-12-22 | 2021-04-02 | 武汉大学 | 基于oco-2数据和wrf-stilt模型反演co2背景场浓度的方法及系统 |
CN113406667A (zh) * | 2021-07-07 | 2021-09-17 | 武汉大学 | 基于星载激光雷达获取城市co2格网化通量的方法及系统 |
CN114547553A (zh) * | 2022-04-27 | 2022-05-27 | 河北先河环保科技股份有限公司 | 二氧化碳排放量的反演方法、装置、设备及存储介质 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20140257782A1 (en) * | 2013-03-08 | 2014-09-11 | The Penn State Research Foundation | System and method for measurement of temporal changes in trace gas fluxes |
-
2022
- 2022-06-07 CN CN202210642075.4A patent/CN114970184B/zh active Active
Patent Citations (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
FR2953021A1 (fr) * | 2009-11-26 | 2011-05-27 | Tanguy Griffon | Methode de mesure des emissions hebdomadaires et annuelles d'un gaz a effet de serre sur une surface donnee |
CN104834823A (zh) * | 2015-05-18 | 2015-08-12 | 中国科学院地理科学与资源研究所 | 基于卫星-地基co2数据联合同化的碳源汇估测方法 |
CN111257241A (zh) * | 2020-01-20 | 2020-06-09 | 中国科学院地理科学与资源研究所 | 一种deei的基于卫星观测的大气二氧化碳浓度反演算法 |
CN111487216A (zh) * | 2020-05-06 | 2020-08-04 | 北京中科锐景科技有限公司 | 一种二氧化碳通量反演方法、系统 |
CN111914388A (zh) * | 2020-06-06 | 2020-11-10 | 中南大学 | 基于时平均最低浓度的城区自然因素co2浓度的计算方法 |
CN111723482A (zh) * | 2020-06-17 | 2020-09-29 | 南京大学 | 一种基于卫星co2柱浓度观测反演地表碳通量的方法 |
CN111881569A (zh) * | 2020-07-24 | 2020-11-03 | 中国科学院大气物理研究所 | 二氧化碳柱浓度的反演方法、装置、存储介质及电子设备 |
CN112597651A (zh) * | 2020-12-22 | 2021-04-02 | 武汉大学 | 基于oco-2数据和wrf-stilt模型反演co2背景场浓度的方法及系统 |
CN113406667A (zh) * | 2021-07-07 | 2021-09-17 | 武汉大学 | 基于星载激光雷达获取城市co2格网化通量的方法及系统 |
CN114547553A (zh) * | 2022-04-27 | 2022-05-27 | 河北先河环保科技股份有限公司 | 二氧化碳排放量的反演方法、装置、设备及存储介质 |
Non-Patent Citations (3)
Title |
---|
A non-linear least squares enhanced POD-4DVar algorithm for data assimilation;Xiangjun Tian等;《Tellus A: Dynamic Meteorology and Oceanography》;20150126;第1-12页 * |
An atmospheric perspective on the carbon budgets of terrestrial ecosystems in China: progress and challenges;Baozhang Chen等;《Science Bulletin》;20210915;第17卷(第15期);第1713-1718页 * |
大气CO2浓度非均匀分布及其对地表升温影响的 研究进展与展望;张帆等;《地球信息科学》;20210831;第23卷(第8期);第1362-1371页 * |
Also Published As
Publication number | Publication date |
---|---|
CN114970184A (zh) | 2022-08-30 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN114970184B (zh) | 同步反演高分辨率人为co2排放与自然co2通量的同化方法及系统 | |
CN111487216B (zh) | 一种二氧化碳通量反演方法、系统 | |
CN106372730B (zh) | 利用机器学习的植被净初级生产力遥感估算方法 | |
CN110046415A (zh) | 一种时空精细化的土壤有机质含量遥感动态反演方法 | |
CN112699959B (zh) | 基于能量泛函模型的多源多尺度降水数据融合方法和装置 | |
CN115186437B (zh) | 碳同位素联合同化模型及甄别人为碳排放与自然碳通量区域同化系统的构建方法 | |
Zhang et al. | Developing a process-based and remote sensing driven crop yield model for maize (PRYM–Maize) and its validation over the Northeast China Plain | |
CN111859054A (zh) | 气象卫星数据的处理方法及装置 | |
Zhang et al. | Improving global gross primary productivity estimation by fusing multi-source data products | |
Lu et al. | Global prediction of gross primary productivity under future climate change | |
Chen et al. | Variation of gross primary productivity dominated by leaf area index in significantly greening area | |
SAITO et al. | Overview of model systems for global carbon dioxide and methane flux estimates using GOSAT and GOSAT-2 observations | |
CN115876948B (zh) | 基于卫星柱浓度和4d-letkf混合同化算法的碳卫星同化系统及其构建方法 | |
Dokoohaki et al. | A novel model–data fusion approach to terrestrial carbon cycle reanalysis across the contiguous US using SIPNET and PEcAn state data assimilation system v. 1.7. 2 | |
Bayarsaikhan et al. | Variations of vegetation net primary productivity and its responses to climate change from 1982 to 2015 in Mongolia | |
Reyes-Muñoz et al. | Inferring global terrestrial carbon fluxes from the synergy of Sentinel 3 & 5P with Gaussian process hybrid models | |
Liu et al. | Study of remote sensing based parameter uncertainty in production efficiency models | |
Gilbert et al. | Machine-Learning Approaches for Assessing Aerosol Optical Depth (AOD) in Ghana, West Africa | |
CN118053519B (zh) | 耦合ResNets和改进CASA模型的区域植被GPP估算方法 | |
Jin et al. | An Efficient Algorithm for Retrieving CO2 in the Atmosphere From Hyperspectral Measurements of Satellites: Application of NLS-4DVar Data Assimilation Method | |
Hu et al. | An Interpolation and Prediction Algorithm for XCO2 Based on Multi-Source Time Series Data | |
Thompson et al. | Top-down approaches | |
CN117496369A (zh) | 一种矿区碳汇数据集构建方法 | |
Fang et al. | Research Article Long-Term Variation of Global GEOV2 and MODIS Leaf Area Index (LAI) and Their Uncertainties: An Insight into the Product Stabilities | |
Bahrami et al. | Developing a Parsimonious Canopy Model (PCM v1. 0) to Predict Forest Gross Primary Productivity and Leaf Area Index |
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 | ||
CB03 | Change of inventor or designer information | ||
CB03 | Change of inventor or designer information |
Inventor after: Chen Baozhang Inventor after: Zhang Huifang Inventor before: Chen Baozhang |
|
GR01 | Patent grant | ||
GR01 | Patent grant |