气候系统翻转点的早期预警方法研究

梅莹 万仕全 谢孝强 顾斌 何文平

(1 南京信息工程大学 物理与光电工程学院，南京 210044;2 南京信息工程大学 空间天气研究所，南京 210044; 3 中山大学 大气科学学院，广东 珠海 519802;4 扬州市气象局，江苏 扬州 225009)

引 言

气候系统的长期演变记录表明，突然的气候变化会严重影响人类社会和生态环境，甚至会导致朝代的更替[1-5]。因此，对气候突变进行预测预警对提前采取减缓和应对措施具有重要意义[6]。然而，气候系统是复杂且非线性的，目前的气候预测理论和方法都难以对其突变进行较为准确的预测[7-11]。对各种复杂系统(如气候系统，生态系统等)临界阈值即临界点的认识[12-15]，使准确预测成为可能。当系统的各子系统的相互作用强度、控制参数或外强迫的强度达到某个阈值时，系统的状态会发生根本性的变化。譬如，每种气态物质都具有一个特定的温度，高于这个温度，无论如何增大压强，气态物质都不会液化，该温度值就是其气液状态转换的临界温度。又如，已有研究表明，大西洋经向翻转环流(Atlantic Meridional Overturning Circulation， AMOC)有多个平衡态，当注入北大西洋高纬度地区的淡水量逐步增加后，AMOC强度随之减弱，而当淡水量增加到一个临界阈值后，AMOC可能会发生突然中断[16]。这将会使区域尺度甚至更大空间尺度的气候状态发生突变。因此，在一个系统趋近其临界阈值的过程中，研究其是否有某种一般性的变化规律，将有助于提前捕捉气候突变来临前的早期预警信号。

近年来的一些研究表明，当一个动力学系统缓慢地趋近其临界点时，系统的一些统计物理量会发生明显变化，可作为突变来临前的早期预警信号，例如减小的恢复速率，变化的偏度和峰度，增大的自相关和方差等[17-23]。然而，各种方法都有一定的适用性。例如，强噪声会导致基于偏度系数和自相关系数的早期预警信号失效[21-22]。因此，发现新的早期预警方法仍是迫切问题。自然界中的许多系统都存在着长程相关性，如水文系统、气候系统、经济系统、生理系统[24-28]等。这意味着气候系统并非是随机演变，而是具有某种自组织性[29-30]。因此，在一个动力学系统趋近其临界翻转点的过程中，其长程相关性是否会有某种一般性的变化特征是一个值得深入研究的问题。若有，则可将其视为一个动力系统发生突变(由临界转换所导致)的早期预警信号。鉴于此，本文利用零维气候模型和两个折叠模型开展了临界转换导致的突变的数值试验，研究了3种不同动力学系统在趋近其自身的临界翻转点的过程中系统长程相关性的一般性变化特征。

1 方法与模型

1.1 去趋势波动分析法

去趋势波动分析法(Detrended Fluctuation Analysis，DFA)是一种能够较好地从时间序列中估计系统长程相关性的方法之一[31-32]，计算步骤如下：

(1)将长为N的时间序列的距平序列X(i)进行累积，得到累积和序列Y(i)：

(1)

(2)将Y(i)划分为长度为s的Ns=[N/s]个不重叠的等长子序列，为保证序列信息的完整性，从序列末端反向再划分一次，共得到2Ns个等长的子序列；

(3)对每个子序列v的数据进行多项式回归拟合，得到相应的局部趋势函数yv(i)，可以是一次、二次或更高次多项式(记为DFA1，DFA2，…)。随后，去除各个子序列内的局部趋势，再计算2Ns个去除子区间局部趋势后的残差序列的平方和的均值，进而开方得到波动函数F(s)，本文试验中采用的是DFA2；

(4)逐步增大s，重复(1)和(3)，可得到一个随窗口s变化的波动函数F(s)，进而获得二者的双对数曲线，若其近似的满足线性分布特征，则表明被分析的时间序列具有分形特征，然后计算其斜率。在非线性动力学中，其斜率通常称为Hurst指数。

