楚雄州青蚕豆种植的气候及变化特征分析

Analysis of Climate and Change Characteristics of Green Broad Bean Cultivation in Chuxiong Prefecture

摘要: 青蚕豆是以采收鲜荚为目的蚕豆的统称。近年楚雄州农科院育成的粮菜兼用型蚕豆品种彝豆3号，在秋种冬末菜荚销售中，增效显著，占楚雄州青蚕豆种植面积的60%。为充分利用当地独特的气候资源，科学种植、降低气候风险，本文用气候倾向分析、Mann-Kendall突变检验、小波周期分析和资料对比法分析了楚雄州早种植主生育期气温、降水、日照，≥10℃积温、产量形成关键期气温、降水和低温时段气温的变化特征。结果表明，主生育期气时段平均温呈0.2℃/10a显著上升趋势变化(P < 0.01)，2005年、2007年发生突变，存在45a左右主周期、6~8a次周期变化；积温呈0.81℃/10a显著上升趋势变化，2008年发生突变；降水量呈−3.683 mm/10a不显著(P > 0.1)减少趋势变化，突变多发生在1960~1970年代和2021年，存在6~8a左右主周期、28a次周期变化；日照时数呈−6.116 h/10a不显著(P > 0.1)减少趋势变化，突变多发生在1980年代中期之前。产质量形成关键生育期时段平均气温呈0.209℃/10a显著上升趋势变化(P < 0.01)，2008年、2010年发生突变；时段降水量呈0.661 mm/10a不显著(P > 0.1)微增加趋势变化，突变与主生育同步发生，存在4a左右的主周期、11~12a的次周期变化。最冷时段平均气温呈0.38℃/10a显著上升趋势变化(P < 0.01)；平均最低气温呈0.65℃/10a呈显著上升趋势变化(P < 0.01)。

Abstract: Green broad bean is a general term for broad bean for the purpose of harvesting fresh pods. n recent years, Yidou No. 3, a variety of both grain and vegetable broad bean developed by Chuxiong Academy of Agricultural Sciences, has significantly increased its efficiency in the sales of vegetable pods at the end of autumn and winter, accounting for 60% of the planting area of green broad bean in Chuxiong Prefecture. In order to make full use of local unique climate resources, scientific planting and climate risk reduction, In this paper, With climate tendency analysis, Mann-Kendall mutation test, wavelet cycle analysis and data comparison methods to analyze the changes of tempera-ture, precipitation, sunshine, temperature ≥10˚C accumulated temperature, key period of yield formation, temperature, precipitation and temperature in low temperature period. The results show that the average gas temperature in the main growth stage showed a significant upward trend of 0.2˚C/10A (P < 0.01), with a sudden change in 2005 and 2007. Here are about 45a main cycle and 6~8a cycle changes. The accumulated temperature showed a significant upward trend of 0.81˚C/10A, with a sudden change in 2008; Precipitation showed an insignificant decreasing trend of −3.683 mm/10a (P > 0.1), and the mutation occurred mostly in 1960~1970s and 2021, There are about 6~8A main cycle and 28a cycle changes; The sunshine duration showed an insignificant decreasing trend of −6.116 h/ 10A (P > 0.1), and the mutation occurred mostly before the mid-1980s. The average temperature during the key growth period of yield and quality formation showed a significant upward trend of 0.209˚C/10A (P < 0.01), and the mutation occurred in 2008 and 2010, The precipitation in the period was 0.661 mm/10a (P > 0.1), and the mutation occurred simultaneously with the main fertility, There is a main cycle of about 4a and a subcycle of 11~12a; The average temperature in the coldest period showed a significant upward trend of 0.38˚C/10A (P < 0.01), The average minimum temperature was 0.65˚C/10A and showed a significant upward trend (P < 0.01).

1. 引言

