用三维曲面插值技术求解水轮机实时运行参数

结少鹏，虞 鸿

(浙江省水利水电勘测设计院，浙江 杭州 310002)

在水利工程规划设计和实际运行调度中一些关键参数或结果的获得常常需要查找特定的已知数据[1- 2]、参数表或曲线图，如实用堰淹没系数曲线[3]、水轮机模型综合特性曲线[4]等；另一方面，一些复杂方程的求解也需要求助于图解法，如矩形断面明渠底流消能水力计算、梯形断面渠道正常水深、临界水深等。虽然一些可以用复杂公式(如隐式公式)表达的问题目前已经能通过程序进行查表和插值[5- 6]，甚至可以通过试算、迭代等数值方法求解，简化了大量设计人员的工作，但还有一些问题因为无法找到对应的数学表达式，必须查找实验或经验图形曲线才能获得对应参数。具体来说，模型试验绘制的曲线，如水轮机转轮综合特性曲线、实用堰淹没系数曲线等，受到众多参数的影响，往往难以用简单函数或其组合进行良好的拟合，每次参数查找基本都需要由设计人员根据曲形进行取值，有较大的随意性，而且往往需要根据多条曲线上的邻近数据进行单向或双向线性插值计算，在种种限制条件下查询的过程不仅效率不高，获得的结果也有较大的随机误差，也难以用程序实现计算过程。

在当前水利数字化和信息化蓬勃发展的背景下[7- 12]，传统的参数曲线处理方法已经不能很好地适应程序化的处理方式和实时精细化控制对参数的要求。目前用于计算的图表通常表现为笛卡尔坐标系中的等值曲线簇，每条线上的各个点都可以认为是分别对应纵、横轴自变量的函数值，因此上述曲线簇可以认为是来自于一个复杂三维曲面的成果；如水电站工程中关键设备水轮机的运转综合特性曲线图，包含横坐标“功率”、纵坐标“水头”及等效率值曲线簇。对水轮机实时工况点效率值的计算，在数学上就构成三维曲面问题，可以通过克里格方法将类似曲面数据网格化，使之能用程序处理。

1 克里格(Kriging)方法简介

克里格法是比较灵活和准确的网格化插值方法之一，它可以通过在整体预测中给予较小的权重来补偿集群数据。

Z(p)为区域Ω上的一个随机过程，p∈Ω。Z(p)在p，p+h两点处的值之差的方差之半为Z(p)在p方向上的变差函数γ(h)。

(1)

Z(p)在p1，p2处的两个随机变量Z(p1)和Z(p2)的二阶混合中心距，即：

