地理信息系统第三章

由于对地理空间中研究对象的表达方法多种多样，如栅格、矢量、DTM、等值线等，而且这些方法的数据表达形式差别较大，所以，对于一个功能齐全的GIS系统，就存在着数据的转换问题。下面将主要介绍这方面的内容。

3．1 坐标系的转换

地理信息系统的空间数据来自各种渠道，它们在输入到正式的 数据库之前，往往采用各自不同的坐标系。在同一个GIS系统中，多种坐标体系并存会给查询、分析带来不便，尤其是不同层的数据需要叠合，不同图幅之间需要拼接时就会有矛盾。因此，正式数据库用的坐标系越少越好，以方便用户使用。为此，需要坐标转换。

最简单的坐标转换是平移、缩放、旋转，其数学公式如下：

Xm(xcosysin)h （3-1）

Ym(xsinycos)k式（3-1）中x,y为转换前的坐标，X,Y为转换后的坐标，为逆时针方向的旋转角度，m分别为X,Y方向的缩放比例，h,k分别为旋转后在X,Y方向的平移距离。

比较复杂的是投影变换。地球是一个近似的椭球体，当用二维的平面坐标系来反映三维的地球表面时，就要用特殊的投影方法，以实现三维到二维的转换。显然，相关的投影方法会不可避免地带来计算结果的误差。这些误差在绘制小范围、大比例尺（如1：500，1：1000）地图时往往忽略不记，而对大范围、小比例尺的地图，在地图的某些部位（一般在边缘）失真就会明显，这样，对距离、面积、形状的量算、分析都可能与实际有较大的误差。为了照顾某些范围的变形、失真，往往需要变换投影，即从一种投影方式转成另一种投影方式，其过程是一系列复杂的数学运算。一般情况下，我国对大比例尺～1：50万的地形图均采用高斯—克吕格投影，而1：100万采用Lambert投影。

3．2 TIN模型的自动生成

61

无论是地理信息系统，还是计算机辅助设计系统，不规则三角网的自动联结(又称三角剖分)占有十分重要的位置。在所有的三角网自动联结算法中，主要遵循以下两条准则： 1.在形成的三角形网格中，扩展点与扩展边两点连线的夹角最大（一般通过余弦定理进行判断）；

2.所形成的三角形的全部边长之和最小。 第一种准则使用最为广泛，几乎所有的文献和商业化软件系统都以它为准以进行三角形的扩展或剖分。因为它所形成的三角形尽可能地接近等边三角形，这将有利于TIN的应用。采用第一种准则在数学上是有其理论依据的，即Delaunay三角剖分。它被国际上确认为“最佳三角剖分方案”。

Delaunay三角剖分的基本概念和性质如下：

设二维平面P上的随机点集为VV1,V2,...,VN。点集V中数据点满足以下条件： 1.N3;

2.数据点面状分布，不都在同一条直线上； 3.没有四或四个以上的点共圆。

令dVi,Vj表示点Vi到点Vj间的欧氏距离，Vi表示平面上距离点Vi较之距离其它任何点更近的点的轨迹，即 ViXP|dX,VidX,Vjj1,2,...,i1,i1,...N

则Vi是一个区域，称为点Vi的Voronoi多边形(此区域也可认为是点Vi的影响区域)。它实际上是点Vi与其它N1个点的连线的垂直平分线所形成的N1个半平面的交，即 ViHVi,Vj

ijN个点的Voronoi多边形构成Voronoi图，如图3-1的虚线所示。构成Voronoi多边形的边称为Voronoi边，它是由两点连线的垂直平分线构成的。如果将对应于每一条Voronoi边的两点相连，则构成一系列的三角形，形成随机点集所构成凸包的一个三角剖分，称为Delaunay三角剖分。Voronoi相邻多边形的公共点称为Voronoi点，它是在没有四个及四个以上点共圆的情况下，每一个三角形外接圆的圆心。此圆心也是每三条Voronoi边的交点，此时，每一Voronoi点被三个多边形所共有。在退化的情况下，共用一个Voronoi点的多边形多余三个，但仍为多边形外接圆的圆心。

62

在进行Delaunay三角剖分时，如果有四个或四个以上的随机点共圆，那么联网过程要作特殊的处理。如扩展边有多个点满足最大夹角时，可按顺时针或逆时针方向直接顺次相连。

图3-1 Delaunay三角剖分示意图

到目前为此，建立TIN的Delaunay三角剖分算法已十分成熟。这些算法的一个最大特点就是把离散随机点所构成的制图区域认为是凸形区域，并且一般不考虑制图区域内的空洞或不连续的突变区域，即使有些文献提出了处理凹形区域的算法，但仍先假设制图区域为凸形的，待联网结束后去掉那些三角形三点都为边界点的三角形。通过研究发现，此算法是具有局限性的，它可能去掉那些合法的三角形(图3-2中的三角形ABC)。

边界CBA

图3-2 凹形区域生成TIN示意图

3．2．1构建TIN的定义和遵循的规则

目前，构建TIN的方法主要有如下三种：三角网生长法，逐点插

63

入法，分割—归并法。在本节中，主要介绍应用最广泛的生长法。

1.遵循Delaunay三角剖分的定义及有关规定 2.几个新的规定 规定 1

制图区域由边界相互封闭的Nii0,1,2,...,n1个子区域组成。 规定 2

第Ni个子区域包含m个数据点，那么，第Ni个子区域三角联网时的数据扩展点只能从Ni边界区域内m个点中寻找，即 PijNi,,,...,n1;j012,,,...,m1 i012规定3：

在生成TIN时，扩展三角形的两边（不包括扩展边）与联网区域内的包括所有诸如山脊线、边界线等特征线不能相交。这样，可以保证联网结果的合理性。 4.正负区判断准则 对于某一直线段AB，它的判别方程如下（式3-2,3-3,图3-3（a）)：

yAxBFx,yx1xx2x1x2x1 (3-2)

Ay2y1/x2x1By2x2y1x1/x2x1 判别点C所处正负区的数学表达式为：

0Fcx,y00C点位于正区C点位于直线上 (3-3) C点位于负区BCCBA A

ED

(a)

(b)

图3-3 三角形扩展准则示意图

64

5.扩展点优先准则

在Voronoi图退化的情况下，三角形外接圆上的数据点数>3。此时，数个点与扩展边连线的夹角相等。从数据点向扩展边作垂线，取其垂足离扩展边中点最近的点为扩展点（如图3-3(b)中的A点）。 6.封闭区域边界点正负区判别值的计算

第Ni个子区域边界点k或k,k1边的正负区判别值等于区域内某点P相对于边界点所在边k,k1的正负区判别值(假设图3-4中的边界点D为k点,边界数据点的顺序为C、D、E、F、G、H、I、C,那么DE边即为k,k1边)。

由于P点在边界内,所以pb_value(i,k)0.0 。对任一边界线段DE,P点的确定方法如下:计算DE的中点坐标并作垂线MN,在区域内的MN线段上任取一点即可。

D(x1,y1)MPCINGFE(x2,y2)H

图3-4 计算封闭区域边界点正负区判别值示意图

3．2．2基于制图边界拓扑结构的正负区判别模型

以Delaunay剖分准则为基础的现有三角剖分算法不适合构建复杂对象TIN的主要原因是它们考虑的数据点属性较为单一，数据在平面上的构型较为简单，边界形状较为规则（以凸形为主）。 通过对原始数据及各类制图边界实体相关关系的研究发现，无论制图区域的数据点构型、属性、边界形状多么复杂，只要先建立原始点、线、面间的各种拓扑关系，再利用与制图原始边界有关的正负区判别准则，就可以彻底地解决复杂对象TIN建立过程中遇到的各种问题。这就是下面将要介绍的基于制图边界拓扑结构的正负区判别剖分模型。

剖分模型的原理及算法主要有以下几点：

