随机物理过程扰动方案在克拉玛依区域集合预报中的应用研究

史永强,张涵斌,刘郁珏,张歆然

强天气现象具有很强的突发性和局地性,其发生发展的动力和热力机制极为复杂,尤其是复杂地形区域的极端降水、低温、大风等极端天气事件预报不确定性较大,导致目前的区域数值模式对此类天气现象的预报准确性较低(任福民等,2014),因此针对此类强天气的预报技术亟待发展。针对数值模式本身的改进(如提高模式分辨率)可以一定程度上改善对此类天气的预报效果(Xue et al.,2020),但是考虑到局地强天气的预报不确定性,需要进一步发展集合概率预报技术(Epstein,1969;Leith,1974,陈静等,2005)来提高极端强天气预报的准确率(Zhang et al.,2021)。

集合预报扰动技术是集合预报研究的一个难点,对于区域集合预报而言,需要考虑初值、模式以及侧边界等多种不确定性来源。集合预报初值扰动方法首先在全球集合预报开始应用,到目前已获得了较大发展并形成了一系列典型的扰动技术,如增长模繁殖法(Breeding Growing Mode,简称BGM;Toth and Kalney,1993,1997),奇异向量法(Singular Vectors,简称SVs;Molteni et al.,1996)以及集合变换卡尔曼滤波法(Ensemble Transform Kalman Filter,简称ETKF;Wang and Bishop,2003;Bowler et al.,2008)等,此类方法也有效地沿用到了区域集合预报中(Stensrud et al.,1994,1999;Du et al.,2003;Walser et al.,2006;Bishop et al.,2009;Bojarova et al.,2011)。

除了初值影响以外,数值预报不准确很大程度上是模式不确定性导致的,因此不管是全球集合预报还是区域集合预报均需要发展能够有效代表模式不确定性的扰动方法 (关吉平和张立凤,2009;王璐和沈学顺,2019)。在中尺度环流的发展过程中,次网格物理过程对于描述相态转换、动量交换、对流触发等起到十分重要的作用(Stensrud and Fritsch,1994;陈静等,2005),因此目前的模式扰动技术主要针对次网格物理过程开展,常用的模式扰动方法包括:1)对不同物理过程选项进行随机组合 (Houtekamer et al.,1996;Stensrud et al.,2000);2)在模式参数化方案中加入随机强迫来体现次网格物理过程的不确定性,如随机物理过程倾向方案(Stochastically Perturbed Parameterization Tendencies,SPPT;Buizza et al.,2007;Lin and Neelin,2000;Li et al.,2008,任志杰等,2011;谭宁等,2013)、随机动能后向散射方法(Stochastic Kinetic Energy Backscatter,SKEB;Shutts,2005;Berner et al.,2009,2011;Bowler et al.,2009);3)针对特定物理过程的随机参数进行扰动(Stochastically Perturbed Parameterization,SPP;Bowler et al.,2009;徐致真等,2019)。其中,SPPT方法由于理论基础较科学且易于实现,在集合预报中应用较为广泛(Buizza et al.,2007)。

SPPT方法自提出以来在国内外开展了一系列的改进研究,如Li et al.(2008)基于一阶自回归模式随机型(Lin and Neelin,2000)进行空间球谐函数展开和傅立叶展开,发展了时空相关的SPPT随机扰动方案;Palmer et al.(2009)对ECMWF集合预报系统中的SPPT随机扰动型进行了更新,构建了单变量的高斯分布随机扰动型,并通过谱展开使随机型在时空上更连续。对于SPPT的应用也从全球集合预报扩展至区域集合预报以及对流尺度集合预报中,如最初该方法于1998年在ECMWF业务全球集合预报系统中获得应用,可以有效增加该系统的集合离散度并降低均方根误差(Buizza et al.,2007);中国气象局地球系统数值预报中心最早基于T213全球集合预报开展了SPPT的应用,表明该方案可有效体现全球集合预报中的模式不确定性,并提高集合离散度(任志杰等,2011;谭宁等,2013);中国气象局地球系统数值预报中心进一步基于CMA-Meso模式开展了SPPT相关工作并应用于CMA区域集合预报系统中,改进了预报技巧(袁月等,2016;陈静和李晓莉,2020)。有学者尝试在更高分辨率的集合预报系统中应用SPPT技术,如李俊等(2015)在对流尺度集合预报中开展了SPPT试验,结果表明SPPT方案对降水量级改变较大,但对降水落区的改善有限;闵锦忠等(2018)基于风暴尺度集合预报对江淮流域一次强降水过程探讨了SPPT的作用,表明天气系统的维持时间与失相关时间尺度关系密切,而天气过程及模式分辨率受失相关空间尺度影响较大;王明欢等(2021)研究了SPPT方案在西部山地对流尺度集合预报中的应用效果,探索了适用于该系统的SPPT参数。综合国内外的研究表明SPPT的引入会对集合预报产生正效果,但也需要针对自身系统特点对SPPT的相关参数进行设置(Qiao et al.,2018;Lupo et al.,2020)。