楚雄州位于云南中北部100˚43'~102˚30'E，24˚13'~26˚30'N之间，地处横断山脉和云贵高原的过渡地带。境内地形复杂，山峦起伏，海拔高差达3101 m；全州辖8县2市，年平均气温 16.3 ℃ ，年平均日照时数2378.9 h，年平均降水量850.0 mm，其中雨季5~10月降水占年降水量的89%。气温年较差小、日较差大、雨热同季，干湿分明、立体性强的中亚热带低纬高原半干旱季风气候特征明显 [1]。秋冬降水适中，温暖多日照的气候，为青蚕豆推广种植提供了基础资源 [2]。蚕豆(Vicia faba L.)，别名胡豆、佛豆、罗汉豆、大豆等，原产于温带，性喜温暖而湿润的气候，属于较耐寒作物 [3]。彝豆3号是楚雄州农科院以天杂30为母本，大姚白花豆为父本经系谱单株选育而成的粮菜兼用型的秋播中早熟型蚕豆新品种 [4]。彝豆3号新品种商品性好，品质优良，中抗锈病，抗赤斑病，抗褐斑病，全生育期180 d左右，鲜食青荚生育期110 d左右。生长特点为前期慢，后期快，历花期长，株型紧凑，茎秆粗壮，株高114 cm左右，叶片卵圆形，叶色淡绿色，黑白花，总茎枝数3.1个，平均单株有效枝2.7个，单株平均实荚12.2荚，单荚粒数1.85粒，鲜荚长9.2 cm宽2.2 cm。鲜籽粒浅绿色，百粒重321.8 g左右，干籽种皮白色、白脐，种皮破裂率为零，百粒重143 g左右。每667 m²产鲜荚1000~1300 kg (干籽粒200~260 kg)，适宜在云南省海拔1400~2200 m的蚕豆种植区种植 [5] [6]。

2. 资料来源

2.1. 资料来源

本文所用气象资料全部来源于云南省气象台开发的(WFOIS)气象服务平台，云南省气象台(MDIS)气象资料管理系统地面气象报表文件及楚雄州自主研发的气候资源管理、咨询与评价系统分析软件平台查询，多年平均值部分来自国家气候中心资料共享。

2.2. 分析因子的选取

根据云南省云豆、楚雄彝豆系列育种项目成果和云南推广青蚕豆特色种植经验 [2] [4]，本文选取对青蚕豆生育期影响权重大的主生育期(10月1日至次年2月)和产质量形成关键期(11月21日至次年2月10日)的时段期间气温、降水、日照及对青蚕豆产量形成有突出制约作用的最冷时段气温等共7个统计统计时段气象因子。

3. 分析方法

3.1. 气候趋势(倾向)分析方法 [7]

气候趋势(倾向)分析法，是分析气象要素随时间变化趋势倾向和变化幅度特征，x_i为样本量为n的气象要素变量，t_i为x_i所对应的时间，用一元一次回归方程描述

$x_i=a+b{t}_i$ ( $i=1,2,3,4,\cdots ,n$ ) (1)

式中：a为回归常数；b为回归系数，a和b可用最小二乖法进行估算，趋势系数用相关系数进行显著性检验。

3.2. K-M检验法(Mann-Kendall) [8]

$S_k={\displaystyle {\sum }_{i=1}^k{\displaystyle {\sum }_j^{j-1}{\alpha }_{ij}}}$ ( $k=1,2,3,4,\cdots ,n$ ) (2)