65

图3-5 复杂地质对象示意图

1.图3-5代表了建立地质对象TIN时可能遇到的所有情况，其中包括正断层(F9)、逆断层(F8)、正断层与正断层的切割(F5与F6)，正断层与逆断层的切割(F1与F2)，逆断层与逆断层切割(F3与F4)，冲刷带，风氧化带以及不同封闭块段之间的紧密相邻和松散相邻。从图中可以看出，在遵循第一小节有关准则及规定的情况下，可以根据地理信息系统拓扑结构的有关概念和算法，确定点状、线状(断层等)，面状(封闭边界等)图形实体的间的拓扑关系。 (1)点的拓扑关系

这里关于点的拓扑关系其实只是指点与封闭区域间的包含关系,即给出点所在封闭区域的序号。建立点的拓扑关系，对于多个相互独立的封闭制图区域的TIN的建立是十分必要的。在处理某一子封闭区域时，只有本区域内的原始数据点才参与三角形扩展点的判别与计算。这样将大大减少计算机的运算量，更为重要的是它将保证扩展点不位于另一封闭子区域。边界数据点太稀疏时，此种情况肯定会出现，见图3-6(a)，扩展AB边时，寻找到的扩展点是D而非C，显然这是错误的。如果出现此种情况，将造成联网过程的混乱，最终无法实现制图区域的网格化。

(2)线状实体间的拓扑关系

线状图形实体间的拓扑关系通过线本身的编号，首尾结点相邻线段的编号以及线段所在的封闭区域的编号确定。

66

1)线段的编号由软件系统根据收索到线段的先后顺序自动确定； 2)建立拓扑结构的图形实体必须为一独立的线段(弧段一般满足要求)。虽然同一属性的线段可能被分为相邻（可能相交，也可能不相交）的几段，但它们的拓扑关系需各自单独建立。图3-6(b)的逆断层F1的下盘被分为两段B--C、E--F，B--C段与其它段的拓扑关系为(2,1,3,0)。2为B--C本身的编号，1为相邻段A--B的编号，3为相邻段C--D的编号，0为封闭圈G的编号。

ADB封闭圈GB断层F1ACCEFD断层F2 (a) (b)

图3-6 点状、线状图形实体拓扑关系处理示意图

(3)面状区域的拓扑关系

面状区域的拓扑关系通过以下三个参数加以确定： 1)封闭区域的编号；

2)相邻内外封闭区域的编号； 3)封闭区域所处的层次编号。

封闭区域的编号即为软件系统在图形数据库中收索到的封闭区域的顺序号i(i=0,1,2...,N)，这些序号以封闭区域从外向里的顺序逐渐增大。规定最外层区域以外的区域编号为-1。在图3-7中，有两个相互独立的制图子区域A和B，A，B内又分别包括C、D、E和F、G、H、I。其中C内又包括J，K。A、B、C、D、E、F、G、H、I、J、K各区域的编号分别为0,1,2,3,4,5,6,7,8,9,10。封闭区域的层次编号与搜索同一级别的封闭区域的顺序有关,以图3-7说明即为:最外层封闭区域A、B为第一层次,层次编号为0,第一层次内的C、D、E、F、G、H、I为第二层次,编号为1,第二层次内的J、K为第三层次,编号为2。首先收索第一层次的封闭区域，然后向内收索第二层次的封闭区域，以次类推，直至无封闭区域为止。

(0，0，-1，2，3，4)表示了A封闭区域与其它子区域间的拓扑

67

关系，其中：第一个0表示层次编号,第二个0表示本身的编号，-1表示它的外层无封闭区域，2，3，4分别表示它包含的子封闭区域的编号；同样，区域C与其它区域的拓扑关系可以通过(1，2，0，9，10)反应出来。

ACDE

FJKHBGI

图3-7 面状区域拓扑关系示意图

图形实体间拓扑关系的建立将有助于正确计算封闭区域边界点的正负区判别值。

(1)计算完所有封闭区域边界点的正负区判别值后，需把层次编号为奇数的封闭区域非逆断层边界点正负区判别值取反。

(2)在建立面状区域的拓扑关系时，各封闭区域的层次编号在此具有十分重要的意义。从图3-7可以看出，对于任意独立的封闭区域，联网数据位于层次编号为偶数和奇数的边界区域之间； 2.边界正负区判别模型的建立

所谓边界正负区判别模型就是根据封闭区域间的拓扑关系并结合二维平面点的正负区判别值计算方法，建立封闭区域内联网数据点与边界间的合法关系。

判别模型的具体描述如下。

如果封闭区域i内的扩展三角形中有边界数据点,那么,使扩展三角形成为合法三角形需满足的条件为:

(1)假设边界点为k,共用k点的两边分别为k1,k和k,k1,相应的正负区判别值为pb_value(i,k1),pb_value(i,k),三角形的另两点为m、n,点m、n相对于k1,k和k,k1边的判别值为pb_value(i,m1),pb_value(i,m2),pb_value(i,n1),pb_value(i,n2),那么, 1)如果k1,k和k,k1两线段在区域内一侧的夹角180.0(图3-8中的角),下面四个不等式就必须同时得到满足:

pb_value(i,m1)*pb_value(i,k1)0.0; (3-4) pb_value(i,m2)*pb_value(i,k)0.0; (3-5)

68

pb_value(i,n1)*pb_value(i,k1)0.0; (3-6) pb_value(i,n2)*pb_value(i,k)0.0; (3-7) 2)如果k1,k和k,k1两线段间在区域内一侧的夹角>180.0(图3-8中的角),(3-4),(3-5)和(3-6),(3-7)两组不等式中各有一个不等式满足条件即可。

(2)如果扩展三角形中边界数据点数大于1,需分别利用上面介绍的算法对三角形顶点进行交叉处理。

正负区判别模型的建立，可以有效且快速的地解决以下问题： 在同一封闭制图区域内寻找三角形的扩展点时，如果扩展点为边界点，那么可以利用此边界点的正负区判别值确定扩展点的合法性。结合图3-8说明即为：

1)假设当前扩展边为AB，利用Delaunay三角剖分准则寻找到的扩展点为C，那么线段EC、CD的正负区判别值分别为pb_value(E)、pb_value(C)。显然,A、B点与EC、CD线段的相对关系满足(3-4)、(3-5)、(3-6)、(3-7)不等式的要求，所以C点合法。这一算法是解决断层问题特别是逆断层的关键算法。逆断层上下盘在空间上有相重的区域，上下盘的数据控制点分别位于对盘的区域，利用现有的三角剖分算法将无法解决这一问题。但利用本文介绍的算法，此一问题将变得较为简单。在图3-8中，按Delaunay三角剖分算法，扩

F2封闭区域PDBCAHIJKαβE

图3-8 正负区判别模型示意图

展边HI的扩展点肯定为J，但因J点为断层F2下盘的控制点，而非当前盘（上盘）的控制点，所以是非法的。根据H、I、J三点的有关正负区判别值就可确定点J的非法。虽然点K在平面上位于封闭区域P内，但它与H、I的关系满足正负区判别模型的算法要求，所以K点才是H--I边的扩展点。

69

建立正负区判别模型是解决逆断层或同一平面坐标有多个高程值联网问题的主要原因之一。

2)解决了任意边界形态的联网问题。边界类型既可以是凹型也可以是凸型，甚至是相邻边界相嵌（多个独立的封闭圈边界）。

3．2．3 三角网的形成过程

1．首先在制图边界寻找相邻的数据点，或者制图区域内最近的两点，把此两点的连线作为扩展边（假设两点分别为第一点、第二点）；

2．利用上面介绍的相关算法得到当前扩展边的扩展点，即三角形的第三点。这样，就形成了第一个三角形；

3．利用第一个三角形的三边（注意，如果从边界开始，就只处理第一点与第三点，第三点与第二点的连线边），使之作为三角形扩展边，得到新的三角形。重复相关步骤，形成新的三角形；