克拉玛依地处准噶尔盆地西北边缘,西有加依尔山,中部、东部为平坦戈壁,其城区主要受背风坡的强下坡风影响,使得大风成为该地区主要气象灾害之一(孙东霞等,2008;辛渝等,2015)。目前克拉玛依气象局业务运行的区域集合预报系统为当地提供了有效的业务支撑,但由于系统仅采用ETKF(Wang and Bishop,2003)初值扰动方案,导致集合成员之间的发散度不够,对极端大风预报的把握能力欠佳。本文尝试针对该集合预报系统进一步增加扰动源,引入SPPT扰动方法,探索适用于克拉玛依区域集合预报的SPPT参数,以期进一步提高集合预报技巧。本研究不仅对克拉玛依区域集合预报的发展具有重要意义,也可为业务集合预报构建提供新方法新思路,具有较好的应用前景。

1.1 克拉玛依区域集合预报系统简介

克拉玛依气象局发展的中尺度区域集合预报系统为日常业务预报工作提供了重要参考,该系统采用WRF(The Weather Research and Forecasting Model)中尺度模式V4.1.2版,模式区域设置为水平分辨率15 km,模拟区域范围共239×180个格点,覆盖中亚及新疆大部分区域(如图1),垂直层次为57层模式面,模式层顶为50 hPa。该系统21个集合成员,每天从00时和12时(世界时,下同)起报两次,预报时效36 h。所有成员物理过程参数化方案设置为Thompson微物理方案,YSU(Yonsei University)边界层方案以及RRTMG(The Rapid Radiative Transfer Model for GCMs)长短波辐射方案,积云对流方案为KF(Kain-Fritsch)方案。

图1 克拉玛依区域集合预报系统范围设置(阴影为地形高度,单位:m;蓝线为克拉玛依区域范围)

目前该系统采用美国国家环境预报中心(National Centers for Environmental Prediction,NCEP)全球预报资料(Global Forecast System,GFS)作为驱动场,为该集合预报系统提供侧边界条件。初值扰动采用ETKF初值扰动方案生成(Wang and Bishop,2003),ETKF扰动循环间隔为12 h,即ETKF每天2次生成初值扰动进行预报。

1.2 随机物理过程倾向方法介绍

根据随机物理过程倾向SPPT扰动方法的原理,模式预报中的不确定性很大程度上来源于物理参数化过程和模式截断误差,且这种误差主要是随机的,为表达这种次网格物理过程的随机误差,在所有物理过程参数化方案计算结束后的总倾向项上施加随机强迫来实现。

根据Berner et al.(2011),WRF中SPPT数学模型:

(1)

右边两项为动力积分过程倾向和物理过程倾向,x为不同变量:

x∈{u,v,T,q}。

(2)

SPPT方案对不同变量采用相同的随机型r(x,y,t)进行扰动,随机型相同可有效避免不同变量计算不协调导致溢出。随机扰动场r(x,y,t)是一个具有时空相关性,取值范围[-1,1]的高斯随机函数。

将随机场r(x,y,t)进行二维傅里叶展开可得:

(3)

其中:K,L分别为纬向和径向波数;t为积分时间;e2πi(kx/X+ly/Y)是由在0

rk,l(t+Δt)=(1-α)rk,l(t)+gk,lεk,l(t)。

(4)

其中:(1-α)为线性回归参数;gk,l为扰动振幅;εk,l(t)为复值高斯噪声过程。

对于扰动振幅gk,l的计算,根据王明欢等(2021),其具有空间自相关特征且受谱空间每个格点的标准差影响,格点标准差越大,扰动振幅gk,l越大。