当Hurst指数H=0.5时，表明系统的演变是随机的，不存在相关性；当H>0.5时，系统的演变存在长期记忆性，即未来的趋势与其过去的演变趋势具有某种相似性；当0

1.2 折叠模型

为了揭示一个动力系统在趋近其翻转点的过程中系统长程相关性的一般性规律，本文着重研究了当一个动力系统的关键控制参数趋近其临界转换阈值的情形。为此，考虑了3个典型的折叠模型，包括两个植被模型和一个简化的零维气候模型。第一个随机动力学模型为改进的单变量Logistic模型，这里主要用于描述植被生物产量与放牧率之间的关系[14]，计算式如下：

(2)

其中：c为放牧率，是植被系统的关键控制参数；系统状态变量V为植被生物量密度；ηV(t)代表外部随机扰动；σV为扰动的强度；植被生长率r，取值为1；环境承载力K，取值为10；参数V0=1，其物理含义为在V=V0时，放牧率为最大放牧率值的一半；系统的临界点由最大放牧率c*决定，c*=r(1+K2/4)/K。

模型二为一个双变量植被模型[14]，包含两个相互耦合的状态变量B(植被生物量)和w(土壤含水量)，计算式如下：

(3)

(4)

其中：R为降水率，为系统的关键控制参数，取值范围为0～3；σ为土壤失水率，其大小为1.0；-αw代表水分蒸发量；λ为植被耗水率，取值为0.12；-λwB则表示植被对水分的吸收量，即植被摄水量；ρ为最大生物增长率，取值为1；μ为最大放牧率，取值为2；Bc含义为生物量承载力，取值为10；B0为放牧率取最大放牧率的二分之一时的生物量，其值为1。公式(3)、(4)的最后一项均代表相应的外部随机扰动及其强度。

模型三是一个简化的零维气候模型[21，33]，描述的是全球气候系统的辐射收支平衡，计算式如下：

(5)

其中：系统状态变量T为全球海洋表面平均温度；相对辐射强迫μ是模型的关键控制参数，在本文试验中的取值范围为0.97～1；t是时间项；热惯性常数c=0.05；εσT4为净输出长波辐射，ε表示有效辐射率(取值为0.69)，Stefan Boltzmann常数σ；太阳常数I0取1 360；行星反照率项为ap=a-bT，其中常数a、b分别取2.8、0.009，对状态变量T添加乘性随机外强迫σTηT，试验中随机外强迫的大小为σT=0.01和0.25。

本文中随机微分方程(组)采用简单的欧拉方法求解，其中模型一和模型二的单个序列长度均为20 000，控制参数以固定步长靠近其临界点，随机外强迫取固定值；考虑到试验过程中的随机性，为了验证算法的鲁棒性，每组试验都重复100次，然后给出了100次试验的DFA2分析结果。试验中，模型三的单个序列长度取为10 000。

2 数值试验

改进的单变量Logistic模型有两个稳定平衡态，分别对应着两个临界参数阈值，即c=1.8和c=2.6。当放牧率参数为1.8

图1 单变量植被模型分岔(实线：系统的稳定平衡态；虚线：不稳定平衡态，下同)Fig.1 The bifurcation diagram of the univariate vegetation model(the solid line represents the stable equilibrium state， and thedashed line denotes unstable one ；same as below)

首先采用去趋势波动分析方法分析了植被密度随时间演变的特征。图2分别给出了放牧率参数在靠近(c=2.2)和远离(c=1.2)临界翻转点(c=2.6)时植被密度的DFA2分析结果。由此可见，在两种情形下，植被密度的波动函数的双对数曲线均近似满足直线分布，其Hurst指数分别约为1.112和0.89。表明单变量植被模型的植被密度的变化呈现出长程相关的特征。此外，对其他参数下的植被密度随时间的演变序列也采用DFA2进行了分析，获得了类似结果。这表明长程相关性特征是单变量植被模型植被密度随时间变化的一个固有的变化特征。此外，靠近翻转点附近时(即c=2.2)的Hurst指数要较远离翻转点时(c=1.2)更大。由此引出：当单变量植被模型的放牧率参数c趋近于其临界转换点的过程中，植被系统的长程相关性是否会持续性增大？如果是，这种长程相关性的持续性增大是否为一般性的变化规律，而与模型具体的数学表达式无关？