4．结束TIN自动生成的条件是无新的三角形生成。

*如果制图区域是一个连续的区域，不存在任何突变信息（如街道、河流、断层），那么，只需利用角度最大准则、正负区判别准则就可解决问题。 3．2．4 实例

利用上面介绍的方法处理了若干不同复杂程度的制图区域,取得了满意的效果。见图3-9，3-10。

图3-9 构建复杂对象TIN模型示意图1

70

图3-10 构建复杂对象TIN模型示意图2 3．2．5 逐点插入法和局部优化算法

下面再简单介绍一下逐点插入法。

算法的基本思路是：先在包含所有数据点的一个多边形中建立 初始三角网，然后将余下的点逐一插入，用LOP算法确保其成为Delaunay三角网。本算法的基本步骤为（如图3-11）： (1)定义一个包含所有数据点的初始多边形；

(2)在初始多边形中建立初始三角网，然后迭代以下步骤，直至所有数据点被处理；

图3-11 逐点插入法示意图

(3)插入一个数据点P，在三角网中找出包含P点的三角形，把P

71

点与三角形的三个顶点相连，生成新的三角形；

(4)用LOP算法优化三角网。Lawson提出的局部优化过程（Local Optimization Procedure）算法是为了生成Delaunay三角网。从理论上说，不论用何种方法生成三角网，只要用LOP算法进行处理，就能使它变成Delaunay三角网。算法的基本含义为：对由两个公共边组成的四边形进行判断，如果其中一个三角形的外接圆包含第四个顶点，则将这个四边形的对角线交换。如图3-12所示。

局部优化算法的含义是：

优化前 优化后 图3-12 LOP算法示意图

3．3 格网或TIN模型的插值

在GIS应用过程中，原始数据点常呈离散分布状态。有时，虽然原始数据点已呈格网或三角网的形式分布，但仍需加密这些数据，故同样存在着数据的内插问题。格网插值的含义如下：设已知一组数据，它们呈离散或规则分布状态（图3-13），现要从这些数据中找出一个函数关系式，使该关系式最好地逼近这些已知的数据，并根据该函数关系式求出区域范围内的任意格网交点或待插值点上的值。根据内插精度的不同，可以将内插方法分为精确的和粗略的两大类，对于每一类，又可根据原始数据点的分布形式分为离散数据和格网数据的插值两小类(图3-14)。鉴于空间数据的特点，这里主要介绍精确类的插值算法。

72

图3-13 数字高程模型插值示意图

图 3-14 插值方法分类图

考虑到数学原理的复杂性，这里只介绍几种较为简单的算法。 3．3．1 利用离散点数据的插值

73

1． 按距离加权平均法 (1) 插值的原理和算法

假设离插值点最近有n个数据点，那么，这n个数据点对插值点 值的影响与距离有关，距离越远，影响越小（图3-15）。

在图3-15中，任意数据点(xi,yi,zi)离插值点A的平面投影距离为di，那么

di[(xixA)(yiyA)] （3-8）

2122在式（3-8）中，i为数据点数。

求出离A点最近的n个数据点的距离di（i1,2,...,n）以求出A点的内插值zA，即：

zi)d1i zAi （3-9） n1()i1di(n(2) 算法的主要步骤

1)输入原始数据点值； 2)确定插值点平面坐标；

3)寻找与插值点最近的n个原始数据点，同时求出它们与插值点的距离；

4)根据式（3-9）计算插值点值；

5)重复2)、3)、4)各个步骤，求出所有插值点的值。

图3-15 按距离加权平均法插值原理示意图

74

2． 按方位取点加权法

按距离加权平均法只要求求得n个最近的数据点，而不考虑数 据点所在的位置，这样就可能造成数据点集中在插值点的某一侧，从而增大插值结果的误差（图3-16）。按方位取点加权法就可克服这一方面的缺点，使插值方法的使用范围更加广泛。

(1) 插值的原理和算法

对A点进行插值时（图3-17），先以A点为中心把平面分成4 个象限，再把每个象限等分成n0份，其结果就把全平面等分成4n0个区域。在每一个区域内寻找一离A最近的点，假设该点的坐标为(xi,yi,zi)，它到A点的距离为di，这样，A点的插值表达式为：

zAcizii14n0d cij1ji4n04n0k1l1lk4n02j （3-10）

2ld22d2j(xjxA)(yjyA)dl2(xlxA)2(ylyA)2式（3-10）中ci之和满足：

cii14n0di1j1ji4n04n02jdk1l1lk4n04n01

2l如果A点坐标刚好与某一原始数据点B的坐标相等，那么A点 的插值就等于B点的值。其原因如下：因A点与B点相重，所以rB0，这样就使除了B点以外的4n01个数据点的ci值都等于零，cB1，即：

75

zAcizi1zB00...0zB

i14n0 图3-16 原始数据点集中一侧示意图 图3-17 按方位取点加权插值原始示意图

(2) 算法的主要步骤 1)输入原始数据点； 2)输入n0；

3)确定插值点平面坐标；

4)在4n0个区域中分别寻找离插值点最近的点，并求出它们之 间的距离；

5)利用（3-10）式计算插值点的值；

6)重复3)、4)、5)各步骤，求出所有的插值点的值。

3．3．2 利用TIN或格网数据的插值

有时，虽然原始数据点以格网或三角网的形式排列，但仍需继续在网格内内插数据。在这种情况下，通常以网格小块为插值区，采用低次项函数拟合已知值。这里只讨论线性内插，双线性多项式内插两种方法。

1． 线性内插

此方法实用于三角网网格内的插值（3-18）。假设ABCD为一平

面，三顶点（x,y,z）坐标已知，现求A点的插值zA。插值函数为： zAa0a1xa2y （3-11） 把B、C、D三点的坐标代入（3-11）式，联立就可以解出a0,a1,a2三个系数，从而求出A点的内插值。

76

图 3-18 线性内插原理示意图 图3-19 双线性内插原理示意图 2． 双线性多项式内插

此方法一般在格网插值时使用（图3-19）。假设插值点A的值在与x轴或y轴平行的方向上分别与坐标y和x成线性比例关系。其表达式为：

 zAa0a1xa2ya3xy （3-12） 利用格网的基本单位长度a,b以及矩形4顶点E、F、C、D的坐标

值就可求得zA。

(1)假设E为坐标原点，EF=a,EC=b。过A点作EC或FD的平行线，它们与EF、CD边分别交于G、H。

(3) 利用E、F、C、D四点的坐标线性内插出zG,zH：

zGzE(zFzE)xaxzHzC(zDzC)ayb （3-13）

(4) 利用zG,zH的值再次进行线性内插，就可得到zA： zAzG(zHzG) （3-14） 联立（3-13）、（3-14）两式，可得：

zAzE(1)(1)zF(1)zC(1)zD （3-15）

如果以格网的行列号为基准，那么，（3-15）式可以变为： zAz(i,j)(1)(1)z(i1,j)(1)z(i,j1)(1)z(i1,j1)

3．4 等值线图的自动生成

77

xaybyxbaxyabxaybxaybyxbaxyabxayb 在三维空间，虽然获取地理研究对象特征点的手段和方法多种 多样，如测绘、钻探、物探、野外调查、航空航天遥感，但获取的数据量仍十分有限，数据采集点在空间上的分布呈离散形式。针对这些特点，GIS工作者的任务就是利用有限的离散空间信息尽量恢复地学变量的分布状态，以研究它们在空间上的变化规律和分布特征，最常用、最普遍的方法就是绘制空间变量的等值线图。