综合以上可知随机型r(x,y,t)与时间和空间相关且扰动振幅可调,因此SPPT扰动随机场主要由失相关时间尺度,失相关空间尺度以及格点标准差决定。较大/较小的失相关时间尺度表示随机型的变化频率较慢/较快;较大/较小的失相关空间表示随机型的水平空间尺度较大/较小;较大/较小的格点标准差表示扰动振幅较大/较小。

1.3 试验方案设计

SPPT随机物理过程扰动方案涉及多个可调参数,因此需要开展相关敏感试验,确定所产生的模式随机扰动对不同参数的敏感性,探索适用于克拉玛依区域集合预报系统的扰动参数。从SPPT方案中三个关键参数入手:失相关时间尺度,失相关空间尺度,格点标准差。WRF模式中给出了一组参数默认值,即失相关时间尺度6 h,失相关空间尺度150 km,格点标准差0.5,敏感性试验即可针对默认值进行增减,表1给出了SPPT的7组敏感性试验的配置,其中P6-150-0.5即为默认参数方案,其他方案均在其基础上对参数进行增减试验,此外单独进行一组无SPPT方案的控制试验(CTRL)作为对比。

通过SPPT敏感试验获得适用于克拉玛依区域集合预报系统的SPPT参数后,基于该集合预报系统设计了两种试验方案,方案一仅采用ETKF初值扰动并作为参照试验,称为ETKF方案,方案二采用ETKF初值扰动与SPPT相结合的方案,称为ETKF-SPPT方案。两组试验设置见表2,试验时段选取2021年2月14日—3月15日连续一个月,试验时段内每组试验每天进行两次起报(00时、12时),每次起报预报时效为36 h。其中2021年2月25日在克拉玛依城区发生了强大风过程,作为批量试验中的重点分析个例。

表1 SPPT引入单模式预报敏感试验方案配置

表2 集合预报对比试验设置

1.4 试验数据

本文集合预报背景场和侧边界条件来自NCEP-GFS全球分析场和预报场,水平分辨率为0.5°×0.5°,该资料每天00、12时获取2次,预报间隔6 h。本文中高空和地面要素检验采用全国探空和地面观测站。

2.1 SPPT参数敏感试验

首先对不同参数设置下SPPT敏感试验的效果进行分析,定义均方根离差为各敏感试验与未加SPPT扰动的控制试验(CTRL)的差异,即

(5)

图2给出了不同参数设置下SPPT试验与控制试验的离差的演变,以评估扰动增长对不同参数的敏感性。由图2首先可以看出所有的扰动试验与控制试验的差异均随着预报时效的增加而增长,但增长程度也存在明显差异,对于500 hPa高空要素(图2a—c)而言,所有方案整体的离差在24 h内均具有增长特征;对于850 hPa低空要素(图2d—f)而言,前18 h的离差处于增长态势,但18 h之后不同方案的离差增长达到饱和,说明SPPT带来的扰动增长在不同高度效果有所差异。对比不同试验参数设置下的扰动试验可以看出,首先对于格点标准差参数,不管是高空还是低空要素,离差最小的是格点标准差降低到0.2的方案,如对于850 hPa纬向风(U850,图2d)24 h离差为0.7 m/s,格点标准差0.55试验的离差明显增长,24 h离差可达1.5 m/s;对于失相关空间尺度敏感试验对比可以看出,不同要素受空间尺度影响程度不同,如500 hPa纬向风(U500)、500 hPa温度(T500)对于50 km和300 km空间尺度,离差增长差异有限,而对于850 hPa相对湿度(RH850)和500 hPa相对湿度(RH500),300 km空间尺度试验离差增长要弱于50 km试验;对于失相关时间尺度敏感试验,失相关时间尺度1 h方案其离差增长最差,如对于18 h的U850,失相关时间尺度1 h试验离差为0.7 m/s,失相关时间尺度9 h试验离差为1.6 m/s。

为了进一步说明SPPT扰动在不同高度的传播特征,图3给出了综合平均7组SPPT扰动试验的不同模式面U风场均方根离差随预报时效的演变。可以看出不同模式面离差随预报时效的增长特征不同,低层(20层模式面以下)离差增长较快,而高层增长较慢,整体呈现出由低层向高层的传播特征,主要是由于SPPT扰动的作用对象是积云对流、微物理过程、边界层过程等物理过程的综合倾向,而这些物理过程在模式低层作用明显;模式中高层主要是动力倾向起作用,受物理过程扰动的影响较为滞后,因此难以快速达到饱和,进而呈现出不断增长特征。

