斜率单变点分析

王建锋

（中国科学院力学研究所，北京，100080）

【摘要】在自然界、社会、经济等领域内，突变现象很常见。若系统的输出序列在某未知时刻起了突然变化，则该时刻即称为变点。变点统计分析的目的是判断变点的存在，确定其位置和个数等。已有的变点分析包括均质变点分析、概率变点分析和模型参数变点分析等。本文提出斜率变点的新概念，它是指曲线斜率加（减）速变化最大的点。本文还结合几个不同类型的实例，提出了寻求单一斜率变点的回归系数二阶差分方法。它可对单调性和凹凸性均单一的曲线求其“转折点”。实例表明，该法具有简单、直观、有效等优点。

【关键词】变点分析　斜率变点　回归系数

1　前言

在自然界、社会、经济等领域内，突变现象很常见，且很重要。研究突变是否发生，何时发生等问题，有助于掌握事件或过程的演化规律，尤其是灾害发生发展规律，从而为灾害的预报、预防和治理提供依据。

系统的输出序列在某未知时刻起了突然变化，该时刻即称为变点。变点统计分析的目的是判断和检验变点的存在、位置、个数，并估计出变点的跃度。变点统计分析对定量分析各种监测数据，研究各种地质灾害的规律等，都是个强有力的工具。

变点分析又分为均值变点分析、概率变点分析以及模型参数变点分析等[1]。某一时刻前、后数据的均值（或概率分布，或某模型参数）发生了显著改变，则该时刻就称为均值（或概率，或某模型参数）变点。但地质问题中还常遇到一条曲线逐渐上升，在某时刻以后突然加速上升；有时遇到一条曲线开始加速下降，到某点以后突然减缓速度下降而逐渐趋于平缓。还有的曲线开始缓速下降，到某点后突然转为快速沿斜线下降。由于这一突然转变的“转折点”是一类重要的特征点，往往与特定问题的特定物理意义相联系，因此准确地确定该变点是重要的。此类“转折点”，本文称其为曲线的斜率变点，简称斜率变点。

本文拟结合岩土体表面蠕变曲线如何寻找第二阶段和第三阶段分界点的问题，求土体特性指标的相关距离δ（或θ）中如何寻找方差折减函数的突然趋缓的“转折点”问题，以及根据e-lgP曲线求前期固结压力P。的问题，来进行斜率变点分析论述。

2　寻找单一斜率变点的基本原理

先研究最简单的斜率变点问题。假定已知一条曲线中有且只有一个斜率变点，问题是如何找到一种方法可简单、定量、准确地确定此变点。

一般的，大多数的测量数据都呈离散的数据点对，并且当测量时间间隔较长时，不易形成能真实反映过程变化的连续曲线。因此，多数观测序列不易用曲线方程表达，因而就无法用微积分求导数的方法求出各点处的斜率，只能用求某点两侧一些相继数据点的线性回归系数的方法近似地求出某时刻点（或距离点）两侧较短时间间隔（或距离）内曲线的斜率。这是因为在较短的时间间隔（或距离）内，曲线是近似直线的。某点两侧曲线斜率之差可反映该点两侧斜率变化的幅度，这可说是用到了一阶差分。但这里拟找的“转折点”并不是斜率变化最大的点，而是斜率局部加（减）速变化最大的点。在监测时间是等时间间隔（或等间距）的条件下，斜率变化幅度的二阶差分就反映了斜率变化的加（减）速度。通过寻找斜率变化的加（减）速度局部极大值的方法，就可找出斜率变点所落在的单位区间。然后，再用类似于求分组数据众数的方法，便可定量地求出斜率变点的估计值了。

3　寻找单一斜率变点的方法和步骤

以下结合3个不同类型的实例，介绍寻找单一斜率变点的方法和步骤。