任一空间变量数据采集点都是以三维坐标(x,y,z)表示的，其中(x,y)标识采集点的平面位置，z表示诸如环境污染参数、地层厚度、标高等测量值。等值线图就是根据(x,y,z)值所作的一种数字正射投影。到目前为此，自动绘制等值线图的方法归纳起来不外乎有两种，即三角网方法和矩形网方法。这两种方法代表的算法和过程差别较大，各自均具有本身的特点。

3．4．1 格网方法

根据等值线追踪算法的不同又可以分为两种方法。一种方法是直接在格网边上作线性内插得到等值点，然后按一定的法则追踪出一条等值线的全部离散点，再利用曲线光滑算法完成等值线的绘制。这种算法的优缺点归纳起来有如下几条：(1)程序设计简单，占用计算机的内存随网格数的增多而增多。就目前计算机软硬件的发展水平来看，它们一般已能完全满足算法对内存的要求；(2)因格网尺寸的差异，有时会造成同一幅图的输出结果出现微小的差别，其原因是如果网格尺寸过大，会忽略一些微小的变化区域；(3)如果格网尺

(a) (b)

图20 格网尺寸大小不同对绘图结果的影响

寸过大，曲线光滑时可能发生曲线相交现象（图3-20(b)）。绘制等

78

值线的另一种方法是利用已有的网格点数据拟合一个光滑曲面，该曲面有一个统一的数学表达式，并且已知网格点都在曲面上。根据函数关系，就可以追踪出当前要生成的等值线，最终绘出等值线图。由于空间变量的随机性和复杂性，实际研究的空间变量的坐标值很难准确满足某一数学函数表达式，故编程时程序较为复杂，曲面一般不能很好地拟合已知空间变量。但该方法也有其自身的优点，即曲线追踪容易、很光滑，且不易出现曲线的相交现象。考虑到两种方法的优缺点，一般在GIS中选取前一种方法，本节也将主要介绍这方面的内容。

纵观格网方法绘制等值线的过程，可以分为5个方面进行论述：(1)内插当前等值点；(2)追踪等值点；(3)开、闭曲线线头、线尾的寻找；(4)其它一些算法；(5)曲线的光滑。在本节中将不介绍曲线光滑算法，相关内容在后面的内容中介绍。

为了便于讨论问题，这里先对格网和有关数组作一说明。

(1)设绘图区域由M*N个网格点组成，x,y方向的分割分别是i1,2,...,M;j1,2,...,N。这样，网格区域内就有(N1)*M条纵边和(M1)*N条横边；矩形区域左小角的坐标为（0，0），x和y方向网格基本单位长度分别为a和b（图3-21）。据此可以得到任一纵横边的坐标值：

xy(i,j,1)xi,j(i1)*axy(i,j,2)yi,j(j1)*bi1,2,...,Mj1,2,...,N

图3-21 与格网有关的标识示意图 图3-22 内插等值点示意图

(2)对于任一网格，以该网格左小角点在x,y方向的分割序号值加

79

以标识；任一网格边，以该边左端点或下部端点在x,y方向的分割序号值表示；网格交点的z值用H(i,j)表示。以图3-21说明即为：

，AB边的标识为（2，3），BDM7,N6，阴影网格的标识为（2，3）

边的标识为（3，3），网格交点F的值表示为H(4,3)100.5。

(3)假设SS(i,j)、HH(i,j)分别标识纵、横边上是否具有当前等值点（i1,2,...,M;j1,2,...,N）；当前追踪等值线的值为HZ。

(5) 对于任一矩形网格，假定网格边的z值在两个相邻网格交点 间是线性变化的 1． 内插当前等值点

图3-22是任一平面投影格网在三维空间上的表现形式。其中矩 形ABCD是空间四边形A'B'C'D'在平面上的投影，A'B'边上的E'点在平面上AB边的投影为E。现求AB边上是否具有HZ值的点。从图中可以看出，只要矩形ABCD的A点坐标确定后，其它3个点的坐标也就确定了。现假设A,B两点的z值分别为H(i,j),H(i1,j)，这样，根据下式就可以确定AB边是否具有当前值。 rHZH(i,j)HZH(i1,j) 或 r

H(i1,j)H(i,j)H(i,j)H(i1,j)如果0r1，就证明AB边上包含有HZ值，即有当前等值线上的

等值点存在。

假设E点为当前等值点，那么

xExy(i,j,1)[(xy(i1,j,1)xy(i,j,1)]*r 或 xExy(i,1j,1)[(xy(i,j,1)xy(i1,j,1)]*r yExy(i,j,2)[(xy(i1,j,2)xy(i,j,2)]*r 或 yExy(i,1j,2)[(xy(i,j,2)xy(i1,j,2)]*r AB代表网格区域的横边。在纵边上求等值点的算法与上述相似，只是所利用的HH(i,j)值和xy数组不同罢了。 为了有利于等值点的自动追踪，需事先计算出网格区域内全部纵横边上的当前等值点，并作如下处理：(1)如果0r1，则令HH(i,j)或SS(i,j)r；(2)如果1r或r0，则令HH(i,j)或SS(i,j)2。

等值线自动追踪时，如果某一边的HH(i,j)或SS(i,j)值大于1，证明此边无当前等值点存在。

80

2． 追踪等值点

等值线的追踪是自动等值线的关键。正确的追踪方法才能够保证追踪出的等值线与实际情况最大程度地保持一致，才能保证等值线的合理联结和不相交。为了达到这些目的，就必须正确地选择等值线的追踪方向和对若干相邻等值点的连接方式。 (1)确定等值线的追踪方向

由于等值线是由一系列位于矩形网格纵横边上的等值点序列组成 的，而这些等值点又是由(x,y)的坐标表达的，所以，在研究等值线的追踪方向时，需从已追踪出的等值点的(x,y)坐标值着手，研究这些点坐标与追踪方向的相关关系，从而找出计算机自动确定追踪方向的算法。由于等值点位于网格边上，而任一网格边又被相邻的两个矩形网格所共有，每一网格周围有四个网格，所以等值线通过相邻网格的走向只有四种可能，即自下向上、自左向右、自上向下、自右向左。

1) 自下向上追踪

见图3-23(a)。假设A、B是当前追踪等值线序列中的最新两点， B在A点之后，它们所在边的标识分别为(i1,j1),(i2,j2)(后面内容中A、B的含义与此相同)。从图中可以看出，A点可以位于Ⅰ号矩形的(i,j),(i1,j)纵边和(i,j)横边，而B点只能位于(i,j1)横边，A点y值与y方向网格基本单位b相除所得商的整数值小于B点y值与b相除所得商的整数值，即j1j2。如果A、B两点所在边满足这一条件，就说明当前等值线的追踪方向是由下向上的。

2) 自左向右追踪

见图3-23(b)。在此种情况下，A点位于Ⅰ号矩形网格的

(i,j),(i,j1)横边和(i,j)纵边，B点位于(i1,j)纵边上。比较A、B两点的横坐标我们不难发现，A点横坐标与x方向网格基本单位长度a相除所得商的取整值一定小于B点横坐标与a值相除所得商的取整值，既int(xA/a)int(xB/a)或i1i2。所以，判断当前追踪等值线是否由左向右追踪的条件是i1i2。

3) 自上向下追踪

见图3-23(c)。A点位于Ⅰ号矩形网格的(i,j),(i1,j)纵边和 (i,j1)横边，而B点则位于(i,j)横边。比较A、B两点的x,y坐标发现，它们不满足上面论述两种情况的判断条件。实际追踪等值线时，先进行自下向上、从左向右两种情况的判断，如不满足各自的条件，

81

(a)自下向上追踪 (b)自左向右追踪

(c)自上向下追踪 (d)自右向左追踪 图3-23 追踪等值线四种情况示意图

那么追踪方向就属于自上向下和由右向左两种情况。区别自上向下、从右向左两个追踪方向的条件是很简单的，只要B点x坐标值满足如下条件，即int(xB/a)*axB,就证明当前的追踪方向是自下向上的。

4) 自右向左追踪

