﻿ Delaunay三角剖分法在降水量插值中的应用
 气象学报  2012, Vol. 70 Issue (6): 1390-1400 PDF
http://dx.doi.org/10.11676/qxxb2012.117

0

#### 文章信息

XIONG Minquan. 2012.
Delaunay三角剖分法在降水量插值中的应用
Delaunay triangulated method with an application to the precipitation interpolation

Acta Meteorologica Sinica, 70(6): 1390-1400.
http://dx.doi.org/10.11676/qxxb2012.117

### 文章历史

Delaunay三角剖分法在降水量插值中的应用

1. 中国科学院大气物理研究所，北京，100029;
2. 中国科学院大学，北京,100049;
3. 国家气象中心，北京，100081

Delaunay triangulated method with an application to the precipitation interpolation
XIONG Minquan1,2,3
1. Institute of Atmospheric Physics, Chinese Acamdemy of Sciences, Beijing 100029, China;
2. University of China Academy of Sciences, Beijing 100049, China;
3. National Meteorological Center, Beijing 100081, China
Abstract: Delaunay triangulated method plays an important role in the spatial analysis. The characteristic of Delaunay triangulation and three kinds of generation algorithm are introduced in this article. The stochastic growth algorithm is also discussed in detail. Based on the principle of the triangulated slice linear interpolation the coefficient is attained and thus the value at an arbitrary point can be calculated through the triangulated slice linear interpolation algorithm. By using the daily precipitation data of the 2400 stations over China in 2008, the interpolation cross verification for the precision of interpolation is made between the triangulated slice linear interpolation algorithm and the inverse distance weighting for the whole China and the eight regions divided in China. The results show that the root mean square error of the triangulated slice linear interpolation algorithm is smaller for the various regions, and it has higher accuracy in the root mean square error and the mean absolute error in the concentrated site. For average error, triangulated slice linear interpolation algorithm has a highlighted high precision and the possible cause is discussed. The value of annual mean error is 0.005 mm and -0.107 mm for the triangulated slice linear interpolation and the inverse distance weighting method, respectively. It is also showed that the rationality of the triangulated slice linear interpolation algorithm. In the meanwhile, the three dimensional structure and the grid field on triangulated slice linear interpolating are showed in graph as well. The results showed that the figures calculated by the triangulated slice linear interpolation algorithm maintain consistency with the observation and have good visual effects. The precipitation distribution grid field by the triangulated slice interpolation gets rid of the inherent defect of the inverse distance weighting, thus improving the precision of interpolation. Finally, the application prospect of Delaunay triangulation in meteorology is discussed.
Key words: Delaunay triangulated networks     Triangulated slice linear interpolation algorithm     Inverse distance weighting     Precipitation
1 引 言

 图 1 降水量三维图(a.泰森多边形，b.Delaunay三角形) Fig. 1 Three-dimensional illustrations of precipitation(a. Thiessen polygon，b.Delaung triangulation)

Delaunay三角剖分法是根据一定规则，建立空间离散点的联系，构建正三角形，解决空间插值的难题。由于Delaunay三角网的良好特性，在地学、测绘、水利、城市规划、交通等多领域被普遍接受并广泛应用(Chris et al，2007;Hassan et al，2011; Sergio et al，2009;赵博等，2007)。Delaunay三角剖分法建立的观测站点间的几何拓扑结构则可以较好地解决邻近点选择问题，以此为基础进行的插值将能提高计算效率和获得更佳的插值精度；随着气象站点资料的网格化发展，Delaunay三角剖分法能发挥其应有的价值。

Delaunay三角剖分法具有严格的数学定义和完备的理论基础。它起源于Voronoi图(简称V图)，1850年Dirichlet讨论了V图概念； Voronoi于1908年在数学上定义了V图，V图主要用于分析空间点集或者其他几何对象和距离有关的问题；1911年，Thiessen在计算面雨量时，根据观测站点的降水量分布，讨论了V图构造及其计算方法；1934年，Delaunay在分析V图的基点间直线嵌入时，首次提出了Delaunay三角网络概念；1975年，Shamos和Hey讨论了使用计算机有效地构建V图，标志着计算几何的诞生，V图和Delaunay三角网是计算几何的重要基石。

