CN116205022A - 基于co2/co比率反演核算高分辨率人为co2排放的方法及系统 - Google Patents
基于co2/co比率反演核算高分辨率人为co2排放的方法及系统 Download PDFInfo
- Publication number
- CN116205022A CN116205022A CN202210640366.XA CN202210640366A CN116205022A CN 116205022 A CN116205022 A CN 116205022A CN 202210640366 A CN202210640366 A CN 202210640366A CN 116205022 A CN116205022 A CN 116205022A
- Authority
- CN
- China
- Prior art keywords
- flux
- inversion
- assimilation
- area
- data
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 313
- 238000007599 discharging Methods 0.000 title claims description 5
- 230000004907 flux Effects 0.000 claims abstract description 742
- 238000011160 research Methods 0.000 claims abstract description 68
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 38
- 230000005540 biological transmission Effects 0.000 claims abstract description 32
- 230000008878 coupling Effects 0.000 claims abstract description 9
- 238000010168 coupling process Methods 0.000 claims abstract description 9
- 238000005859 coupling reaction Methods 0.000 claims abstract description 9
- 238000013022 venting Methods 0.000 claims abstract description 6
- 229910002091 carbon monoxide Inorganic materials 0.000 claims description 1088
- 238000004088 simulation Methods 0.000 claims description 281
- 230000008569 process Effects 0.000 claims description 82
- 239000000126 substance Substances 0.000 claims description 56
- 238000009826 distribution Methods 0.000 claims description 41
- 230000006870 function Effects 0.000 claims description 38
- 238000004458 analytical method Methods 0.000 claims description 28
- 238000006303 photolysis reaction Methods 0.000 claims description 8
- 238000000342 Monte Carlo simulation Methods 0.000 claims description 5
- CBENFWSGALASAD-UHFFFAOYSA-N Ozone Chemical compound [O-][O+]=O CBENFWSGALASAD-UHFFFAOYSA-N 0.000 claims description 4
- 239000000926 atmospheric chemistry Substances 0.000 claims description 4
- 239000000809 air pollutant Substances 0.000 claims description 3
- 231100001243 air pollutant Toxicity 0.000 claims description 3
- 239000003054 catalyst Substances 0.000 claims description 3
- 238000003860 storage Methods 0.000 claims description 3
- 231100000331 toxic Toxicity 0.000 claims description 3
- 230000002588 toxic effect Effects 0.000 claims description 3
- UGFAIRIUMAVXCW-UHFFFAOYSA-N Carbon monoxide Chemical compound [O+]#[C-] UGFAIRIUMAVXCW-UHFFFAOYSA-N 0.000 description 1086
- 208000037516 chromosome inversion disease Diseases 0.000 description 342
- 229910052799 carbon Inorganic materials 0.000 description 175
- 238000002474 experimental method Methods 0.000 description 174
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 170
- 239000010410 layer Substances 0.000 description 69
- 230000002829 reductive effect Effects 0.000 description 64
- 230000008859 change Effects 0.000 description 60
- 230000007246 mechanism Effects 0.000 description 46
- 230000000694 effects Effects 0.000 description 43
- 241000282414 Homo sapiens Species 0.000 description 36
- 238000005457 optimization Methods 0.000 description 26
- 239000011159 matrix material Substances 0.000 description 25
- 238000012544 monitoring process Methods 0.000 description 25
- 239000013598 vector Substances 0.000 description 25
- 239000007789 gas Substances 0.000 description 22
- 239000000203 mixture Substances 0.000 description 22
- 238000006243 chemical reaction Methods 0.000 description 18
- 238000012545 processing Methods 0.000 description 18
- 230000035945 sensitivity Effects 0.000 description 18
- 238000010206 sensitivity analysis Methods 0.000 description 15
- 238000012795 verification Methods 0.000 description 15
- 238000011161 development Methods 0.000 description 14
- 239000002803 fossil fuel Substances 0.000 description 14
- 238000012360 testing method Methods 0.000 description 14
- 230000003321 amplification Effects 0.000 description 13
- 238000003199 nucleic acid amplification method Methods 0.000 description 13
- 230000001932 seasonal effect Effects 0.000 description 13
- 238000000354 decomposition reaction Methods 0.000 description 12
- 238000004519 manufacturing process Methods 0.000 description 12
- 230000009467 reduction Effects 0.000 description 12
- 238000005516 engineering process Methods 0.000 description 11
- 230000007613 environmental effect Effects 0.000 description 11
- 239000000779 smoke Substances 0.000 description 11
- 238000004364 calculation method Methods 0.000 description 10
- 239000005431 greenhouse gas Substances 0.000 description 10
- 238000011835 investigation Methods 0.000 description 10
- 230000006872 improvement Effects 0.000 description 9
- 238000012937 correction Methods 0.000 description 8
- 230000000875 corresponding effect Effects 0.000 description 8
- 238000002156 mixing Methods 0.000 description 8
- 235000011962 puddings Nutrition 0.000 description 8
- PEDCQBHIVMGVHV-UHFFFAOYSA-N Glycerine Chemical compound OCC(O)CO PEDCQBHIVMGVHV-UHFFFAOYSA-N 0.000 description 7
- 230000002159 abnormal effect Effects 0.000 description 7
- MWUXSHHQAYIFBG-UHFFFAOYSA-N Nitric oxide Chemical compound O=[N] MWUXSHHQAYIFBG-UHFFFAOYSA-N 0.000 description 6
- RAHZWNYVWXNFOC-UHFFFAOYSA-N Sulphur dioxide Chemical compound O=S=O RAHZWNYVWXNFOC-UHFFFAOYSA-N 0.000 description 6
- 238000010276 construction Methods 0.000 description 6
- 230000007423 decrease Effects 0.000 description 6
- 239000003344 environmental pollutant Substances 0.000 description 6
- 230000029553 photosynthesis Effects 0.000 description 6
- 238000010672 photosynthesis Methods 0.000 description 6
- 231100000719 pollutant Toxicity 0.000 description 6
- 238000002485 combustion reaction Methods 0.000 description 5
- 239000002356 single layer Substances 0.000 description 5
- 241000894007 species Species 0.000 description 5
- QGZKDVFQNNGYKY-UHFFFAOYSA-N Ammonia Chemical compound N QGZKDVFQNNGYKY-UHFFFAOYSA-N 0.000 description 4
- 238000010521 absorption reaction Methods 0.000 description 4
- 238000001311 chemical methods and process Methods 0.000 description 4
- 239000003245 coal Substances 0.000 description 4
- 238000011156 evaluation Methods 0.000 description 4
- 238000000605 extraction Methods 0.000 description 4
- 238000010438 heat treatment Methods 0.000 description 4
- 230000010355 oscillation Effects 0.000 description 4
- 241000196324 Embryophyta Species 0.000 description 3
- 239000000443 aerosol Substances 0.000 description 3
- 238000013459 approach Methods 0.000 description 3
- 239000000356 contaminant Substances 0.000 description 3
- 230000003247 decreasing effect Effects 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000004927 fusion Effects 0.000 description 3
- 230000003993 interaction Effects 0.000 description 3
- 230000002452 interceptive effect Effects 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 230000008520 organization Effects 0.000 description 3
- 230000002093 peripheral effect Effects 0.000 description 3
- 238000012827 research and development Methods 0.000 description 3
- 238000012216 screening Methods 0.000 description 3
- 238000012876 topography Methods 0.000 description 3
- 238000012546 transfer Methods 0.000 description 3
- 229910018503 SF6 Inorganic materials 0.000 description 2
- 235000005775 Setaria Nutrition 0.000 description 2
- 241000232088 Setaria <nematode> Species 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 2
- 230000009471 action Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 239000004568 cement Substances 0.000 description 2
- 239000012141 concentrate Substances 0.000 description 2
- 238000005520 cutting process Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 230000004069 differentiation Effects 0.000 description 2
- 238000009792 diffusion process Methods 0.000 description 2
- 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 2
- 235000019534 high fructose corn syrup Nutrition 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 238000007689 inspection Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000012423 maintenance Methods 0.000 description 2
- 238000013508 migration Methods 0.000 description 2
- 230000005012 migration Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 229910000069 nitrogen hydride Inorganic materials 0.000 description 2
- 230000036961 partial effect Effects 0.000 description 2
- 238000005192 partition Methods 0.000 description 2
- 239000003208 petroleum Substances 0.000 description 2
- 238000001556 precipitation Methods 0.000 description 2
- 238000002360 preparation method Methods 0.000 description 2
- 238000003908 quality control method Methods 0.000 description 2
- 239000011435 rock Substances 0.000 description 2
- 238000004062 sedimentation Methods 0.000 description 2
- SFZCNBIFKDRMGX-UHFFFAOYSA-N sulfur hexafluoride Chemical compound FS(F)(F)(F)(F)F SFZCNBIFKDRMGX-UHFFFAOYSA-N 0.000 description 2
- 238000010200 validation analysis Methods 0.000 description 2
- 239000012855 volatile organic compound Substances 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- OHMHBGPWCHTMQE-UHFFFAOYSA-N 2,2-dichloro-1,1,1-trifluoroethane Chemical compound FC(F)(F)C(Cl)Cl OHMHBGPWCHTMQE-UHFFFAOYSA-N 0.000 description 1
- 239000002028 Biomass Substances 0.000 description 1
- 101100510615 Caenorhabditis elegans lag-2 gene Proteins 0.000 description 1
- 206010011906 Death Diseases 0.000 description 1
- RPNUMPOLZDHAAY-UHFFFAOYSA-N Diethylenetriamine Chemical compound NCCNCCN RPNUMPOLZDHAAY-UHFFFAOYSA-N 0.000 description 1
- 229910006715 Li—O Inorganic materials 0.000 description 1
- 241001465754 Metazoa Species 0.000 description 1
- 229910002089 NOx Inorganic materials 0.000 description 1
- 241000233855 Orchidaceae Species 0.000 description 1
- 240000006394 Sorghum bicolor Species 0.000 description 1
- 235000015505 Sorghum bicolor subsp. bicolor Nutrition 0.000 description 1
- 235000011684 Sorghum saccharatum Nutrition 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 229910021529 ammonia Inorganic materials 0.000 description 1
- 238000013398 bayesian method Methods 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 150000001721 carbon Chemical class 0.000 description 1
- 238000004177 carbon cycle Methods 0.000 description 1
- KYKAJFCTULSVSH-UHFFFAOYSA-N chloro(fluoro)methane Chemical compound F[C]Cl KYKAJFCTULSVSH-UHFFFAOYSA-N 0.000 description 1
- 238000000658 coextraction Methods 0.000 description 1
- 239000000470 constituent Substances 0.000 description 1
- 238000011109 contamination Methods 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 238000002790 cross-validation Methods 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 238000013480 data collection Methods 0.000 description 1
- 238000013499 data model Methods 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 238000000151 deposition Methods 0.000 description 1
- 230000008021 deposition Effects 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000010304 firing Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 238000007716 flux method Methods 0.000 description 1
- 239000000446 fuel Substances 0.000 description 1
- 239000002440 industrial waste Substances 0.000 description 1
- 239000011261 inert gas Substances 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000012067 mathematical method Methods 0.000 description 1
- 239000012567 medical material Substances 0.000 description 1
- 238000006386 neutralization reaction Methods 0.000 description 1
- 238000007254 oxidation reaction Methods 0.000 description 1
- 239000013618 particulate matter Substances 0.000 description 1
- 230000015843 photosynthesis, light reaction Effects 0.000 description 1
- 150000003254 radicals Chemical class 0.000 description 1
- 230000002285 radioactive effect Effects 0.000 description 1
- 239000000700 radioactive tracer Substances 0.000 description 1
- 230000000241 respiratory effect Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 238000004158 soil respiration Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 229960000909 sulfur hexafluoride Drugs 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- VZGDMQKNWNREIO-UHFFFAOYSA-N tetrachloromethane Chemical compound ClC(Cl)(Cl)Cl VZGDMQKNWNREIO-UHFFFAOYSA-N 0.000 description 1
- TXEYQDLBPFQVAA-UHFFFAOYSA-N tetrafluoromethane Chemical compound FC(F)(F)F TXEYQDLBPFQVAA-UHFFFAOYSA-N 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 239000005436 troposphere Substances 0.000 description 1
- 238000013076 uncertainty analysis Methods 0.000 description 1
- 238000004800 variational method Methods 0.000 description 1
- 239000002912 waste gas Substances 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
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/29—Geographical information databases
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- 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)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Databases & Information Systems (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Remote Sensing (AREA)
- Data Mining & Analysis (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供了一种基于CO2/CO比率反演核算高分辨率人为CO2排放的方法及系统。该方法包括：采用FORTRAN语言结合POD4DVar同化算法，实现正交4维变分POD4DVar同化算法与区域大气传输模型CMAQ的耦合，研发了区域高分辨率同化反演系统TracersTracker；对上甸子区域和瓦里关区域的CO通量进行反演和优化，获取两个研究区域内的最优CO通量；基于排放清单中的CO2/CO排放比及获取的最优CO通量，反演得到上甸子区域和瓦里关区域的人为CO2排放通量；基于地表大气中的ΔCO2/ΔCO浓度比值观测数据及CO地基观测数据，利用本专利研发的TracersTracker系统反演得到上甸子区域和瓦里关区域的高分辨率人为CO2排放通量。本发明提供的基于CO2/CO比率反演核算高分辨率人为CO2排放的方法及系统能够解决大气反演系统无法区分人为CO2排放和自然CO2通量的技术难题。
Description
技术领域
本发明涉及大气科学、全球变化科学及地理科学技术领域，特别是涉及一种基于CO2/CO比率反演核算高分辨率人为CO2排放的方法及系统。
背景技术
碳是自然界中动植物的重要组成元素，它以各种形态存在于大气圈、水圈、生物圈和岩石圈中。碳元素在各个圈层中不停的迁移、转换的过程称为碳循环，碳循环过程涉及面广且过程较为复杂。在自然条件下的碳循环过程中，各个圈层中的碳含量相对稳定，主要受一些剧烈的自然活动影响，如火山喷发、地壳运动等。然而，自工业革命以来，人类活动向大气中排放了大量的CO2，打破了各个圈层中碳含量的基本平衡，尤其以大气圈为甚。联合国政府间气候变化专门委员会(Intergovernmental Panel on Climate Change，IPCC)第五次评估报告指出，人类工业革命进程明显提升了大气中CO2、CH4的水平，CO2浓度是工业革命前的1.4倍，CH4浓度是工业革命前的2.5倍。化石燃料(如煤、石油等)的使用及土地利用变化(如森林砍伐等)是CO2浓度升高的主要原因，因此产生的CO2总量已超过5000亿吨。其中，化石燃料的使用和水泥生产等人为活动排放的CO2占比超过七成(～3350-3950亿吨)，土地利用变化产生的CO2约占三成。工业革命以来，大气中CO2浓度逐年攀升，已从1750年的278ppm上升到2018年的411.2ppm。科学、合理的减少温室气体CO2的排放，首先需要明确全球及区域碳源汇的时空分布特征，进而制定有效的减排计划。虽然碳元素是在各个圈层中迁移转换的，但由于大气圈与人类生产、生活的关系更为紧密，所以碳循环及碳源汇的研究通常以大气圈为参照，并依据各圈层与大气圈碳交互状态描述全球及区域的“碳源”和“碳汇”。“碳源”与“碳汇”的概念源于《京都议定书》，碳源(Carbon Source)是指向大气中释放碳的过程、活动或机制，人类的很多生产活动会向大气中排放CO2，如化石燃料燃烧、水泥生产、森林砍伐等。碳汇(Carbon Sink)是吸收、消耗大气中的CO2，从而使大气中CO2浓度降低的过程、活动或机制，如植物光合作用等。人类活动向大气排放的大量CO2，是地球上重要的“碳源”，海洋和陆地生态系统可以在特定条件下消耗大气中的CO2，起到气候变化缓冲器的作用，是明显的“碳汇”。当前，各国政府及科研机构开展了一系列全球和大洲尺度的碳源汇估算研究，但各个成果之间存在很多不确定性。另外，当前碳源汇估算研究多集中在较大尺度上，较高时空分辨率的区域碳源汇研究成果仍相对匮乏。
区域高分辨率碳源汇研究有助于了解城市化进程在气候变化中的作用。城市是人类生产、生活活动最为集中的区域，其碳排放量约占总排放量的85％左右。但是，当前碳源汇估算研究的成果主要集中在全球和洲际尺度，其时空分辨率较低，难以满足区域尺度应用的需求。区域及城市尺度的碳源汇研究，可以获取更高空间和时间分辨率的碳源汇产品，为更大范围、更大大尺度的碳源汇估算研究提供参考和有力数据支持。在此基础上开展人为碳排放的区分研究，可以精确量化陆地生态系统碳源汇和人为碳源汇，为进一步理解陆地生态系统、城市间相互作用机制提供数据基础，也使我国的城市温室气体减排工作更加有的放矢。
碳源汇估算方法总体上可以分为两类：“自下而上”法和“自上而下”法。“自下而上”法以陆地为研究对象，主要通过调查统计、通量观测、模型模拟等方法对碳源汇进行估算。“自上而下”法则以大气为研究对象，通过大气CO2浓度观测数据、反演算法及大气传输模型组成的反演系统对碳源汇进行估算。
碳源汇估算方法总体上可以分为两类：“自下而上”法和“自上而下”法。“自下而上”法以陆地为研究对象，主要通过调查统计、通量观测、模型模拟等方法对碳源汇进行估算。“自上而下”法则以大气为研究对象，通过大气CO2浓度观测数据、反演算法及大气传输模型组成的反演系统对碳源汇进行估算。
同化反演方法属于“自上而下”的碳源汇估算方法，它利用大气CO2地基及卫星观测数据，结合大气传输模型和数据同化算法，反演得到不同尺度的碳源汇时空分布情况。同化反演方法可以有效规避传统调查统计方法估算碳源汇的弊端，并在很大程度上减少由于参数设置等人为因素引入的不确定性。虽然同化反演方法发展时间不长，但其在碳源汇估算研究中的应用越来越多。在碳同化反演系统中，同化技术或同化算法是核心，其发展经历了若干阶段，最早可以追溯到早期的气象数值模拟。早期的数值模拟方法较为简单且易于实施，如逐步订正法、最优插值法。逐步订正法是一种基于统计估计理论的数据同化方法，这使得数据同化方法有了理论基础，另外逐步订正法引入了背景场的概念，通过确定影响半径来逐个更新网格点的背景值，避免了异常数据震荡的问题。但其在更新网格点背景值时的权重是人为设置的，这在很大程度上依赖操作人员的专业素养及经验。最优插值法也是一种基于统计估计理论的数据同化方法，与逐步订正法不同，最优插值法基于观测资料和背景场提供的先验信息，通过最小方差估计确定最优权重，对背景场进行订正，得到统计意义上的最优解。该方法明显提高了预报场的精度，但是该方法要求观测变量与分析变量之间必须满足线性关系，这大大降低了该方法对多源观测数据(如遥感观测)直接同化的能力。
在碳同化反演研究中，研究人员发展了很多碳源汇模拟模型及同化反演系统。早期的碳源汇研究多采用一维或二维模型模拟碳源汇的时空变化，随着人们对大气物理、化学过程理解的深入，越来越多的三维模型开始投入到实际应用中。Enting等利用GISS(Goddard Institute for Space Studies)大气示踪传输模型估算了1986-1987年、1989-1990年两个阶段的陆地生态系统和大气间的碳通量交换情况。Fan等根据CO2观测数据，利用三维模型GFDL(Geophysical Fluid Dynamics Laboratory)和GCTM(Global ChemicalTransport Model)，同时结合两个海-气先验通量，对陆地生态系统的碳通量进行了优化。Kaminski等在三维模型TM2支持下，结合贝叶斯反演方法和美国国家海洋和大气管理局(National Oceanic and Atmospheric Administration，NOAA)提供的CO2观测数据，对全球的年际碳排放量进行了反演，同时分析了反演的碳通量的季节变化特征。不同模型在模拟大气传输时采用的模拟方法不同，在模拟时存在一定的不确定性，针对此问题，Gurney等提出了大气示踪输送模型比较计划Transcom1～Transcom3，该计划将全球分为23个区，利用若干大气传输模型对所有分区进行碳同化反演并比较了这些模型在模拟时的误差分布情况。该计划反演过程的时间步长为年或月，同时反演过程中分区也较为粗糙，使得反演出来的通量数据时空分辨率都较低。因此，许多研究者在此系统基础上将分区细分并引入嵌套模式，提高了碳源汇反演的时空分辨率。为了更系统、全面的理解全球碳循环过程及碳源汇时空演变机制，美国国家海洋与大气局于2007年研发了全球碳追踪同化系统(CarbonTracker，CT)，该系统以北美为重点研究区域，在世界不同区域提供不同空间分辨率的同化反演结果，其中北美地区分辨率为1°×1°，其他地区分辨率为3°×2°。该系统可以同化地面、飞机、卫星等多源观测数据，已被广泛用于陆地-大气间碳通量的估算研究。相对于Transcom系统，CarbonTracker系统在很多方面做了改进：(1)同化算法的改进。CarbonTracker系统充分考虑了CO2滞后窗口的影响，提出了滞后卡尔曼滤波同化方法，实现了CO2现实浓度与通量之间的动力学关联，使CO2物理、化学过程尽可能的接近真实状态，在很大程度上提高了同化反演精度。(2)支持多源观测数据同化。在Transcom系统中，参与同化的观测数据非常有限且标准不一，这使得同化反演的碳通量存在较大不确定性。CarbonTracker系统对同化反演框架进行了改进，可支持地基、卫星、航测等多源观测数据的同化，降低通量反演中的不确定性。(3)提高反演通量的时空分辨率。在Transcom系统的反演过程中，时空分辨率较低，只粗略的将全球分为23个区，时间分辨率为年或月。CarbonTracker在Transcom系统基础上进行了细化，将全球分为240个区进行CO2通量反演，另外CarbonTracker系统将时间步长设置为天或小时，时间分辨率和空间分辨率均有较大提升。
目前，世界许多国家和地区对Carbon Tracker系统进行了移植和改进，建立了本地化的碳同化反演系统，如CT-Australasia、CT-Europe、CT-China等。同时，基于CarbonTracker及其改进的碳同化反演系统，各国科研人员进行了很多碳源汇的估算研究。Ingrid等改进了CT-Europe系统的核心模块，并用新改进的碳同化反演系统估算了全球的碳通量，揭示了陆地生态系统和海洋生态系统的碳汇分布情况。Laan-Luijkx等利用南美碳追踪系统CarbonTracker South America分析了南美洲亚马逊地区的碳源汇分布情况，指出该地区碳汇在2010年出现了明显的减少。Kim等利用CarbonTracker系统，采用三层嵌套模式，重点估算了亚洲地区的碳源汇，结果表明CarbonTracker系统表现良好。Lim等CarbonTracker系统模拟了韩国首尔都市区的碳源汇变化情况，分析了CO2日间的具体变化规律。
我国在碳同化反演领域的研究起步较晚，但近年来取得了长足的进展。国内较早的碳同化反演系统是南京大学和北京师范大学合作研发的GCAS(Global CarbonAssimilation System)系统，该系统基于GEOS-Chem(Goddard Earth Observing System-Atmospheric Chemistry Model)大气传输模式，采用四维变分同化方法同化观测数据。相较于欧美的全球碳同化系统，该系统对先验碳通量进行了优化，采用了多个陆地生态系统模型的通量估计值，同时还将大气CO2浓度作为一个状态变量进行了优化。与此同时，陈报章等与CarbonTracker系统主要研发者——荷兰Wageningen University大学的WouterPeters教授合作，将CarbonTracker系统移植到中国区域，研发了中国碳追踪同化系统(CarbonTraker-China,CT-China)。CT-China的先验排放清单包括三个部分：人为排放清单、生物质燃烧排放清单和海洋碳源碳汇清单，在先验排放清单基础上，以网格化气象数据为驱动场，同化地表CO2浓度观测数据，反演得到全球网格化的陆地生态系统碳源碳汇时空分布成果。CT-China可以提供中国区域1°×1°，中国以外区域2°×3°分辨率的碳通量数据。该系统得到了NOAA和欧洲CarbonTracker研发团队官方承认，是继北美碳追踪系统(CarbonTracker-North America)和欧洲碳追踪系统(CarbonTracker-Europe)之后的全球第三套完整的全球尺度嵌套式碳同化反演系统。此外，一些研究人员利用CarbonTracker及其改进系统，对中国境内的碳源汇通量结果进行了验证。安兴琴等利用分布在中国的4个大气本底站的观测数据对CarbonTracker的反演结果进行了验证，结果表明CarbonTracker在中国区域具有较好的适应性，系统可用于中国区域碳通量的反演研究。麦博儒等利用瓦里关本底站和番禺观测站的观测数据对改进的CT-2010系统进行了验证，结果显示，CT-2010能较好地反映了近地层CO2浓度的分布状况。李灿等在一个二维大气化学传输模式基础上，研发建立了用于CO2源汇反演的反演模型，利用该模型获取了全球碳源汇的时空分布情况。Zhang等利用CT-China系统，同化反演了中国区域内的碳源汇分布状况，指出中国陆地生态系统是明显的碳汇，在2001-2010十年间吸收了大约0.33Pg C/yr。Jiang等利用贝叶斯方法，采取嵌套模式，在Transcom中全球分区的基础上，对2002-2008年间中国区域和全球其他区域的碳源汇进行了估算，其结果表明虽然全球陆地生态系统碳汇呈增长趋势，中国区域陆地生态系统大部分时间为碳汇。在此基础上，Jiang等利用大气CO2同化反演模型估算出2006-2009年中国陆地碳汇为0.33～0.35Pg C/yr。
经过几十年的发展，同化反演方法在碳源汇估算方面取得了众多成果，且新的反演方法和观测数据也在不断地加入到碳源汇估算研究中来，极大地促进了碳源汇的研究进展，但是当前的碳源汇同化反演研究还存在一些不足，主要表现为：
(1)当前碳源汇同化反演研究多集中在全球或大洲尺度上，同化反演结果分辨率较低，无法满足区域碳源汇研究的精度要求；(2)当前碳源汇估算研究中缺少对人为排放的定量区分，碳源汇估算研究需进一步深入和细化。
近年来，针对上述问题，在对碳排放总量进行同化反演的基础上，研究人员试图定量区分人为碳排放。人为碳排放对大气中CO2浓度起着至关重要的影响，即使在较为偏远的农村和郊区，化石燃料的燃烧产生CO2仍很大程度上影响着大气中CO2的浓度，将人为CO2排放和自然CO2排放进行定量区分，有利于进一步理解碳源汇时空分布格局及其演变机制。区分人为CO2排放的主要手段是借助CO2伴生气体来实现的，例如14C、SF6、CO等。在很多相关研究中，利用以14C为代表的放射性元素对人为排放的CO2进行区分被证明是一种有效的手段。由于化石燃料在地下埋藏了几千万年甚至几亿年，化石燃料中14C的含量已经非常少，所以化石燃料的燃烧等活动几乎不排放14C，研究人员通过分析大气中14CO2的含量占比来确定人为CO2排放量。但是，14CO2观测数据的获取流程复杂且耗费昂贵，从而导致观测数据无法满足持续性研究的需要。因此，研究人员视图利用其他CO2伴生气体进行人为碳排放的定量区分研究，如氯氟碳化物(CFCS)、氢代氯氟碳化物(HFCS)、全氟化碳(PFCS)、六氟化硫(SF6)、一氧化碳(CO)等。借助于CO2伴生气体定量区分人为排放，伴生气体一般需要具有以下几个特性：(1)与CO2相关性高；(2)观测数据易获取；(3)化学机理易于理解。在使用这些CO2伴随排放气体时，应充分了解其物理、化学性质及演变规律，否则在使用过程中会带来更大的不确定性。虽然CFCS、HFCS、PFCS、SF6与CO2排放具有良好相关性，但其观测数据在时空分辨率上难以满足要求，另外这些物质的化学机理较为复杂，在模拟过程中存在很大困难。相较于这些物质，CO与CO2排放的相关性更加紧密，并且在国家环境监测数据库中有较为充足、精度较高的观测数据。但是CO是一种化学性质较为活跃的物质，在大气传输过程中易于与其他物质发生反应，在实际应用时应充分考虑其在空间和时间上的演变规律。
总之，目前基于大气反演模型达到方法估算人为碳排放的结果仍存在着很大的不确定性，而且不能区分自然CO2通量和人为CO2排放，这是评估碳中和行动有效性急需解决的技术问题。
发明内容
本发明要解决的技术问题是提供一种基于CO2/CO比率反演核算高分辨率人为CO2排放的方法及系统，能够解决大气反演系统无法区分人为CO2排放和自然CO2通量的技术难题。
为解决上述技术问题，本发明提供了一种基于CO2/CO比率反演核算高分辨率人为CO2排放的方法，所述方法包括：采用FORTRAN语言和正交4维变分POD4DVar同化算法，实现POD4DVar同化算法与区域大气传输模型CMAQ的耦合，研发了区域高分辨率同化反演系统TracersTracker；在TracersTracker系统的支持下，对上甸子区域和瓦里关区域的CO通量进行反演和优化，获取两个研究区域内的最优CO通量；基于排放清单中的CO2/CO排放比及获取的最优CO通量，反演得到上甸子区域和瓦里关区域的人为CO2排放通量；基于地表大气中的ΔCO2/ΔCO浓度比及CO地基观测数据，利用TracersTracker系统同化CO2观测数据，反演得到上甸子区域和瓦里关区域的高分辨率人为CO2排放的通量；结合自然CO2排放通量，在CO2总通量内区分人为CO2通量，并在此基础上分析两种方法获取的人为CO2排放通量的核算精度及其不确定性。
在一些实施方式中，CMAQ是一个三维欧拉大气化学和传输模拟系统，可以模拟颗粒物、臭氧、有毒空气污染物的80多种物质。
在一些实施方式中，CMAQ包括：化学传输模块、气象化学接口模块、边界场模块、初始场模块、光分解模块。
在一些实施方式中，本专利所研发的TracersTracker CO2同化系统涉及到的第三方软件或函数库包括：IOAPI函数库、NetCDF函数库、Matlab函数库。
在一些实施方式中，TracersTracker，主程序采用FORTRAN90语言编写，编译器为Intel Fortran Complier。
在一些实施方式中，POD4Dvar同化算法采用蒙特卡罗法生成样本集合来表征通量的空间分布格局及时间演变特征。
在一些实施方式中，POD4Dvar同化算法无需积分伴随模式。
在一些实施方式中，同化过程包含若干次同化循环。
在一些实施方式中，采用均方根误差、线性相关系数、绝对误差均值、误差均值对同化系统反演结果的准确性进行评价。
此外，本发明还提供了一种基于CO2/CO比率反演核算高分辨率人为CO2排放的系统，所述系统包括：一个或多个处理器；存储装置，用于存储一个或多个程序；当所述一个或多个程序被所述一个或多个处理器执行，使得所述一个或多个处理器实现根据前文所述的基于CO2/CO比率反演核算高分辨率人为CO2排放的方法。
采用这样的设计后，本发明至少具有以下优点：
针对当前碳同化反演系统分辨率较低的问题，本发明实现了POD4DVar高效同化算法和区域大气传输模型CMAQ的耦合，构建了区域高分辨率碳同化反演系统TracersTracker。针对CO2观测数据匮乏导致区域高分辨率通量无法反演的问题，提出了基于CO2/CO排放比的CO2通量反演方法和基于CO2/CO浓度比的CO2通量反演方法。以地理和气候条件不同、经济发展水平各异的上甸子和瓦里关为研究区域，开展了CO和CO2通量的反演优化研究，对CO2人为通量和自然通量进行了定量区分。主要结论如下：
(1)将本征正交分解技术嵌入传统四维变分同化方法中，开发了POD4DVar同化算法。相较于传统四维变分同化方法，POD4DVar同化算法避免了伴随模式的开发，操作和维护都较为简便，极大的提高了同化反演过程的效率。在此基础上，实现了POD4DVar算法和区域大气传输模型CMAQ的耦合，构建了区域高分辨率碳同化反演系统TracersTracker，该系统可充分吸收观测数据中的信息，有效校正和优化先验通量中的“误差”。基于地基实时观测数据反演CO通量过程中，TracersTracker系统消除了25％的先验模拟误差，与CO观测数据的相关系数均值由0.64提升至0.75。
(2)提出了基于CO2/CO排放比(排放端法)的CO2通量反演方法和基于CO2/CO浓度比(大气端法)的CO2通量反演方法，在研究区域内开展了CO2通量的反演优化研究。排放端法获取的CO2通量极大地校正了先验通量中的“误差”，有效减小了后验模拟误差，明显提升了CO2后验模拟与实际观测的相关系数。与先验模拟相比，后验模拟绝对误差均值在两个区域分别减小了65％和62％，相关系数均值在上甸子区域由0.67提升至0.88，在瓦里关区域由0.69提升至0.89。大气端法获取的CO2通量在少数时段可以优化先验通量，但整体上表现并不稳定。与先验模拟相比，CO2后验模拟误差没有明显减小，在两个区域仅分别减小了2％和3％，相关系数甚至出现了下降，上甸子区域相关系数均值由0.67下降为0.66，瓦里关区域相关系数均值由0.69下降到0.68。大气端法反演CO2通量涉及的步骤较多，在反演同化过程中容易引入新的误差，其反演得到的CO2通量不确定性较大、精度难以有明显提高，在很多时段后验模拟精度甚至会出现明显降低。在当前地表观测分布密度条件下，该方法具有明显的不稳定性。相对于大气端法，排放端法步骤较少、操作简单，可以有效提高CO2通量精度，但要求采用的CO2/CO排放比较为稳定、可靠。总体而言，虽然基于CO2/CO排放比进行CO2通量反演会引入一定程度的不确定性，但在区域CO2观测数据较为匮乏的现状下，该方法仍然是获取区域较高分辨率CO2通量的有效方法。
(3)基于排放端法和大气端法获取的CO2通量，对研究区域内CO2人为通量和自然通量进行了定量区分。其中，上甸子区域是较明显的碳源，并且该区域碳源有逐年变大的趋势，基于CO2/CO排放比和浓度比反演的CO2后验总通量分别为36.5摩尔/秒和38.1摩尔/秒，人为、自然CO2排放比分别为4.9:1和5:1。瓦里关区域是较小的碳汇，并且该区域碳汇有逐年减小的趋势，基于CO2/CO排放比和浓度比反演的CO2后验总通量分别为-4.4摩尔/秒和-4.0摩尔/秒，人为、自然CO2排放比分别为1:2和1:1.8。
附图说明
上述仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，以下结合附图与具体实施方式对本发明作进一步的详细说明。
图1是总技术路线图；
图2是同化TracersTracker系统运行流程图；
图3是同化系统的同化循环示意图；
图4是一月份同化前后CO浓度模拟；
图5是七月份同化前后CO浓度模拟；
图6是一月份同化前后CO浓度模拟误差；
图7是七月份同化前后CO浓度模拟误差；
图8是一月份同化前后CO通量绝对误差(a为同化前，b为同化后)；
图9是七月份同化前后CO通量绝对误差(a为同化前，b为同化后)；
图10是一月份同化前后CO2浓度模拟；
图11是七月份同化前后CO2浓度模拟；
图12是一月份同化前后CO2浓度模拟误差；
图13是七月份同化前后CO2浓度模拟误差；
图14是一月份同化前后CO2通量绝对误差(a为同化前，b为同化后)；
图15是七月份同化前后CO2通量绝对误差(a为同化前，b为同化后)；
图16是主实验CO通量反演结果；
图17是主实验CO2通量反演结果；
图18是扰动样本实验CO通量反演结果；
图19是扰动样本实验CO2通量反演结果；
图20是边界场实验CO通量反演结果；
图21是滞后窗口实验CO通量反演结果；
图22是滞后窗口实验CO2通量反演结果；
图23是化学机制实验CO通量反演结果；
图24是化学机制实验CO2通量反演结果；
图25是空间分辨率实验CO通量反演结果；
图26是空间分辨率实验CO2通量反演结果；
图27是CO通量反演流程图；
图28是排放端法CO2通量反演流程图；
图29是大气端法CO2通量反演流程图；
图30是实施例区域及两层嵌套范围；
图31是上甸子和瓦里关区域内CO观测站点；
图32是上甸子区域CO人为先验通量；
图33是瓦里关区域CO人为先验通量；
图34是上甸子区域CO2人为先验通量；
图35是瓦里关区域CO2人为先验通量；
图36是上甸子区域CO2自然先验通量；
图37是瓦里关区域CO2自然先验通量；
图38是上甸子区域CO观测站点背景浓度；
图39是瓦里关区域CO观测站点背景浓度；
图40是上甸子区域CO后验通量增量；
图41是瓦里关区域CO后验通量增量；
图42是瓦里关区域先验模拟、后验模拟和观测；
图43是上甸子区域CO先验模拟、后验模拟和观测；
图44是北京定陵观测站模拟误差；
图45是保定监测站模拟误差；
图46是廊坊药材公司观测站模拟误差；
图47是廊坊药材公司观测站模拟误差；
图48是唐山十二中观测站模拟误差；
图49是西宁监测站模拟误差；
图50是临夏环保局观测站模拟误差；
图51是上甸子区域不同季节后验通量增比；
图52是瓦里关区域不同季节后验通量增比；
图53是上甸子区域不同季节后验通量增量；
图54是瓦里关区域不同季节后验通量增量；
图55是上甸子区域春季CO模拟误差；
图56是上甸子区域夏季CO模拟误差；
图57是上甸子区域秋季CO模拟误差；
图58是上甸子区域冬季CO模拟误差；
图59是瓦里关区域春季CO模拟误差；
图60是瓦里关区域夏季CO模拟误差；
图61是瓦里关区域秋季CO模拟误差；
图62是瓦里关区域冬季CO模拟误差；
图63是上甸子区域月均CO2/CO排放比；
图64是瓦里关区域月均CO2/CO排放比；
图65是上甸子区域月均ΔCO2/ΔCO浓度比；
图66是瓦里关区域月均ΔCO2/ΔCO浓度比；
图67是上甸子区域排放端法人为CO2后验通量增量；
图68是瓦里关区域排放端法人为CO2后验通量增量；
图69是上甸子区域大气端法人为CO2后验通量增量；
图70是瓦里关区域大气端法人为CO2后验通量增量；
图71是上甸子区域CO2人为通量和自然通量；
图72是瓦里关区域CO2人为通量和自然通量；
图73是CO2观测值和模拟值；
图74是CO2模拟误差分布；
图75是月均CO2/CO排放比相关系数；
图76是上甸子和瓦里关区域月均CO2/CO排放比；
图77是上甸子和瓦里关区域CO观测值。
具体实施方式
以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
技术方案的目标是依据CO2与CO是人类使用化石燃料排到的同源性的特征，研发基于CO2与CO排放比率和浓度比率反演人为CO2排放的方法，实现区域尺度刚分辨率的人为CO2排放的合理估算，解决大气反演系统无法区分人为CO2排放和自然CO2通量的技术难题。核心技术包括如下4个部分：
(1)原创性研发了区域高分辨率区域碳同化反演系统。
将本征正交分解技术(POD)与四维变分同化方法相结合，实现POD4Dvar高效同化算法，提高同化反演系统的运行效率。基于POD4DVar算法和区域大气传输模型CMAQ，研发区域高分辨率同化反演系统TracersTracker。对构建的同化反演系统进行有效性验证和敏感性分析，验证系统的通量反演效果。
(2)开展区域高分辨率CO通量反演研究。
原创性提出区域高分辨率CO通量反演方案，以上甸子和瓦里关为两个靶向研究区域，利用本发明构建的区域同化反演系统同化CO地基观测数据，获取两个研究区域内的高分辨率(3km)CO后验通量。
(3)原创性地提出CO-CO2联合同化的技术方案，研发了反演核算高分辨率人为CO2通量的方法。
利用CO和CO2的伴生相关性，提出两种CO-CO2联合同化方法对CO2通量进行反演，即基于大气端观测数据的反演方法(大气端法)和基于排放端通量数据的反演方法(排放端法)。两种方法分别利用排放端和观测端的CO和CO2相关性信息，实现上甸子区域和瓦里关区域的高分辨率CO2通量的反演和优化，将CO2后验通量的空间分辨率提升至3km。
(4)原创性研发了一种区分CO2人为通量与自然通量的技术。
在已获取的后验CO通量数据和后验CO2通量基础上，定量区分人为CO2排放和自然CO2通量，厘清区域内人为和自然碳排放之间的关系，为深入理解陆地生态系统、城市间相互作用机制提供数据基础。
本发明研发总体技术方案：
首先，使用FORTRAN语言开发了POD4DVar同化算法，在此基础上实现POD4DVar同化算法与区域大气传输模型CMAQ的耦合，研发了区域高分辨率同化反演系统TracersTracker。然后，利用TracersTracker系统同化CO地基观测数据，实现上甸子区域和瓦里关区域CO通量的反演优化。其次，从CO排放清单和浓度观测中获取CO2/CO排放比和ΔCO2/ΔCO浓度比，采用CO2-CO联合同化方法，实现区域高分辨率CO2通量的反演优化。最后，对CO2人为通量和自然通量进行定量区分。发明技术路线如图1所示。
(1)首先，基于FORTRAN语言开发POD4DVar同化算法，实现POD4DVar同化算法与区域大气传输模型CMAQ的耦合，研发区域高分辨率同化反演系统TracersTracker，该系统支持高分辨率CO2和CO通量的反演和优化。
(2)在TracersTracker系统的支持下，对上甸子区域和瓦里关区域的CO通量进行反演和优化，获取两个研究区域内的“最优”CO通量。
(3)基于排放清单中的CO2/CO排放比及(2)中获取的“最优”CO通量，反演得到上甸子区域和瓦里关区域的CO2通量。
(4)基于地表大气中的ΔCO2/ΔCO浓度比及CO地基观测数据，推算得到“CO2观测数据”，利用TracersTracker系统同化“CO2观测数据”，反演得到上甸子区域和瓦里关区域的高分辨率CO2通量。
(5)结合自然CO2排放通量，在CO2总通量内区分人为CO2通量，并在此基础上分析两种方法获取的CO2通量的模拟精度及其不确定性。
所涉及的四个方面的核心方法及关键技术叙述如下。
1同化反演系统构建方案及框架
1.1引言(Preface)
四维变分同化方法可以有效同化历史观测资料，对初始场进行校正和优化。但是，四维变分同化方法在实际应用中也存在很多困难，如切线性伴随模式开发非常困难、高维运算资源消耗大、运算效率低等。POD4DVar是基于本征正交分解的显示四维变分方法，它用较少的基向量捕捉数据的时空演变特征，使高维运算问题可以在维度较低的环境下得到解决，具有较高的运算效率。
本发明研究将POD4DVar高效同化算法和区域大气传输模型CMAQ相耦合，研发了区域碳同化反演系统TracersTracker，该系统可以同化CO、CO2观测数据，反演得到区域高分辨率CO、CO2后验通量。本章详细介绍了构建TracersTracker系统的具体方案和系统框架，并对TracersTracker系统的CO、CO2通量反演结果进行了验证。下面对本征正交分解方法、POD4DVar及TracersTracker同化反演系统的基本原理、构建方案等做详细阐述。
1.2相关理论
1.2.1四维变分同化
变分法是17世纪末发展起来的一门数学分支，主要用来求取泛型函数的极值，最典型的例子是伯努利提出的最速降线问题。变分同化是根据变分法的基本原理，通过同化历史观测数据使分析场达到统计意义上最优的一种数学方法。三维变分同化只能同化特定时刻的观测数据，并且同化结果在时间上也不连续，四维变分同化扩展了三维变分同化的时间维度，可以利用历史观测资料对初始场进行校正和优化。
在四维变分同化方法中，通常先定义代价函数来表示模型模拟和观测之间的“误差”，四维变分同化的任务是寻求该代价函数极小化时所对应的初始场，代价函数可表示如下：
J(x)＝JB+JO (3.1)
这里
JB＝(x-xb)TB-1(x-xb) (3.2)
其中M表示积分预报模型，T表示转置，b表示背景场，k代表tk时刻，w代表各类观测数据，O为观测算子，B和R分别代表背景误差协方差矩阵和观测误差协方差矩阵。
求取式(3.1)最小值本质上是非线性约束条件下的极值求解问题，其基本思路是按照Euler-Lagrange方程求解极值的方式进行，首先需求取式(3.1)的梯度
求取式(3.1)的梯度可分为两部分，即
即
假设A为对称矩阵，
为关于A的函数，则其梯度为
根据上式，计算JB的梯度，得
求取Jo的梯度则较为困难，因为其中含有复杂的模型运算，而模型一般是非线性的。
假设
其中
y＝f(x) (3.9)
则Z是嵌套函数，其梯度为
根据上式，JO部分也是嵌套函数，在求取t0时刻x的最优值时，其梯度为
其中
xi＝Mi[x(t0)] (3.12)
式(3.11)转置可得
则
根据上式，在求取代价函数极值的过程中，需开发模拟模型的伴随模式，以此寻找初始场的梯度来获取函数的极值，即向前积分模型一次，向后积分伴随模式到初始时刻调整初始场，然后再返回向前模式，循环这个过程，直到问题得到解决。但是对于复杂的方程(如CMAQ模型)，其伴随模式的求解是非常困难的。另外，基于伴随模式求解代价函数极值时，需要消耗大量的计算机资源，计算效率较低。
1.2.2本征正交分解
本征正交分解(Proper orthogonal decomposition，POD)在数理统计领域具有广泛的应用，这种方法本身蕴含着降维降阶的思想，可以用较低的维度表征高维空间的时空演变特征，其具体过程可表述如下。
设F表示模型模拟的物理、化学过程，
i＝1,2,…,M表示变量值在
点处模拟的M个状态值的集合，则该集合的平均值为：
则每个模拟值的偏差为
通过上式得到偏差的集合代表着模型模拟值集合的空间特征，假设能够找到该集合的简洁表达方式，那么高维问题可以得到更快速的解决。利用一组基函数及其对应系数可以简洁表达该集合，但要求该组基函数在某种意义下可以代表整个样本集合。假设这组基函数存在，则上述偏差值集合所有数值之和可表示如下：
这里pi为各偏差值的待定系数，各系数确定后可以使得
最大且满足
将式(3.17)带入式(3.19)得特征值问题
其中Rij＝(Ei，Ej)/M，所以Rij为非负对称矩阵，这样问题就得到了简化，只需要求解出M×M矩阵特征值和特征向量即可。取M×M最大特征值和特征向量时，
取得最大值。
离散状态下，同样可以根据K-L展开式最优化表示。假设状态变量的维度为n行×n列那么每个样本向量维度都为n2，则每个样本可表示为：
其中
代表样本向量
的第j个分量，在网格化模型模拟中代表第j个网格值。
的离散协方差为：
这里
是均值向量，E是期望值，经转换得：
和
样本向量的维度为n2，所以协方差矩阵的维度为n2×n2，n2代表网格点数量。
方阵对应若干特征值的特征向量，其中较大的几个特征值对应的特征向量经过正交化后能捕获样本空间的主要特征，可代表整个样本向量空间参与运算。
1.3支持模型及平台
1.3.1 WRF
WRF(Weather Research and Forecast Model)是新一代中尺度数值天气预报模型，广泛应用于大气相关研究和业务预报。WRF模型的研发始于20世纪90年代后期，是美国气象界通力合作的成果，美国大气研究中心(National Center for AtmosphericResearch，NCAR)、美国海洋和大气管理局(National Oceanic and AtmosphericAdministration，NOAA)、美国空军气象局(AirForceWeatherAgency，AFWA)、美国海军研究实验室(United States Naval Research Laboratory，NRL)、俄克拉荷马大学(TheUniversity of Oklahoma，OU)及美国联邦航空管理局(Federal AviationAdministration，FAA)等部门都参与了研发。WRF模型采用FORTRAN90语言编写，它的主要特点是灵活、易维护、可扩展，可适用于各种计算机平台。WRF模型内部采用先进的数据同化算法和嵌套模式，物理模拟过程较为先进，尤其在对流层模拟方面效果较好。WRF模型主要包括两个核心部分：数据同化系统和并行计算模块。该模型可以提供数十米到数千千米范围内的气象服务，所以无论是全球尺度还是局地小尺度，WRF模型都可以为其提供数据支持。
WRF模型主要包括两个模块：WRF模块和WPS模块。WPS模块主要用来将输入数据(主要为基础地理数据和气象数据)转为模型所需格式数据。WRF模型内部含有namelist.wps配置文件，内部配置了起始时间、嵌套坐标及相对位置等信息。WRF运行前需要边界场和初始场数据(主要由real.exe生成)，wrf.exe生成气象场预报数据，配置内容在namelist.input文件中。WRF最终的生成数据为气象场数据，可供其他模型和同化反演系统使用。
在本文中使用的WRF模型版本为3.6.1，其主要作用是处理美国国家环境预报中心(National Centers for Environmental Prediction，NCEP)提供的FNL再分析数据(FinalOperational Global Analysis Data，FNL)得到气象场数据，为CMAQ模型(CommunityMultiscale Air Quality)、SMOKE模型(Sparse Matrix Operator Kernel Emissions)及同化反演系统提供气象场驱动数据。
1.3.2 SMOKE
SMOKE是排放清单处理模型，采用高性能稀疏矩阵算法，可以将先验排放清单处理成特定时空分辨率的清单数据供其他模型使用。SMOKE模型在1996年投入使用，当时该模型设计较为简单，经过若干年的改进和发展，现在SMOKE模型支持区域、生物源、移动(公路和非公路)和点源排放处理标准，可以处理多种污染物，如一氧化碳(CO)、氮氧化物(NOX)、挥发性有机物(VOC)、氨气(NH3)、二氧化硫(SO2)等。通常各机构发布的排放清单数据的时间分辨率较粗，而空气质量模型一般要求排放清单数据为网格化、逐时数据，因此SMOKE的主要作用是对污染物进行化学谱分配、空间分配和时间分配，可用于创建网格化、逐小时及其他特定要求的排放数据文件，用于各种空气质量模型的输入数据，如CMAQ模型、REMSAD模型(The Regional Modeling System for Aerosols and Deposition)、CAMx(TheComprehensive Air quality Model with extensions)和UAM(Urban Airshed Model)等。
在本文中，SMOKE模型版本为3.0，其主要作用是将MIX清单中的CO、CO2通量数据处理成同化反演系统所需的格式。
1.3.3 CMAQ
CMAQ(Community Multiscale Air Quality)是第三代空气质量模型，其设计目的主要是模拟各物质在传输过程中的物理、化学过程。它是一个三维欧拉(即网格)大气化学和传输模拟系统，可以模拟颗粒物(PM)、臭氧(O3)、有毒空气污染物等80多种物质。相对于第一代、第二代空气质量模型单一尺度、单一污染物而言，CMAQ模型采用“一个大气”的基本思想，可以同时处理从局部到半球范围内的多个空间尺度上多种污染物之间的复杂耦合。CMAQ源代码高度透明和模块化，采用FORTRAN90语言编写，通过由空气质量建模社区的所有成员进行扩展和维护。
CMAQ主要包括五个模块：化学传输模块(CMAQ Chemistry-Transport Model，CCTM)、气象化学接口模块(Meteorology-Chemistry Interface Processor，MCIP)、边界场模块(Boundary Conditions Processor，BCON)、初始场模块(Initial ConditionsProcessor，ICON)、光分解模块(Photolysis Rate Processor，JPROC)。气象化学接口模块是化学传输模式与气象模式的接口模块，WRF模型和MM5模型的输出数据无法为CMAQ模型直接使用，需要通过MCIP将这些数据转换为符合要求的数据格式。初始场模块主要作用是为CCTM模块提供初始的网格化物理、化学数据文件，该文件为浓度垂直分布的初始条件文件，可以根据默认文件产生也可从已有模型输出文件生成。边界场模块主要作用是为CCTM模块提供模拟时间窗口内的边界场数据，边界场数据可以从现有CCTM模拟文件生成或从ICON模块生成的三维垂直剖面浓度生成。光分解速率模块的主要作用是为CCTM运行过程提供光分解信息，该模块可以计算每日晴空光分解速率以及气候衍生臭氧柱和光化学厚度数据，根据化学反应机理可产生逐小时不同经纬度的光分解速率。化学传输模块的主要作用是根据化学反应机制的基本原理对空气中的各种物质进行模拟，其模拟结果包括各物质瞬时浓度、逐时平均浓度以及沉降数据等。化学传输模块中主要包括三类模拟过程：化学反应过程、扩散过程以及化学和气象相结合过程。
CMAQ默认支持的痕量气体并不包括CO2，本文在用CMAQ模拟CO2传输过程前对其源码进行了维护，将CO2添加到ic_profile.dat、bc_profile.dat、hrdata_mod.f和hrinit.f等文件中并重新编译。
在本文中，CMAQ模型作为同化系统的“正演”模型，是同化系统的核心模型，主要用来模拟CO、CO2的传输、沉降、化学反应等过程，起到承接大气CO、CO2浓度和陆面通量的作用。
1.3.4 MPICH
MPI(Message Passing Interface)是为并行计算而设计的，它本质上并不是软件或模型，而是一种关于消息传递的规范和标准，世界各大计算机制造商都提供对MPI规范的支持。计算机在解决低维运算问题时，无需进行复杂的任务消息传递，但在解决高维问题时(如污染物传输数值模拟)则需调度多个计算节点进行运算，此时必须使用服务器的高性能运行模式。高性能运行模式通过协调服务器的计算资源，使得作业具有多节点并发进行的能力，而多个节点之间的协作是通过MPI来实现的。MPI通过提供并行运算库来实现作业的并行化，只需对MPI源码进行编译、安装MPI库即可。
基于MPI的相关规范和标准，不同机构和组织推出了MPI的若干个具体实现，其中更新较快、使用较广的是MPICH。随着MPI内容的不断更新，MPICH也不断推出相应版本，这些版本可以从网上免费下载，对MPICH源代码进行编译、安装后即可构建并行计算环境。
本文采用的MPICH为v2版本，主要实现CMAQ模拟CO、CO2浓度时的并行运算。
1.3.5运行平台
本文构建的区域碳同化反演系统TracersTracker，主程序采用FORTRAN90语言编写，编译器为Intel Fortran Complier。系统涉及到的第三方软件或函数库包括：IOAPI函数库、NetCDF函数库、Matlab函数库，这些函数库主要用于数据的读写、处理和分析。Shell脚本负责整个系统各模块的衔接，另外还包括少量的C、C++模块用来数据处理和分析。本系统运行平台为中国科学院地理科学与资源研究所IBM服务器，该服务器操作系统为Linux系统，版本为64位Red Hat Enterprise Linux Server 5.8。该服务器有134个计算节点，每个节点安装10个CPU。
1.4 TracersTracker(TT)同化反演系统构建方案
1.4.1同化算法
传统四维变分同化方法需要开发伴随模式，然而求取非线性模型的伴随模式是非常困难的。为解决传统四维变分同化方法在实际应用中的困难，本文将本征正交分解法和传统四维变分同化方法相结合，开发了POD4DVar降阶四维变分同化算法，该方法可以使通量反演过程在较小维度内得到解决。POD4Dvar算法采用蒙特卡罗法生成样本集合来表征通量的空间分布格局及时间演变特征，同时该方法无需积分伴随模式，在运行和维护上也较为简便。POD4Dvar算法的基本流程可表述如下：
假设在同化时间窗口(0，T)内有S个步长，根据历史预报数据或蒙特卡罗方法产生N个初始场样本x0，n(n＝1，...，N)。在同化窗口内，从N个初始场开始运行模型
可以得到S个时刻的n个状态向量的样本xk，n(k＝1，2，...，S)，这样就可以构建Xn＝(x0，n，x1，n，..，xs，n)T，对该样本集合进行分解得到基向量，进而可得到问题的最优解，具体步骤如下。
(1)根据蒙特卡罗方法的基本原理，对状态变量的初猜场进行N次随机扰动，产生N个扰动样本xn(t0)，n＝1,2,…,N。
(2)假设在一个同化时间窗口内有K个时刻，基于上述扰动产生的N个初猜场运行模型，则可以在K个时刻产生N个状态变量的模拟值，这些模拟值构成一个集合xn(ti)(i＝0,1,2,…,K)。第tk时刻的N个模拟值的均值为
(3)构建距平矩阵A(δX1(ti)，δX2(ti)，...，δXn(ti))，A∈RM×N，其中M＝Mg×Mv×K，Mg和Mv分别是空间格点个数和状态变量个数。第tk时刻的距平矩阵为
(4)根据距平矩阵构建方阵T′＝AAT，这里T′∈RM×M。M表示分析变量的维度，在一般数值模拟过程中M要远大于N，计算量较大。对方阵T′进行转置，得到T＝ATA，T∈RN×N，由于T的维度一般远小于T′，并且T和T′的特征值是相同的，因此可通过求解T的特征值和特征向量来间接求解T′的特征值和特征向量。若λk和Vk分别为T的特征值和特征向量，则满足以下关系：
TVk＝λkVk，k＝1，2，...，N (3.29)
(5)获取本征正交分解基向量。由于T是对称半正定的，因此存在一组(n个)特征值和特征向量，特征值由大到小可表示为λ1≥λ2≥…≥λn≥0及对应特征向量
(k＝1,2,…,n)，V∈RM×N组成的矩阵。为求取T′的特征向量，作变换F＝A×V，F∈RM×N。为实现矩阵的标准正交化，对矩阵F的每一列做如下处理：
其中，Fi为矩阵F的一列，此时，基向量可表示为：
其中
为矩阵F的第k行i列元素。
(6)构建分析变量解的线性表达式。将上述求解得到的基向量按照线性展开，则分析变量的解最终可按下式表示：
这里p为选择的基向量的个数，一般情况下p远小于N，但要求选择的基向量能够最大限度的捕获向量空间的主要特征，即要求下式尽可能接近于1。
选定其中p个基向量后，基向量矩阵可表示如下：
根据代价函数及分析解的表达式，重构代价函数，则：
经转换，得到：
这里，
另外，背景误差协方差矩阵为：
(7)求取最优分析变量解。代价函数在导数为零处取得极值解，所以令
则可以得到各基向量的最优系数α＝(α1，α2，...，αp)T，即
将求解的基向量的最优系数带入分析变量解的表达式，则可以得到“最优”通量。
1.4.2 TracersTracker(TT)系统框架
系统主要基于Linux Shell脚本和FORTRAN90语言实现，Linux Shell脚本主要负责系统内各模块和接口之间的调度和衔接，实现整个系统任务部署、任务监听及数据读写的自动化。FORTRAN90语言主要实现系统的核心同化算法POD4DVar及相关数值运算。同时，系统还包含若干基于C、C++、NCL语言编写的处理程序，主要实现对输入数据(如先验通量、气象场等)的批处理，使所有数据具有统一的分辨率、步长等。总体而言，系统主要包含三个模块，即系统配置模块、并行运算及监听模块、最优求解模块(图1)，各模块的主要功能如下：
(1)系统配置模块
该模块主要实现系统运行所需数据的自动化配置及处理。系统所需数据主要包括：先验通量、气象场、初始场和边界场等。先验通量来自清华大学MEIC团队提供的MIX排放清单，该清单时间分辨率和空间分辨率较低，本文的主要目的是以该清单为先验通量，通过同化系统获取更高分辨率的后验通量。初始场数据为每次同化提供初始浓度，系统初次运行的初始场来自“spin-up”模拟，即在系统正式同化前一段时间内开始模拟浓度变化信息，使研究区域内浓度场达到基本稳定并将其作为初始场供同化系统使用。后期每次同化循环的初始场则来自于上一同化循环获取的浓度模拟场。气象场采用美国国家环境预报中心提供的FNL再分析数据，FNL数据需要经过WRF模型处理才能转换为同化系统支持的格式，在每次同化前系统自动调用WRF模型对FNL进行预处理。除上述数据之外，系统运行前还需配置若干系统参数，如扰动样本数目、同化窗口时长、同化循环次数、同化物种等。系统配置模块在每次同化前自动产生上述输入数据及系统参数，并将这些信息写入sys.config文件中。
(2)并行计算及监听模块
该模块主要实现数值计算任务的自动分配及监听。该模块首先根据sys.config文件内容读取系统配置模块提供的参数配置及数据，然后根据扰动样本数目启动若干线程运行CMAQ模型对同化物种进行模拟。CMAQ运行时采用并发模式，根据服务器空间使用情况自动安排提交的任务数量。同时，监听模块监听每个任务的具体执行情况，任务执行过程中出现问题被迫中断的，系统自动重启任务。任务执行完毕后，系统将模拟结果拷贝至指定文件夹下供后续数值求解。
(3)最优解算模块
该模块包含系统的核心同化算法POD4DVar，通过数值解算求取最优后验通量。首先读取(2)中CMAQ并行运算的模拟结果，然后基于多源观测数据及四维变分代价函数，计算得到最优通量。最优通量既是当前同化窗口内获取的优化后的通量数据，也是下一同化窗口初始时刻的先验通量。
同化系统在实际运行时，具体流程如下(图2)：
(1)配置数据及参数。将先验通量、初始场、气象场、边界场等数据放入指定文件夹下。在配置文件中设置各数据的具体路径信息，同时在配置文件中设置扰动样本数目、同化窗口时长、同化物种等参数。
(2)初始数据批处理。读取先验通量、气象场、初始场、边界场等数据，调用批处理程序将其转换为统一分辨率和步长的初始数据。在初始同化窗口中，初始场从“spin-up”模拟结果读取。在后续同化窗口中，初始场从上一同化窗口的同化结果中读取。
(3)产生通量扰动样本。对通量进行扰动时主要分两种情况，即初次运行同化系统和非初次运行同化系统。初次运行同化系统时，对先验通量进行扰动，非初次运行同化系统时，是对上一同化窗口获取的“最优”通量进行扰动，进而产生通量扰动样本集合，集合中含有的样本的多少根据(1)中设置的扰动样本数目确定。
(4)并行提交模拟运算任务。(3)中产生的所有通量扰动样本，需要在CMAQ模型中进行积分运算，由于进行该运算是需要占用较多计算资源，在提交运算任务时是在不同时间分批进行的。任务执行期间，监听程序实时监测服务器资源占用情况，在已提交任务结束时自动提交下一任务，并将每次任务结束时产生的模拟数据写入指定文件中。基于(3)中产生的N个通量扰动样本，在同化时间窗口内基于CMAQ模型进行积分运算，可以得到N个CMAQ模拟样本集合。
(5)最优求解。读取(4)中基于通量扰动样本获取的模拟样本集合，如果扰动的N个扰动样本合适，则获取的N个CMAQ模拟样本集合可以较好的表征模拟浓度的时空变化特征。利用获取的N个模拟样本构建分析矩阵，经过转换后对其进行正交分解并进行最优通量求解。求解原理是根据四维变分代价函数，结合地基及卫星观测数据，在“误差最小”情况下得到p个基向量及其系数，最后根据获取的基向量及其系数得到最优通量，具体步骤已在1.4.1中说明。
(6)获取下一同化窗口初始场。基于获取的最优通量数据，重新积分CMAQ模型，可获取该同化窗口内的最优模拟值，该模拟值可作为下一同化窗口初始时刻的初始场。
(7)进入下一同化窗口。
1.4.3同化循环
整个同化过程包含若干次同化循环，每次同化循环的时长由同化窗口决定，图3描述了同化系统的同化循环过程，具体可表述如下：
(1)A1→B1，为第1次同化反演过程。时间长度与同化窗口一致，在反演优化CO通量时设置为2天，在反演优化CO2通量时设置为1周。该过程利用TracersTracker反演系统同化观测数据，反演得到该同化窗口内的“最优”通量。
(2)B1→A2→B2，为第2次同化反演过程。时间长度与步骤(1)相同。该过程将步骤(1)得到的后验模拟作为本次同化反演过程的初始场，同化反演本次同化窗口内的“最优”通量。
(3)B2→…→Bn，为第3次到第n此同化反演过程。类似步骤(1)和步骤(2)，依次执行后面若干次同化反演过程，得到对应同化窗口内的反演优化通量，直至覆盖整个所有同化窗口。
1.5评价指标
本文主要采用如下评价指标对同化系统反演结果的准确性进行评价，具体如下。
(1)均方根误差RMSE，用来表征观测值与模拟值之间的误差水平，计算公式如下：
其中，n是观测值的数量，e是模拟值和观测值的差值，o和s分别代表观测值和模拟值。
(2)线性相关系数R，用来表征观测值与模拟值间的相关性，计算公式如下：
其中，n是观测值的数量，e是模拟值和观测值的差值，o和s分别代表观测值和模拟值。
(3)绝对误差均值MAE，用来表征观测值与模拟值之间绝对误差的均值水平，计算公式如下：
在本文中，各变量含义与(1)相同。
(4)误差均值ME，用来表征观测值与模拟值之间误差的均值水平，计算公式如下：
在本文中，各变量含义与(1)相同。
1.6 TracersTracker系统有效性验证
本节通过理想实验对同化系统的有效性进行验证，即在理想状态下检验TracersTracker系统同化观测数据反演通量的准确性，其中CO通量反演有效性检验的具体步骤如下：
(1)构建CO理想通量。将MIX排放清单中CO通量作为理想通量，假设理想通量是真实的、无误差的通量。
(2)构建CO先验通量。在CO理想通量基础上随机添加30％的“误差”作为先验通量，先验通量是需要校正、优化的初始通量。
(3)获取CO理想观测。基于CO理想通量积分CMAQ模型得到CO浓度模拟值，将该模拟值作为“观测数据”，由于该“观测数据”并非真实的观测信息，以下称为理想观测数据。
(4)获取CO后验通量。同化系统同化CO理想观测数据，逐步校正、优化先验通量，反演得到CO后验通量。同化反演过程中，95％的理想观测数据参与同化，其余5％的理想观测数据作为验证数据对后验通量的准确性进行验证。
CO通量反演有效性检验以上甸子区域(图1)为实验区域，实验的参数配置具体为：两层嵌套(图30)，外层(d02)为内层(d01)提供边界场，同化CO观测数据反演得到内层CO后验通量，模型内部采用CB05化学机制，扰动样本数量设置为126，滞后窗口设置为0，空间分辨率设置为3km，上述参数的设置依据第四章的参数敏感性分析结果确定。另外，为深入分析TracersTracker系统在不同季节进行CO通量反演的差异，分别在一月和七月对TracersTracker系统进行了有效性检验。
表3-1CO理想实验结果
实验结果表明(表3-1)，基于添加了“误差”的先验通量，CO模拟误差出现明显增大、相关系数出现明显降低。其中一月份CO先验模拟的绝对误差均值达到0.25毫克/立方米，均方根误差达到0.36毫克/立方米，相关系数也比较低，仅为0.86(图4)。七月份CO先验模拟的绝对误差均值达到0.16毫克/立方米，均方根误差达到0.23毫克/立方米，相关系数也不高，仅为0.87(图5)。
经过同化系统优化后，CO后验通量的模拟误差出现了明显的降低，其中一月份后验模拟绝对误差均值减小至0.11毫克/立方米，均方根误差减小至0.14毫克/立方米，七月份后验模拟的绝对误差均值减小至0.08毫克/立方米，均方根误差减小至0.10毫克/立方米。同时，模拟值与观测值的相关系数也有明显提高，一月份相关系数从0.87提高到了0.98，七月份相关系数从0.87提高到了0.97，模拟值和观测值具有非常好的一致性(表3-1)。由此可见，同化系统对先验通量的校正是非常明显的，在很大程度上消除了先验通量中存在的“误差”，有效提高了CO后验模拟的精度。
另外，无论是在一月还是七月，同化前CO模拟值均明显小于理想观测，该现象在同化后得到校正，同化后CO模拟值基本以0值为中心呈正态分布。另外，同化前后CO误差分布区间也有较大变化。在一月份，同化前CO误差基本分布在[-2.5,2.5]范围内，同化后误差范围降低到[-0.8,0.8]，无论是误差区间还是误差极值均有较明显减小(图6)。同样，在七月份，同化后的误差分布区间也出现了较明显缩小，由同化前的[-2,2]减小到[-0.65,0.65](图7)。
总体而言，上述反演结果表明，TracersTracker系统可以较好的消除先验通量中的“误差”，使后验通量逐步趋近于理想通量。虽然后验通量中仍存在误差，但后验模拟的误差整体上较小、并不明显，后验通量的时空分布上与理想通量一致性较好(图8、图9)。
CO2通量反演的有效性验证实验与CO的基本相同，不同之处是在滞后窗口的设置上。CO和CO2的生命周期有较大差异，在理想实验中CO2的滞后窗口设置为1周，同化系统各参数的具体配置为：两层嵌套(图30)，外层(d02)为内层(d01)提供边界场，同化CO2观测数据反演得到内层CO2后验通量，模型内部采用CB05化学机制，扰动样本数量设置为126，滞后窗口设置为1周，空间分辨率设置为3km。
表3-2 CO2理想实验结果
在CO2理想通量基础上添加“误差”后，CO2模拟值的误差明显增加，其中先验模拟的绝对误差均值在一月份和七月份分别为10.3ppm和9.5ppm，均方根误差分别为13.8ppm和13.5ppm。同时，无论是在一月份还是七月份，CO2先验模拟和理想观测的相关系数也有明显的降低，在一月份两者的相关系数为0.78，七月份两者的相关系数为0.81(图10、图11)。
CO2先验通量在经过同化系统优化后，在一月份和七月份均显著改善了CO2
的模拟精度。在一月份，CO2后验模拟的绝对误差均值由10.3ppm减小至2.31ppm，七月份CO2后验绝对误差均值由9.5ppm减小至2.75ppm，模拟精度有了明显提高。同时，模拟值和理想观测值的相关系数在两个月份也有较大提高，一月份两者的相关系数由0.78提升至0.99，七月份两者的相关系数由0.81提升至0.98，后验模拟与理想观测值保持了较好的一致性，与理想观测值的变化趋势基本一致(表3-2)。
另外，先验通量经过同化系统优化后，CO2后验模拟误差的分布形态也有明显改变。在一月份，同化前CO2误差基本分布在[-45,45]范围内，同化后误差范围降低到[-4.5,4.5]，仅为同化前的10％左右(图12)。在七月份，同化后的误差范围也出现较大缩小，由同化前的[-45,45]减小到[-4.0,4.0](图13)。
从CO2通量的反演结果看，同化系统可以较好的消除先验通量中的“误差”，使后验通量逐步趋近于理想通量。虽然后验通量中仍存在误差，但误差总体上较小、并不明显，后验通量的时空分布上与理想通量一致性非常好(图14、图15)。
上述结果表明，本文构建的TracersTracker同化系统可以有效吸收CO和CO2观测数据中的信息，在很大程度上校正先验通量中的“误差”，获取较高精度的后验通量，切实提高CO、CO2后验模拟的准确性。
1.7 TracersTracker同化反演系统构建方案小结
同化反演系统构建方案包括四部分内容：
第一部分：介绍了本发明系统构建的同化反演系统的相关理论。本文构建的同化反演系统是在传统四维变分同化法和本征正交分解法的基础上实现的，在本章中详细介绍了四维变分同化法和本征正交分解法的基本原理及推导过程。
第二部分：介绍了本发明技术构建的同化反演系统所需的支持模型及运行平台。所构建的同化反演系统是架构在若干模型上的，其中WRF模型为系统提供气象场，SMOKE模型用来处理排放清单，CMAQ模型是同化系统的核心“正演”模型，MPICH提供并行运算服务，本章对这些模型及平台做了详细的介绍并说明了其在同化反演系统中扮演的角色。
第三部分：阐述本发明技术构建的同化反演系统的方案和系统框架。首先介绍了同化反演系统的核心同化算法POD4Dvar及其推导过程，然后详细阐述了该系统的基本框架和功能模块，在此基础上阐述了系统具体的运行流程。
第四部分：对同化反演系统通量反演效果进行了有效性验证。针对CO和CO2，在一月和七月时段内，对系统反演通量的有效性进行了检验，结果表明，本文构建的同化反演系统可以充分吸收CO和CO2观测信息，有效消除先验通量中的“误差”，明显提高模型模拟的精度。
2.TracersTracker同化反演系统敏感性分析
2.1引言(Preface)
在基于发明技术所构建的TracersTracker系统中，输入数据及参数的设置会直接影响系统最终的同化结果，如扰动样本数目、滞后窗口、边界场等，这些输入数据及参数需进行合理优化，以获得同化系统最优的同化效果和同化效率。参数的敏感性分析是进行参数优化的有效手段，可以帮助研究人员理解输入数据和参数对同化结果影响的重要性。参数敏感性分析通过改变输入数据和参数进行测试，定量给出各参数变化对同化结果的具体影响规律。在参数的合理范围内，如果参数值的改变导致同化结果产生明显改变，则说明该参数敏感性较强，反之则说明参数敏感性较弱。最后依据各参数的敏感性强弱程度，可以将各参数设置在合理范围内，既保证同化结果的准确性，同时兼顾同化系统的运行效率。
依据敏感性分析的基本原理，对TracersTracker系统中扰动样本数目、边界场、滞后窗口、化学机制、空间分辨率等输入数据和参数进行了测试，给出了一套较为合理的参数值。
2.2 TracersTracker敏感性实验设计
主要研究CO和CO2的通量反演，所以本节分别对CO和CO2通量反演过程中的输入数据和参数进行了敏感性分析。CO敏感性分析主要包含8个实验，分别对扰动样本、边界场、滞后窗口、化学反应机制、空间分辨率等参数进行敏感性分析，各实验的具体设置如下：
主实验：两层嵌套，外层(d02)为内层(d01)提供边界场，同化CO观测数据反演得到内层CO后验通量，模型内部采用CB05化学机制，扰动样本数量设置为126，滞后窗口设置为0(即不考虑滞后性)，空间分辨率设置为3km。该实验为敏感性分析的主实验，其他实验都以该实验为参照对比分析。
实验1：扰动样本数目实验。扰动样本数目设置为60，其他参数设置与主实验相同。本实验主要研究扰动样本数目对CO通量反演结果的影响。
实验2：扰动样本数目实验。扰动样本数目设置为200，其他参数设置与主实验相同。本实验目的与实验1相同。
实验3：边界场实验。取消两层嵌套模式，采用单层非嵌套模式，其他参数设置与主实验相同。本实验主要研究嵌套模式(即边界场)对CO通量反演结果的影响。
实验4：边界场实验。取消两层嵌套模式，采用三层嵌套模式，其他参数设置与主实验相同。本实验目的与实验3相同。
实验5：滞后窗口实验。滞后窗口设置为2天，其他参数设置与主实验相同。本实验主要研究滞后窗口长短对CO通量反演结果的影响。
实验6：滞后窗口实验。滞后窗口设置为2周，其他参数设置与主实验相同。本实验目的与实验5相同。
实验7：化学机制实验。内部化学机制改为CB06化学机制，其他参数设置与主实验相同。本实验主要研究不同化学机制对CO通量反演结果的影响。
实验8：空间分辨率实验。空间分辨率设置为27km，其他参数设置与主实验相同。本实验主要研究空间分辨率对CO通量反演结果的影响。
CO敏感性实验以上甸子区域为实验区域，TracersTracker系统同化该区域内的17个站点的CO浓度观测数据反演得到该区域内的CO后验通量，利用区域内的5个验证站点的观测数据对后验通量及其模拟值进行验证。
在进行敏感性实验前，首先基于先验通量对CO浓度进行了模拟，先验通量模拟的CO浓度精度明显低于主实验模拟的CO浓度，模拟绝对误差为0.77毫克/立方米，主实验后验通量模拟值的绝对误差为0.59毫克/立方米，模拟精度有明显提升。另外，对先验通量进行优化后，CO模拟浓度相关系数从0.6提高到0.76，CO模拟浓度与观测浓度相关系数也有较大提高(表4-1)。
表4-1主实验CO通量增比及模拟误差
另外，CO后验通量较先验通量出现了较明显的提高，从全年视角看，年均通量从1.14摩尔/秒提升到1.63摩尔/秒，提升了32.9％。从月际变化来看，12个月的CO后验通量较先验通量分别提高了38.7％、34.8％、35.5％、33.2％、29.6％、27.3％、33.2％、26.5％、30.6％、27.9％、36.7％和40.9％，增比较高的为1月和12月(表4-1)。从季节上看，春、冬两季增长比例要明显低于夏、秋两季，所以先验清单中春、冬两季的低估程度要大于夏、秋两季(图16)。
CO2的敏感性实验包含8个实验，分别对扰动样本、边界场、滞后窗口、空间分辨率等参数进行敏感性分析，实验的具体设置如下：
主实验：滞后窗口设置为1周，其他参数配置与CO主实验相同。该实验为CO2敏感性分析的主实验，其他实验都以该实验为参照。
实验1：扰动样本数目实验。参数配置与CO实验1相同。本实验主要研究扰动样本数目对CO2通量反演结果的影响。
实验2：扰动样本数目实验。参数配置与CO实验2相同。本实验主要目的与实验1相同。
实验3：边界场实验。参数配置与CO实验3相同。本实验主要研究嵌套模式(即边界场)对CO2通量反演结果的影响。
实验4：边界场实验。参数配置与CO实验4相同。本实验主要目的与实验3相同。
实验5：滞后窗口实验。参数配置与CO实验5相同。本实验主要研究滞后窗口的长短对CO2通量反演结果的影响。
实验6：滞后窗口实验。参数配置与CO实验6相同。本实验主要目的与实验5相同。
实验7：化学机制实验。参数配置与CO实验7相同。本实验主要研究不同化学机制对CO2通量反演结果的影响。
实验8：空间分辨率实验。参数配置与CO实验8相同。本实验主要研究空间分辨率对CO2通量反演结果的影响。
表4-2主实验CO2通量增比及模拟误差
由于没有充足的CO2浓度观测数据进行敏感性实验，实验中采用的CO2观测数据是将CO观测数据乘以不同比例系数得到的，比例系数是依据大气本底站点中CO、CO2观测数据比例确定的。通过上述方法最终得到上甸子区域内17个站点的“CO2观测数据”和5个站点的“CO2验证数据”，各站点的位置与CO观测站点相同。
CO2先验通量模拟的CO2浓度精度与主实验没有明显区别，先验通量模拟的CO2的绝对误差为11.8ppm，相关系数为0.67。主实验后验模拟值相关系数从先验模拟的0.67变为0.66，绝对误差减小到11.7ppm(表4-2)。
在CO2通量方面，同化反演后的CO2通量大部分出现了明显的升高，整体上从35.6摩尔/秒提升到46.6摩尔/秒，1月-12月分别提高了49.3％、36.9％、20.5％、22.4％、25.9％、16.3％、27.9％、27.2％、36.9％、32.7％、23.2％和46.9％，某些月份的后验通量增长明显，如12月和1月(图17)。
2.3扰动样本实验
在CO扰动样本实验中，采用126个扰动样本与200个扰动样本的CO模拟浓度差异没有明显区别，其相关系数分别为0.76和0.77，两者的绝对误差也很接近，分别为0.59毫克/立方米和0.58毫克/立方米。这说明扰动样本数目增加到一定程度后，同化反演结果的精度很难再有明显提升。相比而言，采用60个扰动样本进行同化反演时，反演结果则存在较大误差，CO模拟浓度的相关系数只有0.42，相关系数出现了明显的降低。同时绝对误差则从0.59毫克/立方米增加到0.78毫克/立方米，模拟精度也出现了明显的下降(表4-3)。
表4-3扰动样本实验CO通量增比及模拟误差
从CO通量反演结果看(图18)，采用200个扰动样本的CO通量反演结果与主实验反演得到的结果较为接近，年均通量由1.63摩尔/秒变为1.65摩尔/秒，几乎没有变化。在全年12个月中，相较于先验通量分别提高了38.2％、38.3％、37.1％、36.7％、31.7％、31.4％、28.6％、30.4％、34.3％、30.8％、34.1％和40.6％，与主实验反演结果的趋势基本保持一致。相比而言，基于60个扰动样本的同化反演结果虽然在总量上没有明显差异、略有下降(1.59摩尔/秒)，但在月际变化中有较为明显的震荡，12个月的变化分别为32.1％、11.8％、9.7％、65.7％、44.4％、10.6％、10.1％、52.7％、39.2％、6.6％、9.2％和54.2％，1月、2月、11月、12月的反演结果均较为异常，出现了很大的升高或降低。
从对比结果看，扰动样本数量在CO通量反演过程中具有重要的作用。扰动样本数量较少时，样本集合无法充分反映CO排放的时空特性，再加上随机产生的不确定性，极有可能造成反演结果的震荡，出现较大偏差。
与CO敏感性实验一样，CO2样本实验中分别在实验1和实验2中分别采用了60和200个扰动样本。实验结果表明，使用126个样本与200个样本得出的后验CO2模拟浓度差异并不明显，其相关系数均分别为0.56和0.58，没有明显变化(表4-4)。
表4-4扰动样本实验CO2通量增比及模拟误差
两者的绝对误差也很接近，分别为11.7ppm和11.4ppm。这说明与CO类似，扰动样本增加到一定程度后，同化反演精度已无法明显提高。相比而言，采用60个扰动样本进行同化反演时，反演结果则存在较大波动和误差，从图6可以看到，CO2浓度的相关系数只有0.51，绝对误差则上升到16.1ppm，误差出现了明显的提高，模拟相关系数出现了明显的降低(表4-4)。
在CO2通量反演结果上(图7)，200样本反演通量全年通量为46.6摩尔/秒，与主实验通量相同。在12个月增长变化为49.9％、37.4％、20.3％、22.7％、25.7％、16.6％、28.5％、27.1％、37.8％、32.4％、23.1％和46.6％，在大部分时段
与主实验通量保持一致，只在少数月份有较小差异。相比而言，60个扰动样本的反演结果则有较大波动和不确定性，其年均通量为47.1摩尔/秒，虽然总体上与主实验结果没有明显变化，但月际变化差异较大，在1月-12月增比分别为26.7％、32.2％、29.3％、13.5％、20.7％、21.4％、39.3％、39.6％、30.7％、28.5％、29.0％和66.6％，与主实验通量月际变化相比，整体波动非常大，存在很大不确定性。
在进行正式的CO和CO2通量同化反演之前，本文针对扰动样本设计并实施了若干次实验，结果表明，扰动样本数量在110-126之间时，同化反演结果皆能取得较为满意的反演精度，当样本数量大于120时，反演结果精度提升较小。理论上，扰动样本数量越多，样本空间对通量的时空变化特性把握越好，同化反演结果也应当越好，但样本数量太多会极大的增加计算成本和时间成本，为兼顾反演精度与效率，在本文的实际反演过程中，扰动样本数量均设置为126。
2.4边界场实验
不设置边界场(单层非嵌套模式)的CO模拟浓度误差值较为明显，绝对误差由主实验中的0.59毫克/立方米增加到0.71毫克/立方米，相关系数由0.76降低为0.54(表4-5)。
重新设置边界场后，绝对误差降到0.59毫克/立方米，相关系数也提高到了0.75，略低于主实验模拟浓度相关系数，在模拟精度上出现了较明显的提高。这说明边界场的有无及准确与否对同化反演结果有着重要的作用(表4-5)。
在通量反演结果看(图8)，采用默认边界场数据的CO通量年均值为1.57摩尔/秒，低于主实验通量。另外，后验通量在月际变化中也有较大区别，1月-12月变化分别为18.9％、6.6％、44.8％、30.0％、56.4％、9.6％、23.4％、45.8％、48.4％、5.9％、41.0％和17.9％，其波动性较大且与主实验结果差异明显，尤其是在2月、5月、6月和10月。重新设置边界场的CO通量年均值为1.64摩尔/秒，接近主实验得到的CO通量结果，月际变化分别为39.8％、34.4％、36.4％、32.1％、30.1％、26.3％、32.8％、27.1％、31.7％、26.8％、36.2％和41.0％，除较小差异外与主实验结果的变化趋势基本保持一致。
表4-5边界场实验CO通量增比及模拟误差
CO2边界场敏感性实验中，单层非嵌套模式后验CO2模拟浓度误差明显升高，相关系数只有0.47，绝对误差达到16.7ppm。采用三层嵌套模式后，相关系数提高到0.68，绝对误差减小到11.6ppm，模拟精度出现了明显提高，基本与主实验通量模拟精度持平。合理设置CO2边界场，可以明显提高CO2同化反演通量的准确性(表4-6)。
表4-6边界场实验CO2通量增比及模拟误差
从通量反演结果上看，单层非嵌套模式后验通量年均值为46.4摩尔/秒，总量上无明显变化，但其月际通量变化明显，默认边界场通量月际变化分别为30.1％、54.9％、16.8％、17.1％、17.7％、20.2％、14.2％、39.2％、49.2％、28.7％、33.5％和40.4％，整体上波动较大，存在较明显不确定性。相比而言，三层嵌套通量无论是整体上还是月际变化，均与主实验通量基本保持一致。三层嵌套通量年均值为46.5摩尔/秒，与主实验通量较为接近，月际变化分别为48.9％、37.6％、20.8％、23.0％、26.4％、16.2％、27.7％、27.3％、36.9％、32.7％、23.5％和46.8％，与主实验通量差异较小、基本一致。
单层非嵌套模式操作简单、计算效率较高，但是会给模拟过程中的边界场带来不确定性，合理设置嵌套层数可以为内层与外层的CO、CO2交互转换提供较为准确的边界场，从而明显提高CO、CO2通量同化反演的准确性。但在实际操作中，嵌套层数越多消耗的计算资源会明显增加，严重影响计算效率。在本文中，实际反演过程中均采用两层嵌套模式。
4.5滞后窗口实验(Experiment for Lag Window)
表4-7滞后窗口实验CO通量增比及模拟误差
实验结果表明(表4-7)，滞后窗口为2天与2周得到的CO后验模拟浓度没有明显差异，与观测值相关性均较好，2天滞后模拟相关性略有提高(R由0.76变为0.77)，2周滞后模拟相关性略有降低(R由0.76变为0.75)，但总体而言变化并不明显。同时两者的绝对误差也未出现明显变化，分别为0.59毫克/立方米和0.61毫克/立方米，滞后2天模拟结果误差有所减小，滞后2周模拟误差有所增大，但都与主实验的误差水平较为接近。
从CO通量反演结果上看，滞后2天的后验通量在总量上、月际变化趋势上均与主实验通量较为接近，其年均通量为1.64摩尔/秒，月际变化分别为38.3％、34.3％、36.4％、32.0％、30.6％、26.1％、33.3％、26.1％、30.6％、27.0％、36.3％和39.9％。滞后2周通量在总量上与主实验通量基本一致(1.63摩尔/秒)。月际变化分别为35.9％、38.6％、33.0％、33.0％、27.6％、27.6％、33.7％、24.3％、31.2％、27.6％、34.4％和42.8％，与主实验结果略有差异但趋势基本一致。
表4-8滞后窗口实验CO2通量增比及模拟误差
CO2滞后窗口实验与CO不同，一般情况下CO2的滞后窗口要比CO长(5-9周)，其主要原因是因为CO2为惰性气体，在大气中较为稳定，需要较长时间与其他区域的大气进行交互融合。由于本文研究区域较小，CO2进行交互融合所需时间也较短，所以本文在对CO观测数据进行同化时并未设置滞后窗口，CO2主实验的滞后窗口设置为1周。同时，本节另外设置了滞后窗口为2天和2周的同化反演实验，以此滞后窗口对同化反演结果的具体影响。
实验结果表明(表4-8)，2天滞后后验模拟与2周滞后后验模拟精度变化不太明显，基本与主实验模拟精度持平，2天滞后后验模拟绝对误差为11.8ppm，与主实验模拟精度基本一样，相关系数为0.66，与主实验模拟相关系数持平。2周滞后后验模拟绝对误差为11.8ppm，与主实验模拟精度相比略有升高，相关系数为0.65，略低于主实验后验模拟相关系数。
反演通量方面，2天滞后后验通量年均值为46.8摩尔/秒，略高于主实验通量年均值，通量月际变化分别为48.9％、37.7％、20.4％、23.4％、25.9％、15.6％、29.0％、28.1％、37.7％、34.0％、24.3％和48.8％，无论是增长幅度还是增长变化趋势都与主实验保持一致。相较于2天滞后，2周滞后后验通量波动则有轻微增加，其通量年均值为46.4摩尔/秒，略低于主实验通量。在通量变化方面，2周滞后窗口通量波动也有增加，其月际通量变化分别为47.4％、36.8％、20.3％、21.8％、25.1％、16.4％、27.3％、27.6％、35.4％、31.8％、22.4％和48.0％，其中9月-12月的通量变化与主实验后验通量差异较明显。
“滞后窗口”是碳同化反演系统中较为重要的参数，尤其是在全球或半球等较大尺度的碳源汇反演研究中。碳元素以不同速度在地球各个圈层中不断迁移转换，在岩石圈中交互转移较慢，在大气圈中较快。所以大气中的CO2浓度是不断变动的，空气监测站点的CO2观测浓度是较长时间内碳源汇不断交互融合、动态演变的结果，当前时间段碳的排放与吸收情况直接决定了一段时间后大气中CO2的浓度。在碳同化反演系统尤其是较大尺度同化反演系统中，“滞后窗口”的设置在很大程度上影响着反演结果的合理性和精确性，因为较大区域内的碳源、碳汇需要在“滞后窗口”时间段内达到稳定状态。滞后窗口设置的长短，很大程度上决定着同化过程吸收的信息的多少，滞后窗口设置的长，则可以抓取更多观测信息，滞后窗口设置的短，则获取的观测信息可能不足。但是，滞后窗口设置过长，会占用过多计算资源，影响同化系统的正常运行。Peters在全球碳同化过程中，通过一系列测试实验，指出5-6周滞后窗口设置可以捕捉80％观测信息，滞后窗口为12周时基本可以捕捉到所有观测信息，但是运算效率较低。所以，在全球尺度的碳同化过程中，一般将滞后窗口设置为5周左右，区域尺度的碳同化过程滞后窗口的设置需要根据区域范围大小进行相应设置。Wouter等通过大量实验证明，5-9周的“滞后窗口”的设置可以使模型较好较快的反演出碳源汇的时空分布状况。另外，不同气体的化学特性不同、在自然界中的转移交换过程各异，滞后窗口的大小也存在很大差异。CO2化学特性较为稳定，不易发生化学反应，在大气中滞留时间较长，所以其滞后窗口也较长。相比而言，CO化学特性要活跃很多，易于发生氧化反应，因此其滞后窗口要小得多。
TracersTracker系统主要应用于较小尺度下(瓦里关区域和上甸子区域)CO2和CO通量的同化反演。相较于大尺度下(如全球尺度)的碳通量反演，区域碳通量反演的滞后窗口要小很多。在本文中，根据滞后窗口敏感性实验结果，在CO2通量反演过程中滞后窗口均设置为1周，由于CO的生命周期较短(大约2个月)，CO通量反演过程中并未考虑滞后窗口，即滞后窗口均设置为0。
2.6化学机制实验
CO主实验采用了CB05化学反应机制来描述各物质化学反应，在CO实验7中，将化学反应机制替换为CB06，以此研究不同化学机制在同化反演过程中的敏感性。
表4-9主实验CO通量增比及模拟误差
如表4-9所示，在模拟误差上，CB05机制模拟误差和CB06机制模拟误差基本持平，绝对误差分别为0.59毫克/立方米和0.58毫克/立方米。在模拟相关性上，采用CB05和CB06化学反应机制对模拟结果几乎没有影响，相关系数均为0.76。
从通量反演结果上看，CB06机制通量年均值为1.64摩尔/秒，略高于主实验通量，月际变化分别为38.3％、35.8％、34.7％、32.5％、30.4％、27.8％、33.2％、27.5％、31.2％、28.3％、36.3％和41.0％，CB05机制通量和CB06机制通量无论是在总量上还是在月际变化趋势上都非常接近。
将化学机制由CB05转为CB06后，CO2与CO的表现类似，在误差和相关性上都没有明显变化。实验结果表明(表4-10)，不同化学机制对CO2后验模拟精度基本没有影响。CB06机制下CO2模拟的绝对误差为11.6ppm，与主实验模拟精度基本一样，相关系数为0.67，与主实验模拟相关系数持平。
表4-10主实验CO2通量增比及模拟误差
从通量反演结果上看，CB06机制下的后验通量年均值为46.6摩尔/秒，与主实验后验通量持平，月际变化分别为49.7％、37.0％、20.6％、22.4％、25.7％、15.6％、27.7％、27.7％、36.5％、32.6％、23.7％和46.8％，CB05机制CO2通量和CB06机制CO2通量无论是在总量上还是在月际变化趋势上都基本一致。
由于不同化学机制对CO和CO2通量反演结果影响甚微，在本文中，实际通量反演过程中均采用CB05化学反应机制。
2.7空间分辨率实验
表4-11空间分辨率实验CO通量增比及模拟误差
以全球尺度为代表的较大尺度下的通量反演结果分辨率较粗，通常以度(°)为基本单位，这种粗分辨率通量结果无法为区域高分辨率通量研究提供支持。本文的同化反演通量数据空间分辨率为3km，相对于全球尺度的通量产品，空间分辨率得到了极大的提升。
在实验8中，对空间分辨率为27km的通量反演结果和模拟精度进行了分析，以此测试在空间分辨率变粗对同化反演结果的具体影响。结果显示(表4-11)，空间分辨率变粗会影响反演通量的模拟精度。在模拟误差方面，27km分辨率通量模拟的绝对误差从0.59毫克/立方米增加到0.63毫克/立方米，误差出现了明显升高。在模拟的相关性上，27km分辨率通量模拟相关系数从3km分辨率的0.76降低到0.66，也出现了降低。总体而言，空间分辨率变粗会给通量反演结果带来不确定性，增大模拟结果的误差。
在通量反演结果上，空间分辨率变粗也会带来较明显的震荡。27km分辨率通量反演结果年均值为1.58摩尔/秒，明显低于主实验通量。在月际变化上也与主实验通量有明显差异，27km分辨率通量12月变化分别为36.0％、30.5％、47.7％、20.2％、16.3％、42.8％、38.1％、6.7％、39.8％、39.4％、16.2％和16.6％，与主实验通量增减有较为明显的差异。
表4-12空间分辨率实验CO2通量增比及模拟误差
CO2主实验中CO2同化反演通量数据空间分辨率为3km，在实验9中，对空间分辨率为27km的通量反演结果和模拟精度进行了分析，以分析空间分辨率变粗对同化反演结果的具体影响。
实验结果表明(表4-12)，将同化反演系统分辨率调大之后，后验通量的模拟精度会出现一定程度的降低，相关系数也会变差。在模拟误差方面，27km分辨通量模拟的绝对误差从11.7ppm增加到13.1ppm，误差出现了升高。在模拟的相关性上，27km分辨通量模拟相关系数从3km分辨的0.66降低到0.62，也出现了降低。总体而言，空间分辨率变粗会给通量反演结果带来不确定性、增加模拟结果的误差。
在通量反演结果上，空间分辨率变粗也会带来较明显的震荡。27km通量反演结果年均值为47.3摩尔/秒，明显高于主实验通量，并且在月际变化上与主实验通量也有明显差异，27km通量12月变化分别为43.8％、42.9％、16.3％、32.2％、15.2％、23.4％、29.6％、15.7％、42.6％、42.4％、28.8％和56.5％，与主实验后验通量有较为明显的差异。
空间分辨率变粗，必然会使系统在反演过程中忽略一些信息，继而引入不同程度的不确定性，所以在区域碳同化反演系统中，在计算资源允许的情况下应尽量提高空间分辨率。在本文中，实际通量反演过程中空间分辨率均设置为3km。
2.8 TracersTracker同化反演系统敏感性分析小结
该部分阐述了同化系统敏感性分析的重要意义，详细说明了本文敏感性实验的具体实施方法，主要包括三部分内容：
第一部分：针对CO设计并实施了若干敏感性实验，通过这些实验对本文研发的TracersTracker同化反演系统进行了测试。CO敏感性实验主要包括主实验、扰动样本实验、边界场实验、滞后窗口实验、化学机制实验和空间分辨率实验。
第二部分：针对CO2设计并实施了若干敏感性实验，通过这些实验对本文研发的TracersTracker同化反演系统进行了测试。CO2敏感性实验包括主实验、扰动样本实验、边界场实验、滞后窗口实验、观测数据实验和空间分辨率实验。
第三部分：根据敏感性实验的结果，确定了同化系统进行实际同化反演时的参数设置，具体为：两层嵌套，CB05化学机制，扰动样本数量为126，CO滞后窗口为0，CO2滞后窗口为1周，空间分辨率为3km。
3 TracersTracker区域高分辨率CO通量反演技术
CO通量同化反演的基本流程是利用TracersTracker系统同化CO地基浓度观测数据，不断对CO先验通量进行校正和优化，最终获取精度更高的CO后验通量。根据第四章参数敏感性分析的结果，TracersTracker系统在CO通量反演过程中的参数配置为：两层嵌套、CB05化学反应机制、扰动样本数目为126，滞后窗口为0、空间分辨率为3km。具体步骤如下：
(1)准备基础数据。
CO通量反演所需数据具体包括：边界场、初始场、气象场及CO先验通量。同化系统采用两层嵌套模式，外层模拟场为内层模拟提供边界场数据。初始场来自spin-up模拟，即在正式开始同化反演前，利用CMAQ模型先模拟一段时间(2015年12月1日-2015年12月31日)CO浓度，使内层CO浓度达到基本稳定，并以最后的CO模拟状态作为初始场。气象场采用FNL再分析数据，利用WRF模型将其转换为同化系统兼容格式。CO先验通量采用MEIC团队提供的MIX排放清单，使用前需利用SMOKE模型将其转为同化系统兼容的格式。
(2)获取CO背景浓度。
即使在没有明显污染源的情况下，大气中也存在一定浓度的CO，称为背景浓度。不同区域的背景浓度存在明显差异，它受当地的经济水平、地形、气候等因素的综合影响，是经过长期的迁移、转换、融合而形成的。因此，各监测站观测的CO浓度主要包括两个部分：背景浓度(CObg)和增量浓度(ΔCO)，其中增量浓度是由当前CO排放活动引起的。CO通量反演的目的是获取当前CO的排放通量，同化的CO观测数据应为增量浓度。所以，在TracersTracker系统同化CO浓度观测数据之前，首先需求出区域内的CO背景浓度，然后将实际观测浓度减去背景浓度，获取当前排放活动所导致的增量浓度，最后将增量浓度参与到通量同化反演中。
求取CO背景浓度的具体步骤主要包括：①.剔除异常数据。异常数据主要是由仪器故障或其他突发情况造成的，这种观测数据无法反应真实的CO浓度，在使用前需将其删除。②.删除“污染”观测。设置阈值(明显含有增量浓度的观测值)，将超过阈值的观测数据删除。③.在剩余的CO浓度观测数据中，以72小时为平滑窗口将观测值超过均值±3倍标准差的观测值移除。④.循环筛选。不断重复步骤③，直至剩下的数据集中的数量趋于稳定为止，拟合剩下的观测数据即为最终的CO背景浓度。
(3)获取CO通量扰动集合。
在CO先验通量的基础上，随机添加30％的扰动，产生指定数目的CO通量扰动样本。若同化系统是首次运行，则扰动是基于CO先验通量进行的。若同化系统不是首次运行，则扰动是对上一同化窗口优化后的CO通量进行的。
(4)获取CO浓度模拟样本。
基于(3)中产生的CO通量扰动样本，运行CMAQ模型模拟CO浓度。CMAQ运行的次数与通量扰动样本数目是相同的，所有模拟任务结束后得到同等数目的CO浓度模拟样本。
(5)最优求解。
读取(4)中CO模拟样本，构建模拟误差矩阵并进行最优通量求解，求解原理是根据四维变分同化方法中的代价函数，结合CO地基站点观测数据，在统计误差最小条件下得到p个基向量及其系数，根据获取的基向量及其系数得到最优通量，具体求解步骤已在3.4.1中给出。
(6)获取下一同化窗口初始场。
基于获取的CO后验通量，重新积分CMAQ模型，可获取该同化窗口内的CO模拟，该模拟作为下一同化窗口初始时刻的初始场。
(7)进入下一同化窗口。
TracersTracker区域高分辨率CO通量反演技术小结
基于TracersTracker系统进行CO通量反演的具体步骤，具体包括六个步骤：
(1)基础数据处理。
(2)计算CO背景浓度。
(3)获取CO通量扰动集合。
(4)获取CO浓度模拟样本集合。
(5)最优通量求解。
(6)获取下一同化窗口初始场。
4.TracersTracker区域高分辨率CO2通量反演技术
总体而言，CO2/CO比例数据可以分为两类：(1)排放过程中释放的CO2和CO的比例数据(CO2/CO排放比)；(2)地表大气中的CO2和CO浓度比例数据(CO2/CO浓度比)。CO2/CO排放比数据获取较为简单，可直接从已发布的排放清单中得到，本发明采用的CO2/CO排放比数据从MIX排放清单中提取。相较于CO2/CO排放比，CO2/CO浓度比提取过程涉及的数据和流程较多。由于地表大气中CO2、CO观测浓度均包含背景浓度，在本文中，CO2/CO浓度比实际为增量浓度比(ΔCO2/ΔCO)，即去除背景浓度后的比值，求取增量浓度比的具体步骤如下：
(1)获取地表CO2和CO浓度观测数据。在本发明中，CO2观测数据采用OCO-2近地表混合比浓度，CO观测数据采用MOPITT地表混合比浓度。由于CO和CO2观测卫星过境时间和路径不同，利用其观测数据求取CO2/CO浓度比时需要进行时间和空间匹配。
(2)在两个研究区域，以上甸子和瓦里关本底站CO2和CO月均观测作为背景浓度，分别求取两个区域内CO2和CO月均增量浓度ΔCO2和ΔCO，继而可得月均ΔCO2/ΔCO浓度比。由于OCO-2和MOPITT空间分辨率不同，在求取ΔCO2/ΔCO浓度比过程中需进行空间匹配。OCO-2数据分辨率较高(2.25km×1.29km)，MOPITT数据分辨率较低(22km×22km)，在相同观测时间内MOPITT单个网格包含多个OCO-2观测值，求取其平均值作为该网格点的CO2观测浓度。
基于上述两种CO2/CO比例数据及当前较为丰富的CO浓度监测和排放清单统计数据，可在很大程度上丰富“CO2观测数据”，克服CO2观测数据缺乏导致无法进行区域CO2通量反演的不足。在本章中，提出了两种基于CO2和CO伴生比例关系对CO2通量进行反演的方法，即基于排放端清单统计数据的反演方法和基于大气端浓度观测数据的反演方法，为方便起见，将两种方法分别简称为“排放端法”和“大气端法”。在本章中，基于“排放端法”和“大气端法”，实现了上甸子区域和瓦里关区域内人为CO2通量的反演优化，对反演结果进行误差分析并讨论了反演过程中的不确定性。
4.1排放端法
在该方法中，首先从MIX排放清单中提取CO2/CO排放比(Rations of CO2 to CO inthe emission inventory，以下简称Remi)，然后在已优化的CO后验通量基础上，按照Remi比例计算得到CO2人为通量。具体步骤如下(图28)：
(1)准备基础数据。
CO2通量反演所需数据包括：边界场、初始场、气象场及CO2先验通量。CO2边界场来自嵌套模拟(图30)，外层模拟为内层模拟提供边界场。初始场主要来自spin-up模拟，在正式开始同化反演前，利用CMAQ模型模拟一段时间(2015年12月1日-2015年12月31日)CO2浓度，使内层CO2浓度达到基本稳定，并以最后的CO2模拟状态作为初始场。气象场来自FNL再分析数据，利用WRF模型将其转换为同化系统支持的格式。CO2先验人为通量来自清华大学MEIC团队提供的MIX排放清单，CO2先验自然通量来自MOD17A2陆地四级标准数据。
(2)获取CO2/CO比例Remi。
本文采用的Remi来自MIX排放清单，该清单中包含覆盖整个中国区域的CO2、CO通量，从中提取两个研究区域内的CO2、CO通量计算可得Remi。
(3)反演CO人为后验通量。
在当前同化窗口中，利用TracersTracker系统同化两个研究区域内的CO地面浓度观测数据，反演得到优化后的人为CO后验通量。该操作具体步骤与第5章相同，人为CO后验通量与第5章反演结果也相同。
(4)获取人为CO2后验通量。
根据反演优化后的人为CO后验通量和Remi，计算人为CO2后验通量，如下：
CO2(anthro)＝CO(inv)×Remi
式中，CO2(anthro)代表人为CO2后验通量，CO(inv)为反演得到的人为CO后验通量，Remi是(2)中提取的CO2/CO排放比。
(5)进入下一同化窗口。
4.2大气端法
在该方法中，首先基于本底站、OCO-2、MOPITT中的CO2和CO观测数据计算ΔCO2/ΔCO浓度比(Rations ofΔCO2 toΔCO in the atmosphere，以下简称Ratm)。然后依据Ratm及CO站点(上甸子区域22个站点，瓦里关区域9个站点)观测数据，计算得到“CO2观测数据”。最后利用TracersTracker系统同化“CO2观测数据”反演得到CO2总通量。具体步骤如下：
(1)准备基础数据。
本方法所需的边界场、初始场、背景场、气象场及CO2先验通量与排放端法相同。
(2)获取ΔCO2/ΔCO浓度比Ratm。具体步骤上面已详细说明，此处不再赘述。
(3)获取“CO2观测数据”。
根据两个研究区域内的Ratm和CO地基观测及其背景浓度，计算得到“CO2观测数据”。在本文中，上甸子区域内有22个CO观测站点，瓦里关区域内有9个CO观测站点，最终在两个区域内可分别获取对应数量的“CO2观测数据”。
(4)获取CO2通量扰动集合。
在CO2先验总通量的基础上，随机添加30％的扰动，产生指定数目的CO2通量扰动样本。若同化系统是首次运行，则扰动是基于CO2先验通量进行的。若同化系统不是首次运行，则扰动是对上一同化窗口获取的优化后的CO2通量进行的。
(5)获取CO2浓度模拟样本。
基于(4)产生的CO2通量扰动样本，积分CMAQ模型模拟CO2浓度。CMAQ积分运算的次数与通量扰动样本数目一致，最终可得到同等数目下的CO2浓度模拟样本。
(6)最优求解。
读取(5)中基于通量样本获取的CO2模拟样本集合，构建模拟误差矩阵并进行最优通量求解，求解原理是根据四维变分同化方法中的代价函数，结合(3)中计算获取的“CO2观测数据”，在统计误差最小条件下得到p个基向量及其系数，
最后根据获取的基向量及其系数得到CO2最优通量。根据得到的CO2后验通量和CO2先验自然通量，可得到人为CO2后验通量。
CO2(anthro)＝CO2(inv)-CO2(nat)
式中，CO2(anthro)代表人为CO2后验通量，CO2(inv)代表大气端法反演的CO2后验总通量，CO2(nat)代表CO2自然通量。
(6)获取下一同化窗口初始场。
基于获取的CO2最优通量数据，重新积分CMAQ模型，可获取该同化窗口内的CO2最优模拟值，该模拟值可作为下一同化窗口初始时刻的初始场。
(7)进入下一同化窗口。
本发明主要创新点：
(1)原创性研发了区域高分辨率同化反演系统TracersTracker。
当前，同化反演系统主要集中在全球和大洲等较大尺度，区域高分辨率同化反演系统较少，缺乏通量反演的精细化研究。本文实现了POD4DVar高效同化算法和区域大气传输模型CMAQ的耦合，研发了区域高分辨率碳同化反演系统TracersTracker。该系统操作简便且无需开发伴随模式，具有较高的运算效率。系统同时支持CO和CO2通量的反演和优化，可有效提高反演通量的空间分辨率，将反演通量空间分辨率提升至千米级。
(2)原创性研发了两种CO2-CO联合同化方法，获取了区域高分辨率CO2通量并对人为通量和自然通量进行了定量区分。
针对当前CO2观测匮乏导致区域小尺度CO2通量反演无法实施的问题，本文提出了基于CO2/CO排放比的CO2通量反演方法和基于CO2/CO浓度比的CO2通量反演方法。前者利用当前已知的排放清单中的CO2/CO伴生比例关系，将CO2/CO排放比作为变量因子参与到同化过程中，最终实现高分辨率(3km)CO2通量的反演和优化。后者从当前较为有限的CO2、CO卫星及大气本底站观测数据中提取地表CO2/CO月均浓度比，依据此浓度比及当前较为丰富的CO地基观测数据推算“CO2观测数据”，利用TracersTracker系统同化“CO2观测数据”反演得到高分辨率(3km)CO2通量。
步骤1设置实施例应用区域，构建区域碳同化系统TacersTracker
2.1设置实施例区域
无论是CO还是CO2，其排放量主要受人为因素(人口、经济水平等)和自然因素(植被种类、碳密度等)影响，本文选取了地理条件、经济发展水平差异较大的两块区域，即瓦里关区域(以瓦里关大气本底站为中心)和上甸子区域(以上甸子大气本底站为中心)，以分析不同人为和自然条件下碳通量的时空分布格局及演变特征。大气本底站是由中国政府与世界气象组织(World Meteorological Organization，WMO)和环境基金(GlobalEnvironment Fund，GEF)合作开展的对全球尺度大气本底污染物浓度进行监测的国际合作项目，可以长期、稳定、连续地获取全球大气本底监测数据，为深入理解温室气体变化规律提供基础数据。
瓦里关区域位于中国西部地区(图30)，以瓦里关本底观测站为中心，长、宽各约270千米，东起东经100.432度，西至东经103.448度，北至北纬37.107度，南至北纬34.781度，主要包含甘肃和青海两省的部分区域，即青海省西宁市、海北藏族自治州、海南藏族自治州、海东地区、黄南藏族自治州以及甘肃省兰州市、临夏回族自治州。该区域属于温带季风气候、高山气候、温带大陆气候的交汇地带，平均海拔3800米，该区域大部分地区为干旱、半干旱荒漠草原及沙洲，平均年降水量370毫米，主要集中在5-9月。瓦里关区域年平均气温零下0.9℃，年最高气温28.3℃，年最低气温-23.3℃，冬寒夏凉、温差较大。据2010年第六次全国人口普查数据，该区域常住人口1082.92万人，人口密度较小。该区域属于经济欠发达地区，工业基础较为薄弱，重工业及污染较为严重的生产、生活活动较少。
上甸子区域位于华北平原北部(图30)，以上甸子本底观测站为中心，长、宽各约260千米，东起东经115.173度，西至东经118.511度，南至北纬38.757度，北至北纬41.273度，主要包含河北省承德市、保定市、廊坊市、唐山市以及北京市、天津市。该区域属于温带季风气候，年降水量400-600毫米，地形主要为平原和山地，是华北平原向蒙古高原的过渡地带，地势西北高、东南低。上甸子区域包含中国经济最发达的工业区之一——京津唐工业区。据2010年第六次全国人口普查数据，该区域常住人口超过7000万，是全国人口密度最大的几个区域之一。该区域工业基础较好、体系门类齐全，石油工业、煤化工业、海洋化工、机械电子工业、冶金工业等都较为发达，是我国北方最大工业聚集区，每年会产生大量的工业废物、废气，大气环境较差。
中国国土面积大，不同地区经济水平、人口密度、气候、地理条件各异，CO和CO2排放强度及变化特征也有较大区别，本文选取的两个研究区域，可以在一定程度上代表上述差异。在本文中，利用构建的区域高分辨率碳同化反演系统TracerTracker，对上甸子和瓦里关区域的CO和CO2通量进行反演，总结不同人文、地理条件下碳排放的时空变化特性。同化系统在进行碳通量反演时采取嵌套模式(图30)，其中d01为外层区域，d02为内层区域，本文主要对内层区域进行碳通量反演优化，外层区域负责为内层区域提供边界场。在同化系统内，模型模拟气体浓度是以网格化的方式进行的，其中上甸子外层区域为78行×78列网格，内层区域为91行×91列网格，瓦里关外层区域为69行×69列网格，内层区域为85行×88列网格。外层区域空间分辨率为27km，内层区域空间分辨率为3km。
一般而言，碳同化反演系统由同化反演算法、大气传输模型、观测数据等部分组成。其中观测数据是进行通量反演优化的核心依据，下面对本文中所使用的数据做详细说明。
步骤2数据收集
2.1 CO观测
表2-1瓦里关区域CO观测站点及验证站点
CO地基观测数据来自环保部门建立的空气质量监测网，该监测网包含国家、省、市、县四级5000余个空气质量监测站点，其中国家级监测站点1436个。各监测站点数据实时对外发布，数据可从全国城市空气质量实时发布平台(http://106.37.208.233:20035/)获取。该平台提供符合GB3095-2012标准的PM2.5，PM10，SO2，NO2，O3，CO六种污染物的逐小时浓度观测数据。在本文中，采集了2016年全年(共366天8784小时)的CO观测数据，即在没有数据缺失的情况下，每个站点每年有8784个观测数据，主要用于CO通量的反演优化。本文中的两个研究区域包含的CO观测站点数量差异较大(图31)，瓦里关区域共有9个观测站点(表2-1)，上甸子区域共有22个观测站点(表2-2)。为方便对反演结果进行交叉验证，分别从两个研究区域选择若干观测站点作为验证站点，其中瓦里关区域2个站点(表2-1)，上甸子区域5个验证站点(表2-2)。
CO卫星观测数据来自美国国家大气研究中心(National Center forAtmospheric Research，NCAR)，由MOPITT扫描仪观测数据解析获得。MOPITT扫描仪搭载在Terra卫星上，该卫星于1999年发射升空，对大气中CO和CH4进行观测。MOPITT可提供CO和CH4的柱总量、6层(850hPa，700hPa，500hPa，350hPa，250hPa，150hPa)及地表CO混合比浓度数据。MOPITT的CO观测数据的水平分辨率为22km×22km，垂直分辨率为3km，经验证精度优于10％。本文采用的是MOPITT的L2级V6版本数据，该数据经过定标和校准，剔除了云污染和天气影响。
表2-2上甸子区域CO观测站点及验证站点
2.2 CO2观测
本文采用的CO2观测数据分为两个部分：地基观测和卫星观测。地基CO2观测数据主要由世界气象组织(World Meteorological Organization，WMO)设立的全球大气观测网(Global Atmosphere Watch Programme，GAW)提供。这种地基实时CO2浓度观测数据具有较高精度、可靠性也较强，但是，地表CO2观测站点稀少，在全球仅设置有200多个大气本底观测站点，并且分布也很不均匀。中国境内仅有5个大气本底观测站(表2-3)并且这些观测数据并未全部对外公开。本文使用的CO2观测数据来自上甸子站点和瓦里关站点，该数据可从世界温室气体数据中心(World Data Centre for Greenhouse Gases，https://gaw.kishou.go.jp/，WDCGG)获取。两个站点的CO2观测作为背景浓度，主要用来获取地表大气中ΔCO2/ΔCO浓度比信息。
表2-3中国境内大气观测本底站
CO2卫星观测数据来自OCO-2(Orbiting Carbon Observatory-2，OCO-2)卫星。NASA(National Aeronautics and Space Administration)计划于2009年发射轨道碳观测者用于获取CO2浓度信息，但发射计划失败。2014年再次发射轨道碳观测者2号(OCO-2)并获得成功，该卫星是全球第二颗碳观测卫星，主要任务是获取地球表面碳源碳汇信息。OCO-2卫星距地高度为705km，每天绕地球14.65圈，每16天就可获取覆盖全球的观测数据，空间分辨率为2.25km×1.29km。OCO-2数据以最优估计理论为基础，采用数据同化方法，对全球CO2浓度数据进行反演和优化，产生了L1-L3等多级数据产品。由于云、雨和气溶胶等天气因素的影响，OCO-2的一部分观测数据质量较低、无法满足精度要求，这部分数据需要剔除。另外，在使用OCO-2卫星数据前，应按照官方发布的操作手册对数据进行质量控制。
OCO-2卫星观测数据可从NASA官方网站(https://search.earthdata.nasa.gov/)下载。本文采用的CO2观测数据是OCO-2的L2版本，观测时间从2016年1月1日到2016年12月31日。另外，由于天气及仪器的问题，部分观测数据是无效的，在CO2卫星观测数据参与到同化反演之前，对其进行了质量控制。
2.3先验通量
CO2和CO人为通量来源于清华大学MEIC(Multi-resolution Emission Inventoryfor China)团队发布的MIX排放清单。MIX清单可以为气象模型和大气化学模型提供最新的、最精确的地表排放通量数据。该清单采用集合开发方法，综合考虑了REAS2、CAPSS、PKU-NH3、MEIC等发布的最新排放清单数据，提供了以2008和2010年为基准年的0.25度分辨率网格化排放清单。目前，MIX清单数据已在许多项目中得到应用，如联合国半球大气污染传输计划(HTAP)、东亚模式比较计划第三期(MICS-AsiaⅢ)等。
表2-4上甸子区域和瓦里关区域CO人为先验通量
本文采用2010年MIX清单v1.0版本，该清单内包含东亚地区SO2、NOx、CO、NH3、NMVOC、PM10、PM2.5、BC、OC、CO2等污染物或温室气体通量数据。MIX清单数据格式为NetCDF，每个文件都包含五个部门的排放数据，即电力行业、工业部门、民用领域、交通领域、农业领域(仅针对NH3)。MIX数据的横向、纵向分辨率都为0.25度，数据覆盖范围四至为:西至东经40.125度、东至东经179.875度,、南至南纬20.125度、北至北纬89.875度。各污染物排放数据文件含有相同的维度，均为560(列)×441(行)×12(月)。MIX网格化数据的单位是不同的，主要分为两种：摩尔/月/网格(SAPRC-99机制和CB05机制下的物种)和吨/月/网格(其他物种)。图32、图33分别为从MIX清单数据中提取的上甸子区域和瓦里关区域的CO先验通量数据。由于在本文的两个研究区域中CO自然源汇很小，在进行同化反演过程中并未考虑此部分，即CO自然通量设置为零。
CO2先验通量主要包括两个部分：自然先验通量和人为先验通量。图34、图35为从MIX清单数据中提取的上甸子区域和瓦里关区域的CO2人为先验通量数据。本文采用的自然先验通量主要来自EOS(Earth Observation System)观测计划的观测数据。EOS由美国国家航空航天局主持实施，是一系列对地观测卫星的简称，Terra和Aqua是其中较为重要的两颗卫星。Terra卫星和Aqua卫星分别发射
表2-5上甸子区域和瓦里关区域CO2人为先验通量
表2-6上甸子区域和瓦里关区域CO2自然先验通量
于1999年和2002年，两颗卫星分别搭载了多种用途广泛的传感器，其中两颗卫星都搭载了中分辨率成像光谱仪MODIS(Moderate Resolution Imaging Spectradiometer)。MODIS用于捕捉数据的波段多达36个，覆盖范围从可见光波段到红外波段，且观测数据更新速度快，平均1-2天即可获取覆盖全球的观测数据。MODIS获取的数据类型也非常广泛，涉及陆地、大气、海洋等，其空间分辨率为250m、500m或1000m，视场宽度为2330km。MODIS标准数据根据其具体内容分为0级、1级数据产品，在1级数据基础上，经过一系列校正和处理过程，最终形成2-4级数据产品，主要分为三种：大气标准数据产品、陆地标准数据产品和海洋标准数据产品。
自然CO2先验通量(图36，图37)从MODIS的MOD17A2产品中提取，MOD17A2产品属于陆地4级标准数据，主要用来获取植被生产力等信息。该数据的空间分辨率为500m，正弦投影，时间分辨率为8天，是8天观测数据的合成产品。该产品可作为数据模型的输入，用于计算陆地能源、碳、水循环过程和植被的生物地球化学。MOD17A2产品中包括关毛初级生产力(GPP)和净光合作用(PSN)等信息，其中PSN为GPP减去呼吸消耗。为与本文的碳同化反演系统保持一致，该数据参与同化反演之前，利用NOAA提供的MRT软件将其转换为3km分辨率的网格数据。
2.4气象数据
气象场数据采用美国国家环境预报中心(National Centers for EnvironmentalPrediction，NCEP)发布的气象全球再分析数据(Final Operational Global AnalysisData，FNL)，该数据采用当前最先进的数据同化技术，对多源观测数据(地基、船舶、飞机、卫星等)进行了质量控制和数据同化处理。FNL再分析数据具有密度大、时次多、分辨率较高、连续性强、内容丰富等特点，数据内含有气温、地面气压、相对湿度、位势高度、海平面气压、海表面温度、土壤参数、风场数据、涡度等27个物理量。FNL再分析数据采用sigma坐标系，垂直方向包括地面层及1000hPa-10hPa共27层，水平分辨率为1°×1°，时间间隔为6h，包含4个时次(世界时0、6、12、18时)的数据。
由于采用的观测数据和处理标准不同，FNL存档数据在不同时间的编码格式及内部参数有所不同，按照存储时间分为三种，具体如下：
(1)1967年7月1日00时到1997年4月12日12时。其编码格式为GRIB1或ON84，空间分辨率为2.5°×2.5°，时间分辨率为12h，大气分为12层。
(2)1997年4月1号00时到2007年6月30号12时。其编码格式为GRIB1，大气分为16层，其他参数与第一种相同。
(3)1999年7月30号18时至今的数据。本文使用的是该时间段内的数据，其编码格式为GRIB1，空间分辨率为1.0°×1.0°，时间分辨率为6h，大气分为26层。
本文采用的为第(3)种FNL数据，时间为2016年1月1日-2016年12月31日，使用前利用WRF模型将其转换为同化系统兼容格式。
步骤3实施例区域模型运行与结果分析
1.TracersTracker区域高分辨率CO通量反演实施例
1.1引言(Preface)
CO是大气中含碳量较高的气体，仅次于CH4和CO2，是全球碳循环的重要中介气体。在适当的气候条件下，CO可与大气中OH自由基等发生化学反应产生CO2，间接影响温室气体浓度。
在TracersTracker区域高分辨率CO通量反演技术中，详细阐述了本发明构建的同化反演系统TracersTracker的核心算法和基本架构，同时验证了该系统在CO通量反演中的有效性。结果表明，TracersTracker系统在CO通量反演中表现较好，基本消除了先验通量中的“误差”，明显提高了后验模拟的精度。但是，上述对同化系统的有效性验证是在理想状态下进行的，在这种状态下“CO观测数据”非常充足且没有“误差”，可以极好的满足同化系统对观测信息的需求。而在实际CO通量反演研究中，CO观测数据是有限的，并且由于仪器、观测手段及其他因素的影响，CO观测数据是存在不同程度误差的。
这里利用真实CO地表站点观测数据，在TracersTracker同化系统支持下，在上甸子区域和瓦里关区域进行了CO通量反演，并对CO后验模拟的准确性进行了验证，同时从季节性变化角度讨论了CO后验通量的变化特征。
1.2数据准备(Data)
1.2.1先验通量
CO先验通量采用2010年MIX排放清单1.0版本，数据单位为摩尔/月/网格，空间分辨率为0.25×0.25度，使用前需利用SMOKE模型将其处理成符合TracersTracker系统要求的网格数据。在上甸子区域，外层先验通量大小为78行×78列，空间分辨率为27km×27km，内层先验通量大小为91行×91列，空间分辨率为3km×3km。在上甸子区域，外层先验通量大小为69行×69列，空间分辨率为27km×27km，内层先验通量大小为85行×88列，空间分辨率为3km×3km。处理后先验通量单位为摩尔/秒/网格。
1.2.2气象场
气象场数据采用美国国家环境预报中心发布的FNL全球再分析数据(GRIB1格式)，同化前需利用WRF模型将其处理成符合系统要求的网格数据。在上甸子区域，外层气象场大小为78行×78列，空间分辨率为27km×27km，内层气象场大小为91行×91列，空间分辨率为3km×3km。在瓦里关区域，外层气象场大小为69行×69列，空间分辨率为27km×27km，内层气象场大小为85行×88列，空间分辨率为3km×3km。
1.2.3 CO观测
CO观测来自环保部地基观测站点，在上甸子区域，共22个CO站点的逐小时观测数据参与同化，观测时间从2016年1月1日到2016年12月31日，CO观测站点分别位于承德市、保定市、廊坊市、唐山市、北京市、天津市(图30)。在瓦里关区域，共9个CO站点的逐小时观测数据参与同化，观测时间从2016年1月1日到2016年12月31日，CO观测站点分别位于西宁市、海北藏族自治州、海南藏族自治州、海东地区、黄南藏族自治州、兰州市、临夏回族自治州(图30)。CO观测数据参与同化前，对其进行了质量控制，具体包括：
(1)剔除“缺测”较多站点。由于仪器故障等问题，CO观测值会出现“缺测”的情况，导致观测数据不连续。在进行同化反演之前，首先剔除“缺测”数据较多的站点。
(2)剔除“异常”观测数据。由于仪器异常和操作失误等原因，观测数据会出现远大于或远小于平均值的现象，将其定义为“异常”数据，在进行同化反演之前，剔除这些“异常”数据。
1.3 TracersTracker区域高分辨率CO通量反演系统运行
见技术方案部分：3.TracersTracker区域高分辨率CO通量反演技术
1.4结果分析
1.4.1背景浓度结果
基于1.3所述的背景浓度求取方法，提取了两个研究区域内各站点的背景浓度，图36、图37分别为上甸子区域和瓦里关区域内部分观测站点背景浓度提取结果。本文从气候学和气象学角度，将每年分为四个季节，3月-5月为春季，6月-8月为夏季，9月-11月为秋季，12月-2月为冬季。
上甸子区域和瓦里关区域经济、地理条件不同，其CO整体浓度水平和变化特征也存在明显差异。整体而言，瓦里关区域的CO浓度要明显远小于上甸子区域(表5-1)，瓦里关区域的CO年均浓度为1.09毫克/立方米，上甸子区域CO年均浓度为1.43毫克/立方米，是瓦里关区域的1.3倍。1月-12月，上甸子区域CO浓度均比瓦里关区域高，但在12月和1月两者差距最大。但从季节性变化而言，两个区域CO浓度变化趋势是相似的，即在夏、秋两季浓度普遍较低，在冬、春两季浓度明显升高。
表5-1上甸子区域和瓦里关区域CO观测及背景浓度
CO背景浓度与站点观测浓度具有高度相关性，其变化趋势与CO浓度观测数据基本一致。上甸子区域背景浓度年均值为1.08毫克/立方米，瓦里关区域背景浓度年均值为0.69毫克/立方米，明显低于上甸子区域。上甸子区域在夏、秋两季的背景浓度为瓦里关区域的1.53倍，冬、春两季是瓦里关区域的1.55倍，就CO增量浓度而言，上甸子区域所占比例更高，其主要原因是上甸子区域经济较为发达，人为生产、生活活动较为活跃，CO排放量较大。
同时可以看到，在同一区域内，CO背景浓度在空间上的变化也是极不平衡的。在上甸子区域，唐山、廊坊、天津的CO背景浓度要远高于其他区域，是排放CO的主要区域。在瓦里关区域，西宁、临夏的CO背景浓度要远高于其他区域，甚至是其他区域的几倍。这些区域主要是工业和人类聚居地，人类活动相对活跃，CO排放量大(图38、图39)。
虽然两个研究区域内CO浓度水平差距较大，但在某些重污染时段(如1月和12月)中，CO浓度均能达到较高水平。这说明即使在CO排放强度不强的区域，不利的气象、地理条件同样可以造成较为严重的污染事件。另外，无论在上甸子区域还是瓦里关区域，CO浓度均在春、冬两季要高于其他两季，主要原因是上甸子区域和瓦里关区域都处于中国北方，在春、冬两季气温较低，居民供暖、出行等活动需要消耗大量煤炭等化石能源，排放的CO明显增多，短时间内推高了大气中CO的浓度。但相较而言，瓦里关区域CO浓度在四季的差距明显小于上甸子区域，说明此区域内春、冬季节供暖等活动强度的提升仅限于少数地区，对整个区域影响有限。整体上，上甸子区域内CO浓度变化较为剧烈，污染较为严重的时段较多，虽然这些时段在春、冬两季较多但不仅限于春、冬两季，在夏、秋两季也会出现较为严重的污染事件。
1.4.2通量反演结果
利用TracersTracker系统同化了两个研究区域的CO观测数据，反演得到了上甸子区域和瓦里关区域的CO后验通量。结果表明(表5-2、表5-3)，CO先验通量在两个区域都存在一定程度的低估，CO后验通量在大多数时段都要比先验通量高。
在上甸子区域，CO后验通量在1月-12月较先验通量分别提高了39.4％、34.3％、34.5％、33.1％、29.0％、28.5％、32.1％、26.3％、31.0％、27.3％、37.5％和39.9％，全年增幅为32.7％。从月均通量变化上看，后验通量增加幅度较大的为1月和12月，增幅较小的为8月和10月。
表5-2上甸子区域CO后验通量增量及增比
从上甸子区域年均通量上看，CO先验通量年均值为1.14摩尔/秒，其中最小值为0.87摩尔/秒，出现在9月，最大值为1.65摩尔/秒，出现在1月。CO后验通量年均为1.62摩尔/秒，约为先验通量的1.3倍，其中最小值为1.12摩尔/秒，出现在9月，最大值为2.24摩尔/秒，出现在1月。整体而言，后验通量的变化趋势与先验通量基本保持一致。
表5-3瓦里关区域CO后验通量增量及增量比例
在瓦里关区域，CO后验通量在大多数时段也比先验通量有所增加，但增幅要远大于上甸子区域。后验通量在1月-12月中较先验通量分别增加了68.0％、66.0％、64.0％、61.4％、59.8％、59.2％、58.7％、58.7％、59.6％、60.9％、64.4％和66.6％，整体增幅约为上甸子区域的1.6倍，其主要原因是瓦里关区域CO先验通量整体有较大幅度偏低。虽然整体增幅较高，但从绝对增量来看，瓦里关区域CO后验通量的净增量(0.05摩尔/秒)要远小于上甸子区域(0.33摩尔/秒)，仅为上甸子区域的七分之一左右。
从瓦里关区域年均通量上看，CO先验通量均值为0.09摩尔/秒，其中最小值为0.05摩尔/秒，出现在5月-9月，最大值为0.15摩尔/秒，出现在12月和1月。CO后验通量均值为0.14摩尔/秒，大约为先验通量的1.6倍，其中最小值为0.07摩尔/秒左右，出现在5-9月，最大值为0.24摩尔/秒，出现在1月和12月，后验通量整体趋势与先验通量也基本保持一致。
从后验通量变化的空间分布特征看，后验通量的增量和增幅具有不同表现。CO后验通量的增量主要集中于城市区域或人口密集区，在上甸子区域中通量增加最多区域在唐山市境内(以唐山市主城区为主)，其次为天津地区和北京地区，承德市区域内通量增加较少。在瓦里关区域，CO后验通量增加较为明显区域主要位于西宁市主城区附近，其次是临夏市城区及周边范围内，其他区域如海北藏族自治州、海南藏族自治州、海东地区等增量较小。无论是在上甸子区域还是瓦里关区域，CO通量增量主要集中在城市主城区及附近区域。由于城区及附近区域集中了大部分的生产、生活活动，所以这些区域内的CO先验通量比其他区域要大的多，在相同增加比例下，这些区域的CO增量也较大。从后验通量增加比例看，两个区域表现不同。在上甸子区域，增幅在空间分布上并无明显趋势。在瓦里关区域，虽然增量较大的区域位于城区及其近郊，但非城区的增幅要明显高于城区及其近郊，其主要原因是瓦里关区域非城区的先验通量较低所致，也是瓦里关区域后验通量增幅明显高于上甸子区域的原因。上述结果表明，不论是在上甸子区域还是瓦里关区域，CO先验通量经过TracersTracker校正、优化后都有较明显的变化，但整体上的通量水平及变化趋势与先验通量仍基本保持一致。
1.4.3模拟误差分析
为检验CO后验通量的准确性，本节利用CO验证站点的观测数据对后验模拟进行了验证。图40、图41分别为瓦里关区域和上甸子区域CO先验模拟、后验模拟及验证站点观测数据。上甸子区域和瓦里关区域CO浓度具有相同变化趋势，即在1月-7月-12月CO浓度为高-低-高，CO先验模拟和后验模拟都可较好的把握这个趋势。但是，CO先验模拟与观测数据之间存在较大误差，尤其是在CO浓度较高的时段，如1月和12月。
在上甸子区域，CO先验模拟与观测数据在北京定陵(BJDL)、保定监测站(BDJCZ)、廊坊药材公司(LFYCGS)、天津空港物流加工园(TJKG)、唐山十二中(TSSEZ)五个验证站点的相关系数分别为0.58、0.66、0.65、0.61和0.52，相关系数较差且存在明显误差，这误差主要是由先验通量中存在明显的不同程度的低估造成的。而CO后验模拟的精度有了较明显的提高，在很大程度上减小了模拟与观测之间的误差，在北京定陵、保定监测站、廊坊药材公司、天津空港物流加工园、唐山十二中五个站点上的相关系数分别提高到了0.82、0.81、0.79、0.78和0.74，CO后验模拟与观测数据的相关系数提升明显。
表5-4 CO先验、后验模拟误差
在瓦里关区域，CO先验模拟与观测数据在西宁监测站(XNJCZ)、临夏环保局(LXHBJ)两个站点上的相关系数分别为0.69和0.64，相关系数比上甸子区域略高，但仍存在较明显的误差。经同化反演优化后，CO后验模拟和观测数据的相关系数分别提高到了0.76和0.76，模拟值与观测值在变化趋势上更加契合。
另外，无论是在上甸子区域还是瓦里关区域，CO后验模拟的误差也出现了明显了下降。在上甸子区域，绝对误差在北京定陵、保定监测站、廊坊药材公司、天津空港物流加工园、唐山十二中观测站点上分别减小至0.34毫克/立方米、0.61毫克/立方米、0.56毫克/立方米、0.45毫克/立方米、1.01毫克/立方米，整体减幅达到31％。在瓦里关区域，绝对误差在西宁监测站、临夏环保局站点上分别减小至0.53毫克/立方米、0.35毫克/立方米，整体减幅达到35％，两个区域的后验模拟误差明显下降。
从CO先验和后验模拟误差分布区间看(图42、图43、图44、图45、图46、图47、图48)，各个站点模拟误差的分布既有相同性也有差异性。整体而言，误差的分布区间在各个站点都呈缩小趋势，误差更为集中(向0值方向)。在上甸子区域，北京定陵、保定监测站、廊坊药材公司、天津空港物流加工园、唐山十二中五个站点的误差均值分别减小至0.07毫克/立方、0.05毫克/立方、0.12毫克/立方、0.06毫克/立方、0.14毫克/立方米。在瓦里关区域，西宁监测站、临夏环保局两个站点的误差均值分别减小至0.06毫克/立方米、0.05毫克/立方米。虽然在各站点的模拟误差均值变化趋势总体上相似，但具体到每个站点，误差变化趋势也不尽相同。有的收缩幅度大(如唐山十二中)，有的收缩幅度小(如北京定陵)，有些误差收缩不太明显(如西宁观测站)，这主要与各站点所处位置及CO观测浓度本身的大小有关。整体上，无论是在上甸子区域还是瓦里关区域，TracersTracker系统很大程度上校正了CO先验通量明显偏低的问题，使后验模拟误差整体上向0趋近。
上述结果表明，TracersTracker系统在实际CO通量反演应用中表现良好，可以在很大程度上校正CO先验通量中的“误差”，有效提高CO通量的准确性。CO后验模拟的误差明显减小，模拟精度有较大提升，模拟结果与实际观测保持了较好的一致性。但同时可以看到，CO通量反演结果在很大程度上依赖于CO观测数据，观测数据较为充足的区域其通量数据可以得到有效且较为准确的优化，观测数据较为匮乏的区域通量的优化在很大程度上依赖其周边区域的观测信息，这在一定程度上会提高通量反演的不确定性。
1.4.4季节性变化分析
在中国北方，人为CO排放有明显的季节性变化，如冬季的CO排放强度是明显高于夏季的。在CO后验通量中，不同季节CO增量和增幅也有明显差异。在本节中，从季节性变化角度分析了CO后验通量的时空变化特征及其后验模拟精度的不确定性。
从后验通量增量看(表5-5)，在上甸子区域，四个季节中冬季增量最大，达到0.43摩尔/秒，其他三个季节无明显差异，分别为0.31摩尔/秒、0.31摩尔/秒和0.27摩尔/秒，只有冬季增量的三分之二左右。瓦里关区域和上甸子区域类似，CO后验通量在冬季的增量最大，达到0.12摩尔/秒，春季和夏季通量增量最小，为0.03摩尔/秒，只有冬季增量的四分之一左右。
表5-5不同季节CO后验通量
在上甸子区域，不同季节CO后验通量增幅是不同的，冬季增幅最高，达到37.8％，夏季增幅最低，仅为28.9％，春、秋两季增幅居中，约为32％(表5-5)，CO净增量与增幅表现完全一致。在瓦里关区域，所有季节增幅都要远高于上甸子区域，但整体变化趋势与上甸子区域基本一致。瓦里关区域冬季增加比例最大，达到66.8％，略高于其他季节，夏季增幅最低，只有59％左右，春、秋两季增幅基本持平，约为62％。
整体而言，CO后验通量中变化较大的季节主要是冬季，冬季居民供暖活动集中，短时间内推高了CO的排放量。另外，随着居民生活水平的提高，冬季取暖对燃料的需求的增长速度明显高于其他季节，CO排放量在冬季增长比例较其他季节也明显偏高。
在两个研究区域，不同季节先验模拟和后验模拟的表现有相似之处，也存在差异(表5-6)。在上甸子区域，CO先验模拟与观测数据在春季、夏季、秋季、冬季的相关系数分别为0.55、0.62、0.61和0.61，后验模拟的相关系数分别0.72、0.71、0.78和0.84，后验模拟与观测数据的契合程度有明显提高。在瓦里关区域，先验模拟与观测数据在春季、夏季、秋季、冬季的相关系数分别为0.55、0.58、0.62和0.63，后验模拟的相关系数则分别提升至0.69、0.74、0.78和0.79。整体上，无论是在上甸子区域还是瓦里关区域，后验模拟与观测数据的相关系数都有明显的提升，上甸子区域提升幅度更大。另外，无论是先验模拟还是后验模拟，秋、冬两季CO模拟的变化趋势与观测数据具有更好的一致性。同时可以看到，无论是在上甸子区域还是瓦里关区域，后验模拟精度与先验通量具有一定联系，先验通量会在一定程度上影响后验模拟，其主要原因时后验通量的获取是以先验通量为基础的。
表5-6季节性CO后验通量增量及增比
先验模拟和后验模拟的绝对误差与相关性变化趋势大致相同。在上甸子区域，春、夏、秋、冬四季先验模拟的绝对误差分别为0.62毫克/立方米、0.59毫克/立方米、0.72毫克/立方米、1.17毫克/立方米，后验模拟的绝对误差分别减小至0.49毫克/立方米、0.45毫克/立方米、0.58毫克/立方米、0.89毫克/立方米，误差有明显减小，尤其是在秋、冬两季。在瓦里关区域，春、夏、秋、冬四季的绝对误差分别为0.42毫克/立方米、0.39毫克/立方米、0.69毫克/立方米、0.84毫克/立方米，在后验模拟中分别减小至0.34毫克/立方米、0.29毫克/立方米、0.54毫克/立方米、0.61毫克/立方米，误差减小趋势与上甸子区域基本一致。
从误差分布上看，不同季节误差分布及其变化趋势有明显不同。总体上，后验模拟误差在不同季节的分布区间都呈缩小趋势。在上甸子区域，春、夏、秋、冬四季的误差均值分别减小至0.10毫克/立方、0.08毫克/立方、0.12毫克/立方和0.07毫克/立方米。在瓦里关区域，春、夏、秋、冬四季的误差均值分别减小至0.08毫克/立方米、0.08毫克/立方米、0.03毫克/立方和0.06毫克/立方。相比先验模拟，两个区域的误差明显向0值集中。
在误差分布区间上，不同季节变化趋势相同但存在差异。在上甸子区域，先验模拟在春、夏、秋、冬四个季节误差分布区间(累积百分比10％到90％)分别为[-1.13,0.57]、[-0.85,0.40]、[-1.35,0.75]、[-2.70,1.19]，后验模拟在春、夏、秋、冬四个季节误差范围分别为[-0.69,0.88]、[-0.47,0.69]、[-0.80,0.99]、[-1.63,1.42]，误差范围在春、夏、秋、冬四季分别减小了8.0％、7.8％、14.8％、21.6％。在瓦里关区域，先验模拟在春、夏、秋、冬四个季节误差范围分别为[-0.86,0.33]、[-0.65,0.27]、[-1.55,0.72]、[-1.79,0.86]，后验模拟在春、夏、秋、冬四个季节误差范围分别为[-0.47,0.51]、[-0.35,0.51]、[-0.96,0.92]、[-1.17,1.13]，误差范围在春、夏、秋、冬四季分别减小了17.7％、6.6％、17.2％、13.3％。无论是在上甸子区域还是瓦里关区域，秋、冬两季误差区间收缩较大，春、夏两季误差区间收缩较小，秋、冬季CO后验模拟精度提升更加明显。
在中国北方的秋、冬季，气温降低会使以石油、煤炭为代表的化石燃料消耗大量增加，进而在短时间内使空气中的CO浓度急剧升高。另外在局部气候变化条件作用下，在较小区域内，可以造成CO的迅速扩散或累积，进而影响空气中CO的浓度。一般情况下，对于变化较为频繁的CO浓度，大气传输模型或空气质量模型在模拟时难度会增加，但从本文的两个研究区域的反演结果看，后验模拟可以较好的反映冬季CO浓度的起伏变化，这说明本文的同化系统较好的优化了CO先验通量中的“偏差”，有效消除了明显的模拟误差，提高了后验模拟精度。
1.5小结
本实施例利用TracersTracker同化反演系统对本文的两个研究区域(上甸子区域和瓦里关区域)的CO通量进行了反演和优化，对反演优化后的后验通量进行了阐述并分析和讨论了后验通量的时空变化特征，主要包括三部分内容：
(1)阐述了TracersTracker系统进行CO通量反演所需的数据，详细说明了数据的处理过程。
(2)对CO通量反演结果进行分析和讨论。首先，展示了CO后验通量反演结果，上甸子区域CO后验通量较先验通量增加了32.7％，瓦里关区域CO后验通量较先验通量增加了62.3％。然后，对比分析了CO观测值、先验模拟值、后验模拟值的误差分布情况。从后验模拟误差分布来看，TracersTracker同化反演系统可以在很大程度上优化先验清单，有效提高模型的CO模拟精度。最后，分析阐述了后验CO通量在不同季节的具体变化情况。结果表明，两个区域的CO后验通量变化情况是不同的。在通量增量来看，在上甸子区域，四个季节中冬季增量最大，达0.55摩尔/秒，其他三个季节增量较为接近无明显差异。瓦里关区域和上甸子区域类似，CO后验通量在冬季的增量最大，达0.08摩尔/秒，夏季通量增幅、增量最小。
2 TracersTracker区域高分辨率CO2通量反演实施例
2.1引言(Preface)
CO2是大气中重要的温室气体，对全球环境变化有着重要的影响，确定CO2排放强度及其时空演变特征具有重要意义。但是，当前的CO2通量反演系统主要集中在全球或大洲尺度上，其反演结果分辨率较低，无法满足区域小尺度碳通量研究的精细化要求。
在区域同化反演系统中，观测数据的质量直接决定着反演过程的可行性和最终反演结果的准确性，较高分辨率的通量反演成果通常需要较高时空分辨率的观测信息，观测数据的不足会给反演结果带来不确定性甚至错误。在理想实验中(第3章)，基于充足的“CO2观测”，同化系统可以获取较高精度的CO2通量，有效提高后验模拟的准确性。在对CO通量进行反演优化时(第5章)，两个研究区域也具有较为充足的CO地基观测站点，并且各站点的观测数据具有较高的时间分辨率，也可以较好的满足区域小尺度范围的CO通量反演研究。
但是，当前CO2地基观测站点稀少，时间分辨率和空间分辨率较低，无法满足区域小尺度范围内的CO2通量反演。虽然碳卫星在一定程度上丰富了CO2观测信息，但碳卫星的观测时间(依赖于卫星过境时间)和观测位置(依赖于卫星过境路径)存在严重的不连续性，远远不能满足较小区域内通量反演对观测数据的需求。
在本章中，研究利用CO2和CO的时空演变相关性，在现有较为丰富的CO浓度观测、排放统计数据的支持下，对上甸子区域和瓦里关区域的CO2通量进行反演优化，对反演得到的CO2通量进行准确性验证和不确定性分析。
2.2数据准备(Data)
2.2.1先验通量
CO2先验通量包括两部分，即人为先验通量和自然先验通量，总通量为两者之和，其空间分辨率为0.25×0.25度。使用前需利用SMOKE模型将CO2总通量处理成符合TracersTracker系统要求的网格数据。在上甸子区域，外层先验通量大小为78行×78列，空间分辨率为27km×27km，内层先验通量大小为91行×91列，空间分辨率为3km×3km。在上甸子区域，外层先验通量大小为69行×69列，空间分辨率为27km×27km，内层先验通量大小为85行×88列，空间分辨率为3km×3km。处理后先验通量单位为摩尔/秒/网格。CO先验通量及处理方式与5.2.1所述相同。
2.2.2气象场
气象场数据与5.2.2所述相同。
2.2.3CO观测
CO观测数据主要包括地基观测和卫星观测，地基观测数据与5.2.3所述相同。CO卫星观测采用MOPITT的L2级V6版本数据，该数据可提供CO柱总量及地表混合比浓度数据。本文采用CO地表混合比浓度数据，观测时间为2016年1月1日-2016年12月31日，其主要作用是结合CO2卫星近地表混合比浓度求取CO2/CO浓度比。
2.2.4 CO2观测
CO2观测采用OCO-2的V7r版本数据，观测时间为2016年1月1日-2016年12月31日。由于云、雨和气溶胶等因素的影响，CO2卫星观测数据质量有较大差异，OCO-2数据中除包含浓度数据外，还包含数据质量标识信息，质量标记为0代表质量好，1代表质量差。按照官方发布的数据说明手册，对CO2观测数据进行了筛选，剔除质量较差的数据，具体筛选指标如表6-1所示。
表6-1 OCO-2观测数据筛选指标
在本文中，CO2卫星观测主要有两方面的作用：(1)结合CO卫星近地表混合比浓度数据求取CO2/CO浓度比，具体步骤在6.3中阐述；(2)对后验CO2模拟进行精度验证及误差分析。
2.2.5本底站观测
本文采用的大气本底站观测数据来自上甸子站点和瓦里关站点，两个站点分别位于两个研究区域的中心附近。本文采用两个站点的月均浓度数据，分别作为上甸子区域和瓦里关区域的CO、CO2的背景浓度，结合CO和CO2浓度数据，求取CO2/CO浓度比。
2.3 TracersTracker区域高分辨率CO2通量反演系统运行
见技术方案：4.TracersTracker区域高分辨率CO2通量反演技术
2.4结果分析
2.4.1 CO2/CO提取结果
根据技术方案中“4.TracersTracker区域高分辨率CO2通量反演技术”所述方法，Remi提取结果如图60、图61所示。在上甸子区域，Remi绝大多数集中在10到40之间(约占85％)，全年均值为27.89。Remi最大值为184.1，主要分布在内东南沿海区域，由于该区域靠近海洋，其先验清单中的CO含量较低，导致CO2/CO比值偏高。Remi最低值为10.13，出现在西北方向的1、2月间，该部分区域城市密度低、人口较少，人为活动强度明显低于其他区域(表6-1)。这说明在远离城市的较偏远的农村地区，化石燃料燃烧效率低，容易出现燃烧不充分的现象，所以人为活动释放的CO比例要比城市及其附近区域要高。从季节性变化角度看，夏、秋两季的Remi明显比其他两季高，约为春、冬两季的1.3倍。其主要原因是夏、秋两季产生CO的人类活动的频率和强度明显比春、冬两季要小，春、冬两季的化石燃料燃烧等活动明显增多、变强，排放CO2的同时释放了更多的CO，继而减小了CO2/CO的比例。另外，Remi较高的区域主要分布在城区及其附近区域，尤其是在京-津-唐工业区范围内，该区域是中国工业最为发达、经济较为活跃的几个区域之一，频繁的工业活动产生了大量的CO2和CO。
表6-2上甸子和瓦里关区域月均CO2/CO排放比
在瓦里关区域，Remi主要分布在8.95和33(约占85％)之间，全年均值为26.68，与上甸子区域CO2/CO排放比基本持平，Remi月际变化趋势在两个区域内也基本保持一致。Remi在瓦里关区域的最大值为200.39，主要分布在研究区域的正北方向，该区域集中了瓦里关区域的大多数的工业和人口，人为活动强度大，CO2和CO排放量大。Remi最小值为8.95，主要位于研究区域的正和西南方向，该区域人口密度较低，工业活动强度较弱。从季节性变化角度看，瓦里关区域和上甸子区域略有不同，夏季的Remi要明显高于其他季节(表6-1)。
总体而言，无论是在上甸子区域还是瓦里关区域，较高的Remi一般出现在夏季或秋季，其主要原因是这两季的CO排放量较低，很大程度上提高了CO2/CO比例。另外，在两个区域内，较高的Remi主要分布在城市及其近郊，其原因可能是化石燃料等燃烧工艺提高首先出现在城市及附近区域，在相同条件下燃烧等量的化石燃料排放的CO较低，从而提高了城市及附近区域的CO2/CO比例。
表6-3上甸子区域月均ΔCO2/ΔCO浓度比
无论是在上甸子区域还是瓦里关区域，地表大气中的CO2和CO浓度的变化趋势较为相似，即在年初浓度普遍较高，由春天开始逐步递减到夏天达到最低，继而逐步升高至冬季达到最高(图62、图63)。但CO2与CO浓度具体变化情况是不同的，在两个区域内变化也有差异。CO2全年浓度变化较为平缓，月均浓度差异并不明显，变化较大的主要在冬、春两季出现。CO浓度变化幅度则较大，不同季节浓度差异可达150ppb。另外，上甸子区域CO浓度明显高于瓦里关区域，其变化幅度也要大于瓦里关区域。
根据6.3.2所述方法，Ratm提取结果如表6-3、表6-4所示。两个区域内的Ratm平均水平较接近，Ratm年均值分别为41.2和40.7，但各月的变化幅度不同。实际上，Ratm主要由CO浓度水平决定，因为CO2浓度全年变化平缓，远小于CO的变化幅度。其主要原因是CO浓度变化较易受其他因素影响，CO可与其他物质发生反应而被消耗掉，其他物质也可以被氧化而产生CO(如CH4等)。但整体而言，变化趋势明显，即春、夏两季较低，秋、冬两季较高。
表6-4瓦里关区域月均ΔCO2/ΔCO浓度比
2.4.2 CO2通量反演结果
如前所述，在CO2通量反演过程中，CO2自然通量时间分辨率为8天，在8天同化窗口内，CO2自然通量是设置为恒定不变的，所以这里仅讨论人为CO2后验通量的时空变化。图67、图68分别为两个研究区域基于排放端法反演得到的人为CO2通量增量结果，图69、图70分别为两个研究区域基于大气端法反演得到的人为CO2通量增量结果。
在每个同化窗口内，排放端法CO2人为后验通量是CO人为后验通量与相应比例相乘得到的，因此该方法获取的CO2人为后验通量增量比例与第5章获取的CO人为后验通量增量比例保持一致。在上甸子区域(表6-5)，排放端法CO2人为后验通量全年均值为47.34摩尔/秒，较人为先验通量的35.57摩尔/秒增幅达32.7％。经反演优化后，CO2排放最弱时段为5月，通量为40.53摩尔/秒，排放最强时段为12月，通量为66.12摩尔/秒，约为最弱时段的1.6倍，与人为先验通量排放强度的变化趋势一致。在瓦里关区域(表6-5)，CO2人为后验通量全年均值为4.76摩尔/秒，约为人为先验通量的1.6倍，通量的增长比例明显高于上甸子区域。其中最小值为3.46摩尔/秒，出现在8月，最大值为6.08摩尔/秒，出现在1月和12月，与人为先验通量基本一致。虽然上甸子区域和瓦里关区域的CO2人为后验通量在全年不同月份增长比例差异很大，但从全年增长趋势看，两个区域的CO2人为后验通量变化趋势基本一致，即在冬季增长比例最高，在夏季增长比例最低。
与排放端法CO2人为后验通量类似，大气端法CO2人为后验通量在冬季增长比例也高于其他季节，但全年月均增量及增量比例差异非常明显、波动较大。
在上甸子区域(图66、图68、表6-5)，CO2人为后验通量年均值为46.35摩尔/秒，大概为人为先验通量的1.3倍，与排放端法人为后验通量基本持平。其中最小值为37.85摩尔/秒，出现在4月，最大值为66.12摩尔/秒，出现在12月。在瓦里关区域(图69、图70、表6-6)，大气端法CO2人为后验通量在大多数时段也都要比人为先验通量有所提高，在1月-12月分别提高了62.6％、43.3％、29.9％、44.2％、39.8％、34.2％、42.8％、51.2％、53.2％、49.7％、46.2％和66.1％。与排放端法人为后验通量相比，虽然大气端法人为后验通量在变化趋势上较为类似，但变化幅度差异较大，大气端法人为后验通量增加比例明显小于排放端法。在总量上来讲，大气端法CO2人为后验通量均值为4.31摩尔/秒，仅为排放端法的80％左右，差距较大。在1月-12月中，最小值为3.11摩尔/秒，出现在6月，最大值为6.06摩尔/秒，出现在12月。
从季节变化角度看，在上甸子区域中(表6-6)，排放端法CO2人为后验通量在春、夏、秋、冬四季分别提高了32.2％、28.9％、31.9％和37.8％，整体上提高了32.7％。在瓦里关区域(表6-6)，排放端法CO2人为后验通量在春季、夏季、秋季、冬季分别提高了61.7％、58.9％、61.6％和66.8％，整体上提高了62.3％。在上甸子区域中(表6-6)，大气端法CO2人为后验通量在春季、夏季、秋季、冬季分别提高了22.7％、23.7％、30.2％和44.7％，整体上提高了30.3％。在瓦里关区域(表6-6)，大气端法CO2人为后验通量在春季、夏季、秋季、冬季分别提高了38.0％、42.7％、53.1％和57.3％，整体上提高了46.9％。
从两个研究区域的CO2通量反演结果看，排放端法和大气端法得到的CO2人为后验通量在时空分布上有相似之处，但差异也非常明显。从人为后验通量增量上看，两种方法反演的CO2人为后验通量在冬季的增量均比其他季节要高，尤其比夏季要高，主要原因是夏季的CO2人为先验通量非常低。但是，从人为后验通量增长比例角度看，大气端法人为后验通量与排放端法人为后验通量具有非常大的差异，其增长比例与排放端法人为后验通量增长比例完全不同。大气端法人为后验通量在冬季较大，在夏季较小，与排放端法人为后验通量的变化趋势相反。另外，1月-12月的通量增量和增长比例的波动也比排放端法人为后验通量大的多。
表6-6排放端法和大气端法CO2后验通量增比
另外，从两个研究区域内CO2人为后验通量增量来看，通量增量主要集中于城市区域或人口密集区，在上甸子区域中通量增长最多区域在唐山市境内，以唐山市主城区为主，其次为天津地区和北京地区，承德市区域内通量增长较低。但从通量增长比例角度看，增长率较高的主要区域并不位于主城区，一般情况下位于城区附近区域或郊区，这主要得益于主城区污染企业的外移以及城郊区域自身经济的发展。在瓦里关区域，CO2通量增量较为明显区域主要位于西宁市主城区附近，其次是临夏市城区及周边范围内，其他区域如海北藏族自治州、海南藏族自治州、海东地区等增量较小。但从后验通量增比看，非城区区域增比要明显高于城区及其近郊，其主要原因是这些区域的先验通量偏低较多所致。
2.4.3人为CO2通量区分
大气中CO2的浓度不仅受人类活动影响，还在很大程度上受自然环境的制约，在人类不断向大气排放CO2的同时，自然界也在进行着CO2吸收和排放活动，如“土壤呼吸”会释放CO2、植物光合作用会吸收CO2，其他生态系统如稻田、沼泽、草原等在不同季节扮演不同角色。在本文中，CO2通量包括人为CO2通量和自然CO2通量，总通量为两者之和。在上一节中，基于两种反演方法分别获取了人为CO2后验通量和CO2后验总通量，本节主要对人为CO2通量和自然CO2通量进行对比分析，图68、图69分别为上甸子和瓦里关两个区域内人为和自然CO2通量的具体组成，其中橙色为人为CO2通量，绿色为自然CO2通量。
总体来看，上甸子区域CO2先验总通量、排放端法后验总通量、大气端法后验总通量年均值分别为25.5摩尔/秒、38.1摩尔/秒和36.5摩尔/秒，整个区域每年都向大气中排放CO2，是明显的碳源(表6-7)。瓦里关区域CO2先验总通量、排放端法后验总通量、大气端法后验总通量年均值分别为-5.8摩尔/秒、-4.0摩尔/秒和-4.4摩尔/秒，整个区域每年都从大气中吸收CO2，是明显的碳汇(表6-8)。
但是，从月际CO2通量变化情况看，两个区域表现有相同之处，但也存在较大差异(表6-7、表6-8)。在人为CO2通量方面，无论是在先验通量还是后验通量中，均是由春季较低值逐步增长，在冬季达到最大值。自然CO2通量与人为CO2通量的月际变化正好相反，一般是夏季较高、冬季较低。两个区域CO2自然通量年均值分别为-10.1摩尔/秒和-9.8摩尔/秒，整体水平相当，且自然CO2通量最高值一般出现在7、8月间。虽然两个区域在CO2人为通量和自然通量变化趋势均较为相似，但在月际CO2总通量上是有明显差异的。在上甸子区域，先验CO2总通量、排放端法CO2总通量、大气端法CO2总通量在1月-12月均为碳源，先验总通量中最小值在9月，最大值在12月，排放端法总通量中最小值在9月，最大值在12月，大气端法总通量中最小值在5月，最大值在12月，三者整体变化趋势基本一致。瓦里关区域则不同，先验CO2总通量、排放端法CO2总通量、大气端法CO2总通量在不同月份表现有明显差异，5月-10月基本为碳汇，其他时间为碳源，整个区域随着时间迁移在碳源、碳汇角色中切换。
表6-7上甸子区域CO2先验和后验总通量
表6-8瓦里关区域CO2先验和后验总通量
从季节角度看，两个区域CO2总通量也有不同表现。在上甸子区域，先验CO2总通量、排放端法CO2总通量、大气端法CO2总通量在各季节均表现为碳源且整体变化趋势基本一致。在瓦里关区域，先验CO2总通量中春、夏、秋三季均为碳汇，仅有冬季为碳源。在排放端法CO2总通量和大气端法CO2总通量中，夏、秋两季为碳汇，春、冬两季为碳源。
整体而言，在上甸子区域，人为活动排放的CO2的量要大于自然吸收CO2的能力，使该区域在特定时段地表空气中的CO2浓度迅速升高，尤其是在人口密集、工业活动强度大的区域，如唐山、天津等，使其空气质量迅速变差甚至出现严重污染事件。在瓦里关区域，整个区域内自然吸收CO2的能力大于人为活动排放的CO2的量，但是在局部区域，如西宁、临夏等，人为活动排放的CO2的量超出自然吸收CO2的速度，仍然可以在某些时段造成污染事件。
上甸子区域是中国经济最发达的区域之一，高强度的生产、生活活动产生了大量的温室气体(如CO2、CH4等)，这在很大程度上提升了大气中CO2的浓度。相比而言，瓦里关区域人口较少，经济欠发达，工业基础较为薄弱，在日常生产、生活中产生的温室气体相对较少。同时可以看到，虽然两个区域内由于人类活动产生的CO2在总量上相差较大，但季节性变化趋势较为相似，即在冬季(12月到次年2月)一般达到最高值，在夏季(6月-8月)降到最低值。主要的原因是上甸子区域和瓦里关区域都位于中国北方，冬季居民供暖等消耗了大量的化石燃料等能源，排放了大量的CO2到大气中。不同之处在于，在低强度的人类生产生活活动下，瓦里关区域的人为CO2通量的季节性变化不是很明显。与人为CO2排放源时空变化特征不同的是，自然CO2排放通量一般在冬季较小，随着温度升高在夏季达到最大，其具体的时空变化主要受当地生态系统影响。同时，瓦里关区域的自然CO2通量(绝对值)一般要比上甸子区域的大，这也说明了瓦里关区域的陆地生态系统发育更好、光合作用更强。
2.4.4误差分析及讨论
本节利用CO2卫星观测数据对先验模拟、排放端法后验模拟、大气端法后验模拟的准确性进行验证。CO2后验模拟采用的通量为CO2总通量，即人为通量和自然通量之和。除CO2通量外，先验模拟和后验模拟的参数配置均相同。另外，CO2模拟输出为大气垂直方向23层，空间分辨率为3km×3km的网格浓度，为与卫星观测的柱浓度数据进行对比，首先需将CO2模拟浓度转为柱浓度数据，转换公式如下：
其中，H(x)代表CMAQ模拟的CO2插值到GOSAT和OCO-2卫星观测的各气压层，
为CMAQ模拟的柱浓度，
是卫星先验柱浓度，a是卫星柱平均核，h是气压权重函数，ya是卫星先验廓线。图42、图43为先验模拟、后验模拟结果，图44为先验模拟、后验模拟的误差分布情况。
表6-9 CO2先验和后验模拟误差
在上甸子区域，先验模拟可以反映CO2浓度变化的基本趋势，但其相关系数仅为0.67，对CO2浓度变化细节的把握仍存在不足。另外，先验模拟误差均值为-4.91ppm，模拟有明显低估。在瓦里关区域，先验模拟表现也较差，相关系数仅为0.69，但基本能反映出CO2浓度变化的趋势(图42)。同时，瓦里关区域后验模拟误差水平比上甸子区域略小，但整体上还是有明显低估(误差均值-3.34ppm)。整体而言，无论是在上甸子区域还是瓦里关区域，先验模拟与观测数据的相关系数均较低并且有明显的误差。
与先验模拟相比，大气端法后验模拟值的准确性在某些时段和区域都有不同程度的提高，但整体上模拟精度并没有明显提升甚至出现了下降(表6-9)。在上甸子区域，后验模拟的误差绝对值由11.8ppm减小到11.6ppm，绝对误差水平基本没有变化。误差均值从-4.91ppm减小为-0.54ppm，模拟低估表现有明显减小。但是，后验模拟与观测之间的相关性表现较差，相关系数从0.67变为0.66，没有提升反而出现了下降。在瓦里关区域，后验模拟的绝对误差有所减小，由11.2ppm减小到10.9ppm。误差均值从-3.49ppm减小到-2.67ppm，模拟低估表现有所缓解。但是，后验模拟与观测之间的相关性出现小幅度下降，相关系数从0.69变为0.68。显然，大气端法反演得到的CO2后验通量存在较大不确定性，无法提高后验模拟的精度。
与大气端法不同，无论是在上甸子区域还是瓦里关区域，排放端法后验模拟精度均有明显的提升，尤其在污染较为严重的时段，如1月和12月(表6-9)。在上甸子区域，后验模拟的误差绝对值由11.8ppm减小到4.15ppm，误差整体水平有大幅度的减小。误差均值从-4.91ppm提升至-1.03ppm，模拟低估表现有明显缓解。另外，后验模拟与观测之间的相关性也有较大提升，相关系数由0.67提升到0.88。在瓦里关区域，后验模拟的绝对误差也有明显的减小，由11.2ppm减小到4.26ppm。误差均值从-3.49ppm提升到了-0.64ppm，模拟低估表现明显缓解。此外，后验模拟与观测之间的相关性也有较大提升，相关系数从0.69变为0.89。
在误差分布上，先验模拟、大气端法后验模拟和排放端法后验模拟的误差分布情况有较大不同(表6-9)。无论是大气端法后验模拟还是排放端法后验模拟，都在很大程度上改变了CO2模拟误差的分布形态。在先验模拟中，上甸子区域和瓦里关区域的误差区间(累积百分比10％到90％)分别为[-22.8,12.8]和[-21.6,12.5]。在排放端法后验模拟中，误差区间分别变为[-10.5,7.5]、[-10.8,6.2]，误差范围明显向0值方向收缩，说明排放端法后验模拟在很大程度上减小了CO2模拟误差。在大气端法后验模拟中，误差区间为[-19.5,17.5]、[-21.5,13.8]，虽然某些时段模拟误差有所减小(负值部分)，但也有误差增大的现象(正值部分)，整体上误差并没有明显改变。这说明大气端法后验通量存在较大程度的不确定性，无法有效提高后验模拟的精度。
上述结果表明，无论是在上甸子区域还是瓦里关区域，排放端法均比大气端法表现要好，它极大的消除了CO2先验通量中的“误差”，CO2后验模拟准确性明显提升。而大气端法表现较不稳定，虽然CO2后验模拟误差有所减小，但相关系数有所下降，整体上CO2后验通量的准确性没有明显提升。
理论上，无论是排放端法还是大气端法，基于CO2/CO相关性数据进行CO2通量反演时均会引入新的不确定性，其主要原因是参与同化的观测数据是间接推算得到的，而非直接的观测或统计数据。但不确定性在两种方法中的影响程度是不同的，由此得到的CO2通量反演结果也存在明显差异。
在排放端法中，反演结果的不确定性主要来自排放清单中的CO2/CO比例数据。本文使用的排放清单为清华大学MEIC团队发布的MIX清单，该清单已被应用到诸多研究中，如东亚模式比较计划第三期模拟研究、联合国半球大气污染传输计划，其准确性是经过多方验证的，因此从清单中提取的CO2/CO排放比数据也较为可靠。从数据本身看，无论是在上甸子区域还是瓦里关区域，各月的CO2/CO排放比相关性非常高，上甸子区域12个月排放相关系数均值为0.94，瓦里关区域12个月排放相关系数均值为0.95(图72)，这说明两个区域内CO2/CO排放比在不同时段是整体变化趋势基本一致，相对波动幅度也较小。另外，两个区域间的CO2/CO排放比在不同月份变化趋势相似度也较高(图76)，夏、秋两季较高，冬、春两季较小，并且两个区域CO2/CO排放比大小也相近。这说明即使在不同区域，CO2/CO排放比也具有较好的稳定性。综上所述，无论是在上甸区域还是瓦里关区域，CO2/CO排放比整体上表现较为稳定、可靠，在已优化的CO后验通量基础上，依据此排放比获取的CO2后验通量也是“优化”的。
相较于排放端法，大气端法涉及的数据和步骤较多，不确定性来源也较多，主要包括以下几个方面：
(1)观测数据较少且空间分布不均匀。在CO通量反演中，两个研究区域内共有31个地基观测站点，CO浓度观测数据具有较高空间分辨率。与CO观测站点相比，CO2观测站点及卫星观测稀少，在两个研究区域内仅各有1个大气本底观测站点，卫星观测数据严重依赖卫星过境路径，空间分布不连续且分布很不均匀。利用这些有限的观测数据的均值代表整个区域CO2浓度变化，一定程度上忽略了CO2浓度空间变化特征，从而引起新的误差和不确定性。这也是当前利用卫星数据进行区域CO2通量反演普遍存在的问题。
(2)观测时间连续性较差，月均浓度比掩盖了细节变化。在CO通量反演中，地基观测站点可提供逐小时CO浓度观测数据。由于OCO-2观测数据和MOPITT观测数据在时间上连续性较差，实际反演过程中采用月均CO2/CO浓度比。但是，CO2/CO浓度比变化幅度较大，其主要原因是排放到大气中的CO化学特性较为活跃，可与其他物质发生反应而被消耗掉，其他物质也可被氧化产生CO。图74为大气端法采用的MOPITT获取的CO地表混合比浓度数据，从中可以看出，CO浓度在短时间内可以出现剧烈变化，最高浓度和最低浓度的差值可达1500ppb(表6-10)。在CO观测数据充足的情况下，这种变化可通过多次观测取均值的方法抵消，但是在CO观测数据不足的情况下，CO的浓度变化会明显增加CO2/CO浓度比的“偏差”，给最终的CO2通量反演结果带来了较大不确定性。
(3)多源观测数据标准不一。在求取CO2/CO浓度比过程中，采用的数据包括地基观测、OCO-2卫星观测、MOPITT卫星观测，各数据间的观测仪器、方法有很大不同，同时各数据的时空分辨率也有明显差异，虽然进行了时间和空间匹配，但仍存在不同程度的不确定性。
虽然排放端法和大气端法在进行CO2通量反演是会引入新的不确定性，但本文的CO2通量反演结果表明，依据优化后的CO通量及较为可靠的CO2/CO排放比，排放端法明显提高了CO2后验通量的准确性，CO2后验模拟误差明显减小，在当前CO2观测数据无法满足区域高分辨率通量反演的情况下，排放端法仍不失为较为可靠的CO2通量反演方案，可有效提高CO2通量的时空分辨率。大气端法涉及的数据和流程较多，引入的不确定性也较多，同时CO2、CO观测数据在时间、空间、标准上都有不同程度的兼容问题，最终导致CO2通量反演结果的准确性没有明显提升，在当前CO2观测数据较为缺乏的情况下，利用该方法进行区域CO2通量反演仍不太稳定。
表6-10 CO卫星观测数据
2.5小结
为克服区域CO2观测数据不足导致通量无法反演的问题，本发明提出了两种基于CO数据的CO2通量反演方法，即基于CO2/CO排放比的通量反演方法(排放端法)和基于CO2/CO浓度比的通量反演方法(大气端法)，充分利用当前区域尺度上较为充足的CO观测数据和排放清单数据，反演得到了区域高分辨率CO2通量，并在此基础上区分了CO2人为通量和自然通量，最后对两种方法获取的CO2通量的准确性进行了验证并讨论了其中的不确定性。主要包括四部分内容：
第一部分：阐述了TracersTracker系统进行CO2通量反演所需的数据并介绍了数据的处理过程。
第二部分：详细阐述了两种基于CO数据的CO2通量反演方法的基本思路和实施步骤。排放端法利用MIX清单中的CO2/CO排放比，结合第5章已获取的CO后验通量，反演得到了区域CO2后验通量。大气端法利用有限的卫星、地基浓度观测数据，获取地表大气CO2/CO浓度比，依据此浓度比及当前较为丰富的CO地基观测数据推算“CO2观测数据”，利用TracersTracker系统同化“CO2观测数据”反演得到了区域高分辨率(3km)CO2通量。
第三部分：CO2通量反演结果及讨论。结果显示，两个研究区域内CO2通量增长量主要集中于城市区域或人口密集区，在上甸子区域中通量增长最多区域在唐山市境内，以唐山市主城区为主，其次为天津地区和北京地区，承德市区域内通量增长较低。在瓦里关区域，CO2通量增长较为明显区域主要位于西宁市主城区附近，其次是临夏市城区及周边范围内，其他区域如海北藏族自治州、海南藏族自治州、海东地区等增长量较小。
上甸子区域在多数时间扮演碳源的角色，CO2先验通量、排放端法后验通量、大气端法后验通量年均值分别为25.5摩尔/秒、38.1摩尔/秒和36.5摩尔/秒，整个区域碳源在逐步增加。瓦里关区域在多数时间为明显的碳汇，在CO2先验通量、排放端法后验通量、大气端法后验通量中的年均值分别为-5.8摩尔/秒，-4.0摩尔/秒和-4.4摩尔/秒，整个区域碳汇呈减小趋势
第四部分：CO2后验通量误差分析及不确定性讨论。无论是在模拟误差还是模拟相关性上，排放端法获取的CO2后验通量都要优于大气端法，排放端法后验模拟误差明显减小，模拟精度、相关系数有较大提升。大气端法后验模拟在某些时段和区域内精度有所提高，但后验模拟中仍存在较大误差，模拟相关性也未出现明显提升。
以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰，均落在本发明的保护范围内。
Claims (10)
1.一种基于CO2/CO比率反演核算高分辨率人为CO2排放的方法，其特征在于，包括：
采用FORTRAN语言和正交4维变分POD4DVar同化算法，实现POD4DVar同化算法与区域大气传输模型CMAQ的耦合，研发了区域高分辨率同化反演系统TracersTracker；
在TracersTracker系统的支持下，对上甸子区域和瓦里关区域的CO通量进行反演和优化，获取两个研究区域内的最优CO通量；
基于排放清单中的CO2/CO排放比及获取的最优CO通量，反演得到上甸子区域和瓦里关区域的人为CO2排放通量；
基于地表大气中的ΔCO2/ΔCO浓度比及CO地基观测数据，利用TracersTracker系统同化CO2观测数据，反演得到上甸子区域和瓦里关区域的高分辨率人为CO2排放的通量；
结合自然CO2排放通量，在CO2总通量内区分人为CO2通量，并在此基础上分析两种方法获取的人为CO2排放通量的核算精度及其不确定性。
2.根据权利要求1所述的基于CO2/CO比率反演核算高分辨率人为CO2排放的方法，其特征在于，CMAQ是一个三维欧拉大气化学和传输模拟系统，可以模拟颗粒物、臭氧、有毒空气污染物的80多种物质。
3.根据权利要求1所述的基于CO2/CO比率反演核算高分辨率人为CO2排放的方法，其特征在于，CMAQ包括：化学传输模块、气象化学接口模块、边界场模块、初始场模块、光分解模块。
4.根据权利要求1所述的基于CO2/CO比率反演核算高分辨率人为CO2排放的方法，其特征在于，本专利所研发的TracersTracker CO2同化系统涉及到的第三方软件或函数库包括：IOAPI函数库、NetCDF函数库、Matlab函数库。
5.根据权利要求1所述的基于CO2/CO比率反演核算高分辨率人为CO2排放的方法，其特征在于，TracersTracker，主程序采用FORTRAN90语言编写，编译器为Intel FortranComplier。
6.根据权利要求1所述的基于CO2/CO比率反演核算高分辨率人为CO2排放的方法，其特征在于，POD4Dvar同化算法采用蒙特卡罗法生成样本集合来表征通量的空间分布格局及时间演变特征。
7.根据权利要求6所述的基于CO2/CO比率反演核算高分辨率人为CO2排放的方法，其特征在于，POD4Dvar同化算法无需积分伴随模式。
8.根据权利要求1所述的基于CO2/CO比率反演核算高分辨率人为CO2排放的方法，其特征在于，同化过程包含若干次同化循环。
9.根据权利要求1所述的基于CO2/CO比率反演核算高分辨率人为CO2排放的方法，其特征在于，采用均方根误差、线性相关系数、绝对误差均值、误差均值对同化系统反演结果的准确性进行评价。
10.一种基于CO2/CO比率反演核算高分辨率人为CO2排放的系统，其特征在于，包括：
一个或多个处理器；
存储装置，用于存储一个或多个程序；
当所述一个或多个程序被所述一个或多个处理器执行，使得所述一个或多个处理器实现根据权利要求1至9任意一项所述的基于CO2/CO比率反演核算高分辨率人为CO2排放的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210640366.XA CN116205022A (zh) | 2022-06-07 | 2022-06-07 | 基于co2/co比率反演核算高分辨率人为co2排放的方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210640366.XA CN116205022A (zh) | 2022-06-07 | 2022-06-07 | 基于co2/co比率反演核算高分辨率人为co2排放的方法及系统 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116205022A true CN116205022A (zh) | 2023-06-02 |
Family
ID=86510097
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210640366.XA Pending CN116205022A (zh) | 2022-06-07 | 2022-06-07 | 基于co2/co比率反演核算高分辨率人为co2排放的方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116205022A (zh) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101326268B1 (ko) * | 2012-12-12 | 2013-11-11 | 사단법인대기환경모델링센터 | 씨엠에이큐용 배출량 입력자료 생성 방법 |
| CN111487216A (zh) * | 2020-05-06 | 2020-08-04 | 北京中科锐景科技有限公司 | 一种二氧化碳通量反演方法、系统 |
| CN114564841A (zh) * | 2022-03-03 | 2022-05-31 | 上海市环境科学研究院 | 城市大气排放清单反演方法、系统、设备及存储介质 |
-
2022
- 2022-06-07 CN CN202210640366.XA patent/CN116205022A/zh active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101326268B1 (ko) * | 2012-12-12 | 2013-11-11 | 사단법인대기환경모델링센터 | 씨엠에이큐용 배출량 입력자료 생성 방법 |
| CN111487216A (zh) * | 2020-05-06 | 2020-08-04 | 北京中科锐景科技有限公司 | 一种二氧化碳通量反演方法、系统 |
| CN114564841A (zh) * | 2022-03-03 | 2022-05-31 | 上海市环境科学研究院 | 城市大气排放清单反演方法、系统、设备及存储介质 |
Non-Patent Citations (1)
| Title |
|---|
| 鲁立江: "区域高分辨率碳同化系统研发及人为碳排放估算研究", 中国博士学位论文全文数据库（电子期刊）, no. 2021, 15 March 2021 (2021-03-15), pages 009 - 2 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111401634B (zh) | 一种获取气候信息处理方法、系统、存储介质 | |
| Engardt et al. | Deposition of sulphur and nitrogen in Europe 1900–2050. Model calculations and comparison to historical observations | |
| Ding et al. | Intercomparison of NO x emission inventories over East Asia | |
| Nielsen et al. | Chemical mechanisms and their applications in the Goddard Earth Observing System (GEOS) earth system model | |
| Carmichael et al. | Evaluating regional emission estimates using the TRACE‐P observations | |
| Sillman et al. | Reactive mercury in the troposphere: Model formation and results for Florida, the northeastern United States, and the Atlantic Ocean | |
| CN111487216B (zh) | 一种二氧化碳通量反演方法、系统 | |
| Palmer et al. | A measurement-based verification framework for UK greenhouse gas emissions: an overview of the Greenhouse gAs Uk and Global Emissions (GAUGE) project | |
| Ukhov et al. | Study of SO pollution in the Middle East using MERRA‐2, CAMS data assimilation products, and high‐resolution WRF‐Chem Simulations | |
| Lian et al. | Sensitivity to the sources of uncertainties in the modeling of atmospheric CO 2 concentration within and in the vicinity of Paris | |
| Xiang et al. | Nitrous oxide (N2O) emissions from California based on 2010 CalNex airborne measurements | |
| CN115186437B (zh) | 碳同位素联合同化模型及甄别人为碳排放与自然碳通量区域同化系统的构建方法 | |
| CN116070500B (zh) | 一种基于深度学习空气质量浓度场模拟仿真器 | |
| Park et al. | Numerical simulation of atmospheric CO2 concentration and flux over the Korean Peninsula using WRF-VPRM model during Korus-AQ 2016 campaign | |
| Chen et al. | Evaluation of regional CO2 mole fractions in the ECMWF CAMS real‐time atmospheric analysis and NOAA carbontracker near‐real‐time reanalysis with airborne observations from ACT‐America field campaigns | |
| Turnbull et al. | Radiocarbon in the atmosphere | |
| He et al. | China's terrestrial carbon sink over 2010–2015 constrained by satellite observations of atmospheric CO2 and land surface variables | |
| Keller et al. | Description of the NASA GEOS composition forecast modeling system GEOS-CF v1. 0 | |
| Wu et al. | Three-dimensional spatiotemporal variability of CO2 in suburban and urban areas of Shaoxing City in the Yangtze River Delta, China | |
| Banos et al. | Assessment of the data assimilation framework for the Rapid Refresh Forecast System v0. 1 and impacts on forecasts of a convective storm case study | |
| Sekiya et al. | Impacts of horizontal resolution on global data assimilation of satellite measurements for tropospheric chemistry analysis | |
| Fu et al. | MICS-Asia II: Modeling gaseous pollutants and evaluating an advanced modeling system over East Asia | |
| Nabavi et al. | Spatiotemporal variation of radionuclide dispersion from nuclear power plant accidents using FLEXPART mini-ensemble modeling | |
| CN116205022A (zh) | 基于co2/co比率反演核算高分辨率人为co2排放的方法及系统 | |
| Brasseur et al. | Ensemble forecasts of air quality in eastern China-Part 1: Model description and implementation of the MarcoPolo-Panda prediction system, version 1 |
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 |