其中， ${\alpha }_{ij}=\{\begin{array}{l}1x_i>x_j\\ 0x_i\le x_j\end{array}$ $1\le j\le I$ 定义统计变量

$U{F}_k=\frac{\left|s_k-E\left(s_k\right)\right|}{\left|\sqrt{V_{\alpha \gamma }\left(S_k\right)}\right|}$ ( $k=1,2,3,4,\cdots ,n$ ) (3)

式中： $E\left(S_k\right)=k\left(k+\text{1}\right)/4$ ； $V_{ar}\left(S_k\right)=k\left(k-1\right)\left(2k+5\right)/72$。

UF_k为标准正态分布，给定显著水平а ( $а=0.05$, $u_{0.05}=1.96$ )，若 $\left|U{F}_1\right|\ge U_{a/2}$，则表明系存在明显的趋势变化。将时间系列按逆系排列，再按照上式计算，同时使

$\{\begin{array}{l}U{B}_k=-U{F}_k\\ k=n+1-k\end{array}$ ( $k=1,2,3,4,\cdots ,n$ ) (4)

通过分析统计量UF_k和UB_k可进一步分析系列X的趋势变化，并明确突变时间，指出突变的区域。若UF_k值大于0则表示系列呈上升趋势，小于0则表示呈下降趋势；当它们超过临界值时，表示上升或下降趋势显著。若UF_k和UB_k这两条曲线出现交点，且交点在临界值线之间，则交点对应的时刻就是突变开始时刻。

3.3. 小波分析基本原理 [8] [9]

小波分析的基本思想是用一簇小波函数系来表示或逼近某一信号或函数。因此，小波函数是小波分析的关键，它是指具有震荡性、能够迅速衰减到零的一类函数，即小波函数 $\psi \left(t\right)\in L^2\left(R\right)$ 且满足：

${\displaystyle {\int }_{-\infty }^{+\infty }\psi \left(t\right)\text{d}t=0}$ (5)

式中： $\psi \left(t\right)$ 为基小波函数，它可通过尺度的伸缩和时间轴上的平移构成一簇函数系：

${\psi }_{a,b}\left(t\right)={\left|a\right|}^{-1/2}\psi \left(\frac{t-b}{a}\right)$ 其中， $a,b\in R$， $a\ne 0$ (6)

式中： ${\psi }_{a,b}\left(t\right)$ 为子小波；a为尺度因子，反映小波的周期长度；b为平移因子，反应时间上的平移。若 ${\psi }_{a,b}\left(t\right)$ 是由式(2)给出的子小波，对于给定的能量有限信号 $f\left(t\right)\in L^2\left(R\right)$，其连续小波变换(Continue Wavelet Transform，简写为CWT)为：

$W_f\left(a,b\right)={\left|a\right|}^{-1/2}{\displaystyle {\int }_{R}f\left(t\right)\bar{\psi }\left(\frac{t-b}{a}\right)\text{d}t}$ (7)

式中： $W_f\left(a,b\right)$ 为小波变换系数； $f\left(t\right)$ 为一个信号或平方可积函数；a为伸缩尺度；b为平移参数； $\bar{\psi }\left(\frac{x-b}{a}\right)$ 为 $\psi \left(\frac{x-b}{a}\right)$ 的复共轭函数。地学中观测到的时间序列数据大多是离散的，设函数 $f\left(k\Delta t\right)$，( $k=1,2,\cdots ,N$ ； $\Delta t$ 为取样间隔)，则式(3)的离散小波变换形式为：

$W_f\left(a,b\right)={\left|a\right|}^{-1/2}\Delta t{\displaystyle \underset{k=1}{\overset{N}{\sum }}f\left(k\Delta t\right)\bar{\psi }\left(\frac{k\Delta t-b}{a}\right)}$ (8)

由式(3)或式(4)可知小波分析的基本原理，即通过增加或减小伸缩尺度a来得到信号的低频或高频信息，然后分析信号的概貌或细节，实现对信号不同时间尺度和空间局部特征的分析。将小波系数的平方值在b域上积分，就可得到小波方差，即

$Var\left(a\right)={\displaystyle {\int }_{-\infty }^{\infty }{\left|W_f\left(a,b\right)\right|}^2\text{d}b}$ (9)

小波方差随尺度a的变化过程，称为小波方差图。由式(5)可知，它能反映信号波动的能量随尺度a的分布。因此，小波方差图可用来确定信号中不同种尺度扰动的相对强度和存在的主要时间尺度，即主周期。

4. 结果与分析

4.1. 楚雄州青蚕豆主生育期气温、降水和日照倾向变化及突变分析

青蚕豆种植是近来年云南省开发利用低纬高原秋冬降水适中，冬春温暖、光照丰富的独特气候资源，发展特色种植，促农业增收的一项新、特型种植模式。该种植模式对气候资源的利用特点是，光合有效性高，制约因素是产量形成关键期和最冷时段的气温、降水和日照时长变化。在全球气候气候变暖背景影响下，青蚕豆生育期几个主要影响时段统计结果表明，各统计时段平均气温和积温呈上升趋势变化、主生育期时段平均降水量呈减少趋势变化显著，日照时数变化趋势不显著。61年来，最远10年(1961~1970年)与最近10年(2011~2020年)青蚕豆主生育期时段平均气温上升1.0℃、积温增加300℃∙d、降水量29 mm、日照时数21小时，产量形成关键期时段平均气温1.0℃、降水量增加4 mm，最冷时段平均气温上升1.7℃、平均最低气温上升3.2℃。见表1。

Table 1. Period average of temperature, precipitation and sunshine during the growth period of green broad bean in Chuxiong Prefecture (unit: ˚C, ˚C∙d, mm, h)

表1. 楚雄州青蚕豆生育期气温降水日照年代纪平均(单位：℃、℃∙d、mm、h)

4.1.1. 主生育期气温、降水、日照倾向变化及突变分析

楚雄州青蚕豆主要生长期(10月至次年2月)气温、降水和日照的年际变化趋势倾向。图1左a，近61年来，青蚕豆主要生育期时段平均气温以0.2℃/10a的倾向率呈显著上升趋势变化，趋势系数0.565，通过了99%的信度检验(P < 0.01)；图1左b，时段平均降水量呈减少趋势变化，变化不显著，变化倾向率−3.683 mm/10a，趋势系数0.1253，未通过90%的信度检验(P > 0.1)；图1左c，时段平均日照时数呈微减少趋势变化，变化不显著，变化倾向率−6.116 h/10a，趋势系数0.1170，未通过90%的信度检验(P > 0.1)。

UF曲线分析，图1右a近61年来，青蚕豆主要生育期时段平均气温1961~1965年上升变化趋势显著，1966~1992年呈缓慢下降变化，1993~2021年呈显著上升，其中2005年和2007年发生了突变，曲线相交于信度 ± 1.96 (α = 0.05)；图1右b时段降水量1962~1972年降水量呈微上升变化，1973年~1988年降呈微减少变化，1989年~1992年呈增加变化，1993年~2013年呈微减少变化，2014年~2021年呈微增加变化，其中1967年、1969年、1971年、1973年、1975年和2021年发生突变，曲线相交于信度 ± 1.96 (α = 0.05)；图1右c时段日照时数1961年~1969年呈跳跃式增加变化，1970年~1977年呈阶梯式减少变化，1978年~1980年呈增加变化，1981年~1988年呈“V”字型变化，1989年~1996年呈减少变化，1997年~2003年呈微增加变化，2004年~2008年呈微减少变化，2019年~2014年呈增加变化，2015年~2021年呈“V”字型变化，其中1975年、1977年、1981年，1985年发生了突变，曲线相交于信度 ± 1.96 (α = 0.05)。

Figure 1. Left a, b, c are the evolution of temperature, precipitation and sunshine during the growth period of broad bean, respectively; Right d, e, f, show the UF curve analysis of temperature, precipitation, and sunshine mutation during the growth period of broad bean, respectively

图1. a，b，c，分别为蚕豆生育期气温、降水量和日照变化演变；d，e，f，分别为蚕豆生育期气温、降水量和日照突变的UF曲线分析

4.1.2. 主生育期气温、降水周期变化分析

楚雄州青蚕豆主生育期(10月至次年2月)气温、降水和日照时数的时间尺度上振荡周期分析。图2(a)，图2(b)主生育期时段平均气温存在2~4a、6~8a、11a、16~17a和45a左右的周期振荡，以45a左右的周期振荡最强，为第一主周期，其它依次为6~8a、2~4a、11a和16~17a，45a左右周期在2000~2021年与16~17a周期叠加，16~17a左右周期在1962~1974年与6~8a周期叠加，6~8a周期振荡稳定少变，16~17a左右周期振荡在1975年至今最明显。图2(c)，图2(d)主生育期时段降水存在2~3a、5a、6~8a和28a左右的周期振荡，以6~8a的周期振荡最强，为第一主周期，其它依次为28a、2~3a、5a，除第一主周期外，28a左右的周期振荡在61年的分析期间稳定少变。

Figure 2. a, b Analysis of temperature cycle changes during the main growth period; c, d Analysis of precipitation cycle changes during the growth period of green broad bean

图2. a，b为主生育期气温周期变化分析；c，d青蚕豆生育期降水周期变化分析

4.2. 楚雄州青蚕豆产质量形成关键期气温、降水倾向变化及突变分析

4.2.1. 产质量形成关键期气温、降水、日照倾向变化及突变分析

楚雄州青蚕豆产质量形成关键生育期(11月21日至2月10日)气温、降水和日照的年际变化趋势倾向。图3(a)近61年来，青蚕豆产质量形成关键生育期时段平均气温以0.209℃/10a的倾向率呈显著上升趋势变化，趋势系数0.4922，通过了99%的信度检验(P < 0.01)；图3(b)时段平均降水量呈微增加趋势变化，变化不显著，变化倾向率0.661 mm/10a，趋势系数0.057，未通过90%的信度检验(P > 0.1)。

UF曲线分析，图右上b近61年来，青蚕豆产量形成关键期平均气温1961~1966年上升变化趋势显著，1967~1983年呈缓慢下降变化，1984~1992呈“n”字型微变化，1993~2021年呈显著上升，其中2008年和2010年发生了突变，曲线相交于信度 ± 1.96 (α = 0.05)；图右下b时段降水量1962~1972年降水量呈微上升变化，1973年~1988年降呈微减少变化，1989年~1994年呈增加变化，1995年~2013年呈阶梯式微减少变化，2014年~2021年呈微增加变化，其中1967年、1969年、1971年、1973年、1975年和2021年发生突变，曲线相交于信度 ± 1.96 (α = 0.05)。

4.2.2. 产量形成关键期气温、降水周期变化分析

楚雄州青蚕豆产量形成关键期(11月21日至2月10日) (10月至次年2月)气温、降水时间尺度上的振荡周期分析。图4(a)，图4(b)，时段平均气温存在2a、7a、10~11a、17a和45a左右的周期振荡，以7a左右的周期振荡最强，为第一主周期，其它依次为45a，2a、10~11a和17a，61年间，7a和45a周期振荡稳定少变，17a周期振荡在1975年至今最明显；图4(c)，图4(d)时段降水存在2a、4a、8a、11~12a和21~23a左右的周期振荡，以4a左右的周期振荡最强，为第一主周期，其它依次为11~12a、2a、21~23a、8a，61年间，21~23a周期振荡在稳定少变，11~12a周期振荡在1990~2007a表现最明显。

Figure 3. a, b show the evolution of temperature and precipitation during the critical growth period; c, d show the abrupt change of temperature and precipitation during the critical growth period, respectively

图3. a，b为生长关键期气温、降水演变；c，d为生长关键期气温、降水突变

Figure 4. a, b Green broad bean yield formation key period temperature cycle changes; c, d Changes in the precipitation cycle during the critical period of green broad bean yield formation

图4. a，b青蚕豆产量形成关键期气温周期变化；c，d青蚕豆产量形成关键期降水周期变化

4.3. 楚雄州青蚕豆主生育期积温和最冷时段气温倾向及突变分析

楚雄州青蚕豆主生育期(10月至次年2月)时段积温及其变化是青蚕豆种植热量适宜度的权衡指标，与产质量形成成正比，最冷时段(12月15日至次年1月15日)气温对青蚕豆产量的形成影响权重最大，低温天气危害是青蚕豆种植中最大气象风险，严重的低温持续天气可致绝收 [2]。

4.3.1. 青蚕豆主生育期积温变化特征分析

楚雄州青蚕豆主要生长期(10月至次年2月)积温的年际变化趋势倾向。图5(a)近61年来，青蚕豆主要生育期时段积温以0.81℃/10a的倾向率呈显著上升趋势变化，趋势系数0.611，通过了99%的信度检验(P < 0.01)。

UF曲线分析，近61年来，青蚕豆主要生育期时段积温1961~1965年上升变化趋势显著，1966~1982年呈缓慢下降变化，1983~1992年小幅微升降变化，1993~2021年呈显著上升，其中2005年和2007年发生了突变，曲线相交于信度 ± 1.96 (α = 0.05)，图5(b)。

Figure 5. a Evolution of accumulated temperature during the growth period of green broad bean; b Mutation of accumulated temperature during growth period of green faba bean

图5. a青蚕豆生育期积温演变；b青蚕豆生育期积温突变

4.3.2. 青蚕豆生育期最冷时段气温变化及检验分析

楚雄州青蚕豆生育期最冷时段(12月12日至次年1月15日)平均气温以0.38℃/10a的倾向率呈显著上升趋势变化，趋势系数0.587，通过了99%的信度检验(P < 0.01)，图6(a)；平均最低气温以0.65℃/10a的倾向率呈显著上升趋势变化，趋势系数0.742，通过了99%的信度检验(P < 0.01)，图6(b)。

Figure 6. a Evolution of average temperature at the coldest time during the growth period of green faba bean; b Evolution of average minimum temperature during the coldest period of the growth period of green faba bean

图6. a青蚕豆生育期最冷时段平均气温演变；b青蚕豆生育期最冷时段平均最低气温演变

5. 讨论与结论

5.1. 讨论

1) 本文用气象资料统计对比、气候变化趋势、要素突变和小波分析方法，定量和定性分析了1961~2021年间，彝豆3号青蚕豆种植生主育期(10月至次年2月)，产量形成关键期(11月21日至2月10日)和最冷时段(12月11日至次年1月15日)时段的平均气温、最低气温、积温、降水和日照时数的变化特征。分析结论为青蚕豆种植规划、栽培技术和自然灾害防御方案制定等提供了决策依据。

2) 楚雄州青蚕豆种植是充分利用秋粮后茬残余肥力、适宜的土壤墒情、丰富的冬春光热资源及适中的秋冬降水，推广以提高农业农、农村和农民收入为宗旨的创新性特色种植新模式。彝豆3号是楚雄州农科院自主新近选育的优质本地青蚕豆新品种，气候适应性强、商品性好，抗病性强，产量稳定，推广潜力大，是青种植蚕豆类当家推广品种。

