2010, Vol. 53
2. 北京大学湍流与复杂系统研究国家重点实验室, 北京 100871;
3. 中国科学院大气物理研究所LASG, 北京 100029
2. State Key Laboratory for Turbulence and Complex Systems, Peking University, Beijing 100871, China;
3. LASG, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China
气候系统是一个复杂的非线性耗散系统.气候变量的时间序列是不同时间尺度的内部震荡及外部强迫因子叠加而成的结果，所含的时间尺度跨越几个量级.因此气候过程可视为多元随机过程，具有非周期、非线性的特点.但同时，气候变量的时间序列也存在相似性和确定性的特征.应用混沌和分形的方法定量地分析气候变量的非线性演化过程，揭示了大气的变化呈现不同时间尺度间自相似的分形结构[1~3].简而言之，大尺度上是高低值的震荡，而每个峰值和谷值区间，又可以分为更小尺度的高低值震荡.
分形理论解释了复杂气候系统中有序与无序的统一，以及确定性与随机性的统一.除相似性之外，分形结构的另一个重要特征是气候系统具有无特征尺度性，同时对应的分形结构具有标度不变性，存在自相似的标度律.标度理论的发展给分形集一个定量的描述途径[4, 5].
通过分形维数，可以定量地描述气候系统变量的自相似结构特征.但研究发现，应用单一的分形维，如容量维或者信息维，可能无法全面描述非均匀的多分形结构[5].多重分形（Multi-fractals），又称为多标度分形，或称为多重分形测度（Multi-fractal Measure），已经在大气科学领域广泛发现并研究[6~9]，用以刻画气候系统的非均匀和各向异性特征.
分析多分维结构，主要有两种方法.其一是采用广义分维（Generalized dimension）[6~8]，另一种是采用多分形谱（奇异谱Spectrum of Singularities）[9]的方法来描述多分形现象.传统的估计多分形谱的方法难以克服非平稳的影响，在本文中，应用去趋势波动分析（Detrended Fluctuation Analysis，简称DFA）的一种扩展方法，多分形去趋势波动分析方法（Multi-fractal Detrended Fluctuation Analysis，简称MF-DFA）分析北京市不同气候变量的时间序列，分析其多分形特征，并得出各自的多分形谱.
DFA方法最早由Peng等[10]在研究DNA序列的分子链结构的幂函数相关特性时提出，可以有效地去除各阶趋势项，消除其中的虚假相关，得到非平稳时间序列中真实的长程相关性.目前DFA方法已经广泛地应用于气象及地球物理领域中[11~16]. Kantelhardt等[17]在2002年基于DFA方法，与标准配分函数基础上的多分形公式体系相联系，发展了多分形去趋势波动分析方法，用以探测非平稳时间序列中的多分形特征.MF-DFA方法在很多领域都有着广泛的应用[18~21].
2 研究方法在本文中，去掉原始序列的季节循环之后，考虑一个给定的等间距逐日时间序列（i=1，2，…，N），MF-DFA方法的计算步骤[17]如下：
第1步，计算距平的累加值，得到一个新的廓线序列：
|
(1) |
x表示原始序列的平均值.廓线序列保留原始序列变量特性的同时，降低了噪声的水平.
第2步，在应用波动分析前，将廓线序列Y（i）等分为大小为s的两组互不重叠的等时间段Ns=N/s，一组顺序操作，一组逆序操作，以保证在N不能被s整除的情况下所有的数据信息都能够得以利用.
第3步，在每个时间段v内，应用不同阶的多项式拟合pvk来去掉局部的趋势，多项式拟合的阶数即对应DFA方法的阶数（k=i，DFA（i））.然后得到的是去掉趋势的时间序列：
|
(2) |
第4步，计算每个时间区间内去趋势序列的二次方波动，即Ys（i）序列的方差：
|
(3) |
第5步，不同时间段的均方根偏差给出了q阶DFA波动函数：
|
(4) |
第6步，对不同的时间段长度s重复上述计算.
第7步，在s和Fq（s）的双对数坐标图上，如果原序列xi是长程幂律相关的，则Fq（s）与s成幂律关系，即Fq（s）-s（hq），h（q）描述了q阶波动函数的标度行为.
h（q）也称为标度指数，反映了不同时间尺度上的自相似分形行为.当q > 0的情况下，如公式（4）所示，h（q）突出的是大振幅脉动的标度行为，而在q < 0的情况下，反映的则是小振幅脉动的标度行为.对于平稳时间序列而言，h（2）就是广泛应用的Hurst指数，因此我们称h（q）为广义的Hurst指数（generalized Hurstexponent）.对于单分形的时间序列，所有时间尺度上都用一个标度指数来表征，h（q）独立于q的变化；而对于多分形的时间序列，h（q）则随着q的变化而变化.一般情况下，表征大振幅脉动的h（q）值要小于表征小振幅脉动的h（q）值.
计算多分形奇异谱f（α）是描述多分形时间序列的一种传统的定量方法[22].它是通过各自的奇异性标度指数α（singular exponent，也称为Holder指数[22, 23]）来描述的.Holder指数的大小反映时间序列中某一点上的奇异性的强度：Holder指数越大，表示序列越平滑和规则，奇异性的强度越小；反过来，Holder指数越小，奇异性的强度则越大.而所有奇异指数的变化程度就体现了多分形时间序列波动的复杂程度.α是q的函数，也可以写为α（q），再根据质量指数谱τ（q）与h（q）之间的关系，多分形奇异谱f（α）与广义Hurst指数h（q）之间可以通过Legendre变换联系在一起[17, 24]：
|
(5) |
|
(6) |
|
(7) |
因此通过计算广义Hurst指数，我们就可以估计多分形时间序列的奇异谱，后者反映了时间序列中奇异性（如极值事件，及突变等奇异行为）的含量.
3 数据说明本文所用的数据来自于中国气象局信息中心发布的中国地面气候资料日值数据集.选用了北京站（站号54511）9种常规数据，包括平均气温（TMEAN）、最高气温（TMAX）、最低气温（TMIN）、平均气压（PMEAN）、最高气压（PMAX）、最低气压（PMIN）、相对湿度（R-H）、平均风速（WMEAN）和日照时数（SUN）；以及由最高气温和最低气温计算得到的气温日较差（DTR）数据.选取的资料长度为1956年至2005年，共50年数据.
4 结果及讨论图 1给出了日平均气温和气温日较差去掉年循环后截取的同时间段距平序列，如图所示，去掉年循环后，平均气温的距平序列仍然包含多尺度的振荡特征，存在不同尺度间的层次结构，这种非规则的动力特性反映了不同尺度的振荡间可能存在不同的标度行为.而气温日较差的序列则基本只包含高频的振荡，表现为序列值的剧烈变化及快速涨落.
|
图 1 平均气温和气温日较差去掉年循环后截取的同时间段距平序列（a）平均气温；（b）气温日较差. Fig. 1 Time series of anomaly daily mean temperature record (a) and diurnal temperature range (b) for the station Beijing |
图 2显示的是平均气温和气温日较差分别取不同阶数q（q=-6，-5，-4，-3，-2，-1，1，2，3，4，5，6）值时的lg（Fq（s））和lgs函数关系图，标度区域取10~1000天，（a，b）为平均气温（TMEAN）的情况；（c，d）为气温日较差（DTR）的情况.当阶数q固定时，在双对数图中分析波动函数Fq（s）与窗口s的关系，可以看到lg（Fq（s））和lgs之间存在显著的线性拟合关系，即Fq（s）与s存在幂律关系Fq（s）-sh（q），h（q）描述了q阶波动函数的标度行为，说明对于不同时间尺度的气候振荡，平均气温序列和气温日较差序列都是长程幂律相关的.
|
图 2 不同阶数q（q=-6，-5，-4，-3，-2，-1，1，2，3，4，5，6）值时的lg（Fq（s））和lgs函数关系图 由下至上对应阶数q=-6到q=+6.（a）平均气温；（c）气温日较差.不同阶数的序列在垂直方向作一定的移动以保证清楚可视.（b）和（d）分别是平均气温和气温日较差，对应-6阶到+6阶最小二乘法拟合的结果. Fig. 2 Log-Log plots of the MF-DFA curves of daily mean temperature (a) and diurnal temperature range (c) for the station Beijing From the top to the bottom curves correspond to different q (from q=-6 to q=+ 6) and are shifted vertically for clarity. Corresponding least square fit results for TMEAN and DTR are shown in (b) and (d). |
但平均气温序列与气温日较差序列的标度行为存在显著的差异，图 2a中平均气温序列对不同阶数q，由负的q值到正的q值，广义Hurst指数h（q）的值越来越小，显示大振幅振荡的标度指数值要小于小振幅振荡的标度指数值.而图 2c中气温日较差序列的h（q）随着阶数q值的变化，则基本不发生变化，近似为一个常数.从对应的线性拟合图 2b中可以清楚地看到斜率，也即指数h（q）的值随阶数q值的变化，这种广义Hurst指数对阶数值的依赖性，清楚地表明平均气温的序列具有多重分形结构.而气温日较差序列对应不同阶数，拟合的斜率基本不变，表现出单分形的特征.
图 3中给出了北京地区不同气候变量序列取不同阶数q（q=-6，-5，-4，-3，-2，-1，1，2，3，4，5，6）值时最小二乘法拟合Fq（s）与s函数关系的广义Hurst指数h（q）值.结果表明平均气温（TMEAN）、最高气温（TMAX）、最低气温（TMIN）、平均气压（PMEAN）、最高气压（PMAX）、最低气压（PMIN）、相对湿度（R-H）和平均风速（WMEAN）呈现多分形的特征.而气温日较差（DTR）和日照时数（SUN）序列则呈现单分形的特征.这说明它们之间所受制约的物理机制可能并不同.日照时数的单分形特征，说明温压湿风等常规气象场的奇异性，并不是由于辐射的直接强迫作用，而是可能来源于背景环流场的变化，锋面气旋系统带来的非对称及突变性是可能的原因[2, 25].而气温日较差的单分形特征，也说明了环流场中这种奇异性在最高气温和最低气温中具有位相的一致性.
|
图 3
q取-6到+6时，不同气候变量的广义Hurst指数h（q）与阶数q的对应关系
黑色空心方块对应MF-DFA拟合的结果，黑色实心圆对应二项式串级模式估计的结果. （a）平均气温；（b）气温日较差；（c）平均气压；（d）相对湿度；（e）平均风速；（f）日照时数. Fig. 3 h{q) versus g plots for different climate variables Open black squares: estimated from MF-DFA results like in Fig. 2; solid black circles: obtained by fits of the two-parameter binomial model. (a) Daily mean temperature (TMEAN); (b) Diurnal temperature range (DTR); (c) Daily mean pressure (PMEAN); (d) Relative Humidity (R-II); (e) Daily mean wind speed (WMEAN); () Sunshine duration (SUN). |
通常广义Hurst指数可以用如下公式来估计：
|
(8) |
这里引入一个扩展的二项式串级模式[22]，该方法是描述多分形标度行为的规范方法之一，细节内容见文献[17].最后，f=0时多分形奇异谱f（α）的谱宽可以由下式给出：
|
(9) |
由公式（8）结合前面公式（6）和公式（7），得αmax=
图 3中也给出了取不同阶数q（q=-6，-5，-4，-3，-2，-1，1，2，3，4，5，6）值时二项式串级模式估计得出的广义Hurst指数h（q）值，对比拟合的h（q）可以看出这个扩展的二项式串级模式可以很好地估计h（q）值与q值之间的关系.由公式（6）（7）（8）得出拟合参数a和b的值，就可以代入公式（9）中用二项式串级模式来估计时间序列多分形谱的谱宽.图 4中给出了北京地区各种气候变量用扩展二项式串级模式拟合的h（q）值与q值关系图.呈现多分形特征的平均气温（TMEAN）、最高气温（TMAX）、最低气温（TMIN）、平均气压（PMEAN）、相对湿度（R-H）和平均风速（WMEAN）均对应较宽的多分形谱，并且多分形谱宽的值在量级上比较接近.近似相近的值也说明，温压湿风场中表现出的奇异性特征，源于同一背景环流场的作用.呈现单分形特征的气温日较差（DTR）和日照时数（SUN）则对应较窄的多分形谱.参数a、b及多分形谱宽Δα的具体值在表 1中给出.
|
图 4 广义Hurst指数与阶数q的对应关系 （a）平均气温、最高气温、最低气温；（b）气温日较差；（c）平均气压、相对湿度及平均风速；（d）日照时数. Fig. 4 Generalized Hurst exponent h (q), as a function of for different climate variables (a) Daily mean temperature (TMEAN), daily maximum temperature (TMAX), and daily minimum temperature (TMIN); (b) Diurnai temperature range (DTR); (c) Daily mean pressure (PMEAN), Relative Humidity (R-II), and daily mean wind speed (WMEAN); (d) Sunshine duration (SUN). |
|
表 1 不同气候变量基于二次多项式串级模式拟合的参数a、b及多分形谱宽Δα的具体值 Table 1 Specific value of a, 6 and the width of the multi-fractal spectrum Aa, by fits of the two-parameter binomial model |
一般来说，时间序列中观测到的多分形标度行为主要是由两个原因造成的.其一是由于较宽的概率密度函数分布造成的，其二是由于大振幅脉动与小振幅脉动间不同的长程幂律相关特性造成的.通过分析随机洗牌的代用数据，我们可以很好地区分以上两种原因造成的多分形标度行为，这是由于随机洗牌后，原始序列的长程相关性质会被破坏掉，而原有的概率密度函数分布则得以保留.如果多分形标度行为是由长程相关性带来的，则随机洗牌后，所有拟合的h（q）值都将变成h（q）=0.5.如果两种原因皆存在，随机洗牌后的拟合的h（q）值比原值偏小.如图 5所示，平均气温和气温日较差序列经随机洗牌后，拟合的h（q）值都变成0.5，对应白噪声（其他变量的结果也与它们一致），这说明，北京地区气候变量序列的多分形标度行为，是由内在的长程相关性引起的，而这种长程幂律相关的性质，被认为是气候可预报性的重要依据.
|
图 5 随机洗牌后的lg（Fq（s））和lgs函数关系图 （a）平均气温；（b）气温日较差.不同阶数的序列在垂直方向作一定的移动以保证清楚可视. Fig. 5 Log-Log plots of the MF-DFA curves for the shuffled data of daily mean temperature (TMEAN) (a), diurnal temperature range (DTR) (b) From the top to the bottom curves correspond to different q (from q=-6 to q=+ 6) and are shifted vertically for clarity. |
本文通过多分形去趋势波动分析方法，和引入一个扩展的二项式串级模式，定量地分析了北京地区不同气候变量的多分形结构特征.使得用两个参数a、b就可以定量估计多分形谱宽Δα的具体值，从而很好地描述多分形结构特征.得到如下结论：
（1）多分形标度行为，在除气温日较差、日照时数以外的气候变量中广泛存在，对应地，代表序列中奇异性含量的多分形谱，表现出较宽的分布.而气温日较差和日照时数则表现为单分形的特征，对应地，多分形谱为窄的分布.日照时数的单分形特征，说明温压湿风等常规气象场的奇异性，并不是由于辐射的直接强迫，而是源于环流场的变化，锋面气旋系统带来的非对称及突变性是可能的原因[2, 25].而气温日较差的单分形特征，也说明了环流场中这种奇异性在最高气温和最低气温中具有位相的一致性.
（2）通过引入一个扩展的二项式串级模式，可以很好地模拟不同气候变量的多分形标度行为.同时，通过引入这个模式，可以定量地计算描述多分形行为的奇异谱.结果表明，具有多分形特征的变量序列对应较宽的奇异谱，而具有单分形特征的变量序列对应较窄的奇异谱.
（3）通过对原始序列洗牌，揭示了所分析气候变量的多分形行为，是由大振幅脉动与小振幅脉动间不同的长程相关性引起的，因此定量地描述多分形标度行为，为气候可预报性建立了理论基础.
| [1] | Dmowska R, Saltzman B. Advance in Geophysics:Long-range Persistence in Geophysical Time Series. San Diego, London, Boston, New York, Sydney, Tokyo, Toronto: Academic Press, 1999 : 14 -16. |
| [2] | Feng T, Fu Z T, Deng X, et al. A brief description to different multi-fractal behavior of daily wind speed records over China. Physics Letters A , 2009, 373: 4134-4141. DOI:10.1016/j.physleta.2009.09.032 |
| [3] | 时少英, 刘式达, 付遵涛, 等. 天气和气候的时间序列特征分析. 地球物理学报 , 2005, 48(2): 259–264. Shi S Y, Liu S D, Fu Z T, et al. The characteristic analysis of weather and climate time series. Chinese J. Geophys. (in Chinese) , 2005, 48(2): 259-264. |
| [4] | 荣平平, 刘式达. 不同时间尺度下气候变化基本特征的探索. 气候与环境研究 , 1997, 2(1): 77–82. Rong P P, Liu S D. An investigation of basic characteristics of climate change on different time scales. Climatic and Environmental Research (in Chinese) , 1997, 2(1): 77-82. |
| [5] | 刘式达. 地球系统模拟和混沌时间序列. 地球物理学报 , 1990, 33(2): 144–153. Liu S D. Earth system modeling and chaotic time series. Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese) , 1990, 33(2): 144-153. |
| [6] | 陈辉, 郭世昌. 昆明地区气候变化的多分形特征. 气候与环境研究 , 1997, 2(4): 361–368. Chen H, Guo S C. The multifractal characteristics of climate change at Kunming. Climatic and Environmental Research (in Chinese) , 1997, 2(4): 361-368. |
| [7] | 郭世昌, 常有礼, 陈辉, 等. 中国云南地区气候变化的多分维分析. 云南大学学报(自然科学版) , 2004, 26(4): 325–330. Guo S C, Chang Y L, Chen H, et al. The multifractal characteristics of climate change in Yunnan, China. Journal of Yunnan University (in Chinese) , 2004, 26(4): 325-330. |
| [8] | 严绍瑾, 彭永清, 张运刚. 一维气温时间序列的多重分形研究. 热带气象学报 , 1996, 12(3): 207–211. Yan S J, Peng Y Q, Zhang Y G. Study on multifractal of 1-D temperature time series. Journal of Tropical Meteorology (in Chinese) , 1996, 12(3): 207-211. |
| [9] | 江田汉, 邓莲堂. 全球气温变化的多分形谱. 热带气象学报 , 2004, 20(6): 673–678. Jiang T H, Deng L T. Multifractal spectra of global temperature. Journal of Tropical Meteorology (in Chinese) , 2004, 20(6): 673-678. |
| [10] | Peng C K, Buldyrev S V, Havlin S, et al. Mosaic organization of DNA nucleotides. Phys. Rev. E , 1994, 49(2): 1685-1689. DOI:10.1103/PhysRevE.49.1685 |
| [11] | Kurnaz M L. Application of detrended fluctuation analysis to monthly average of the maximum daily temperatures to resolve different climates. Fractals , 2004, 12: 365-373. DOI:10.1142/S0218348X04002665 |
| [12] | Chen X, Lin G X, Fu Z T. Long-range correlations in daily relative humidity fluctuations:a new indexes to characterize the climate regions over China. Geophysical Research Letters , 2007, 34: L07804. |
| [13] | Lin G X, Chen X, Fu Z T. Temporal-spatial diversities of long-range correlation for relative humidity over China. Physica A , 2007, 383: 585-594. DOI:10.1016/j.physa.2007.04.059 |
| [14] | Govindan R B, Kantz H. Long-term correlations and multifractality in surface wind speed. Europhys. Lett. , 2004, 68: 184-190. DOI:10.1209/epl/i2004-10188-3 |
| [15] | Kavasseri R G, Nagarajan R. Evidence of crossover phenomena in wind-speed data. IEEE Transaction on Circuits and Systems I , 2004, 51: 2255-2262. DOI:10.1109/TCSI.2004.836846 |
| [16] | Ivanova K, Ausloos M. Application of the detrended fluctuation analysis (DFA) method for describing cloud breaking. Physica A , 1999, 274: 349-354. DOI:10.1016/S0378-4371(99)00312-X |
| [17] | Kantelhardt J W, Zschiegner S A, Koscielny-Bunde E, et al. Multi-fractal detrended fluctuation analysis of nonstationary time series. Physica A , 2002, 316: 87-114. DOI:10.1016/S0378-4371(02)01383-3 |
| [18] | Kantelhardt J W, Rybski D, Zschiegner S A, et al. Multi-fractality of river runoff and precipitation:comparison of fluctuation analysis and wavelet methods. Physica A , 2003, 330: 240-245. DOI:10.1016/j.physa.2003.08.019 |
| [19] | Vitanov N K, Yankulova E D. Multi-fractal analysis of the long-range correlations in the cardiac dynamics of Drosophila melanogaster. Chaos Solitons Fractals , 2006, 28: 768-775. DOI:10.1016/j.chaos.2005.08.082 |
| [20] | Telesca L, Lapenna V, Macchiato M. Multi-fractal fluctuations in seismic interspike series. Physica A , 2005, 354: 629-640. DOI:10.1016/j.physa.2005.02.053 |
| [21] | Matia K, Ashkenazy Y, Stanley H E. Multi-fractal properties of price fluctuations of stocks and commodities. Europhys. Lett. , 2003, 61: 422-428. DOI:10.1209/epl/i2003-00194-y |
| [22] | 苟学强, 张义军, 董万胜, 等. 基于小波的地闪首次回击辐射场的多重分形分析. 地球物理学报 , 2007, 50(1): 101–105. Gou X Q, Zhang Y J, Dong W S, et al. Wavelet-based multifractal analysis of the radiation field of first return stroke in cloud-to-ground discharge. Chinese J. Geophys. (in Chinese) , 2007, 50(1): 101-105. |
| [23] | 李春峰, ChristopherLiner. 基于小波多尺度分析的奇性指数:一种新地震属性. 地球物理学报 , 2005, 48(4): 882–888. Li C F, Christopher Liner. Singularity exponent from wavelet-based multiscale analysis:a new seismic attribute. Chinese J. Geophys. (in Chinese) , 2005, 48(4): 882-888. |
| [24] | Feder J. Fractals. New York: Plenum Press, 1988 . |
| [25] | Ashkenazy Y, Feliks Y, Gildor H, et al. Asymmetry of daily temperature records. J. Atmos. Sci. , 2008, 65: 3327-3336. DOI:10.1175/2008JAS2662.1 |