多测度加权克里金法在高边坡变形稳定性分析中的应用

2009年6月水 利 学 报SH UI LI X UE BAO第40卷　第6期收稿日期:2008204229基金项目:国家自然科学基金项目(50579013);国家自然科学基金重点项目(50539110);江苏省研究生创计划项目(CX07B -137z )作者简介:刘志平(1982-),男,江西九江人,博士生,主要从事卫星导航定位与形变监测数据处理研究。

E 2mail :zhpliu @文章编号:055929350(2009)0620709207多测度加权克里金法在高边坡变形稳定性分析中的应用刘志平1,何秀凤1,张淑辉2(11河海大学卫星及空间信息应用研究所,江苏南京　210098;21武警工程学院信息技术教研室,陕西西安　710086)摘要:针对传统的克里金法仅以欧式距离测度描述空间结构相关性的缺点,基于最小方差准则,提出了综合曼氏、欧式和切氏三种距离测度的多测度加权克里金方法,同时给出了一种简便的替代最小方差准则的改进定权准则。

利用小湾水电站2号山梁高边坡监测点不同季节的变形速率数据进行实例计算,结果表明多测度加权克里金法提高了插值精度与可靠性,且改进定权准则的估计值接近理论最优值。

应用所提出方法得到的2号山梁高边坡2004年春夏秋冬多时相变形速率场,有效地揭示了高边坡自下而上牵引式变形逐步得到控制并趋于平稳的演化特征,为边坡变形稳定性评估与防治提供依据。

关键词:地质统计学;克里金法;高边坡;变形速率场;稳定性分析中图分类号:T U454文献标识码:A1　研究背景由于我国独特而复杂的地形地质条件和社会发展需求,高边坡稳定性问题作为大型水电工程建设中亟待解决的关键科学与工程技术难题[1]日益凸显,并成为工程地质、岩石力学和大地测量领域的热点课题。

多时相变形速率场可以有效地揭示高边坡稳定性及其演化规律,因而变形监测新技术的引入对边坡工程的设计与施工在安全稳定和经济合理的协调统一中具有重要作用[2-3]。

但是,边坡变形监测点位置的选择既受工程地质条件控制,又受工程施工情况约束。

此外,昂贵的监测设备也不允许无限地增加监测点,从而导致边坡工程区监测点分布欠合理且测点数目相对稀少。

因此,基于有限的监测点变形速率数据,利用空间插值方法估计工程区未知点的变形速率并结合规则格网技术实现可视化,是完整地反映边坡变形稳定性空间分布的重要手段。

空间插值方法主要有几何方法、统计方法、地质统计方法、函数方法、随机模拟方法、确定性模拟方法以及综合方法[4],其中以克里金法(K riging )为代表的地质统计方法提供了揭示空间变量结构相关性与随机性的理论[5-6],并己从创立之初的地质采矿学[7]逐渐渗透到农林[8]、环境[9]、气象[10]、测绘[11-13]等众多学科。

然而,上述应用中较少涉及克里金法的估计原理或准则。

运用克里金法的核心工作是建立变差函数模型[14],传统的克里金法仅采用欧几里德距离测度计算变差函数。

但是,不同的距离测度对理论变差函数模型的选择及模型参数拟合结果有重要影响。

此外,曼哈顿、欧几里德和切比雪夫距离测度主要分别用来检测特征平面中菱形、圆形及正方形结构的数据子集[15]。

鉴此,本文利用以上三种常用的距离测度检测数据子集结构的互补特性,提出基于最小方差准则的多测度加权克里金方法,同时给出一种简便的替代最小方差准则的改进定权准则,并将克里金法应用到小湾水电站2号山梁高边坡2004年春、夏、秋和冬4季节监测点变形速率的规则格网插值计算,以期利用该方法获取的春夏秋冬多时相变形速率场较真实地揭示边坡变形稳定性及其演化规律,说明规则格网化和所提出的多测度克里金法相结合对探索边坡变形机理与开展边坡治理工作有重要的应用价值。

2　多测度加权克里金法211　基本原理　克里金法由法国巴黎国立高等矿业学院教授G 1Matheron 和南非矿业工程师D.G.K rige 创立,其数学工具是变差函数。

假设位置x i 处的随机变量为Z i ,若其满足二阶平稳条件,则样本变差函数γ(d )的计算式为[14]γ(d )=12N d∑Ndi =1|Z i +d -Z i |2(1)式中:d 表示分离距离;N d 表示对应分离距离为d 时样本对的个数。

由式(1)可知,样本变差函数值与位置间的距离相关。

若对分离距离d 采用不同的距离测度,则最终拟合得到的理论变差函数模型也会有差异。

传统的克里金法中分离距离d 均是采用欧几里德单一距离测度进行计算。

虽然欧式距离测度在检测圆形结构数据子集具有明显优势,但对菱形和正方形结构数据子集检测效果却较差[15]。

因此,在样本数据子集结构不清楚或变化较大的情况下,分离距离d 均按欧式距离测度计算是不够准确和不切合实际的。

为充分利用不同距离测度在检测相应数据子集结构的优势以便更加真实地刻画空间结构相关关系,本文提出多测度加权克里金方法。

该方法利用最小方差准则综合曼氏、欧式和切氏三种常用距离测度下的普通克里金法估值,得出更高精度和更为可靠的结果。

为书写方便,将曼氏、欧式和切氏距离测度及其相关量分别用上标1、2、3区别,则多测度加权克里金法估计式可表达如下^Z 0=∑3k =1p k ・^Z k(2)式中:^Z 0表示位于x 0处的多测度加权克里金法估计值;p k表示距离测度为d k时的权重;^Z k0表示距离测度为d k时的普通克里金法估计值。