3) 本文用统计对比、气候变化趋势和要素突变分析法分析了楚雄州青蚕豆种植的主生育期、产量形成关键期和最冷时段气温、积温的变化，结论得出：各统计时段热量呈显著上升趋势变化，降水量和日照时数呈微减少趋势变化，验证了楚雄州秋冬气候变暖趋势显著与全球气候变暖一致 [10]，区域气候变化对楚雄州推广青蚕豆种植影响呈正相关。气候变化事实和分析结论预示着未来青蚕豆稳产增收的气象风险随年际变化下降，资源基础更优、种植增效和推广潜力大。

5.2. 结论

1) 10年统计时段分析对比，61年中，青蚕豆主生育期各要素前40年变幅小，近20年变幅大，其中平均气温1961~1970年比2011~2020年上升1.0℃、积温增加300℃∙d、降水量减少29 mm、日照时数减少21小时，产量形成关键期时段平均气温上升1.0℃、降水量增加4 mm，最冷时段平均气温上升1.7℃、平均最低气温上升3.2℃。

2) 气候倾向分析和Mann-Kendall突变检验分析结论，主生育期(10月至次年2月)时段平均气温呈0.2℃/10a显著上升趋势变化(P < 0.01)，2005年、2007年发生突变；积温以0.81℃/10a呈显著上升趋势变化，趋势系数0.611 (P < 0.01)，2005年和2007年发生了突变；降水呈−3.683 mm/10a不显著减少趋势变化(P > 0.1)，1967年、1969年、1971年、1973年、1975年和2021年发生突变；日照时数呈−6.116h/10a不显著减少趋势变化(P > 0.1)，1975年、1977年、1981年，1985年发生突变。产量形成关键期(11月21日至2月10日)平均气温呈0.209℃/10a显著上升趋势变化(P < 0.01)，2008年、2010年发生突变；降水量呈微增加趋势不显著变化，变化倾向率0.661 mm/10a，趋势系数0.057 (P > 0.1)，1967年、1969年、1971年、1973年、1975年和2021年发生突变。

3) 小波周期变化分析，主生育期(10月至次年2月)时段平均气温存在45a左右主振荡周期及6~8a、2~4a、11a、16~17a的次振荡周期变化，降水量存在6~8a的主振荡周期及28a、2~3a、5a的次振荡周期变化；产量形成关键期(11月21日至2月10日)平均气温存在7a左右主振荡周期及45a，2a、10~11a和17a次振荡周期变化，降水量存在4a左右的主振荡周期及11~12a、2a、21~23a和8a的次振荡周期变化。

4) 气候资料对比、气候趋势变和突变分析结论均验证了楚雄州青蚕豆各生育期气温、积温呈上升趋势变化，降水量呈不显著微下降趋势变化，突变均发生在2008年前后，且突变次数多，特别是1960年代至1970年代前中期，日照时数呈不显著微下降趋势变化，突变次数较多，主要发生在1970年代中至1980年代中期，与全球气候变暖趋势一致。

项目来源

云南省陈新专家工作站；项目编号：202105AF150008；云南省气象局基层台站气象科技创新与能力提升计划项目(STIAP)：项目编号：STIAP202222。

参考文献

NOTES

^*通讯作者。