Cov{Z(p1)，Z(p2)}=E[Z(p1)·Z(p2)]-E[Z(p1]·E[Z(p2)]

(2)

Ω上有n个样本点，Zi=Z(pi)在pi处的值，则p0处的最优线性估计值为

(3)

(4)

求解方程组，可得

(5)

2 算例分析

2.1 工程概况

某水库工程以供水、防洪为主，兼顾灌溉、发电及改善河道水环境等综合利用。水库配套水电站型式为地面引水式，发电厂房位于大坝左岸下游侧。电站装设3台卧式混流式水轮发电机组，单机容量1200kW。水轮机运转综合特性曲线如图1所示。

工程正在建设水库智慧运行管理平台，基于设备在线监测系统及平台计算分析功能，精确掌握各类机电设备的运行状态数据，以便进行智能运行管理。为了满足电站精细化管理的需求，需要实时计算水轮机运行工况点效率，同时计算水轮机的过流量，并将工况点实时展示在运转综合特性曲线图中。由于机组容量较小，厂房空间有限，机组未装设流量和效率监测设备，因此水轮机流量和效率数据的采集是智慧管理平台开发的一个难点。

在电站的自动化监测系统中，可以实时读取的参数有：发电机功率Pg，蜗壳进水口压力PW，下游尾水位ZT等。

水轮机的过流量QT可通过式(6)计算：

(6)

(7)

式中，Pg—发电机功率，kW;ηg—发电机效率，可近似为常数。

H=PW+ZW-ZT

(8)

式中，PW—蜗壳进水口压力，m;ZW—蜗壳进水口高程，为常数，m；ZT—电站尾水位，m。

经以上初步分析可知，通过电站现有的监测系统参数，可以间接计算得到水轮机运转综合特性曲线图中工况点的横、纵坐标值(即：水轮机出力PT和净水头H)，只要能计算出水轮机的实时效率值，平台开发的数据问题就可迎刃而解。而水轮机效率值的计算必须采用程序化的处理方法，以满足流量数据的实时计算和实时显示。

2.2 等值线的处理及选点

水轮机运转综合特性曲线是由模型试验数值制成的。对于模型试验难以模拟的过渡过程工况，曲线图上有些区域是空白的。传统的人工计算方法下，水轮机效率值是根据出力PT和净水头H从运转综合特性曲线近似查得，遇到的也是稳态常规运行工况，故不需要对等效率线进行特殊处理。

而智慧管理平台在计算水轮机效率时，会涉及水轮机的开机、停机等过渡过程工况，如果在这个区域内没有曲线数据，则会显示错误。因此，需要对运转综合特性曲线图进行处理，用专业经验按照曲线簇的变化趋势，将外围的等效率曲线补全，保证水轮机运行范围内的效率数据是完整的，以便计算机能够正常计算处理。处理后的运转综合特性曲线图如图2所示。

图1 水轮机运转综合特性曲线

图2 处理后的水轮机运转综合特性曲线

为了将等效率曲线簇生成三维曲面，需要将等效率曲线的坐标点数据导出至Civil 3D软件中，因此需要先进行选点。每条等效率线上的点的选择应该能反映曲线的走势，并保证曲线的精度。在曲线的平直段，点的设置可以适当稀疏一些，而在曲率半径小的位置，点的设置应当加密。每条曲线上的点应将其Z坐标设置为效率值。水轮机运转综合特性曲线的点的设置如图3所示。

图3 等效率线的点设置

通过list命令将CAD中综合特性曲线上点的x(对应PT)、y(对应H)、z(对应ηT)坐标值导出提取至Excel表中，进行适当的坐标变化，形成(PT，H，ηT)数据表，见表1。

表1 从CAD中提取的综合特性曲线数据

2.3 三维曲面生成及栅格数据的输出

将表1中的特性曲线等值线数据导入至Civil 3D软件中，以便生成三维曲面。曲面生成过程中，选择自然邻近内插法或克里格方法对曲面进行光滑处理以提高曲面精度。由此，综合特性曲线图中离散的等值线簇转变为连续的三维曲面。如图4所示。

图4 克里格方法处理的三维曲面图(运转综合特性曲线)

三维曲面生成后利用指定的坐标生成点，可以直接取得曲面Z坐标的数值，就完成了一次数据查询，完全免去了插值计算的过程。

利用Civil 3D软件的按指定栅格处生成点命令，将离散成满足使用精度要求的网格化数据(如图5所示)，并按指定的格式输出“X-Y-Z”的数据文件供软件开发使用(实际文件是按地形数据“东距-北距-高程”等格式导出的)。原有的三维曲线数据点是粗线条的、不均匀的、不规律的，经过克里格方法转换以后，数据点变为精细的且均匀分布的网格化数据，便于计算机进行程序化的查询计算。Civil 3D软件导出的数据见表2。

图5 水轮机运行范围内的栅格设置

表2 Civil 3D导出的网格化数据表

将Civil 3D软件导出的特性曲线网格化数据上传至智慧管理平台，作为水轮机实时效率计算的基础数据。

2.4 水轮机实时参数计算

以表2中的一系列栅格数据为基础，根据水轮机的实时工况参数(出力PT(坐标值x)、实时净水头H(坐标值y))，计算水轮机实时工况点p与每个格栅点i的距离Li：

(9)

式中，Li—工况点p与格栅点i的距离；xp—实时工况点p的横坐标；yp—实时工况点p的纵坐标；xi—格栅点i的横坐标；yi—格栅点i的纵坐标。

通过式(9)的排序计算，找到离实时工况点p最近的4个格栅点A、B、C和D，然后选择最近4个格栅点的效率平均值，作为实时工况点p的效率值ηp。适当提高格栅矩阵的密度会提高计算精度。

ηp=(ηA+ηB+ηC+ηD)/4

(10)

得到水轮机实时效率值ηp后，即可按式(1)计算水轮机实时过流量QT。如此可得到水轮机的所有实时运行参数，满足工程运行管理的需求。

3 结语

(1)水利工程的建设中涉及较多的复杂曲线图，这些参数曲线难以找到对应的数学表达式或拟合函数。本文提出的程序化的插值计算方法，解决了传统人工读数和插值计算精度不高且效率低下的问题。

(2)在综合特性曲线处理时，补全等效率曲线簇需要专业经验，依赖执行人的判断力，精度难以控制。但考虑到水轮机在此区域内运行的时间很短，处理误差对整体的监控效果影响十分有限。

(3)现有的水轮机综合特性曲线都是根据试验数据绘制而成的。而实际上，模型试验时工况点的位置很难刚好处于等效率线上，等效率线也是由试验人员借助经验描绘，本身就存在一定的误差。若采用本文所阐述的方法，可以直接将所有试验工况点的原始数据导入Civil 3D软件，不需绘制等效率线，可以避免人为因素导致的随机误差。