图2 不同SPPT扰动试验与控制试验的预报均方根离差随预报时效的演变:(a)500 hPa纬向风;(b)500 hPa温度;(c)500 hPa相对湿度;(d)850 hPa纬向风;(e)850 hPa温度;(f)850 hPa相对湿度

图3 SPPT扰动试验与控制试验各模式面纬向风均方根离差随预报时效的演变(单位:m/s)

综上所述,SPPT扰动整体呈现出由低层向高层的传播特征,且大/小的格点标准差利于/不利于离差增长,但是试验过程中也发现如果格点标准差过大会导致计算崩溃,因此通过增加格点标准差来提高集合离散度的前提是保证计算稳定;此外小的失相关空间尺度有利于低空要素尤其是水汽的离散度增长,试验中50 km的失相关空间尺度已接近模式最小可分辨尺度(15 km分辨率可解析的最小尺度为30 km),因此在本文集合系统中保留50 km失相关空间尺度的设置;敏感试验中大的失相关时间尺度有利于离散度增长,说明失相关时间尺度不宜过短(即随机强迫型不宜频繁变化),通过进一步试验表明,继续提高失相关时间尺度,其离差并未获得明显改进(图略),说明9 h的失相关时间尺度已接近饱和状态。因此通过敏感试验,对克拉玛依区域集合预报系统SPPT方案关键参数设置为:格点标准差0.55,失相关空间尺度50 km,失相关时间尺度9 h。

2.2 集合扰动特征

采用2.1确定的SPPT参数开展ETKF和ETKF-SPPT两套集合的对比试验。首先对集合预报扰动特征进行具体分析。对于一个Mi×Mj的二维格点场,引入平均绝对扰动:

(6)

其中:Ik(i，j)为第k个成员在格点(i,j)的初值,Ictl(i，j)为对应的控制预报初值,N为扰动的集合成员数(本研究中为21个)。

图4 ETKF(a—c)和ETKF-SPPT(d—f)两种方案在2021年2月24日00时的500 hPa的u分量风集合扰动场演变(单位:m/s):(a、d)00 h;(b、e)12 h;(c、f)24 h

图4给出了两种方案2021年2月24日起报的500 hPa的u分量风的平均绝对扰动(简称集合扰动)随预报时效的演变,图4a可以看出,经过ETKF方案在初始时刻会产生一定的扰动分布,且扰动形态具有一定的流形特征,12 h预报时效扰动范围有所增加,如模式中部巴尔喀什湖附近扰动量级基本达到2 m/s,但随着预报时效进一步延长,扰动增长达到饱和甚至缩小,如24 h基本上在模式区域内虽然也存在离散度大值中心,但是整体范围有明显缩小,尤其是模式范围内上游区域,离散度逐渐消散,仅在下游地区保持一定的离散度,这是由于单一侧边界在上游区域外侧施加作用,使得不同成员趋向一致,这种现象随着预报时效的延长表现得尤为显著。而对于ETKF-SPPT方案,其初始离散度与ETKF扰动分布一致,在加入了SPPT扰动之后,随着模式积分过程其整体离散度有明显增加,如24 h预报时效ETKF-SPPT比ETKF扰动范围和量级更大,说明引入模式扰动可以有效改善原系统的扰动受限问题。

图5 所有成员u扰动绝对值平均的垂直分布(不同线型代表分别表示不同预报时效):(a)ETKF方案;(b)ETKF-SPPT方案

从平均绝对扰动随时间演变特征(图5)来看,对于ETKF方案,初始时刻具有一定的扰动幅度,但扰动增长不够明显,其中0—3 h在各个层次均存在一定增长,但是3 h之后在高层增长已经饱和,在模式低层扰动量随着预报时效的延长还有所衰减,12 h扰动在1 m/s以下,甚至小于初始扰动的量值;ETKF-SPPT方案扰动能量在各个层次相对于ETKF方案有明显增大,12 h扰动最大值超过2 m/s。说明模式扰动在模式积分过程中不断施加作用,使扰动保持增长态势。

