随着对大气环境质量研究的深入，大气中氧气含量变化也越来越受到关注。研究表明，越往前追溯，地球大气含氧量就越高，大多数科学家都认同史前大气平均含氧量在30%至35%之间，而目前只有21%。在一些污染严重、人口密度较大的城市和工业园区，大气含氧量只有15%甚至更低^{[1, 2, 3, 4]}。测量数据表明，低层大气氧气含量水平正在以平均每年2 mg/m³的速率缓慢减少^{[5, 6]}，北半球与南半球大气中的氧气不但逐年减少，而且会随季节变化^[7]。大气中含氧量低于19.5%时，人就会产生缺氧症状，人类所有的慢性疼痛以及疾病均由细胞缺氧引起^{[8, 9]}。

对于大气中氧气含量的研究，国内外学者除了通过测定化石气泡中的氧气含量^[10]、直接监测^{[6, 7]}、及建立生物化学模型估算大气中氧气含量^{[11, 12, 13]}的方法外，运用碳氧平衡法，依据植物生产力估算区域氧气释放量，间接揭示大气氧气含量变化，此类研究在国内外也有较多报道^{[14, 15, 16, 17, 18, 19]}。其研究方法可概括为两种：一是通过估算区域氧气消耗量，折算出区域氧气释放量^[20]；二是通过估算区域内绿色植物生产的有机碳量，按有机碳与氧气的质量比，计算区域氧气释放量。比较而言，第一种方法主观性较强，所以运用第二种方法的研究成果较多。所不同的是，有的学者依据植被初级生产力(GPP)^{[21, 22]}推算，有的学者依据净初级生产力(NPP)^{[23, 24, 25]} 推算。GPP是植物通过光合作用生产的初始有机碳量，NPP是考虑了植物自养呼吸后剩余的有机碳量。由于依据不同，估算出的氧气释放量差异较大。相比GPP和NPP，NEP是净初级生产力减去异养呼吸消耗的同化产物，是考虑了呼吸消耗后的净光合产量^{[26, 27]}，因此依据NEP估算得出的为氧气净释放量，更能准确指示区域大气中氧气含量的变化。而目前相关研究还较少，已有研究也很少结合遥感及GIS技术，分析区域绿色植被氧气释放量空间分布及变化特征。

绿色植被面积减少会导致氧气释放量减少^[1]。森林是陆地绿色生态系统的重要组成部分，具有固碳释氧的重要生态功能^[25]，地球上60%以上的氧气来自森林生态系统^[20]。黑龙江省是中国森林资源最丰富的省份，森林面积约占全国森林总面积的1/7。其广袤的森林是氧气生产源，二氧化碳汇，对调节区域生态环境起着重要作用。但20世纪百年间，由于人类活动影响，黑龙江省森林面积减少了37.163%，森林面积减少所引起的区域氧气释放量变化，此类研究还未见报道。

本文以1900年、1949年、1986年和2009年黑龙江省森林覆盖信息为依据，采用C-FIX模型模拟森林净生态系统生产力(NEP)，基于碳氧平衡原理，运用GIS空间分析方法，模拟了黑龙江省区域20世纪森林氧气释放量的时空变化，为探索人类活动对生态环境的影响提供依据。

黑龙江省(121°11′—135°05′ E，43°26′—53°33′ N)是全国最大的林业省份之一，林业生态地位十分重要。森林类型主要有：北温带落叶针叶林，分布在大兴安岭山地北部，树种以兴安落叶松(Larix gmelini)为主；中温带针阔叶混交林，分布于小兴安岭和东部山地。黑龙江省属大陆性季风气候，年平均气温在-4—5 ℃。气温由东南向西北逐渐降低，南北差近10 ℃。全省年平均降水量多介于400—650 mm。太阳辐射资源比较丰富，年太阳辐射总量在46×10⁸—50×10⁸ J/m²之间。其中，5—9月的太阳辐射总量占全年的54%—60%。全省日照时数在2300—2800 h，其中生长季日照时数占总量的44%—48%。土壤有机质含量高于全国其它地区，黑土、黑钙土和草甸土等占耕地的60%以上，是世界著名的三大黑土带之一^[28]。