图2 植被密度的DFA2分析结果(α为Hurst指数，下同)：(a)靠近临界点(c=2.2)；(b)远离临界点(c=1.2)Fig.2 The DFA2 results for the vegetation biomass in univariate vegetation model (α: Hurst exponent; same as below): (a) close to the bifurcation point (c=2.2); (b) far from the bifurcation point (c=1.2)

为了回答上述问题，本文利用不同模型进行了较为系统的研究，主要考虑了仅仅通过改变动力系统的关键控制参数，使得系统缓慢地趋近其临界转换点。在单变量植被模型的试验中，考虑了在两种不同随机强迫强度下，仅仅沿着其分岔图的上支(图1)逐步线性地增大放牧率参数c，使得系统缓慢地趋近其临界翻转点(c=2.6)。在此过程中，系统的其他参数和外强迫的强度保持为常数。DFA2的分析结果表明，无论随机强迫强度相对小还是大的情形，在单变量植被系统不断靠近临界点的过程中，系统的Hurst指数也在不断增大，且越接近临界点，其增加速率越快(图3)。为了定量衡量Hurst指数变化趋势是否显著，本文引入肯德尔秩相关系数(τ)。100次重复试验的结果表明，其平均Hurst指数的肯德尔秩相关系数分别为τ=0.981 1(随机强迫强度为σ=0.05时)和τ=1(随机强迫强度为σ=0.25时)，增加趋势均是统计显著的(显著性水平均为α=0.05)。且100次重复试验的结果的离散度也较小，这表明DFA2分析获得的长程相关性的增大趋势可以作为单变量植被模型趋近其系统翻转点的早期预警指示信号。

图3 单变量植被系统的长程相关性随着放牧率参数变化：(a)相对较小的随机强迫强度σ=0.05；(b)相对较大的随机强迫σ=0.25Fig.3 The long-range correlation as a function of the grazing rate c in univariate vegetation model: (a) relatively small magnitude of the stochastic forcing(σ=0.05); (b) relatively large magnitude of the stochastic forcing(σ=0.25)

图4给出了双变量植被模型的分岔。其中，系统的关键控制参数为降水率R，实线代表系统的两种稳定状态：有植被状态(上支)和荒漠化状态(下支)，点线表示不稳定状态。系统状态变量为植被生物量B和土壤含水量w，此处以植被生物量B为研究对象。图5给出了双变量植被模型在降水率与随机外强迫强度取不同值时植被生物量的DFA2分析结果，发现其均能够近似地满足直线分布。在相同随机强迫强度下，随着靠近系统的临界翻转点，系统的长程相关性明显增强(Hurst指数变大)，是否这是双变量植被模型中的一般性变化规律？

图4 双变量植被模型分岔(两个临界降水率阈值R=1.06和R=2.0)Fig.4 The bifurcation diagram of the bivariate vegetation modelwith two critical points:R=1.06 and R=2.0

图5 当降水率参数远离和靠近临界翻转点(R=1.06)时，植被密度的DFA2分析结果：(a)远离临界点，R=2.0，外强迫强度σB=0.05；(b)靠近临界点， R=1.3，外强迫强度σB=0.05；(c)远离临界点，R=2.0，外强迫强度σB=0.25；(d)靠近临界点， R=1.3，外强迫强度σB=0.25Fig.5 The DFA2 results for the vegetation biomass in bivariate vegetation model when the rainfall rate R is (a，c; R=2.0) far from or (b，d;R=1.3) close to the bifurcation point (R=1.06)