由于等值线的追踪方向有四种情况所以排除上面已讲的三个方向后，就必然属于自右向左方向。见图3-23(d)，A点位于Ⅰ号矩形网格的(i,j),(i,j1)横边和(i1,j)纵边，B点位于(i,j)纵边，B点的y坐标满足:int(yB/b)*byB。

尽管等值线的追踪只有四个方向，但如果对各方向追踪的先后顺序不同，必然也造成判断条件的差异，自下向上、由左向右、自上向下、由右向左只是一种较为常用的判断顺序，上面叙述的判断条件也是针对这一顺序而言的。

(2)追踪等值点

82

确定等值线追踪方向的目的，是为了正确地追踪出下一等值点C。等值线的追踪方向一旦确定，就可以在Ⅱ号矩形网格的除了包含B点以外的另三边（图3-24）上寻找C点的位置。由于这三边都有可能有等值点存在，所以等值线上B和C点的连接就可能有三种形式。图3-24中，C1,C2,C3为Ⅱ号矩形三边上可能存在的等值点，ABC（C等于C1或C2或C3）走向为ABC1,ABC2,ABC3之一。图3-25是图3-24

图3-24 等值线的 (a) (b) (c) 三种可能连结方式 图3-25 Ⅱ号网格等值点的连接方式

Ⅱ号矩形网格待同值等值点追踪完毕后的三种可能表现形式。显然，图3-25(c)是错误的，因为等值线发生了相交现象，这在等值线图上

83

是不允许出现的，必须加以排除。正确地选择ABC的连结方式，这在手工作图时比较容易解决，制图人员可以根据肉眼的观察、作图的经验和所研究空间变量的总体变化趋势等综合因素确定出等值线的正确走向，以避免等值线的相交或不适当的连结。但是，计算机自动完成此类工作就变得比较复杂，具有相同等值点数的数据构型，就可能连出不同形状的等值线，如图3-26。为了避免图3-26所示的多解性，就必须规定一些法则，使之作为追踪等值点的算法基础。综合手工作图的经验，结合计算机自动成图的特点，可以得出保证等值点联结唯一性的法则：先考虑等值线原来的前进方向，再考虑C1,C2,C3的分布情况以及它们与B点距离的远近。由于前面的内容已论述了确定等值线前进方向的算法，故这里的研究重点就是如何选择C点，即如何从C1,C2,C3中选择一点作为将要追踪出的等值点。为了避免出现图3-25(c)所示的错误，除非C1,C2,C3中只有C2点存在，否则就必须从C1,C3中选择一点，使之作为C点。从C1,C3中选取等值线时，以它们与B点距离最近的一点作为C点。这样做的原因有以下两个方面：1)两点距离最近，说明它们的相关性最大，两点的连线相对较为合理；2)可以避免图3-27中ABC3连线所带来的不良后果。图3-27中，因E点为奇异点（它的值与当前等值线相等）等原因，经数据预处理后，等值点B,C1,F可能离E点最近，显然ABC3连线不合实际情况，最佳连线为GABC1FD，这样就不会使A-B、FC1段分别位于不同的等值线上或在同一等值线上但不相邻，计算机绘出后产生等值线相交的假象。为了简化算法，以C1,C2点到B点所在矩形边的垂直距离近似代替它们之间的直线距离，一般情况下，此种方法的处理结果与实际情况完全一致。

综上所述，计算机自动追踪等值点时，需遵循以下三条规则： 1)C1,C3都存在时，选取与B点所在矩形边垂直距离最短的点作为C点。与C1,C2点相应的HH、SS数组值是判断距离远近的基础数据；

2)只有C1或C2或C3时，它们即为C点。

根据这些规定，就可以完成不同方向对等值点的追踪。

1)当j1j2时，表明当前追踪等值点的方向是由下向上，C点位于(i2j2),(i21,j2)纵边和(i2,j21)横边。根据追踪等值点的三原则，可得选择等值点C的顺序：如果两条纵边都有等值点，选取SS(i2,j2),SS(i21,j2)最小的点作为C点；如果两条纵边只有一个等值点，此点即为C点；如果两纵边都没有等值点，那么C点一定位于

84

(i2,j21)横边上；

2)当i1i2时，表明是由左向右追踪等值点，C点只能位于(i2j2),(i2,j21)横边和(i21,j2)纵边上。追踪等值点C的算法同样遵循三条规则；

3)排除以上两种情况，如果int(xB/a)*axB，表明当前追踪方向是自上向下，此时，等值点C只能位于(i2,j21),(i21,j21)纵边和(i2,j21)横边。同理，追踪等值点C的三条规则仍是算法的基础；

4)如果上述三种条件都得不到满足，那么目前追踪的方向就只能是由右向左。此时，等值点C位于(i21,j21),(i21,j2)横边和(i21,j2)纵边，利用三条规则就可正确地追踪出等值点C。

实际追踪等值线时，A、B、C三点处于动态变化之中，需不断更换A、B值。另外，因每一矩形网格纵横边最多只有一个等值点，所以，需要注销掉那些已追踪出等值点的边，以避免等值点的重复使用，造成等值线的混乱或计算机的死循环。相应的算法很简单，只要令相应的SS(i,j)或HH(i,j)=2即可。

3.开、闭曲线线头线尾的寻找

在任何一幅等值线图上，曲线都可分为开曲线和闭曲线两种类 型。因为制图区域总是有范围的，对于矩形网而言，它的四边即为制图边界。开、闭曲线可作如下定义：1)开曲线：曲线的首尾（线头或线尾）同时位于制图边界上。如图3-28(a)中的等值线A；2)闭

(a)寻找开曲线线头 (b)寻找闭曲线线头 图3-28 寻找开、闭曲线线头示意图

曲线：曲线上所有的点都位于制图边界内，且首尾点的x,y坐标值相等。见图3-28(b)中的等值线B。根据开、闭曲线的特点，就可以确定寻找线头、线尾的算法。

85

(1)开曲线线头的寻找

因线头位于制图区域的边界，故无需制图区域的所有纵横边都 参加运算和判断，只讨论位于边界上的纵横边，其顺序为：底边、左边、上边、右边（图3-28(a)）。

1)从底边上寻找开曲线线头。令i1,2,...,M1,在此循环过程中，只要0HH(i,j)1，就说明(i,1)横边包含有当前等值点，计算出此点的坐标，使之作为开曲线的线头。为了同一追踪算法，这里需假设A点的存在，它所在边的标识为(i,0)，即i1i,j10。这样，就可根据从下向上的追踪算法进行等值点的追踪；

2)从左边界寻找开曲线线头。令j1,2,...,N1，在此循环过程中，主要0SS(1,j)1，就可以确定开曲线线头所在的位置。找到线头后，令i10,j1j；

3)从上边界寻找开曲线线头。令i1,2,...,M1，在此循环过程中，只要0HH(i,N)1，就找到了等值线线头，此时，令i1i,j1N1；

4)从右边界寻找等值线线头。令j1,2,...,N1，在此循环过程中，只要0SS(M,j)1，就证明网格(M,j)纵边具有等值点，此等值点就是开曲线的线头。找到线头后，令i1M1,j1j。

(2)闭曲线线头的寻找

因闭曲线位于制图区域范围内，没有等值点位于网格边界上， 线头线尾重合，所以线头的寻找方法较开曲线简单，只需要从一个方向寻找即可。这里规定寻找方向是从左向右。具体算法如下：令i2,3,...,M1,j1,2,...,N1,只要0SS(i,j)1，就证明(i,j)纵边上有等值点存在，即找到了闭曲线的线头（图3-28(b)）。 (3)开闭曲线线尾的寻找或曲线终结的条件

根据开闭曲线的性质，可以把曲线追踪的结束条件归纳为如下两点：

1)追踪开曲线时，如果当前追踪出的等值线点C位于网格区域边界上，C点就是线尾点，曲线追踪完毕；