(1)1900、1949、1986年黑龙江省森林空间分布。李景文等1993年编辑出版的《黑龙江省森林》中，附有清末民初(1896年)、1949年、1986年1 ∶ 300万比例尺的黑龙江省森林分布纸质图^[29]。该系列图件描绘了20世纪80年代前黑龙江省森林类型及其空间分布状况。本文将森林分布纸质图扫描，与具有地理坐标系统的图件进行坐标配准，通过定义地理数据库的属性表及点、线、面矢量化跟踪，完成此系列纸质图的矢量化(图 1—图 3)。

图 1 1900年黑龙江省森林覆盖空间分布 Fig.1 Spatial distribution of forest covers of the Heilongjiang Province in 1900

图选项

图 2 1949年黑龙江省森林覆盖空间分布 Fig.2 Spatial distribution of forest covers of the Heilongjiang Province in 1949

图选项

图 3 1986年黑龙江省森林覆盖空间分布 Fig.3 Spatial distribution of forest covers of the Heilongjiang Province in 1986

图选项

(2)2009年黑龙江省森林空间分布。依据2009年9月ETM⁺遥感影像(覆盖黑龙江省共39景影像，条带号为113—123)，采用监督分类法将土地利用类型分为6类：耕地、有林地、草地、水域、建设用地和未利用地，提取有林地覆盖图(图 4)。ETM⁺遥感影像从中国科学院遥感研究所购买。

图 4 2009年黑龙江省森林覆盖空间分布 Fig.4 Spatial distribution of forest covers of the Heilongjiang Province in 2009

图选项

为保证模拟结果具有可比性，并减少异常年份波动影响，月平均气温、地面辐射通量分别采用附近5年月均值代替。如1900年7月的平均气温、辐射采用1900—1905年7月均值代替。

(1)月平均气温 资料来源于英国东英吉利大学气候研究所(CRU)提供的0.5°×0.5°逐月网格数据。CRU气候要素资料是根据陆面实际观测资料插值到0.5°×0.5°的经纬度网格数据。CRU资料在中国区域尺度上，与实际观测资料序列吻合得很好，通过了显著性检验^[30]。本文提取1900—1905年、1949—1953年、1986—1990年、2005—2009年覆盖黑龙江省及附近的网格资料共256个点，采用Kriging插值方法，将各年、各月图层叠加，并进行栅格运算，得到1900年、1949年、1986年、2009年黑龙江省3—11月气温1 km×1 km栅格图共36幅。

(2)月地面辐射通量 因CRU资料中没有地面辐射要素，需要依据云量和天文辐射估算^[31]：

点击浏览公式

式中，S_g,d为地面辐射日总量(MJ m^-2 d^-1)；C_a为地面总辐射与天文辐射的比例^[31]；S_0,d为天文辐射日总量(MJ m^-2 d^-1)。

首先计算各纬度天文辐射日总量(S_g,d)^[32]，然后提取1900—1905年、1949—1953年、1986—1990年、2005—2009年上述256个网格点的云量资料，按式(1)计算各网格地面太阳辐射通量，用Kriging插值方法，将图叠加平均后，得到1900年、1949年、1986年、2009年3—11月地面辐射通量1 km×1 km栅格图共36幅。

(1)1986年、2009年森林植被指数(NDVI)获取。1984—1988年、2005—2009年3—11月植被指数(NDVI)从美国航天局NASA网站下载，分辨率1 km×1 km。图层叠加平均后，与黑龙江森林分布图叠加，得到1986年、2009年黑龙江省3—11月NDVI分布图共18幅。

