2. 北京航天飞行控制中心, 北京 100094;
3. 中国科学院大学, 北京 100049;
4. 航天飞行动力学技术重点实验室, 北京 100094
2. Beijing Aerospace Control Center, Beijing 100094, China;
3. University of Chinese Academy of Sciences, Beijing 100049, China;
4. Science and Technology on Aerospace Flight Dynamics Laboratory, Beijing 100094, China
磁暴期间,太阳风能量以增强的对流电场、皮特森耗散电流和高能粒子沉降等形式进入高层大气,使得极区热层大气加热膨胀,并通过大气环流向中低纬区域扩散,引起全球大气密度扰动[1-2].扰动的一个重要表现为大气密度增大,导致中低轨卫星拖曳力增加,卫星轨道降低,提前陨落,甚至失踪;另外,扰动也会出现强的大气密度波以及局部低大气密度“空洞”,造成卫星轨道不稳定,给精密定轨和轨道预报带来困难.
近半个世纪以来,在大量探测数据和高层大气理论知识的基础上发展了各种半经验的高层大气模型,如CIRA系列、Jacchia系列、DTM系列、MSIS系列等.2002年Picone在MSIS模式的基础上提出了NRLMSISE00模式,首次将反常氧和热氧原子计入总的大气密度中,是当前被广泛应用的最先进的大气模式[3].但是进一步的数据分析表明,在磁暴期间,该模式不能刻画热层大气密度的时空分布特征,并且实测密度值与模式预报值的差异达到100%[4].因此研究磁暴期间热层大气密度时空分布规律,对提高现有热层大气模式的空间天气扰动预报能力具有重要意义,对保障中低轨道航天器安全运行具有重要应用价值.
随着星载加速度计技术的发展,利用高精度三轴加速度仪测量卫星飞行方向的非保守力加速度,从而反演出高精度、高分辨率的实测大气密度,为研究热层大气密度的时空分布规律提供了新途径. 2000年发射的CHAMP卫星和2002年发射的GRACE卫星由于纬度覆盖面广,在轨运行时间长,所携带的加速度仪测量精度高且稳定,因此在暴时大气密度观测方面取得了巨大成功.Bruinsma等将2000-2002年期间CHAMP卫星数据反演的密度与DTM-94、DTM-2000以及MSIS-86模型预测的结果相比较,发现模型对密度平均低估10%~20%,方差达到16%~20%[5].Sutton等选择了2004年7月20-29日发生的三个连续磁暴,利用CHAMP卫星数据计算了暴时大气密度对焦耳加热指数的响应时间,低纬度大气密度的响应时间为3~4 h,而在中高纬度的响应时间小于2 h[6]. Sutton对2003年10月29-30日超级磁暴期间的CHAMP数据分析,认为白天北半球的扰动大于南半球,夜间密度扰动传播速度比白天更大[7].Liu以2003年10月29-30日,30-31日,11月20-22日的三个超级磁暴为研究对象,分析了400km高度上大气密度扰动的半球不对称性以及由高纬到低纬的传播速度,认为白天夏季半球的密度扰动大于冬季半球,这与Sutton的结论基本一致;在磁暴的传播速度方面,Liu得出的结论与Sutton不同,他认为无论夜间还是日间,夏季半球的传播速度都大于冬季半球,冬季半球的日间传播速度大于夜间[8-9].
上述研究大气密度时空规律采用的方法一般是分别绘制白天和夜间的基于时间-纬度二维分布的密度图,用不同的颜色表示密度变化.这种方法的优点是简单直观,适合研究大尺度的密度变化,但是由于缺乏量化分析,得出的结论带有作者的主观性,比如Sutton和Liu使用了相同的数据但在磁暴传播速度方面得出了不同结论.另外这种方法不能把时间和纬度两个影响因素区分开.
经验正交分析法(Empirical Orthogonal Function analysis,简称EOF)应用于大气密度的研究是一种新的思路,可以把时间因素和纬度因素分离出来.周旭等对2003-2007年CHAMP卫星的热层大气密度数据进行了EOF分析,提取出了太阳活动周变化与年变化等分量,还得到了赤道密度异常(EMA)结构[10].Matsuo和Forbes对2001-2008年的CHAMP卫星密度数据进行了EOF分解,阐述了前4项正交函数对应的物理意义[11].Lei对2002-2010年地磁平静期的CHAMP卫星和GRACE卫星的观测数据进行联合EOF分析,揭示了热层大气密度的年变化和季节变化,认为热层密度的季节变化与纬度有很强的相关性,春分和秋分时刻往往出现最大值,在中高纬的冬季出现最小值;南半球的年变化比北半球更为明显[12].这些研究不但验证了热层大气密度已知的一些规律,如赤道密度异常,还揭示了一些新的细节,这说明EOF分析应用于大气密度研究有其独到的优势.
关于EOF在大气密度方面的研究主要侧重于长期的、平静期的密度变化[10-12],针对磁暴期间的研究鲜有报道.本文对2002-2006年多个磁暴期间的400km高度上的CHAMP卫星密度数据进行了经验正交分解,分析了暴时热层大气密度的全球分布特征,以及第一阶EOF级数与ap和Dst指数的关系,为进一步建立暴时热层大气密度扰动模式奠定基础.
2 CHAMP卫星数据处理与经验正交分析方法 2.1 加速度仪数据处理CHAMP卫星是德国2000年发射的一颗科学小卫星,主要进行地球重力场、磁场、电离层和大气层的观测.在轨运行10年,其轨道为近圆轨道,倾角87.3°,轨道高度从初始454km下降到350km左右,轨道周期约90 min,每3个多月可覆盖所有地方时[13].
CHAMP卫星上搭载的加速度仪具有精度高、分辨率高的特点.在卫星飞行方向上,精度达到10-9m·s-2,换算成大气密度相当于10-15~10-14 kg·m-3 [14].在CHAMP卫星高度上,大气密度约为10-13~10-11 kg·m-3,因此对于热层大气密度研究来说,加速度仪数据具有足够高的精度.德国地球科学研究中心负责管理和发布CHAMP的科学数据,本文使用的加速度仪数据属于Level-2级,每10s得到一个大气密度测量值,在日地空间环境瞬时剧变的磁暴期间,这样的高分辨率对于揭示热层大气的局部细节有很大帮助.
加速度仪测得的是卫星飞行方向的非保守力加速度总和,在300~500km高度上,主要的非保守力有大气阻力、太阳辐射压力和地球反照,如(1)式:
(1) |
地球反照加速度aalb量级约为10-10m·s-2,相对于加速度仪10-9m·s-2的测量精度,可忽略不计.太阳光压加速度asolar量级约10-8m·s-2,利用Box-Wing光压模型扣除[15],如(2)式:
(2) |
其中,e是地影因子,ps是距离太阳1个天文单位处的光压强度,N是光照面个数,Ai是每个光照面的面积,m是卫星质量,θi是每个光照面法向与光线的夹角,ni是每个光照面的单位法向量,s是卫星到太阳的单位矢量,Crdi是每个光照面的散射系数,Crsi是每个光照面的镜面反射系数.
大气密度ρ的反演由(3)式给出:
(3) |
其中,CD为拖曳系数,一般取2.2,s是卫星迎风面积,v是沿飞行方向卫星与大气的相对速度. CHAMP卫星有13个外表面,卫星姿态不同,迎风面积也不同,对于任意时刻的迎风面积,由(4)式计算:
(4) |
βi是卫星相对大气的运动方向与各表面的夹角,si是每个表面的面积.
为了消除卫星轨道高度变化引起的大气密度变化,利用NRLMSISE00模式[3]将观测资料标准化到400km高度,如(5)式.ρ(400)为修正后400km高度上的实测密度,z为卫星所在高度,ρ(z)为真实高度上的实测密度,ρM(400)、ρM(z)分别为NRLMSISE00模式计算的400km高度和真实高度上的大气密度.
(5) |
在大气密度反演中存在一定误差,其误差源主要有以下三个方面[5, 7].第一,拖曳系数CD的不准确会引起密度5%~10%的误差;第二,加速度仪标定参数会引起2%的系统误差;第三,上述计算过程中未考虑中性风的影响,地磁平静期,中性风在极区和赤道区引起的密度误差分别小于7.5%和1%,磁暴期间中性风速度增大,在极区引起的误差可达到20%,但由于暴时大气密度增大200%以上,所以中性风引起的误差在本文中可以忽略.
2.2 经验正交分析法经验正交分析法的核心思想是对原始数据引入一种能使其简化,并剔除冗余信息的线性变换.这类方法被广泛应用于大气科学及地球科学各领域,是一个重要的数值分析工具.在数学上,它具有较为完整的理论框架;在物理上,它能简洁而巧妙地揭示出不同学科领域内各自的物理内涵;在计算方法和描述方式上,它又能适应于其它数学的、动力的、统计学的技巧.
将标准化到400km的实测大气密度平滑处理成时间和纬度的网格点,记为x(t,φ),时间间隔为1h,纬度间隔为5°.由于白天和夜间的热层大气具有明显不同的特性,故把地方时为0800-1700和2000-0600的数据分开来处理.经验正交分析就是把密度分解成时间函数Z和空间函数V,如(6)式,也可写成矩阵形式X=VZ.
(6) |
其中,i=1,2,…,m,j=1,2,…,n;t和φ分别表示每个网格点对应的时间和纬度;vk为第k阶EOF基函数,表征大气密度的纬度变化;zk为第k阶EOF系数,表征大气密度的时间变化;p为EOF级数的项数.V和Z都是正交的,本文采用文献[16]的方法计算EOF系数和基函数.
令A=XXT=VZZTVT,A为实对称矩阵,根据实对称矩阵分解定理,A=VΛVT,Λ是A的特征值组成的对角矩阵.把A的特征值从大到小排列,λ1≥λ2≥λ3≥…≥λm,对应的特征向量v1,v2,…,vm组成V的列向量,V求出后,由Z=VTX可以求出Z.
经验正交函数具有快速收敛的特点,只要取前p阶(
按照Gonzalez提出的标准,磁暴强度按照Dst大小分类,Dst<-100的磁暴为强磁暴[17].按此标准,我们选取了2002-2006年期间的36个强磁暴为研究对象,表 1列出了这些磁暴的最大ap值、最小Dst值和发生日期.有些磁暴事件,虽然Dst没有达到-100,但是由于地磁扰动较大,表现为ap指数较大,所以也将其作为研究对象.两种指数对磁暴的刻画是有差异的,比如以ap指数来看,13号磁暴的强度比14号大,而从Dst指数来看,结果刚好相反.
Dst指数在暴时往往呈现先增大后减小再慢慢恢复的趋势,而ap指数的变化较简单,因此磁暴的开始和结束时间的截取主要依据ap指数.当ap达到暴前日均值ap的2倍或ap>48时,认为磁暴开始;当ap恢复到暴前日均值附近或ap<27时,认为磁暴结束.磁暴选取时间过短,导致数据不完整,可能会漏掉某些暴时密度特征;选取时间过长,会使结果中包含平静期数据,弱化暴时的密度特征.
3 分析结果 3.1 暴时实测大气密度EOF分解结果按照2.2节的方法,把400km暴时实测大气密度按照地方时分为白天和夜间两组,再把纬度和时间的网格点,分解为时间函数Z和空间函数V.以1号磁暴白天大气密度为例,图 1a、1c和1e给出了特征值λ按降序排列所对应的前三阶时间系数,图 1b、1d和1f为相应的EOF基函数.可以看出时间系数Z1到Z3迅速降低,其对应的EOF基函数从大尺度过渡到中小尺度.磁暴发生时,全球几乎各个纬度上的密度都增大,因此密度绝对值的增大主要反映在第一阶时间系数上,第一阶EOF基函数反映了大尺度上密度随纬度的变化.另外,磁暴期间还会出现强的大气密度波以及短时的、局部的低大气密度空洞,这主要反映在第二、三阶以及更低阶的特征向量场上.
使用与图 1相同观测资料,图 2a给出了CHAMP卫星实测密度的分布图,图 2b-2d给出了由不同阶数的EOF特征向量场绘制的大气密度随地理纬度和UT的分布图.第一阶EOF特征向量场描述的是大气密度在大尺度上的粗略变化,随着阶数的增大,短时的、局部的扰动被揭示出来,阶数越高的向量场与实测密度的分布图越相似.高阶向量场描述的多是小幅度的、短暂的、小尺度的密度变化,往往具有偶发性,甚至是白噪声的特征,而低阶向量场描述的是大幅度的、持续时间较长的、大尺度的密度变化,有比较明确的物理意义.本文对所有磁暴第一阶特征向量场的贡献率Q进行了统计,发现所有的Q值都大于95%,平均值在98%左右,因此本文后续研究采用第一阶特征向量场
图 3和图 4所示分别是白天和夜间第一阶EOF函数随纬度分布的曲线.可以看出,大气密度的纬度分布与季节密切相关,6-8月的白天,处于夏天的北半球的大气密度大于处于冬季的南半球,如图 3a,11-1月的白天,处于夏天的南半球的大气密度大于处于冬天的北半球,如图 3c,这与Liu和Sutton对2003年10-11月三个超级磁暴的分析结论一致[7-9].除此之外,本文对更多的磁暴事件分析发现,夜间也有类似特征,如图 3b和图 3d.总之,无论白天还是夜间,夏季半球的密度大于冬季半球.暴前平静期也呈现类似的不对称结构,如图 5,这说明太阳辐射对南北半球的不同是产生季节不对称性的主要原因.另外,夏季半球的焦耳加热率是冬季半球的两倍[18],当太阳风携带高能粒子进入极区时,不同的焦耳加热率导致大气升温幅度不同,这也是造成季节不对称结构的一个可能机制.
春秋季节南北半球的大气密度几乎对称分布,如图 4a和图 4b.大部分磁暴,尤其是春秋季节,白天EOF函数曲线呈现出赤道密度异常(EMA)结构,即低纬地区密度最大值在南北纬20°~30°附近,赤道附近略低,高纬度地区至两极逐渐降低;夜间EOF曲线呈现倒抛物线的变化趋势,即赤道和低纬度地区密度值较小,高纬度地区至两极密度缓慢增大.无论白天还是夜间,大部分磁暴在中低纬地区的密度变化幅度较小,而在高纬度和极区附近的密度变化各不相同.比如2002年4月和10月发生的磁暴,白天密度分布从高纬至极区呈平缓下降趋势,南北半球对称分布,如图 4a-4d所示;而2005年6月和8月发生的磁暴,白天密度在北半球随纬度增高而增大,呈现出非常明显的南北半球不对称结构,如图 3a-3b.这些现象反映了暴时极区复杂的物理环境对密度的影响.
从磁暴发生的时间来看,对于发生时间相近的磁暴事件(间隔20天以内),EOF函数的变化曲线非常相似,时间间隔越短,相似度越高,如图 4c-4d所示,EOF曲线几乎重合.图 5以2002年5月23日、2003年6月18日和2004年2月11日的三个磁暴为例,给出了暴时和暴前三天的EOF曲线,可以看出暴前平静期三天的EOF曲线几乎重合,暴时的EOF曲线比平静期略有变化,但形状基本相同.以图 5a为例,磁暴发生之前,密度最大值在15°S和30°N附近,赤道附近密度值相对较小,密度最小值在两极区,南极区最小;磁暴发生期间,密度绝对值增大,但随纬度分布的相对幅度和暴前变化不大,密度的最大值和最小值与暴前所在的纬度区域一样.
从上述相邻磁暴以及磁暴前后EOF曲线的分布可以得出结论,大气密度的纬度变化在一段时期内非常稳定,尤其是平静期,即使磁暴对这种稳定性的破坏也不大,暴时密度的增大和扰动主要体现在时间函数上,而纬度分布上的变化很小.这是因为大气密度的变化与大气温度变化密切相关,暴时皮特森耗散电流在极区对大气产生焦耳加热,富含较重分子成分的大气向上升并全球对流,造成大气密度在空间分布上的变化,从大尺度上来看,全球对流和大气温度的变化都是缓慢的,因此在较长时间范围内,大气密度的纬度分布特征比较稳定.
3.3 暴时大气密度与ap指数和Dst指数的关系按照2.2节的方法分解得到的EOF系数矩阵Z第一列是第一阶EOF函数的权系数,记为Z1,表征了大尺度上大气密度随时间变化的趋势,时间间隔为1h.对时间间隔为3h的ap指数进行插值,处理成每小时一个采样点.ap指数是对磁场扰动等级的划分,不是对物理量的准确测量,而且也不连续,插值后的ap指数不会改变原来的磁场扰动等级,因此插值误差不会影响后续分析.Dst指数时间间隔为1h,故无需处理.利用互相关分析法,计算白天和夜间Z1与ap指数和Dst指数的互相关函数,得到Z1相对ap指数和Dst指数的响应时间,以及最大相关系数.
3.3.1 暴时大气密度相对ap指数和Dst指数的响应时间大气密度与ap指数和Dst指数具有最大相关系数时对应的延迟时间如图 6所示,不同磁暴的延迟时间差异较大.大气密度滞后ap指数的时间主要集中在2~6h的范围内,其中,白天平均滞后3.2 h,夜间平均滞后4.3h.大约50%的磁暴事件中,大气密度滞后Dst指数0~1h,在其余事件中,大气密度响应比Dst指数反而提前0~2h,其中,白天大气密度响应平均提前1h,夜间响应几乎同步.
上述统计说明,ap指数比Dst指数能够更快速地对磁暴事件做出响应,这是因为ap指数由全球13个地磁台站测量得到的地磁扰动强度的指数,是行星际磁情指数,而Dst是磁暴时赤道环电流的强度,反映了低纬度和赤道地区地磁扰动[19].由于磁暴总是在高纬度和极区首先发生,需要一定时间才能形成赤道环电流,故大气密度扰动滞后ap指数的时间比Dst指数要长.
3.3.2 暴时大气密度与ap指数和Dst指数的相关性表征时间变化的权重系数Z1与ap指数、Dst指数具有较好的相关性,见表 2.定义相关系数0.8以上为强相关,白天大气密度与ap指数强相关的磁暴事件占81%,与Dst指数强相关的磁暴事件占75%;夜间大气密度与ap指数强相关的磁暴事件占92%,与Dst指数强相关的磁暴事件占78%.
总的来看,夜间大气密度与两类指数的相关性强于白天,这是因为白天大气密度除了受地磁影响外,还受到太阳辐射的影响.大气密度与ap指数的相关性比与Dst指数的相关性强,这与Dst指数的测量手段和物理意义相关,Dst指数是由5个中低纬地磁台H变化均值归算到赤道得到的,主要描述低纬度地区的地磁变化.周云良等对2003年11月20日的磁暴进行了分析,认为环电流过程对中低纬度区大气密度扰动影响较大[20].本文研究的大气密度范围覆盖±85°,因此与反映全球地磁扰动的ap指数相关性更好.
4 结论本文利用CHAMP卫星加速度仪数据对36个强磁暴期间的大气密度进行EOF分解,对第一阶EOF基函数和系数统计分析,结果表明:
(1)大部分磁暴的白天大气密度分布呈现出赤道密度异常结构,夜间则呈现倒抛物线的形状.
(2)大气密度的半球不对称性与季节相关,无论白天还是夜间,夏季半球的密度大于冬季半球,在春秋季节南北半球的大气密度几乎对称分布.
(3)大气密度的纬度分布特征在若干天内具有良好的稳定性,即使磁暴对这种稳定性的破坏也不大.
(4)大气密度的时间函数与ap指数、Dst指数具有较好的相关性,其中夜间的相关性大于白天,与ap指数的相关性大于Dst指数.
(5)大气密度对磁暴的响应速度在日照区比在阴影区快,ap指数对磁暴的响应速度比Dst指数快.
本文将EOF方法应用于磁暴期间热层密度的分析在国内尚属首次,并发现了暴时大气密度纬度分布的稳定性特征.以上这些结论对于提高现有热层大气模式在暴时的预报能力具有重要参考价值,如何把定性的规律转化成定量的数学模型,建立暴时热层大气密度模式,是下一步的研究工作.
[1] | 涂传诒. 日地空间物理学. 北京: 科学出版社, 1988 . Tu C Y. Solar-Terrestrial Space Physics (in Chinese). Beijing: Solar-Terrestrial Space Physics: Science Press, 1988 . |
[2] | Tsurutani B T, Gonzalez W D, Kamide Y, et al. How does the thermosphere and ionosphere react to a geomagnetic storm?. Geophysical Monograph Series , 1997, 98: 203-225. |
[3] | Picone J M, Hedin A E, Drob D P, et al. NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues. J. Geophys. Res. , 2002, 107(A12): 1468-1483. DOI:10.1029/2002JA009430 |
[4] | Rhoden E A, Forbes J M, Marcos F A. The influence of geomagnetic and solar variabilities on lower thermosphere density. Journal of Atmospheric and Solar-Terrestrial Physics , 2000, 62(11): 999-1013. DOI:10.1016/S1364-6826(00)00066-3 |
[5] | Bruinsma S, Tamagnan D, Biancale R. Atmospheric densities derived from CHAMP/STAR accelerometer observations. Planet. Space Sci. , 2004, 52(4): 297-312. DOI:10.1016/j.pss.2003.11.004 |
[6] | Sutton E K, Forbes J M, Knipp D J. Rapid response of the thermosphere to variations in Joule heating. J. Geophys. Res. , 2009, 114(A4): A04319. DOI:10.1029/2008JA013667 |
[7] | Sutton E K, Forbes J M, Nerem R S. Global thermospheric neutral density and wind response to the severe 2003 geomagnetic storms from CHAMP accelerometer data. J. Geophys. Res. , 2005, 110(A9): A09S40. DOI:10.1029/2004JA010985 |
[8] | Liu H, Lühr H, Henize V, et al. Global distribution of the thermospheric total mass density derived from CHAMP. J. Geophys. Res. , 2005, 110: A04301. |
[9] | Liu H, Lühr H. Strong disturbance of the upper thermospheric density due to magnetic storms: CHAMP observations. J. Geophys. Res. , 2005, 110(A9): A09S29. DOI:10.11029/2004JA010908 |
[10] | 周旭, 万卫星, 赵必强, 等. 基于CHAMP卫星观测数据对热层大气密度的经验正交分析. 空间科学学报 , 2010, 30(3): 228–234. Zhou X, Wan W X, Zhao B Q, et al. Empirical orthogonal function (EOF) analysis on the thermospheric total mass density retrieved from CHAMP observation. Chin. J. Space Sci. (in Chinese) , 2010, 30(3): 228-234. |
[11] | Matsuo T, Forbes J M. Principal modes of thermospheric density variability: Empirical orthogonal function analysis of CHAMP 2001-2008 data. J. Geophys. Res. , 2009, 115(A7): A07309. DOI:10.1029/2009JA015109 |
[12] | Lei J H, Matsuo T, Dou X K, et al. Annual and semiannual variations of thermospheric density: EOF analysis of CHAMP and GRACE data. J. Geophys. Res. , 2012, 117(A1): A01310. DOI:10.1029/2011JA017324 |
[13] | Reigber C, Lühr H, Schwintzer P. CHAMP mission status. Advances in Space Research , 2002, 30(2): 129-134. DOI:10.1016/S0273-1177(02)00276-4 |
[14] | Sutton E K, Nerem R S, Forbes J F. Density and winds in the thermosphere deduced from accelerometer data. Journal of Spacecraft and Rockets , 2006, 44(6): 1210-1219. DOI:10.2514/1.28641 |
[15] | Milani A, Nobili A M, Farinella P. Non-gravitational Perturbations and Satellite Geodesy. Bristol: Taylor and Francis, 1987 . |
[16] | Storch H V, Zwiers F W. Statistical Analysis in Climate Research. Cambridge: Cambridge University Press, 2002 . |
[17] | Gonzalez W D, Joselyn J A, Kamide Y, et al. What is a geomagnetic storm?. J.Geophys. Res. , 1994, 99(A7): 5771-5792. |
[18] | Forbes J M, Gonzalez R, Marcos F A, et al. Magnetic storm response of lower thermosphere density. J. Geophys. Res. , 1996, 101(A2): 2313-2319. DOI:10.1029/95JA02721 |
[19] | 徐文耀. 地磁活动指数的过去、现在和未来. 地球物理学进展 , 2009, 24(3): 830–841. Xu W Y. Yesterday, today and tomorrow of geomagnetic indices. Progress in Geophys (in Chinese) , 2009, 24(3): 830-841. |
[20] | 周云良, 马淑英, LührH, 等. 2003年11月超强磁暴热层大气密度扰动及其与焦耳加热和环电流指数的关系--CHAMP卫星观测. 地球物理学报 , 2007, 50(4): 986–994. Zhou Y L, Ma S Y, Lühr H, et al. Changes of thermospheric mass density and their relations with Joule heating and ring current index during Nov. 2003 superstorm-CHAMP observations. Chinese J. Geophys (in Chinese) , 2007, 50(4): 986-994. |