例1.日本Tohoku铁道线上Asamushi滑坡的累计位移监测数据[21如表1所示（数据由图1量出）。

表1　日本Tohoku铁路线 Asamushi滑坡累计位移

图1　日本Tohoku铁路线 Asamushi滑坡时间—累计位移曲线变点分析结果

这是由滑坡表面测得的位移—时间曲线，它明显包含了第二、第三蠕变阶段的数据。第二阶段的曲线段近似直线，反映出等速蠕变特征；但第三阶段则是加速蠕变阶段，其曲线明显加速上扬（图1）。这就是说，在这个例子中，已知包含且只包含一个斜率变点，只要能找到斜率局部加速变化最大的点，那么就意味着找到了这个斜率变点了。以下分步介绍寻找这个变点的回归系数二阶差分法。

（1）取定探索点：由于监测数据是等时间间隔的，为此选取两个相邻观测时间点的中点为探索点，构成探索点序列。如例1中的探索点序列为t_i=21.25,21.75,22.25,22.75,23.25，…，26.75。

（2）以各探索点为中心构造滑动窗口，以便计算出探索点前后附近曲线的斜率（即探索点前后各若干数据点的线性回归系数）。由于参加回归的数据点数 n的多少会影响回归系数的取值，故在探索点前、后各取一样多（n）数据点构成滑动窗口，这样可在同等条件下进行前、后斜率的对比。又由于曲线只在较短的时间（或较小的距离）内才近似直线，故 n也不能取得太大。这里分别取 n=2,3,4构成三套滑动窗口。

（3）在以 t_i为滑动窗口中，对探索点 t_i前（或左）的n个数据点作线性回归，求出回归系数，记为

（t_i）；同样，对探索点 t_i后（或右）的n个数据点作线性回归，求出回归系数，记为

（t_i）。显然，n=2时算出的回归系数反映局部的斜率性态较多，反映整体的斜率性态则不够，又由于点数太少（只是两点，可连一直线），随机性较大，统计意义不够。相反，n=4时算出的回归系数反映整体的斜率性态较好，反映局部的的斜率性态则较差，又由于点数较多，统计意义较强，较少随机性；n=3则介乎二者之间。因此，应当更着重 n=4时的结果。因此，将n=2,3,4时计算出的

进行加权平均，权取 n²。于是，当对某个 t_i，

均存在时，加权平均数为：

地质灾害调查与监测技术方*文集

（4）对每个探索点 t_i，计算

之差，并记为∆S(t_i），即

地质灾害调查与监测技术方*文集

∆S(t_i）可说是t_i点前、后曲线斜率的增量（或改变量），也可理解为是t_i点后、前曲线斜率的一阶差分，它的大小反映了 t_i点处斜率增加的幅度。

（5）对∆S(t_i）的序列再计算二阶差分，即：

地质灾害调查与监测技术方*文集

这个二阶差分值的大小就反映了在区间（t_i-1,t_i）内曲线斜率加速度变化的大小。△²S(t_i）也构成一序列。

（6）沿着 t_i从小到大的序列，寻找²S(t_i）序列中出现的最大值（它比前、后两个值均大）。设它所对应的区间为（t_i-1,t_i），它就应是斜率变点所在的区间。再把它前、后两个相邻区间（t_i-2,t_i-1）及（t_i,t_i+1）和它们对应的△²S(t_i-1）、△²S(t_i+1）值利用起来，用与分组数据众数类似的方法，进行线性内插，即可求出斜率变点 t＊的精确值，计算公式如下：

地质灾害调查与监测技术方*文集

对例1数据用回归系数二阶差分法进行计算，将计算的中间结果及最终结果列于图2中。

从图2中可以看出，取得最大的二阶差分为²S(t_i)=23.63，它对应的区间为（23.25, 23.75），它的前、后两个相邻区间为（22.75,23.25）和（23.75,24.25）。它们对应的二阶差分值为△²S(t_i-1)＝18.80，△²S(t_i+1)=11.47。

按公式（4）可算得：

地质灾害调查与监测技术方*文集

从图1直观来看，t^＊作为该曲线的斜率变点是合适的，它可作为蠕变曲线第二阶段与第三阶段的分界点。找出这个分界点有重要实际意义，它可提醒观测者何时开始就应对该滑坡动态进行加密监测，也可使现场工程师有可能只用第三阶段的数据点，对滑坡发生时间作出更为精确的预报。因此，斜率变点的确定有助于解决滑坡蠕变曲线第二、第三阶段的划分问题。此外，它还可以帮助寻找曲线由快速下降到突然转缓，以及曲线由平缓下降突然转为快速沿斜线下降的“转折点”问题。

图2　日本Tohoku铁路线 Asamushi滑坡时间—累计位移曲线变点分析Spreadsheet配置

4　寻找曲线由快速下降到突然转缓趋于平稳的“转折点”

地质学中遇到这类问题较多。但过去往往只凭肉眼观察来人为确定一个“转折点”，没有客观定量分析的方法。而这种问题的解决不但对理论研究有帮助，而且具有重要的实际经济效益，如确定经济合理的取样间距、勘探网密度等。也就是俗称的寻找“平衡点”的问题。下面将结合求土体相关距离的实际例子来说明斜率变点分析在其中的应用。

例2.已知从加拿大温哥华及近郊Haney处土体静力触探试验数据计算得出滞后距 T与方差折减函数Г²(T）如表2、图3所示。

表2　Haney处 T与Γ²(T）数据表

例1中的曲线是开始等速上升，经过斜率变点后加速上扬。而例2中的曲线则是开始快速下降，到变点时突然减速下降，并逐渐趋缓，甚至趋于一条水平的渐近线。这里需要寻找的是突然减速下降的“转折点”。由于例2中曲线虽与例1中曲线在直观上很不相同，但都是凹曲线。因此，随着 T_i（或 t_i）的增加，曲线的斜率总是增大的。因此，一阶差分还是用公式：

但到计算二阶差分²S(T_i）时就与例1时有所不同。例1中序列S(T_i）一般是随着 ti的增大而递增的，故²S(T_i)＝S(t_i)－ S(t_i-1）；而例2中序列∆S（T_i）一般是随着 T_i的增大而递降的，故此时计算二阶差分就应倒过来减：

图3　Haney处方差折减函数Γ²（T）图

地质灾害调查与监测技术方*文集

其余求 T^＊的方法、公式与方程（4）类似。

对例2中数据的计算结果列于图4。从图4可以看出，取得最大值的二阶差分为²S（T_i）=0.1822，它对应的区间为（1.1,1.3）。它的前、后两个相邻区间为（0.9,1.1）和（1.1, 1.3），它们对应的二阶差分值为²S（T_i-1）＝0.1679和²S（T_i+1）=0.1314，按照（4）可得：

地质灾害调查与监测技术方*文集

这就是“转折点”的T的估计值。再根据表2中 T与Γ²（T）的数据，用线性内插法求出Γ²（T^*）来：

图4　Haney处方差折减函数Γ₂(T）变点分析结果

地质灾害调查与监测技术方*文集

∴Γ²（T^*）=0.3378-0.054×0.7195=0.2989

于是，Haney处土体相关距离δ≈T^＊·Г²（T^*）=1.1439×0.2989=0.3419（m）。这与原文中计算出的结果δ＝0.324m非常接近，相对误差只有5.5%。原文中 T^*=1.2，如果按本文量出的数据，则Г²（1.2）=0.2838，于是δ≈1.2×0.2838=0.34056（m），此与本文结果更为接近，相对误差仅为0.39%。原文中时 T^*=1.2时，Γ²（T^＊）＝0.27，故算出δ≈1.2×0.27=0.324（m）。可见，误差主要是本文量测数据造成的，与方法本身无关。可见，斜率变点方法提供了计算相关距离的另一有效途径。

5寻找曲线由慢速平缓下降到突然转为快速沿斜线下降的“转折点”

土力学中根据e-lgP高压固结实验曲线求取前期固结压力P_c，即是此类问题的典型代表。

例3.垂直压力P的对数与隙比的试验数据^[4]见表3、图5。首先运用线性内插使试验数据等间隔化。变点分析结果见图6。

表3　e-lgP数据表

图5　高压固结实验e-lgP曲线

图6　高压固结实验e-lgP曲线变点分析结果

由于曲线是凹的，对每个_Li算出的

总大于

因此由公式（2）可推：

地质灾害调查与监测技术方*文集

由于算出的∆S（L_i）序列基本上是递增的，故计算 △²S（L_i）的公式仍可沿用公式（3），即：

地质灾害调查与监测技术方*文集

从图6的最后一栏看出，最大的△²S（L_i）为0.0413，所对应的区间为（2.45,2.55）。它的前后两个相邻区间为（2.35,2.45）和（2.55,2.65），它们对应的二阶差分值为△²S(L_i-1)＝0.0361和△²S（Lⁱ⁺¹）=0.0137，按照（4）可得：

地质灾害调查与监测技术方*文集

即斜率变点的估计值为（lgP)^*=2.4658。因此，前期固结压力 P_c的估计值为

=102.4658=292.28kPa。原文P_c的计算结果为313.91kPa，两者极为接近，相对误差仅为6.89%。

6　结论

在自然界、社会、经济等问题的科学研究中，常遇到要寻求一条曲线的“转折点”问题。当曲线是单调递增（或单调递减），且是凹的（或凸的）曲线时，就意味着曲线中有且只有一个斜率变点存在，可用回归系数二阶差分法确定此变点。

当曲线不是单调且凹凸性有变化时，可以将曲线先划分为单调且凹凸性单一的曲线段，然后分段求取斜率变点。

回归系数二阶差分法求得的斜率变点是斜率最大加速的点，而不是斜率变化量最大的点。该点往往是在沿着近似直线的变化后，开始突然加速转向变为沿着另一条直线或曲线变化的“转折点”，它多在斜率变化量最大的点（或曲率最大的点）之前发生。

应用回归系数二阶差分法需要是等间距数列，且间隔要相对较密，即数据点对足够多，最好大于20或30对，至少也要有13对。这是因为运用滑动窗口时存在边缘效应，并且算出的回归系数不能很好地代表曲线的斜率。若原始数据是不等间距的，可以采用线性内插的方法使之转化为等间隔的较稠密的数据。

曲线单调增、减及凹、凸不同时，计算斜率变点的公式稍有不同：①单调增、凹性时，用公式（2）、（3），如例1；②单调降、凹性时，用公式（2）、（5），如例2；③单调降、凸性时，用公式（6）、（3），如例3；④单调增、凸性时，用公式（6）、（5）。

参考文献

[1]项静恬，史九恩著.非线性系统中数据处理的统计方法.北京：科学出版社，1997:1～43

[2]M Saito.Forecasting time of slope failure by tertiary creep.Proceedings of the Seventh International Conference on Soil Mechanics and Foundation Engineering,Mexico,1969,Vo1.2:677～683

[3]Damika Wickremesinghe and R G Campanella.Scale of fluctuation as a descriptor of soil variability.Probabilistic Methods in Geotechnical Engineering,Li＆Lo(eds),Balkema,Rotterdam,1983:233～239

[4]Jean-Pierre Bardet.Experimental soil mechanics.Prentice-Hall,Inc.,1997:297～306

[5]陈希孺.变点统计分析简介.数理统计与管理，vol.No.1～4∶1～43

[6]Chen X R.Inference in a simple change-point model.Science in China,Ser.A,1988,(6)

[7]Keyue Ding.A lower confidence bound for the change point after a sequential CUSUM test.Journal of Statistical Planning and Inference,2003:115,311～326

[8]Douglas M Hawkins.Fitting multiple change-point models to data,Computational Statistics＆Data Analysis,37,2001:323～341

[9]Wayne A Taylor.Change-point analysis:a powerful new tool for detecting changes.2000,http://www.variation.com/cpa/tech/changepoint.html

[10]Wayne A Taylor.A pattern test for distinguishing between autoregressive and mean-shift data.2000,http://www.variation.com/cpa/tech/changepoint.html

本文如未解决您的问题请添加抖音号：51dongshi（抖音搜索懂视），直接咨询即可。