(2)1900年、1949年森林植被指数(NDVI)获取。将1986年黑龙江省森林分布图与1986年各月森林NDVI图进行叠加，提取各纬度不同林型3—11月份NDVI(表 1，以48°—49°N为例)，再按纬度按林型赋值到1900年、1949年森林分布图中，得到1900年、1949年黑龙江省3—11月1 km×1 km森林NDVI栅格图18幅(以7月份为例，图 5—图 6)。

图 5 1900年7月NDVI分布图 Fig.5 NDVI distribution map of the Heilongjiang Province in July 1900

图选项

图 6 1949年7月NDVI分布图 Fig.6 NDVI distribution map of the Heilongjiang Province in July 1949

图选项

C-FIX模型，是Veroustraete在2002年提出的网格化光能利用率模型，作者通过对模型进行多次验证和参数修正，使之可以在区域尺度上估算2个基本的碳循环分量(GPP、NPP、NEP)。该模型的优点是要求输入的参数少、计算效率高、输出结果的时空分辨率高^[26]。对于每一给定的网格，C-FIX模型使用以下公式来逐步计算每天的GPP、NPP和NEP值，单位为gC m^-2 d^-1，下标d代表日值：

点击浏览公式

式中，p(T_atm)为归一化气温依赖因子，取值[0, 1]；CO₂fert为归一化施肥效应因子；fAPAR为植被可吸收的光合作用有效辐射比例(系数)，取值[0, 1]。

点击浏览公式

点击浏览公式

点击浏览公式

式中，K_m为核糖1，5二磷酸羧化氧化酶与CO₂的结合能力; K₀为O₂的阻止力常数；E_a为气体常数(KJ/mol)；A为常数。

如果T≥15 ℃(288.13 K)，那么E_a1 = 59.4 kJ/mol，且A₁ = 2.419×10¹³。

如果T<15 ℃(288.13 K)，那么E_a2 = 109.6 kJ/mol，且A₂ = 1.976×10²²。

计算K₀时，A₀ = 8240，E_a0 = 13.9135 kJ/mol。

点击浏览公式

其中，A_t= 7.87×10^-5，E_at = -42.8969 kJ/mol。

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

式中，A_d为植被自养呼吸率；R_h,d为土壤异养呼吸通量(gC m^-2 d^-1)；K_s,y为全年日平均异养呼吸效率(gC m^-2 d^-1)。

上述式中还未说明的参数的定义和单位如表 2所示。

除表 2中已赋值的参数外，日平均气温、日地面辐射通量、大气CO₂混合浓度、NDVI是需要输入的参数，K_m、K₀等为中间计算参数。

植物光合作用是植物叶绿素在阳光照射下吸收空气中二氧化碳(CO₂)和水(H₂O)，并将其转化成有机物(C₆H₁₂O₆)和氧气。同时，植物还通过呼吸作用消耗有机物和氧气，去掉呼吸消耗后剩余的有机碳量，称为净生态系统生产力(NEP，gC m^-2 d^-1)。NEP所对应的氧气量即为植物所释放的净氧气量，NEP与氧气质量比为1 ∶ 2.667。

采用ArcGIS空间叠加、栅格计算、空间统计等分析方法，计算黑龙江省及各行政区内森林变化及各月、年森林氧气释放量。

为研究森林变化对森林氧气释放量的影响，设计控制实验：假设百年森林覆盖性质不变，采用1900年森林覆盖信息及NDVI、2009年气候要素驱动C-FIX模式，得到控制实验下的森林氧气释放量模拟值。将此模拟值减去2009年实际森林氧气释放量模拟值，所得到的差值为百年森林变化引起的森林氧气释放量变化值，表达式为：

点击浏览公式

式中，Δ_forest表示森林变化引起的氧气释放变化量，O_2,contro1为控制实验氧气释放量模拟值，O_2,2009为2009年实际森林氧气释放量模拟值。

黑龙江省1900—2009年森林面积呈持续减少趋势。森林总面积减少了10.667万km²，森林覆盖率减少了23.568%，大庆市森林面积减少最多，约减少了94.673%(表 3)。空间上以松嫩平原和三江平原减少显著(图 7)。