Delaunay三角网是V图的对偶图，具有两个重要特性：(1)空外接圆性质：在点集P构建的Delaunay三角网中，每个三角形的外接圆中均不包含点集P中其他任意点，即任意4点不共圆。(2)最大化最小角性质：三角网中每个三角形的角度尽可能最大，即在两个相邻三角形构成的四边行中，当变换四边形的对角线时，6个内角的最小角不再增大，并且使三角形各边尽量接近等边。

 图 2 新增点Pr位于ΔPiPjPk内部 Fig. 2 Newly created point Pr within the triangulated slice ΔPiPjPk

S1、S2、S3的计算结果均大于0，表示点Pr位于ΔPiPjPk中；若S1、S2、S3至少有一个负数，那么PrΔPiPjPk之外；若S1、S2、S3有零值，那么PrΔPiPjPk某条边或顶点上。

3种生成算法各有优缺点，随机增长法较容易实现，占用内存小。本文将采用随机增长法进行Delaunay三角剖分(实现的代码见附录1)。3 Delaunay三角分片线性插值

 图 3 中国区域划分 Fig. 3 Division of China into the regions

 图 4 三角分片线性插值法和反距离权重法交叉检验结果(a. 全中国均方根误差对比，b. 全中国平均绝对误差对比，c. 全中国平均误差对比，d. 长江中下游地区均方根误差对比，e. 长江中下游地区平均绝对误差对比，f. 长江中下游地区平均误差对比) Fig. 4 Cross validations between the triangulated slice linear interpolation algorithm and the inverse distance weighting method(a. root mean square error for the whole nation，b. mean absolute error for the whole nation，c. mean error for the whole nation，d. root mean square error for the middle and lower reaches of the Yangtze River，e. mean absolute error for the middle and lower reaches of the Yangtze River，and f. mean error for the middle and lower reaches of the Yangtze River)

 图 5 2008年6月22日08时长江中下游地区24 h降水量(mm)图形显示(a. 降水量实况和观测站点分布，b. 降水量Delaunay三角网的三维结构，c. 三角分片线性插值后的降水量格点场，d. 反距离权重法插值后的降水量格点场；等值线起始10 mm，间隔20 mm) Fig. 5 24 h accumulative precipitation graphic display for the middle and lower reaches of the Yangtze River at 08：00 BT 22 June 2008(a. the precipitation observation and the station distribution，b. the three-dimensional precipitation field based on Delaunay triangulated networks，c. the grid field from the triangulated slice linear interpolation algorithm，and d. the grid field from the inverse distance weighting algorithm)

(1)三角分片线性插值是精确的估值方法，根据观测点进行三角剖分，再对待估格点进行三角分片线性插值，使用交叉检验方法分析插值精度，相对于反距离权重法的插值误差，三角分片线性插值的均方根误差有一定的优势；尤其在平均误差比较中，三角分片线性插值优势突出；相对于全中国范围，在站点密集地区(例如长江中下游地区)，三角分片线性插值法优势有增大的趋势；在降水量格点场图形个例对比中，三角分片线性插值方法消除了反距离权重法的“牛眼”现象，在格点场上更为精确地描述了降水量分布。

(2)Delaunay三角剖分能构建多面连续的降水面，可以立体地展现降水空间变化趋势，达到较好的视觉效果；由于是以观测点为三角形顶点，三角形的倾斜和标高都由3个数据点决定，任意点的估值都要受到该三角形表面的限制，因此，多面连续的Delaunay三角面是精确的插值面，能准确描述降水的三维结构。