2)追踪闭曲线时，如果C点与线头的坐标值相等，说明闭曲线追踪完毕。

4.其它一些算法

(1)同值等值线分支的识别

86

在制图区域内，无论是开曲线还是闭曲线，某一数值的等值线可 能不止一条，它有数个分支，如图3-29中的100这条等值线。正识别分支的算法较为简单：待某条等值线追踪完毕后，继续寻找同值的线头并追踪出等值线，直到没有线头为此。只有同值等值线的所有分支追踪完毕后，才可进行其它数值等值线的追踪和识别。

为了避免开、闭曲线等值点的相互混用，造成等值线的混乱，具体编程时需先进行开曲线的追踪，再进行闭曲线的追踪。 (2)网格交点奇异z值的消除

自动追踪等值线时，有时会遇到如图3-30(a)所示的情况，网 格交点的z值与等值线值相等。如果不对这种情况进行处理，就会

造成等值线追踪的混乱。由于B点为Ⅱ号矩形网格的两边所共有，所以不能确定追踪方向是由下向上还是从左到右，但只要在100这个值的基础上加或减一个如0.01的很小的数，就可消除奇异点的影响（图3-30(b)）。这样，就把等值线的追踪方向变为由下往上。由于0.01这个数很小（如果等值线值很小，如0.1,那么加减的数应更小。一般情况下，假设等值线线的等间距为H0，那么，加或减的数为H0/5000.0），屏幕显示或其它输出设备输出等值线后，肉眼分辨不出因它而引起的差异，不会造成视觉上曲线位置的变化。对这种算法归

87

纳即为，追踪某一数值的等值线时，如某一网格交点的特征值与等值线值相等，那么在特征值上加或减一个如0.01的校正数。

3．4．2 三角网方法

虽然矩形网方法绘制等值线图具有程序设计简单、等值线不易相交等优点，但它要求原始数据点分布较为均匀或均匀，被研究的空间变量连续分布。实际上，这些条件是很难或不易满足的。针对这些情况，人们又发展了绘制等值线的不规则三角网（TIN）方法。这种方法不仅适合于规则分布的原始数据点，也适合于不规则、甚至畸形分布的原始数据点。通过三角网方法绘制等值线的过程是建立在TIN模型的基础上的。实际上，手工绘制等值线的过程也是采用三角网的方法：首先把相邻数据点联成三角网，再在三角形边内插当前等值点，通过肉眼的分析和观察，把相同的等值点连接起来就勾绘出了等值线。在建立TIN模型的基础上，自动生成等值线图的步骤如下：

1． 内插等值点的平面位置

假设三角形三顶点的特征值分别为Z1,Z2,Z3，那么，三角形边等 值点的分布情况如下：

(1) Z1,Z2,Z3不相等，其最大最小值区间内包含当前等值线值。只要每边两端点特征值满足：

(HZZ1)*(HZZ2)0,则该边无等值点，否则必有等值点； (HZZ1)*(HZZ3)0,则该边无等值点，否则必有等值点； (HZZ2)*(HZZ3)0,则该边无等值点，否则必有等值点。 图3-31给出了此种情况的示意图,HZ为当前等值线的等值线值。 (2) Z1,Z2,Z3都大于或小于当前等值线值，三边上都没有等值线点存在，见图3-32。

(3) Z1,Z2,Z3至少有两个相等，其最小值区间内包含当前等值线值，见图3-33。

88

(a) (b)

图3-31 等值线分布情况类型一示意图

从图3-31、图3-32、图3-33可以看出，任一具有等值点的三角形，等值点只位于三边中的两边，并且每边的个数为1，这说明三角形所包含的等值点数为2。根据这一特点，如果在三角形的某

边找到等值点后，那么另两边之一必定包含另一等值点。

确定三角形等值点所在边后，利用线性内插方法就可以求出等

89

图 3-34 计算三角形边等值点坐标示意图

值点坐标，其数学计算公式如下（图3-34）：

xCxA(xBxA)HZZA

ZBZAHZZA

ZBZA yCyA(yByA)式中HZ是当前等值线值。

在三角网生成等值线过程中，同样会遇到处理奇异点的问题， 处理方法与矩形网方法相同。

2． 追踪等值线

与矩形网方法追踪等值线一样，三角网方法同样存在着线头的 寻找、等值点的追踪、寻找结束标志等过程。

(1)寻找等值线的线头 1) 寻找开曲线线头

开曲线线头一定位于联网区域的边界上。对于现成的TIN，线 头位于那些只被一个三角形所有的三角形边上。

2) 寻找闭曲线线头

根据闭曲线的特点，线头的寻找只能在三角网的内部进行，即 闭曲线的线头位于被两个三角形所共有的三角形边上。

(2)等值点的追踪

找到开、闭曲线的线头后，就需要顺序地追踪出一条等值线的全部等值点。在三角网中，任一等值点都具有如下特点：它既是等值线

90

进入另一三角形的入口点，又是等值线走出当前三角形的出口点。这说明，除了边界等值点外，等值线上任一等值点都被两个相邻的三角形所共有。利用这一原理，就可以实现等值点的追踪。

1)得到共用当前等值点的三角形及该包含等值点的三角形边； 2)在新三角形的另两边上寻找等值点。

重复1)、2)，就能够得到等值线上的全部控制点。 (3)线尾的寻找

线尾的寻找与矩形网中的相应算法完全相同。

1)开曲线：如果追踪出的等值线点（不包括等值线点序列的第一点），位于TIN边界上，该点即为线尾点；

2)闭曲线：如果追踪出的等值线点与等值线点序列的第一点完全相等，该等值线追踪结束。

无论是开曲线、还是闭曲线，同一数值的曲线可能有很多分支，所以追踪完一条等值线后，还需继续追踪与之同值的分支。待同一数值的所有分支追踪完毕后，才能进行下一数值等值线的追踪。

3．6 曲线的光滑

这里我们主要介绍张力样条函数。 主要原理如下：

在GIS系统中，相当多的线划图形实体都是光滑的曲线（如等值线，等高线，剖面线等）。所以，正确地选择满足应用要求的曲线光滑算法具有十分重要的意义。在现有的商业化GIS或CAD系统中，用于曲线光滑的算法主要有如下几种：线形迭代法；Bezier函数法；五点光滑法；三次样条函数法等。通过对相关算法数学原理的分析研究可知，前两种方法拟合的光滑曲线并不通过已知的结点（图3-

图3-35 曲线光滑示意图

91

35(a)中的A、B、C）。显然，对于许多应用这是不允许的；后两种方法拟合的光滑曲线在大多数情况下满足实际需求，但如果变量在空间的分布变化较大，那么可能造成相邻曲线的相交（图3-35(b)中的L1、L2）。在实际应用中，我们认为，如果曲线的光滑算法满足以下两个条件，即既可以形成光滑的曲线，也可以形成分段的线性函数，那么此种算法将是十分优越的，这样，既可以保证曲线的光滑，又可自动地尽量避免曲线的相交。通过研究发现，张力样条函数满足以上条件。

1. 张力样条函数光滑曲线的数学原理

张力样条函数是Schweikert为消除三次样条插值函数有时会出 现多余的拐点而引入的。其基本构思是分段插值函数为直线插值和两个双曲函数shx和chx的线性组合： f(x)c1c2xc3shxc4chx (3-16)

在（3-16）式中，为张力系数。

（3-16）式整体具有连续二阶导数，的作用是控制拐点的位置和曲线的形态。 通过有关复杂的数学推导，可以得到张力样条函数的函数表达式：

fx1fxishxi1xfxi1shxxi2shhifxixi1xfxi1xxiyiyi122hihi

（3-17）

式中xixxi1hixi1xii0,1,2,...,n12. 张力样条函数的特点

(1)如果0，那么