图 7 黑龙江省1900—2009年百年森林面积变化 Fig.7 Change of forest area in the Heilongjiang Province from 1900 to 2009

图选项

黑龙江省1900—1986年除樟子松林、针阔混交林面积有所增加外，其它林型面积均减少，其中云冷杉林及红松林面积减少超过80%(表 4)。

将1900年、1949年、1986年及2009年日平均气温、日平均辐射(均用月平均值代替)、CO₂浓度^[33](用年平均浓度代替)及其它参数输入C-FIX模型，计算NEP，换算成氧气释放量。采用ArcGIS空间统计功能计算出各月及年森林氧气释放量(表 5)，年氧气释放量空间分布如图 8—图 11。

图 8 黑龙江省1900年森林氧气释放量年总量空间分布 Fig.8 Spatial distribution of total annual oxygen release of the Heilongjiang Province in 1900

图选项

图 9 黑龙江省1949年森林氧气释放量年总量空间分布 Fig.9 Spatial distribution of total annual oxygen release of the Heilongjiang Province in 1949

图选项

图 10 黑龙江省1986年森林氧气释放量年总量空间分布 Fig.10 Spatial distribution of total annual oxygen release of the Heilongjiang Province in 1986

图选项

图 11 黑龙江省2009年森林氧气释放量年总量空间分布 Fig.11 Spatial distribution of total annual oxygen release of the Heilongjiang Province in 2009

图选项

可见，20世纪百年间黑龙江省全年森林氧气释放量呈减少趋势，2009年较1900年减少了5621.564 万t，占1900年的26.805%。全年森林氧气释放量空间分布趋势基本一致，均呈现南高北低趋势；但空间变化明显，1900年全省范围内均有森林覆盖，尤其是西南部有面积广泛的森林能释放氧气，其氧气释放量也最高，而到1949年后其森林基本消失，变为全省最低区域，基本为0；同样特征也出现在三江平原。

黑龙江省各行政区百年间森林氧气释放量除伊春微弱增加外，其余行政区均呈减少趋势(表 6)。其中大庆市、齐齐哈尔市减少显著，超过90%。将1900年及2009年各行政区氧气释放量进行显著性检验，其概率检验值为0.000，说明百年来黑龙江省森林氧气释放量减少极显著。

黑龙江省近百年间(1900—1986年)，由云冷杉林释放的氧气量减少了91.855%，红松林减少了78.905%，落叶混交林、落叶松林分别减少了25.043%、23.936%；针阔混交林、樟子松林氧气生产量有所增加。单位面积云冷杉林、红松林释放氧气量最多，落叶松林、樟子松林稍低(表 7)。

控制实验结果表明，O_2,contro1为25952.407 万t，O_2,2009为15350.343 万t，Δ_forest为-10602.060万t，相当于1900年氧气释放量的50.554%，说明百年来森林面积变化导致黑龙江省区域年氧气释放量减少了约50%。

将2009年氧气释放量分布图与控制试验模拟值分布图叠加相减，得到森林面积变化影响下森林氧气释放量空间变化图(图 12)。对比图 7，氧气释放量增加区域与森林面积不变或增加区域吻合，减小区域与森林面积减少区域吻合；而1900—2009年黑龙江省年平均温度提高了1.15 ℃，年平均辐射量增加0.0186 MJ m^-2 a^-1，大气中CO₂浓度增加了84 mg/m³，气候变化有利于森林氧气释放量增加，表 7中结果也说明各林型单位面积氧气释放量呈增加趋势)。由此可见，20世纪由于人类活动引起的森林面积减少，是影响黑龙江省森林氧气释放量减少的主要因素。

图 12 百年森林面积变化对森林氧气释放量的影响模拟 Fig.12 Simulated effects of forest area change during a hundred years on quantity of oxygen release

图选项