(3)Delaunay三角剖分法为估值中样本点采集提供了有效、可靠的几何拓扑结构，也可作为各类气象要素插值的重要工具之一。气象数据有范围广、随时空变化强的特点，是多元散乱数据，当需要获得有序的格点数据时，插值结果的好坏依赖于观测点的分布。因此，合理的样本点采集方法是保证插值算法效率和精度的首要条件。根据Delaunay三角剖分的两个重要性质，剖分的结果是尽可能地接近正三角形。在和反距离权重法的比较分析中，验证了Delaunay三角剖分法在站点选择上的合理性。

(4)气象数据格点化应用中还存在估值精度和计算效率的平衡问题。一般认为，随着采样点数量的增大，估值准确率也会增大；但是样本点数量在增加的过程中，会出现两种情况，一是计算耗时，需要处理的数据量越多，计算量也必然会上升，在一些较复杂的插值算法中，其对计算机的存储空间、计算速度都会有更高的要求；二是样本点超过一定的数量，插值误差可能会增大。Delaunay三角剖分为解决此类问题提供了有价值的实验分析手段，即可以使用Delaunay三角剖分法进行邻近点搜索。

(5)面雨量的计算是Delaunay三角剖分的主要应用之一。目前，中国面雨量的计算缺乏客观、合理的技术方法，Delaunay三角剖分是具有严格的数学定义，完备的理论基础和成熟的算法。构建Delaunay三角网是自动完成的，无需人工干预或其他主观分析，并能充分利用已知站点数据，在客观分析基础上得到的面雨量将更具说服力。

(6)三角分片线性插值是无参数插值方法，由于Delaunay三角网具有唯一性，三角分片线性插值结果也是确定的。无需通过大量的试验获得合理的参数，三角分片线性插值是客观、高效的算法。

1. 在所有点的坐标中，寻找xy坐标绝对值的最大值，M=Max({|x1|，|x2|…|xn|}∪{|y1|，|y2|…|yn|})，确定P-1、P-2、P-3，建立ΔP-1P-2P-3，P-1=(0，3×M)，P-2=(3×M，0)，P-3=(-3×M，-3×M);ΔP-1P-2P-3将包含P中所有点。

2. 将点集P初始化为一个单独的三角形ΔP-1P-2P-3，初始网络用D表示。

3. 在点集P中随机排列，并对每个点分别赋予序号：P1、P2、P3…Pn

4. for i=1 to n

5. do(将点Pr插入到D中)

6. 在D中找到与Pr落点相对应的ΔPiPjPk

7. if(Pr位于ΔPiPjPk内部)

8. then 将ΔPiPjPk的3个顶点Pi、Pj、Pk都和Pr相连，生成新的3条边，ΔPiPjPk将划分成3个子三角形：ΔPiPjPr、ΔPrPjPk、ΔPiPrPk

9. 删除ΔPiPjPk的3条边，并将新生成的3个子三角形并入到D中。

10. EDGE(Pr，Pi，Pj，D)

11. EDGE(Pr，Pk，Pj，D)

12. EDGE(Pr，pi，Pk，D)

13. else(Pr恰好落在ΔPiPjPk的某一条边上，假设：Pr落在PiPj上)

14．连接PrPk，将PiPj对应的另一个三角形ΔPiPjPm上的点PmPr相连。ΔPiPjPkΔPiPjPm将被划分成4个三角形：ΔPrPjPk、ΔPrPjPi、ΔPrPjPm、ΔPiPrPm

15. 删除ΔPiPjPk、ΔPiPjPm，将ΔPrPjPk、ΔPrPjPi、ΔPrPjPm、ΔPiPrPm并入到D

16. EDGE(Pr，Pi，Pk，D)

17. EDGE(Pr，Pk，Pj，D)

18. EDGE(Pr，Pi，Pm，D)

19. EDGE(Pr，Pj，Pm，D)

20. 删除点P-1、P-2、P-3和与之相关联的边。

1. 搜索与Pr位置相对应的ΔPiPjPk

2. if(Pr在圆PiPjPk中)

3. then删除PiPj，连接PrPkPiPrPrPj，得到新的三角形ΔPiPrPkΔPjPrPk，并将新的三角形并入到D中。