图6 集合平均均方根误差(RMSE)与离散度(SPREAD)随时间演变特征:(a)T850;(b)U850;(c)T2 m;(d)U10 m

2.3 集合预报检验结果

对批量试验时段内的结果进行评分统计,进一步定量分析引入SPPT模式扰动的预报效果。采用的集合预报检验方法为集合平均的均方根误差(Root Mean Square Error,RMSE)、集合离散度以及Talagrand分布。批量试验检验主要包括等压面要素以及地面要素。

集合平均的均方根误差(RMSE)和集合离散度(SPREAD)是集合预报最常用的检验方法。本文首先对低空和近地面温度、纬向风场进行检验。图6给出了试验时段内统计的850 hPa温度和纬向风(T850,U850)、2 m温度(T2m)以及10 m纬向风(U10m)的检验结果。从图中可知ETKF-SPPT方案的RMSE在36 h短时效内均小于ETKF方案,尤其是短预报时效最为明显,如12 h预报时效,ETKF集合的RMSE为1.3 K,而ETKF-SPPT集合的RMSE为1.55 K,说明物理过程扰动方法有效降低了集合预报短预报时效的误差;从集合离散度来看,ETKF方案的离散度在短时效内略优于方案,T850的36 h集合离散度为0.61 K,而ETKF-SPPT方案的36 h集合离散度为0.52 K,对于其他层次和要素的检验结果也较为类似,这里不再赘述。以上结果说明SPPT方法可以有效提高集合预报的准确性,减小预报误差,增加集合离散度。

Talagrand分布(Hamill,2001)是一种衡量集合预报的可靠性的方法,该方法通过将集合成员排列成若干区间,统计观测落在各区间的频率来获得柱状分布。Talagrand分布平缓表示观测落在各个成员区间的频率相当,集合预报系统较为可靠。“U”型分布说明集合成员不够发散,“L”型或“J”型分布表明集合成员存在预报负偏差或正偏差。图7给出了两套集合500、850 hPa纬向风和温度24 h预报时效的Talagrand分布。总体上ETKF的Talagrand呈U型分布,表现出一定的离散度不足;ETKF-SPPT方案Talagrand分布更平缓,离散程度有了明显改善,其他要素检验结果也类似(图略),说明ETKF-SPPT方案的概率分布可靠性更高。

图7 ETKF和ETKF-SPPT方案24 h预报Talagrand分布:(a)U850;(b)T850

2.4 大风个例预报效果

为了探索克拉玛依集合预报系统对当地强大风天气的预报能力,研究了试验时段内2月25日一次典型的强大风个例,本次个例中克拉玛依站瞬时风力可达21.3 m/s,给城区造成了一定破坏。图8分别给出了两种集合对克拉玛依站2021年2月24日00时起报的30 h预报(至25日06时)各集合成员以及集合平均预报的瞬时风速演变,同时给出了观测实况。图中的实况风速演变可以看出该大风天气过程在克拉玛依站具有明显的强风时段,在25日03时风速较小,仅为6 m/s左右,但在25日06时起风,达20 m/s以上。从集合成员预报可以看出,ETKF集合各成员预报不够发散,且对起风时段的变化描述不够准确,24 h、27 h预报时效显著高估了实况风速,24 h集合平均预报风速已接近10 m/s;ETKF-SPPT成员来看,首先各个成员预报较为发散,27 h预报有若干成员非常接近实况量级,尚未起风,集合平均风速在10 m/s以下,在30 h预报时效集合平均风速才达20 m/s左右,各个成员预报相对于ETKF集合要更加接近实况量级。

为了进一步分析ETKF和ETKF-SPPT两套集合方案在克拉玛依大风预报中的效果,研究了25日03时两套集合的预报情况。图9a给出了观测的2021年2月25日03时克拉玛依站极其周边站点的风速分布,可以看出该时刻克拉玛依站尚未起风,风速值在2～5 m/s,而上游的几个站点风速均达到5～10 m/s或10 m/s以上,说明大风尚未达到山脚下的克拉玛依站。图9b、c给出了ETKF和ETKF-SPPT两种集合方法2021年2月25日03时的时刻集合平均预报风速分布。ETKF集合(图9b)可以看出其集合平均预报10 m/s覆盖范围较大,说明集合整体预报相对于实况偏强,大风已影响到了克拉玛依站;对于ETKF-SPPT集合(图 9c),10 m/s风速区域尚未影响到克拉玛依城区,克拉玛依站风速在5～10 m/s,与实况对应较好,改善了ETKF集合起风过早的现象。