(1)较多学者采用碳氧平衡法对不同区域森林年氧气释放量进行了估算。如张颖^[22]等根据GPP估算出郑州市2010年森林年氧气释放量为24.07 t/hm²；彭江颖等、管东生等按NPP估算森林年释放氧气12.70、10.59 t/hm^{2[23, 24]}；苏继申等利用NPP估算阔叶混交林、针叶混交林年氧气释放量分别为12.52、13.51 t/hm^2[34]；邓毅等估算出常绿阔叶林、落叶阔叶林、北方针叶林年释放氧气为18.53、17.11、11.39 t/hm^2[35]；Nakazawa等认为阔叶林年放出氧气18.25 t/hm^2[36]；Berbigier等估算法国森林年氧气释放量为11.468 t/hm^2[37]；本文依据NEP模拟的黑龙江省森林年氧气释放量为8.055 t/hm²，不同林型介于6.049—8.697 t/hm²之间。由于研究区域及依据不同，以上结果在数量上虽有所差异，但均处于同一量级；基于GPP或NPP估算结果，由于没有充分考虑呼吸消耗，因此数量上比本文模拟结果偏大；不同学者采用不同方法模拟NPP(如采用森林蓄积量、生物量、文献查阅及经验估算等)，是产生氧气释放量差异的原因。本文模拟结果表明，云冷杉林、红松林年氧气释放量最多，这一结果与苏继申研究结果相吻合^[34]。

(2)近几年来，C-FIX模型较多地被国内外学者应用。如Veroustraete等使用NOAA/AVHRR的NDVI数据分别对比利时、全欧洲植被NPP和NEP进行了模拟^{[38, 39]}；卢玲使用SPOT/VEGETATION遥感数据估算了中国西部地区NPP和NEP^[26]；张冬有利用C-FIX模型较好地模拟了黑龙江省森林NPP^[28]。结果表明，C-FIX模型能较好地模拟区域植被各级生产力。本文采用C-FIX模型模拟了黑龙江省森林NEP，其精确性与森林氧气释放量模拟值的精度密切相关。因此，对模拟精度进行了检验。利用伊春1986年森林资源清查的蓄积量数据，在伊春区域随机抽取37个样区，根据蓄积量—生产力关系模型^[40]，计算出NPP，再采用东北地区土壤异养呼吸模型^[41]，计算出相应样区的NEP，并将单位换算为gC m^-2 a^-1。与本文相应样区NEP模拟值进行对比，均方差根(RMSE)为16.33，为原始值的3.8%—5.0%，模型拟合度较高。另外，由于黑龙江省年平均气温较低，本文模拟的黑龙江省森林NEP年均值为316.727 gC m^-2 a^-1，略低于国内模拟的各种林型NEP平均值(578 gC m^-2 a^-1)^[42]，结果合理。

(3)1900年、1949年和1986年黑龙江森林分布图是本文模拟森林生产力的基础图件，因此需对其矢量化精度进行检验。①20世纪80年代末黑龙江省政府及省测绘局对全省森林资源野外普查，将其分为96种林型并对其空间分布进行了详细记录。因此，以1986年为例，使用ArcGIS中Fishnet工具创建该年森林分布10 km×10 km格网，全省共生成10692个格网，其中含有林地格网共2088个，随机选择100个格网的森林类型与实际普查进行对比，吻合率近100%；②1949年、1986年矢量化的森林面积，与历史文献中^[43]森林清查数据分别相差-2.093万hm²、1.921万hm²，占总量的10%左右，误差较小。但历史文献中关于清末民初森林面积记载数据相差较大，如1917年满铁调查机构根据地面典型调查和实地调查，认为黑龙江省森林面积为3301.73万hm^2[44]，与本文矢量化结果相差4.315万hm²，而文献^[29]记载黑龙江省清末民初森林面积为4692.9万hm²，其面积大于黑龙江省行政区面积，也远大于1917年调查数据。分析其原因可能是历史文献中记载的清末民初森林面积，包括了山地坡地上的森林面积，而本文矢量化数据统计过程中，对坡地森林面积统计有所遗漏。③本文对1949年、1986年各行政区矢量化森林面积与文献^[29]记载数据也进行了对比，除大兴安岭和黑河市误差较大外，其他行政区误差均在10%之内，而历史文献中记载的大兴安岭及小兴安岭的森林面积均是估计值^[43]。综合以上结论，本文对黑龙江省历史时期森林分布矢量化精度是可信的。

