斜率的单变点分析
中国科学院力学研究所，北京，100080)
巨灾在自然、社会和经济领域非常普遍。如果系统的输出序列在某个未知时刻突然发生变化，则该时刻称为变点。变化点统计分析的目的是判断变化点的存在，确定其位置和数量。现有的变点分析包括均值变点分析、概率变点分析和模型参数变点分析。本文提出了坡度变化点的新概念，它是指曲线坡度的加速度(减速度)变化最大的点。结合几个不同类型的例子，本文还提出了用回归系数的二阶差分法来寻找单个斜率变化点。它可以找到具有单调性和凹度的曲线的“转折点”。实例表明，该方法简单、直观、有效。
变点分析；斜率；变点；回归系数
1前言
在自然、社会和经济领域，突变是非常普遍和重要的。研究突变是否发生以及何时发生，有助于把握事件或过程的演化规律，特别是灾害的发生发展规律，从而为灾害预测、预防和管理提供依据。
系统的输出序列在未知时刻突然发生变化，称为变点。变点统计分析的目的是判断和检验变点的存在、位置和数量，估计变点的跳变。变点统计分析是定量分析各种监测数据、研究各种地质灾害规律的有力工具。
变点分析分为均值变点分析、概率变点分析和模型参数变点分析[1]。某一时刻前后数据的均值(或概率分布，或模型参数)发生了显著变化，那么该时刻称为均值(或概率，或模型参数)变化点。但在地质问题中，往往有一条曲线是逐渐上升的，在某一时刻后突然加速；有时候遇到曲线开始加速，过了某个点后突然减速，逐渐趋于平缓。其他曲线开始缓慢下降，到某一点后突然转为沿对角线快速下降。因为这种突变的“转折点”是一个重要的特征点，它往往与具体问题的具体物理意义联系在一起，所以准确确定变化点很重要。这种“转折点”在本文中称为曲线的斜率变化点。
本文提出了如何根据岩土表面的蠕变曲线寻找第二阶段和第三阶段的分界点，如何在土性指标的相关距离δ(或θ)中寻找方差折减函数的“转折点”，如何根据e-lgP曲线寻找前期固结压力P。分析和讨论斜率变化点。
2寻找单一斜率变化点的基本原理
首先，研究了最简单的斜率变点问题。假设一条曲线只有一个斜率变化点，问题是如何简单、定量、准确地找到确定这个变化点的方法。
一般来说，测量数据大多是离散的数据点对，当测量时间间隔较长时，很难形成能够真实反映过程变化的连续曲线。因此，大部分观测序列不容易用曲线方程来表示，所以不能用微分学的方法得到每一点的斜率，而在短时间间隔(或距离)内某一时间点(或距离点)两侧曲线的斜率只能通过求某一点两侧一些连续数据点的线性回归系数来近似得到。这是因为在短时间间隔(或距离)内，曲线近似是直的。一点两侧曲线的斜率之差可以反映该点两侧斜率变化的幅度，可以说是利用了一阶差分。但这里要找的“转折点”并不是坡度变化最大的点，而是坡度局部加速(减速)变化最大的点。在监测时间为等时间间隔(或等间隔)的条件下，斜率变化幅度的二阶差反映了斜率变化的加速度(减速度)。通过找到斜率变化的加速度(减速度)的局部最大值，可以找到斜率变化点所在的单位区间。然后，利用类似于寻找分组数据的模式的方法，可以定量地得到斜率变化点的估计值。
3寻找单一斜率变化点的方法和步骤
通过三个不同类型的例子介绍了寻找单一斜率变化点的方法和步骤。
示例1。日本东北铁路线上浅木石滑坡[21]累计位移监测数据见表1(数据由图1测得)。
表1日本东北铁路浅沼滑坡累积位移
图1日本东北铁路浅沼滑坡时间-累积位移曲线变点分析结果。
这是从滑坡面测得的位移-时间曲线，显然包含了第二、三蠕变阶段的数据。第二阶段曲线段近似为直线，体现等速爬行的特点；而第三阶段是加速蠕变阶段，其曲线明显加速(图1)。也就是说，在这个例子中，已知只有一个斜率变化点。只要能找到坡度局部加速度变化最大的点，就说明找到了这个坡度变化点。逐步介绍了寻找这一变化点的回归系数二阶差分法。
(1)选取勘探点:由于监测数据是等时间间隔的，所以选取相邻两个观测时间点的中点作为勘探点，形成勘探点序列。比如1中的勘探点顺序是Ti = 21.25，21.75，22.25，22.75，23.25，…，26.75。
(2)以每个探索点为中心构造滑动窗口，从而计算出探索点周围曲线的斜率(即探索点前后几个数据点的线性回归系数)。因为参与回归的数据点数n会影响回归系数的取值，所以在勘探点前后取相同数量(n)个数据点形成滑动窗口，这样就可以在相同条件下比较前后的斜率。而且由于曲线只是在短时间内(或者短距离内)近似为直线，所以n不能太大。这里，分别取n = 2、3和4形成三组滑动窗口。
(3)在有ti的滑动窗口中，对探索点ti之前(或左侧)的N个数据点进行线性回归，得到回归系数，记为(Ti)；同样，对勘探点ti后面(或右边)的N个数据点进行线性回归，得到回归系数，记为(ti)。显然，当n=2时，计算出的回归系数反映了更多的局部斜坡行为，而不是足够的整体斜坡行为，并且由于点太少(只有两个点可以用直线连接)，所以它是随机的，没有足够的统计意义。相反，n=4时计算的回归系数较好地反映了整体边坡性状，但反映局部边坡性状较差，且由于点数较多，具有更强的统计意义，随机性较小；N=3介于两者之间。所以我们更应该关注n=4时的结果。因此，对n = 2、3、4时计算的值进行加权平均，权重为n2。因此，当某个ti的所有值都存在时，加权平均值为:
地质灾害调查与监测技术方法论文集
(4)对于每个勘探点ti，计算差值并记录为？S(ti)，即
地质灾害调查与监测技术方法论文集
S(ti)可以说是ti点前后曲线斜率的增量(或变化量)，也可以理解为ti点前后曲线斜率的一阶差，其大小反映了ti点处斜率的增加量。
(5)对吗？S(ti)序列，然后计算二阶差分，即:
地质灾害调查与监测技术方法论文集
这个二阶微分值的大小反映了区间(ti-1，ti)内曲线斜率加速度变化的大小。△2S(ti)也构成一个序列。
(6)沿着ti从小到大的顺序，在2S (ti)序列中寻找最大值(大于前两个值和后两个值)。设其对应区间为(ti-1，ti)，应该是斜率变化点所在的区间。然后利用相邻的两个区间(ti-2，ti-1)和(ti，ti+1)及其对应的值△2S(ti-1)，采用分组数据的方式。
地质灾害调查与监测技术方法论文集
实例1的数据采用回归系数二阶差分法计算，计算的中间和最终结果见图2。
从图2可以看出，最大二阶差分为2S(ti)=23.63，其对应的区间为(23.25，23.75)，其相邻的两个区间为(22.75，23.25)和(23.75，24.25)。它们对应的二阶微分值为△ 2s (Ti-1) = 18.80，△2s(Ti+1)= 11.47。
根据公式(4)，可以计算如下:
地质灾害调查与监测技术方法论文集
直观上，t *适合作为曲线的斜率变化点，它可以作为蠕变曲线第二阶段和第三阶段的分界点。找出这个分界点具有重要的现实意义，它可以提醒观测者何时开始对滑坡动态进行密集监测，也可以使现场工程师仅利用第三阶段的数据点就能对滑坡发生时间做出更准确的预测。因此，斜率变化点的确定有助于解决滑坡蠕变曲线第二、三阶段的划分问题。此外，还可以帮助寻找曲线从快速下降到突然放缓的“拐点”，以及沿对角线从平缓下降到快速下降的曲线。
图2日本东北铁路浅沼滑坡时间-累积位移曲线变点分析
4.寻找曲线从快速下降到突然放缓并趋于稳定的“拐点”。
地质学中有很多这样的问题。但在过去，一个“转折点”仅通过肉眼观察来人为确定，没有客观的定量分析方法。该问题的解决不仅有助于理论研究，而且具有重要的实际经济效益，如确定经济合理的采样间距和勘探网密度。也就是俗称的寻找“平衡点”的问题。下面将通过一个寻找土壤相关距离的实例来说明斜率变化点分析的应用。
例2。已知从加拿大温哥华和哈尼的土的静力触探试验数据计算的滞后距离t和方差折减函数г 2 (t)见表2和图3。
表2 Haney的t和γ 2 (t)数据表
例1中的曲线开始匀速上升，然后在变坡点后加速。而例2中的曲线开始快速下降，然后在到达变化点时突然变慢，并逐渐变慢，甚至趋于水平渐近线。这里需要寻找的是突然减速下降的“拐点”。虽然例2中的曲线与例1中的曲线差别很大，但都是凹曲线。所以，随着ti(或Ti)的增大，曲线的斜率总是增大的。所以一阶差分还是用公式:但是在计算二阶差分2S(Ti)时，和例子1不一样。在例1中，序列S(ti)一般随着ti的增加而增加，所以2S(Ti)= S(Ti)-S(Ti-1)；以及示例2中序列？S(Ti)一般随着Ti的增加而减小，所以二阶差分要反过来计算:
图3 Haney方差约简函数γ 2 (t)图
地质灾害调查与监测技术方法论文集
其他求t *的方法和公式和方程(4)类似。
实施例2中数据的计算结果如图4所示。从图4可以看出，最大值的二阶差分为2S(Ti)=0.1822，其对应的区间为(1.1.1.3)。它的两个相邻区间为(0.9，1.1)和(1.1，1.3)，它们对应的二阶差值为2s (Ti-1) = 0.10。
地质灾害调查与监测技术方法论文集
这是t在“转折点”的估计值。然后，根据表2中t和γ 2 (t)的数据，通过线性插值得到γ 2 (t * ):
图4 Haney处方差异缩减函数γ 2 (t)的变点分析结果
地质灾害调查与监测技术方法论文集
∴γ2(t*)=0.3378-0.054×0.7195=0.2989
因此，Haney处土壤的相对距离为δгt *г2(t *)= 1.1439×0.2989 = 0.3419(m)。这与原文计算的结果δ = 0.324 m非常接近，相对误差仅为5.5%。原文中T*=1.2。如果用本文的数据来衡量，г 2 (1.2) = 0.2838，那么δ≈1.2×0.2838 = 0.34056(m)，更接近本文的结果，相对误差只有0。原文中，当T*=1.2时，γ 2 (t *) = 0.27，所以计算出δ≈1.2×0.27=0.324(m)。可以看出，误差主要是本文的实测数据造成的，与方法本身无关。可以看出，斜率变点法提供了另一种计算相关距离的有效方法。
5寻找曲线缓慢平缓下降到沿对角线突然转为快速下降的“转折点”。
根据土力学中的e-lgP高压固结实验曲线计算前期固结压力Pc就是这类问题的典型代表。
例3。垂直压力P的对数和间隙比的实验数据[4]见表3和图5。首先，使用线性插值使测试数据等距。变化点分析结果如图6所示。
表3 e-lgP数据表
图5高压固结试验的e-lgP曲线
图6高压固结试验e-lgP曲线变化点分析结果
因为曲线是凹的，所以为每个Li计算的总和大于该值，所以可以从公式(2)中推断出:
地质灾害调查与监测技术方法论文集
因为算计？S(Li)序列基本上是递增的，因此计算△2S(Li)的公式仍可遵循公式(3)，即:
地质灾害调查与监测技术方法论文集
从图6最后一列可以看出，最大△2S(Li)为0.0413，对应的区间为(2.45，2.55)。其相邻的两个区间为(2.35，2.45)和(2.55，2.65)，对应的二阶微分值为△ 2s (Li-1) = 0.0361和△ 2s (Li+1) = 0.065438。
地质灾害调查与监测技术方法论文集
即，斜率变化点的估计值是(lgP)*=2.4658。因此，预固结压力Pc的估算值为=102.4658=292.28kPa，原Pc的计算结果为313.91kPa，非常接近，相对误差仅为6.89%。
6结论
在自然、社会、经济等问题的科学研究中，我们经常会遇到寻找一条曲线的“转折点”的问题。当曲线单调递增(或递减)且凹(或凸)时，说明曲线只有一个斜率变化点，可用回归系数二阶差分法确定。
当曲线不单调且凹凸变化时，可将曲线分割成单一凹凸的单调曲线段，然后分段求斜率变化点。
回归系数二阶差分法得到的斜率变化点是斜率加速最大的点，而不是斜率变化最大的点。这个点往往是沿一条近似直线变化后突然加速转向沿另一条直线或曲线变化的“转折点”，通常发生在斜率变化最大的点(或曲率最大的点)之前。
回归系数二阶差分法的应用，需要间隔相对较近的等间距序列，即有足够多的数据点对，最好在20或30对以上，最少13对。这是因为使用滑动窗口时存在边缘效应，计算出的回归系数不能很好地代表曲线的斜率。如果原始数据不是等距的，可以使用线性插值将其转换为等间隔的更密集的数据。
曲线单调递增或递减和凹或凸时计算斜率变化点的公式略有不同:①曲线单调递增或递减时，用公式(2)和(3)，如1；(2)当单调递减且下凹时，使用公式(2)和(5)，如例2；(3)当单调递减凸性时，使用公式(6)和(3)，如在示例3中；④单调递增凸性时，使用公式(6)和(5)。
参考
[1]项，史久恩。非线性系统中数据处理的统计方法。北京:科学出版社，1997: 1 ~ 43。
[2]斋藤先生。用第三纪蠕变预测边坡破坏时间。第七届国际土力学和基础工程会议录，墨西哥，1969，Vo1.2:677~683
[3]达米卡·维克拉马辛哈和R . G .坎帕内拉。作为土壤可变性描述符的波动尺度。岩土工程中的概率方法，李、罗主编，鹿特丹巴尔凯马，1983:233~239
[4]让-皮埃尔·巴尔代。实验土力学。普伦蒂斯霍尔公司，1997:297~306
陈喜儒。变化点统计分析导论。数理统计与管理，第65438卷+0 ~ 4: 1 ~ 43。
[6]陈×R .简单变点模型中的推论.中国科学系列。一，1988，(6)
[7]丁。序贯累积和检验后变化点的置信下限。统计规划与推理学报，2003:115，311~326
[8]道格拉斯·霍金斯。数据的多变点模型拟合，计算统计与数据分析，37，2001:323~341
9韦恩·泰勒。变化点分析:检测变化的强有力的新工具，2000，/cpa/tech/changepoint.html
韦恩·泰勒。区分自回归数据和均值漂移数据的模式检验，2000年，/cpa/tech/changepoint.html