fxfxixi1xxxi （3-18） fxi1hihi 上式为一标准的三次样条函数表达式。三次样条函数可以保证曲

线既光滑，又通过已知数据点。 (2)当时

fxyixi1xxxi （3-17） yi1hihi 92

此时张力样条函数退化为分段线性函数。分段线性函数可以绝对

保证线段不相交。

从张力样条函数优良的数学特性可以看出，对某一地质变量分布区域，只要选择合适的张力系数，就可以尽量避免合乎要求的光滑曲线相交。

图

3-36 张力样条函数光滑曲线示意图1

3. 确定张力系数的算法研究

合适的张力系数可以避免相邻光滑曲线的相交。根据研究，在 不考虑相邻线段相交而且得到光滑曲线的前提下，任一光滑曲线的的计算方法都是一样的，既如（3-18）式所示。

.n/s[n] （3-18） 15

（3-18）式中，n,s[n]分别为线段的原始数据点数和累计长度。

为了使张力样条函数适应不同复杂程度地质变量曲线的拟合和光滑的要求，本文在经过大量实验的基础上提出了逐点计算的算法。算法的主要步骤如下：

(1)通过（3-18）式计算任一曲线的基本张力样条函数，并赋予任一原始数据点的属性为0；

(2)判断相邻线段是否相交，如果相交，记录下相交线段的首尾点编号，这些点的属性值为1；

(3)逐点赋予张力样条函数值，相应的规则如下：属性值为0的数

93

据点的张力系数值为，属性值为1的张力系数值为7.0。

图3-36、图3-37是利用上述算法处理的两幅图形。虽然与这两幅图形相应的空间变量前者简单，后者复杂，但它们都满足实际需求。

图3-37 张力样条函数光滑曲线示意图2

/* ------------------------- 张力样条函数拟合曲线-------------------------*/ #include //头文件 #include #include #include #include #include

extern double map_xmi,map_xma;//图形区域x的最大最小值 extern int id_cur;

extern long map_scale;//图形的比例尺 //计算由两点构成线的角度

extern double angle(double x1,double y1,double x0,double y0); //n:曲线的原始数据点数

//id:是否重新计算拟合曲线的参数 //x,y,z：曲线的三维坐标