^Z k 0及其方差(^σk )2的表达式为^Z k0=∑ni =1λk i ・Z k i(^σk )2=∑ni =1λki・γk i 0+μk(3)式中:Z ki表示距离测度为d k时处于x k i的观测值;λki 、γki 0和μk分别表示距离测度为d k时的普通克里金法系数、样本变差函数以及拉格朗日乘数。

在式(2)～式(3)的基础上,根据误差传播定理可得多测度加权克里金法方差估计式^σ2=[p 1　p 2　p 3]D p 1p2p3(4)式中:D ∈R3×3表示单测度克里金法估计值方差阵,其元素由下式确定D kl =R kl ・^σk ^σl =r kl ・cor (d k ,d l )・^σk ^σl(5)式中:R kl 为测度相关矩阵R 的元素,表示各距离测度用于特征平面检测的相关系数;r kl 为最大测度相关矩阵r 的元素,表示各距离测度用于特征平面检测的最大相关系数;cor (・,・)为相关系数算子;d k表示距离测度为d k时检测的特征向量。

r 、d k 表达式如下:r =12π122π1π412π41d k =[(d k10)2　(d k20)2　…　(d k n 0)2]T (6)式中:d ki 0表示距离测度为d k时位置x ki 与x 0之间的距离,且满足d ki 0≤d k(i +1)0。

212　权的估计方法　由式(2)—(6)可以看出,多测度加权克里金法利用了3种单测度克里金估值进行加权计算,并通过引入测度相关矩阵考虑了不同距离测度下所得估计的相关性以合理评价计算结果。

此外,权重p k的估计体现了多测度加权克里金法本质和意义,以下利用最小方差准则推导p k估计方法。

考虑∑3k =1p k=1这一限制条件,在求多测度加权克里金方差的极小值时须引入拉格朗日乘数κ,按条件极值法构造函数f =(σk )2+κ∑3k =1p k-1(7)由式(2)可知,(σk)2是权重p k的显函数,为使多测度加权克里金法估计方差最小,只需将函数f 对权重p k求偏导数并令其为零,则有9f9pk =0　(k =1,2,3)(8)将式(8)与限制条件联立求解,即得最小方差准则下的权重,本文称之为理论值。

然而,上述方法得到的权重可能会出现负值情况。

虽可以借助现代全局优化搜索算法加以解决,但计算较为复杂,耗机时较多。

因此,本文从表达不同距离测度下的变差函数在描述空间结构相关性的物理意义出发,探讨相对简便的改进定权准则,即基于不同距离测度下估值方差的定权准则,其权重取值如下p k =(σk )-2∑3k =1(σk )-2(9) 式(9)中,将p k取为单测度克里金法方差(σk)2。

倒数的归一化值,隐式地表达了多测度变差函数在综合描述空间结构相关性的作用。

分析式(7)—(9)可得,多测度加权克里金法发挥了不同距离测度的结构检测能力,其估计方差依赖于测度相关矩阵R ,且上限为最大单测度克里金方差,较传统的单测度克里金法应更为真实可靠。

此外,基于最小方差准则的估计结果为理论最优,而计算相对简便的改进定权准则可获得近似理论最优值。

3　基于克里金法高边坡变形稳定性分析311　边坡稳定性的多时相变形速率场描述　边坡滑动的发生大多具有明显的前兆现象,而滑体的变形是最为直接、最易获取的信息,因此,在工程实践中可采用变形速率场刻画边坡稳定性状态。

若监测点在时间尺度T s 后三维变形为{Δx (T s ),Δy (T s ),Δz (T s )},则三维合变形为Δ(T s )=Δ2x (T s )+Δ2y (T s )+Δ2z (T s )(10) 边坡变形演化与季节有着密切的关系,由此本文以春、夏、秋和冬4季为时间周期T c ,以天为时间尺度T s ,将时间周期T c 内各时间尺度T s 的三维合变形取平均以表征各时间周期内的变形速率,即Δ(T s )=1T c ∑i ∈TcΔi(T s)(11) 由式(10)可知,变形后的测点坐标为{x+Δx (Tx),y+Δy(Ts),z+Δz(Ts)}。

但考虑到实际工程中变形值与坐标值相比一般为微小量,将变形后的测点坐标近似为{x,y,z}。

此外,结合式(11)可建立平面坐标下的测点在周期Tc内时间尺度为T s的变形速率的状态量表示为X={x,y, Δ(T s)}(12) 若边坡工程区全部监测点变形速率均由式(12)表达,则构成边坡变形稳定性的完整空间分布,综合春、夏、秋和冬季便得到揭露边坡变形稳定性的季节2多时相变形速率场。

312　小湾水电站高边坡变形速率克里金插值结果与分析　小湾水电站地处云南省南涧县与凤庆县交界的澜沧江中游,主坝高29210m,是当今世界上在建坝中的最高拱坝。

电站枢纽区山高坡陡,河谷呈“V”字型。

因此边坡变形稳定控制成为电站建设中的重大工程问题,其中左岸2号山梁饮水沟堆积体边坡长约70010m,平均宽19010m,前缘高程113010m,后缘高程159010m,平均坡度为30°～35°,是电站边坡监测的重点。

由于地表变形监测点可以大范围覆盖边坡,掌握边坡整体变形,并且可以在边坡开挖先期埋点监测,反映边坡变形过程,为此,在坡面上埋设了较多的地表变形监测点。

整个2号山梁高边坡先后布置了113个全站仪测点和16个G PS测点。

但受边坡开挖施工等诸多因素影响,许多监测点遭破坏而被废弃。