图8 2021年2月14日00时—15日06时的克拉玛依站的逐3 h的10 m风速演变(单位:m/s;红线为实况,蓝线为集合平均预报,灰线为集合成员预报):(a)ETKF方案;(b)ETKF-SPPT方案

图9 2021年2月25日03时风速观测及集合平均预报:(a)观测(点为站点风速,单位:m/s;阴影为地形高度,单位:m),(b)ETKF方案;(c)ETKF-SPPT方案预报(阴影为风速,单位:m/s;等值线为地形高度,单位:m)

该大风个例表明,在引入SPPT模式扰动方案后,集合成员对强天气的捕捉能力更强,有效改善了原有集合系统仅采用初值扰动导致的离散度不足,改进了克拉玛依集合预报系统对强大风天气的定时定量预报。

本文基于克拉玛依集合预报系统探索了通过SPPT模式扰动改进集合效果的方法,尝试通过敏感试验获得适用于该系统的SPPT扰动参数,并对比了仅采用ETKF初值扰动以及ETKF初值扰动加随机物理过程模式扰动(ETKF-SPPT)两组集合试验,得出以下结论:

1)随机物理过程倾向扰动参数敏感试验表明,大/小的格点标准差,利于/不利于离差增长;小的失相关空间尺度有利于低空要素尤其是水汽的离散度增长,而大的失相关时间尺度有利于离散度增长;SPPT扰动具有从低层向高层的传播特征。通过敏感试验,获得了适用于克拉玛依区域集合预报系统SPPT方案关键参数。

2)对比了ETKF和ETKF-SPPT两种集合方案,通过扰动分析表明仅采用ETKF初值扰动方案,受单一侧边界条件的影响,各成员之间的离散度发展受到约束,模式范围内上游地区扰动能量逐渐消散,下游地区的扰动能量难以保持;引入随机物理过程模式扰动后,可在模式积分过程中不断产生新的扰动能量,使得模式区域内集合扰动增长能力显著提高。

3)集合预报检验结果表明相对于ETKF方案,引入模式扰动可以大幅增加集合离散度,并减少预报误差;对于集合预报的概率技巧评分具有较明显提升。

4)大风个例试验结果表明ETKF-SPPT集合能够有效改善原ETKF集合的起风时间,各成员更为发散,对局地强大风的量级和时段具有更准确的预报能力,尤其是对原方案起风时间过早的现象有所改进。

以上结论表明仅采用初值扰动的克拉玛依集合预报系统效果有限,各个成员预报不够发散,难以有效表达局地强大风的预报不确定性,通过增加模式扰动可以有效改进集合预报离散度,若干成员能够获得更好的定时定量预报,因此具有较好的实际应用价值。考虑到一些业务部门资料获取及存储条件的限制,并未通过全球集合来为区域集合提供侧边界扰动,如本文中的克拉玛依区域集合预报系统,今后随着系统的升级引入侧边界扰动后会对集合预报性能产生进一步提升。

猜你喜欢格点时间尺度克拉玛依带有超二次位势无限格点上的基态行波解数学物理学报(2022年5期)2022-10-09“引水工程”对克拉玛依局地小气候的改善农业灾害研究(2022年1期)2022-05-07CaputoΔ型分数阶时间尺度Noether 定理1)力学学报(2021年7期)2021-11-09一种电离层TEC格点预测模型空间科学学报(2021年4期)2021-08-30交直流混合微电网多时间尺度协同控制能源工程(2021年1期)2021-04-13时间尺度上非迁移完整力学系统的Lagrange 方程与Nielsen 方程苏州科技大学学报(自然科学版)(2021年1期)2021-03-24格点计算器数理化解题研究·综合版(2020年3期)2020-09-10克拉玛依稠油MOA基质沥青改性及应用石油沥青(2020年1期)2020-05-25时间尺度上完整非保守力学系统的Noether定理苏州科技大学学报(自然科学版)(2020年1期)2020-04-13Love Remains in the Human World校园英语·上旬(2018年5期)2018-06-30