(4)在模拟过程中，NDVI是输入参数之一，而1900年、1949年没有NDVI数据源，因此如何获取成为关键。陈云浩等指出，虽然NDVI是多种因素(气象、水文、地貌)综合作用的结果，但在中国东北，地形变化和植被类型是影响NDVI变化的主要因素，气候变化对中国东北大面积地带性植被NDVI变化没有明显影响^[44]。本文随机选择落叶松林、红松林、落叶混交林为代表，在1982—2009年8月NDVI图层中，随机各提取10个样本进行平均，分析得出近30年来森林NDVI年变化很小，变异系数分别为0.0028、0.0021和0.0028。因此，本文依据近代森林与NDVI数据源，提取不同纬度各林型NDVI，然后按同纬度同林型赋值到历史时期森林分布图中，获取历史时期森林NDVI。具体操作见2.3。

用两种方法对历史时期NDVI进行检验。①采用赋值法重建2009年森林NDVI，与网站下载同期NDVI比较。在每个纬度内随机选取同一林型NDVI样本共50个，其平均值相差很小，落叶松林、红松林、落叶混交林、樟子松林、云冷杉林、针阔混交林分别为-0.00019、-0.00125、-0.0001、-0.001333、0.0002、0.00347，其样本间标准差分别为0.002456、0.002174、0.000604、0.00265、0.005293、0.00264，说明样本间差异很小。②以1986年为例，分别按实际NDVI和赋值NDVI，计算了黑龙江省全年森林氧气释放量，仅相差12.369万t，占总量的0.084%。以上说明本文所提出的对历史时期森林NDVI赋值方法，对森林氧气释放量模拟值影响很小。

(5)本文没有考虑森林林分状况及林龄等森林性质对氧气释放量的影响，而森林处于不同林龄期，NDVI可能相同，在同样气候条件下，其碳积累速率及积累量有可能不同，最终会导致氧气释放量不同。这些生理因素对森林氧气释放量的影响需要今后更深入的研究加以说明。此外，本文仍然具有较多的不确定性。比如历史时期黑龙江省森林分布图，是目前能获取到的描述历史时期黑龙江省森林分布最详细的图件，但是该图的精度可能值得商榷；纸质图件矢量化及遥感影像解译过程中，会出现图形纠正、边界跟踪、误判等人为误差；在图层转化、叠加、裁切、拼接等过程中，也会出现界线对接、空间统计等误差。以上因素都会给模拟结果带来误差。

本文利用历史文献及遥感影像资料，重建了20世纪黑龙江省4个时间断面的森林空间分布图，基于C-FIX模型及碳氧平衡方法，模拟并分析了20世纪黑龙江省森林氧气释放量及其时空变化趋势，为研究过去百年黑龙江省森林变化及由此引起的生态环境效应提供了基础数据，为更进一步研究人类活动对区域环境变化的影响提供了依据。

黑龙江省百年森林面积呈持续减少趋势，森林总面积减少了10.667万km²，森林覆盖率减少了23.568%；各行政区森林面积均呈减少趋势，大庆市森林面积减少最多，约减少了94.673%，空间上以松嫩平原和三江平原减少显著；除樟子松林、针阔混交林面积有所增加外，其它林型面积均减少，其中云冷杉林、红松林减少率在80%以上；黑龙江省百年来森林年氧气释放量呈减少趋势，共减少5621.564万t，约占1900年的26.805%；百年来各行政区森林年氧气释放量除伊春微弱增加外，其余行政区均呈减少趋势，其中大庆市、齐齐哈尔市减少显著，减少了90%以上；人类活动引起的森林面积减少是影响黑龙江省森林氧气释放量减少的主要因素。