4. EDGE(Pr，Pi，Pk，D)

5. EDGE(Pj，Pr，Pk，D)

 李庆祥,李伟. 2007. 近半个世纪中国区域历史气温网格数据集的建立. 气象学报,65(2):293-300 田付友,程明虎,张亚萍等. 2010. 校准雨量计密度对雷达联合雨量计估测流域平均面雨量的影响. 气象学报,68(5):717-730 吴宗敏. 2008. 散乱数据拟合的模型、方法和理论. 北京:科学出版社,1-40 杨贤为,何素兰. 1987. 江淮平原二类气象站网的设计. 气象学报, 45(1):104-110 赵博,孟小红,侯建全. 2007. 三维地震建模与可视化. 地球科学, 32(4):549-533 赵瑞霞,李伟,王玉彬等. 2007. 空间结构函数在北京地区气象观测站网设计中的应用. 应用气象学报,18(1):94-101 Bowyer A. 1981. Computing dirichlet tessellations. Computer J,24:162-166 Chris G, Maciej D. 2007. Dynamic cartography using Voronoi/delaunay methods. ISPRS workshop on updating geo-spatial database with imagery & the 5th ISPRS workshop on DMGISs, 41-47 de Berg M, van Kreveld M, Overmars M, et al. 2000. Computational Geometry Algorithms and Applications (Second Edition). Springer, 340pp Drozdov O A, Shepelevskij A A. 1946. The interpolation in a stochastic field of meteorological elements and its application to meteorological map and network rationalisation problem. Trudy Niu Gugms Series,1:13 Green P J,Sibson R. 1978 .Computing dirichlet tesselations in the Plane. Computer J, 21(2): 168-173 Hassan C, Majid H, Timothy P, et al. 2011. Delaunay triangulation as a new coverage measurement method in wireless sensor network. Sensors, 11(3):3162-3176 Hulme M. 1992. A 1951-80 global land precipitation climatology for the evaluation of general circulation models. Climate Dyn,7(1):57-72 Leonidas G, Jorge S. 1985. Primitives for the manipulation of general subdivisions and the computation of Voronoi diagrams. ACM Transactions on Graphics,4(2):74-123 Lewis B A, Robinson J S. 1978. Triangulation of planar regions with applications. Computer J, 21(4): 324-332 Ling T, Faisal H, George J H. 2010. Transfer of satellite rainfall uncertainty from gauged to ungauged regional and seasonal time scales. J Hydrometeor, 11(6):1263-1274 Lingas A. 1986. The greedy and delaunay triangulations are not bad in the average case. Information Processing Letters, 22(1): 25-31 Mccullagh M J, Ross C G T. 1980. Delaunay triangulation of a random data set for is arithmic mapping. Cartographic J, 17: 93-99 Miles R E. 1969. Solution to problem (probability distribution of a network of triangles). SIAM, 11:399-402 Schuurmans J M, Bierkens M F P, Pebesma E J R, et al. 2007. Automatic prediction of high-resolution daily rainfall fields for multiple extents: the potential of operational radar. J Hydrometeor,8:1204-1224 Shamos M I, Hoey D. 1975. Closest-point problems//Proceedings of the 16th annual symposium on the foundations of computer science,151-162 Sibson R. 1978. Locally equiangular triangulations. Computer J, 21 (3):243-245 Sergio O, Damon S, Laura R, et al. 2009. Three-dimensional ray tracing on delaunay-based reconstructed surface. Applied Optics, 48(20):3886-3893 Till H M V, Steve W L, Hoshin V G, et al. 2010. Multicriteria design of rain gauge networks for flash flood prediction in semiarid catchments with complex terrain. Water Resources Res, 46(11):554-670 Tsai V J D. 1993. Delaunay triangulations in TIN creation: an overview and a linear-time algorithm. Int J GIS, 7(6):501-524 Watson D F. 1981. Computing the n-dimension delaunay tessellation with application to Voronoi polytopes. Computer J,24:167-172