int splin(int n,int id,double *x,double *y,double *z) {

94

double *xx,*yy,*zz;

double *xp,*yp,*s,*a,*b,*c,*hp;

double dtx1,dtx2,ds1,dty1,dty2,six1,siy1,sixn,siyn; double dx1,dy1,ds2,dx2,dy2,sma,ds,dlt;

double d1,d2,d3,q,ss,e1,e2,e3,x1,y1,alf,zk; int i,i1,k,k1,n1,j,p_number,mao; if(n<=2) return n; mao=n;

dlt=(map_xma-map_xmi)/600.0; if(dlt<0.001) dlt=1.0;

xx=(double *)malloc((n+2)*sizeof(double)); if(xx==NULL) exit(1);

yy=(double *)malloc((n+2)*sizeof(double)); if(yy==NULL) exit(1);

zz=(double *)malloc((n+2)*sizeof(double)); if(zz==NULL) exit(1);

xp=(double *)malloc((n+2)*sizeof(double)); if(xp==NULL) exit(1);

yp=(double *)malloc((n+2)*sizeof(double)); if(yp==NULL) exit(1);

s=(double *)malloc((n+2)*sizeof(double)); if(s==NULL) exit(1);

a=(double *)malloc((n+2)*sizeof(double)); if(a==NULL) exit(1);

b=(double *)malloc((n+2)*sizeof(double)); if(b==NULL) exit(1);

c=(double *)malloc((n+2)*sizeof(double)); if(c==NULL) exit(1);

hp=(double *)malloc((n+2)*sizeof(double)); if(hp==NULL) exit(1); memcpy(xx,x,(n+2)*8); memcpy(yy,y,(n+2)*8); n--; n1=n-1; if(id==0){ alf=angle(xx[1],yy[1],xx[0],yy[0]); alf=alf*3.14159/180.0; *z=alf; }

else alf=*z;

dtx1=xx[1]-xx[0]; dty1=yy[1]-yy[0]; ds1=hypot(dtx1,dty1); six1=cos(alf);

95

siy1=sin(alf);

if((xx[0]==xx[n])&&(yy[0]==yy[n])){ sixn=six1; siyn=siy1; *(z+n)=alf; } else{ if(id==0){

alf=angle(xx[n],yy[n],xx[n1],yy[n1]); alf=alf*3.14159/180.0; *(z+n)=alf; }

else alf=*(z+n); sixn=cos(alf); siyn=sin(alf); }

dx1=six1; dy1=siy1;

xp[0]=dx1-six1; yp[0]=dy1-siy1; hp[0]=ds1; s[0]=0.0; s[1]=ds1; d3=ds1;

for(i=1;i<=n1;i++){ dtx2=xx[i+1]-xx[i]; dty2=yy[i+1]-yy[i]; ds2=hypot(dtx2,dty2); if((id2==0)||(*(z+i)<0.0)){ alf=angle(xx[i+1],yy[i+1],xx[i],yy[i]); alf=alf*3.14159/180.0; *(z+i)=alf; } else alf=*(z+i); dx2=cos(alf); dy2=sin(alf); xp[i]=dx2-dx1; yp[i]=dy2-dy1; hp[i]=ds2; dx1=dx2; dy1=dy2; i1=i+1; s[i1]=s[i]+ds2; d3=d3+ds2;

96

}

d3=d3*(1000.0/map_scale);

if(d3<200.0) dlt=1.5*(map_scale/1000.0); if(d3<100.0) dlt=0.75*(map_scale/1000.0); if(d3<20.0) dlt=0.4*(map_scale/1000.0); if(d3<10.0) dlt=0.15*(map_scale/1000.0); if(d3<5.0) dlt=0.1*(map_scale/1000.0); xp[n]=sixn-dx1; yp[n]=siyn-dy1; sma=1.5*n/s[n]; ds=sma*hp[0];

d1=sma*(((exp(ds)+exp(-ds))/2)/sinh(ds))-1/hp[0]; b[0]=d1;

c[0]=1/hp[0]-sma/sinh(ds); a[0]=c[0];

for(i=1;i<=n1;i++){ ds=sma*hp[i]; c[i]=1/hp[i]-sma/sinh(ds); a[i]=c[i-1]; d2=sma*(((exp(ds)+exp(-ds))/2)/sinh(ds))-1/hp[i]; b[i]=d1+d2; d1=d2; }

b[n]=d1;

xp[0]=xp[0]/b[0]; yp[0]=yp[0]/b[0]; q=b[0];

for(k=0;k<=n1;k++){ b[k]=c[k]/q; q=b[k+1]-a[k]*b[k]; k1=k+1; xp[k1]=(xp[k1]-a[k]*xp[k])/q; yp[k1]=(yp[k1]-a[k]*yp[k])/q; }

for(i=0;i<=n1;i++){ k=n-i-1; xp[k]=xp[k]-b[k]*xp[k+1]; yp[k]=yp[k]-b[k]*yp[k+1]; }

memcpy(zz,z,(mao+2)*8);

x[0]=xx[0];y[0]=yy[0];z[0]=zz[0]; p_number=1;

for(i=0;i<=n1;i++){ zk=(int)hp[i];

97

}

k=(int)(zk/dlt); for(j=1;jp_number=p_number+1; }

free(xp);free(yp);free(hp);free(s);free(a);free(b);free(c); free(xx);free(yy);free(zz); return(p_number);

3．7 栅格、矢量数据的相互转换

地理信息系统空间数据类型主要有矢量和栅格结构。矢量结构包含有拓扑信息，通常应用于空间关系的分析；栅格数据则易于表示面状要素，主要应用于空间分析和图象处理。由于栅格和矢量数据在GIS应用过程中各有其优缺点，所以，一般情况下，同一个GIS系统能够处理、存储栅格和矢量数据。对同一研究区域而言，有时为了分析处理问题的方便，需要实现栅格和矢量数据间的转换（如扫描图象的矢量化，地形图的栅格化）。

98

3．5．1 矢量向栅格的转换

图3-38 栅格单元属性值的确定

从矢量向栅格转换过程中，应尽量保持矢量图形的精度。在决定属性值时尽可能保持空间变量的真实性和最大信息量。在图3-38中，格网单元对应几种不同的属性值，而每一单元只能取一个值。在这种情况下，有如下一些取值方法。

(1)中心点法：用处于格网单元0处的地物类型或空间特征决定属性值。此时，该单元属性值为C。此法常用于连续分布的地理要素，如降雨量分布、大气污染等；

(2)面积占优法：以占单元面积最大的地物类型和空间特征决定格网单元的属性值。此时，栅格单元的属性值为B。面积占优法适合分类较细、地物类别斑块较小的情况；

(3)重要性法：根据格网单元内不同地物的重要性，选取最重要的地物类型代表相应的格网单元的属性值。这种方法对于特别重要的地理实体，尽管其面积很小或不在格网的中心，也采取保留的原则。重要性法常用于具有特殊意义而面积较小的地理要素，特别是具有点、线状分布的地理要素，如城镇、交通枢纽、河流水系等。

在进行弧段或多边形的矢量化时，可以利用上述三种方法确定格网的属性值。

为了逼近原图或原始数据精度，除了采用上述几种取值方法外，还可以采用提高图象分辨率的方法。这种方法可以提高转换的精度，更接近真实状态，表现更细小的地物类型。当然，图象分辨率的提高，将大大增加数据量。

1． 基于弧段数据的栅格化方法

99

首先计算所有弧段结点或中间点所在的格网位置，并赋予该结 点正确的属性，然后根据下面的算法完成弧段的栅格化。

图 3-39 基于弧段数据的栅格化方法示意图

如图3-39所示。利用弧段的数据列与格网的行列线相交，以得 到正确的栅格化结果。实际计算时，需逐段处理弧段中的局部直线段，待处理完某一局部线段后，再进行下一局部线段的处理（如处理完AB后，再处理BC段），直至完成整条弧段的处理。局部线段与行列线求交后，存储交点坐标，并对x或y从小到大或从大到小排序。根据排序结果，相邻交点所构成线段通过的格网需赋予属性值。

2． 基于多边形数据的栅格化方法 如图3-40。可以采用的算法如下；

(1)得到构成多边形的数据列（如ABCDEFGHA）；

图 3-40 基于多边形数据栅格化示意图

(2)根据拓扑包含关系得到多边形内的比当前多边形低一个层次 的所有多边形的数据（如JKLMNOPJ，QRSTUQ）；

(3)逐行或逐列与(1)、(2)得到的数据列求交，得到所有的交点。

100

在进行下一步骤时，需对交点等于结点或中间点的情况进行处理。 如果与交点相邻的两点位于行或列的同一侧，那么，不删除与当前交点相等的另一交点（如图3-41中A点）,反之，需删除另一交点(图3-41中B点)；

(4)以x或y对交点进行由小到大的排序；

(5)奇数点至偶数点间的格网需赋予多边形的属性值，而偶数至 奇数间的格网则不处理。

图 3-41 多边形栅格化奇异点处理示意图

3．5．2 栅格向矢量的转换

把栅格单元中的空间信息转换为几何图形的过程叫矢量化。矢量化的过程要保证以下两点：

(1) 拓扑转换，即保持栅格表示出的连通性和邻接性。否则，转 换出的图形是杂乱无章的，没有任何实用价值的；

(2) 转换空间对象正确的外形。 栅格向矢量转换的主要步骤为： (1) 二值化：一般情况下，栅格数据是按0～255的不同灰度值 表达的。为了简化追踪算法，需把256个灰阶压缩为2个灰阶，即0和1两级。为此，假设任一格网的灰度值为G(i,j),阈值为T，那么，根据下式就可以得到二值图。 B(i,j)10如果G(i,j)T如果G(i,j)T

(2) 细化：细化是消除线划横断面栅格数的差异，使得每一条线 只保留代表其轴线或周围轮廓线（对多边形而言）位置的单个栅格的

101

宽度。对于栅格线划的细化方法，一般采用“剥皮法”。剥皮法的实质是剥掉等于一个栅格宽的一层，直到最后留下彼此连通的由单个栅格点组成的图形。

(3) 跟踪：跟踪的目的是把细化后的栅格数据整理为从结点出发 的线段或闭合的线条，并以矢量形式加以存储。跟踪时，根据人为规定的搜索方向（如沿图幅边界的顺时针或逆时针方向），从起始点开始，在保证趋势的情况下对八个邻域进行搜索，依次得到相邻点，最终得到完整的弧段或多边形。

(4) 去除多余点及曲线光滑：由于搜索是逐个栅格进行的，所以， 弧段或多边形的数据列十分密集。为了减少存储量，在保证线段精度的情况下可以删除部分数据点。可以采用如下算法删除多余点：图3-42，计算当前点A与相邻点B、C组成的线段BA、AC的夹角BAC，如果BAC大于某一固定值（如175度），就删除B点。当然，如果删除B点后，AC距离太长，那么，删除B点是不合适的。

(5) 拓扑关系的生成：判断弧段与多边形间的空间关系，以形成 完整的拓扑结构并建立与属性数据的关系。

图3-42 删除多余点示意图

3．8 不同图象、图形格式数据的转换

3．6．1 不同图象格式间的转换

常见的图象格式有BMP、GIF(Graphics Interchange Format)、TIFF（标记图象文件格式，Tagged Image File Format）等。由于同一种类型图象格式在世界范围内是完全统一的，所以，它们在不同的GIS系统间的数据交流应该不成问题。这里所讲的数据转换是指不同格式图象数据间的转换。

每一种图象文件的格式都具有本身的特点和规定，只要知道了这些规定，就可以方便地完成数据间的转换。实际上，BMP等图象文件

102

的格式是已知的。

3．6．2 不同GIS矢量结构数据的转换

不管是GIS的通用开发平台，还是专业应用系统，相应的数据库（图形或属性）文件结构和图形实体的数据结构都具有或多或少的差别，特别是专业系统更是如此。为了实现数据的共享，就必须实现系统间的数据转换。一般情况下，可以利用如下方法完成数据的转换工作。

1．如果某一系统具备把数据库内容（包括文件参数，图形实体的数据）转换为文本文件(*.dat,*.txt)的功能，那么，可以方便地对这些数据进行处理，通过编程把这些数据读入自己的系统中；

2．如果已知某一系统的数据文件结构、图形实体的数据结构，而通过这些结构存入外部设备（如硬盘）的数据没有经过加密处理，那么，可以编程直接读取数据库数据。

3．如果已知某一系统的数据文件结构、图形实体的数据结构，而通过这些结构存入外部设备（如硬盘）的数据已经经过加密处理，那么，必须先得到加密的具体算法，再通过编程直接读取数据库数据。

4．到目前为此，不同的GIS系统都可以把矢量图形数据转换为DXF(Drawing Interchange Infomat)格式或读入DXF格式的数据文件，而DXF格式在全世界是统一的，而且是公开发布的。所以，不同系统间可以通过DXF格式进行数据交流。

利用GIS技术设计一个北京大学实时学生上课情况查询系统（假设学生的位置能够精确定位）。系统的功能如下：

1． 查询某一学生所在的教室； 2． 查询某一学生相邻的学生姓名； 3． 统计某一教室的学生总人数； 4． 统计全校的学生上课人数。

假设教室、学生的分布如下图所示，请提出数据模型和数据结 构及其它处理方法。

课堂讨论

103

图 3-43教室及学生位置分布图

104