类似于单变量植被模型，这里仍然考虑了在两种不同随机强迫强度下，仅仅沿着双变量植被模型分岔图的上支(图4)逐步线性地减小降水率参数R，使得系统趋近其临界翻转点(R=1.06)。在此过程中，系统的其他参数和外强迫的强度保持为常数。数值试验的结果显示，无论对于随机强迫强度相对小或大的情形，在双变量植被系统不断靠近临界点的过程中，系统的Hurst指数也在不断地增大，且越接近临界点，Hurst指数增加速率越快(图6)。双变量植被系统长程相关性的增大趋势也是显著的(显著性水平为α=0.05)。这与由单变量植被模型试验获得的结果一致，证实了长程相关性的增大趋势可以作为双变量植被模型趋近其系统翻转点的早期预警指示信号。

图6 双变量植被系统的长程相关性随着降水率变化：(a)相对较小的随机强迫强度，σB=0.05；(b)相对较大的随机强迫σB=0.25Fig.6 The long-range correlation as a function of the rainfall rate c in bivariate vegetation model: (a) relatively small magnitude of thestochastic forcing (σB=0.05); (b) relatively large magnitude of the stochastic forcing (σB=0.25)

太阳辐射相对强度的减小会使得全球海表温度降低，当其减小到某个临界阈值后，大气温度将会发生一次突然的转变。零维气候模型中，如持续减小相对辐射强迫参数μ，当其约减小到系统的临界翻转点(μ*≈0.97)附近时，温度将会由一个相对暖的状态突然转变到一个相对寒冷的状态。图7给出了零维气候模型在相对辐射强迫参数μ取不同值时海表平均温度T的DFA2分析结果，可以发现平均温度的演化具有很好的分形特征，且在靠近临界翻转点时，系统的Hurst指数明显较远离翻转点时更大。本文试验逐步减小相对辐射强迫强度使得零维气候模型趋近其临界翻转点，借此考察了系统的长程相关性特征有何种变化。研究显示，零维气候模型的DFA2分析结果与前两个模型试验结果基本一致，即随着系统越来越靠近临界点，系统的Hurst指数逐步变大，且越靠近临界翻转点这种增大的趋势越明显(图8)。值得指出的是，与前两个模型相比，零维气候模型在远离临界点时，Hurst指数的变化更为平缓些。两种情况下Hurst指数的均值序列的肯德尔秩相关系数相同(均为0.966 7)，而且Hurst指数增加是统计显著的(显著性水平为α=0.05)。

图7 零维气候模型平均温度的DFA2分析结果：(a)远离临界点(μ=0.99)；(b)靠近临界点(μ=0.97)Fig.7 The DFA2 results for a simple zero-dimensional climate model when the relatively radiative forcing is (a，μ≈0.99) far from or (b，μ≈0.97) close to the bifurcation point (μ*≈0.97)

图8 零维气候模型中平均温度的Hurst指数随相对辐射强迫变化:(a)随机外强迫；(b)随机外强迫 Fig.8 The long-range correlation as a function of the relatively radiative forcing in a simple zero-dimensional climate model：(a) relatively small magnitude of the stochastic forcing(σ=0.01);(b) relatively large magnitude of the stochastic forcing (σ=0.25)

3 结论

自然界中的许多系统都存在着长程相关性，这已被大量的科学研究所证实，如气候系统[30]。因此，在一个系统趋近其自身临界翻转点的过程中，这种长程相关性是否有某种一般性的变化特征是一个值得研究的科学问题。为此，本文利用3个典型的折叠模型，研究了动力系统在趋近其临界翻转点的过程中，系统长程相关性的变化特征。结果显示，在动力系统的控制参数缓慢的趋近其临界翻转点(越过该点，系统的状态将发生突然的、截然不同的转换)的过程中，表征系统长程相关性的Hurst指数呈现持续增加的趋势，且这种趋势具有统计显著意义(显著性水平α=0.05)。无论是对于单变量植被模型还是双变量植被模型，或者零维气候模型，这一结论均成立。这意味着系统长程相关性的增大，可能是一个动力系统趋近其自身临界翻转点的一般性变化规律，可将其用作各种动力学系统趋近其临界翻转点的早期预警信号。需要指出的是，对于那些快速趋近系统翻转点的情况，由于系统的控制参数等的变化时间较短，本文所提出的早期预警方法无法给出有效的预警，且现有其他的一些预警理论和方法也对此无能为力。