基于球面QTM的大范围有源淹没算法研究

邢华桥，侯妙乐，王　磊，姜晓轶

(1. 中国矿业大学(北京)地球科学与测绘工程学院，北京 100083； 2. 北京建筑大学，

北京 100044； 3. 国家海洋信息中心，天津 300171)

Submergence Analysis Algorithm for Large Area Based on Spherical QTM

XING Huaqiao，HOU Miaole，WANG Lei，JIANG Xiaoyi



基于球面QTM的大范围有源淹没算法研究

邢华桥1，侯妙乐2，王磊1，姜晓轶3

(1. 中国矿业大学(北京)地球科学与测绘工程学院，北京 100083； 2. 北京建筑大学，

北京 100044； 3. 国家海洋信息中心，天津 300171)

Submergence Analysis Algorithm for Large Area Based on Spherical QTM

XING Huaqiao，HOU Miaole，WANG Lei，JIANG Xiaoyi

摘要：面对海平面上升后沿海区域大范围淹没模拟的应用需求，提出了一种基于球面QTM的大范围有源淹没算法，首先构建了研究区域的QTM三角格网集，并利用海陆分离的思想确定了初始淹没单元；然后基于QTM的12邻近搜索设计了面向多目标三角格网的淹没范围扩张方法，将初始淹没单元递归扩张得到最终的淹没范围；最后应用Visual C#语言及Microsoft DirectX三维图形接口设计开发了相应的试验原型系统，并利用天地图高分辨率影像和SRTM地形数据对试验区域进行了海平面上升后的大范围淹没模拟试验。结果表明，该算法具有较好的实用性，可以有效地对海平面上升后的大区域范围进行动态地淹没模拟，对辅助防灾减灾具有一定意义。

关键词：球面QTM；海平面上升；海陆分离；12邻近搜索；有源淹没

一、引言

随着全球气候变暖等极端天气的频发，海平面异常出现的次数明显增多，异常值明显增大，近30年来全球海平面上升速率有加大的趋势[1]，而同期中国沿海海平面的上升速率高于全球的平均值[2-3]。海平面上升带来的最为直接和严重的影响是沿海地区的淹没问题。因此，准确、科学地预测海平面上升后的淹没范围，对海洋灾害的预测预报、灾害过程的动态模拟，以及促进海洋灾害机理的研究、辅助防灾减灾等都具有重要意义。

近年来，国内外研究学者对淹没分析的研究取得了一定的进展，淹没分析主要分为无源淹没和有源淹没[4-6]，无源淹没将所有高程值低于给定水位的区域均记为淹没区，有源淹没除了考虑上述情况外，还考虑到了“流通”的淹没情形，即只有与淹没区域相连通且高程值低于给定水位的区域才被视为淹没区。目前的淹没分析算法多以研究区域的DEM数据为基础[7]，利用种子蔓延算法进行淹没区域的模拟。如张东华等提出一种基于DEM的洪水有源淹没算法，在GIS技术的基础上应用数字高程模型(DEM)的格网模型进行洪水淹没分析[8]；丁志雄等提出基于格网模型的洪水淹没分析算法，在给定洪水水位和洪量两种条件下，分别基于三角形格网模型和任意多边形格网模型进行了水淹模拟[9]；聂汉江等提出了基于库群联合调度和DEM的区域洪水淹没范围模拟方法，以联合调度模型得到的累积成灾水量为淹没水量，以研究区域的DEM为数据基础，运用体积法对一定体量的洪水淹没进行了模拟[10]。但这些研究大都是针对小范围区域，采用传统的平面数据模型进行的洪水淹没分析，而海平面上升带来的淹没通常范围较大，传统的淹没分析方法并不能满足这种大范围淹没分析的需求。

针对这一问题，本文基于球面四元三角网(quaternary triangular mesh，QTM)的12邻近搜索，提出一种面向大范围区域的有源淹没算法，以实现海平面上升后研究区域的淹没动态模拟。

二、QTM理论知识

QTM是由内接于球体的正八面体经多次递归剖分而成的球面拟合三角格网集[11]，具有符合球面的特性和本质是栅格的特点，可以用QTM模拟地球表面；另外，QTM形状和面积的近似相等能保证格网之间具有简洁的邻近关系，便于格网间的邻近查询、分析等处理[12]。

1. QTM剖分与编码

在QTM初始化剖分时，内接于球体的正八面体顶点占据球面主要点(包括两极)，而边的投影与赤道、主子午线，以及90°、180°、270°子午线重合，如图1所示。在对初始的QTM进一步细化剖分时，用大弧平分法找到三角形边的中点，通过连接它们将球面三角形分成4个小的三角形，以此类推进行递归，对整个球面进行近似均匀的剖分[13]，形成了全球的QTM离散格网集。

图1　QTM初始剖分单元

QTM编码不仅暗含着格网的空间位置，而且表示了格网的剖分层次，因此QTM编码是进行邻近搜索及连通性分析的基础。本文采用的编码方案为固定方向编码[12]，与其他编码方案相比，该方案的所有格网都具有固定方向，更利于邻近搜索。如图2所示，编码的第1位为八分体号，从第2位开始，0表示中间三角格网，1表示顶(底)三角格网，2表示左三角格网，3表示右三角格网，即可以用一个四进制的Morton码来标示。因此，QTM的编码M由八分码M0加上Morton码组成，即M=M0q1q2…qi。

图2　固定方向编码方案

2. QTM的邻近关系

邻近关系是GIS空间分析中必不可少的一类空间关系，是球面实体扩张和空间索引的基础[14-15]，在实际应用中(如水淹分析、最短路径分析等)具有重要的意义。如图3所示，QTM三角格网的邻近三角格网分为具有公共边的边邻近(edge-adjacent)三角格网(编号为E)和具有公共点的角邻近(vertex-adjacent)三角格网(编号为V)。将球面按QTM剖分后，每个球面三角格网都有3个边邻近三角格网，位于南、北极的8个顶角三角形和赤道上经度为0、90°、180°(-180°)、-90°的16个顶角三角形均具有7个角邻近三角形(如图3(a)所示)，其他位置的非顶角三角形均具有9个角邻近三角形(如图3(b)所示)。根据QTM的邻近关系，其邻近搜索分为三邻近搜索和12邻近搜索，由于海水淹没不具有方向性，即海水可能向所有方向蔓延，因此，利用12邻近搜索更符合实际情况。本文主要考虑QTM格网的12邻近搜索[16]。

图3　QTM的邻近三角格网

三、基于QTM 12邻近搜索的有源淹没算法

1. 总体思路

首先构建研究区域的一定层次的QTM格网集，基于海陆分层的思想确定初始淹没单元，并以此为基础进行12邻近搜索获取其邻近格网，如果邻近格网的高程低于给定的上升水位，即HQTM

图4　基于QTM 12邻近搜索的有源淹没算法流程

2. 基于海陆分离思想的初始淹没单元判定

确定初始淹没区域是海平面上升分析与模拟的首要问题，海平面上升后的淹没过程是从陆地与海洋相接处开始的。因此，本文利用海陆分离的思想，将研究区域的QTM按照3个顶点的高程划分为陆地和海洋，并确定陆地边缘格网。具体方法为：对于每个格网，如果其3个顶点均位于海洋(顶点的高程值均为-9999 m，如图5(a)所示)则为海洋格网；如果其3个顶点均位于陆地(顶点的高程值均大于-9999 m，如图5(b)所示)则为陆地格网；如果其3个顶点并不完全位于陆地或海洋，那么将该格网视为陆地边缘格网(如图5(c)所示)，所有陆地边缘格网构成的集合即为初始淹没单元。

图5　QTM三角格网陆地与海洋的判别

3. 面向多目标三角格网的淹没范围扩张

由淹没单元通过邻近搜索扩张得到其邻近格网时，淹没单元往往不是单独的一个QTM格网，而是多个连续的格网，这就对邻近搜索提出了特殊的需求，即需要对多个三角格网同时进行扩张，而面向单目标三角格网的扩张方法会造成格网的重复扩张，从而增加了系统的计算量。本文在淹没范围的扩展过程中，剔除扩张得到的邻近格网中的重复格网，同时为进一步减小工作量，每一次扩张只将上一次扩张得到的被淹没格网向外扩张，而不再处理已经判断过的格网。面向多目标三角格网的淹没范围扩张过程如图6所示。该扩张方法的过程具体如下：

1) 利用12邻近搜索算法依次计算初始淹没单元S中所有格网的邻近格网，并从中剔除重复格网，将不重复的邻近格网集合记作G。

2) 利用双线性插值算法计算格网集G中每个格网中心点的高程。

3) 将格网集G中每个格网的高程与当前水位数据进行对比，高程低于当前水位的格网为被淹没的格网，记作G′。

4) 以G′代替步骤1)中的初始淹没单元S，执行步骤1)—步骤3)，直至G中所有格网的高程均高于当前水位(即G′为空)或整个试验区的格网搜索完毕。

图6　面向多目标三角格网的淹没范围扩张过程

四、试验与分析

为了对算法进行验证，本文以C# 3.0为开发语言，DirectX为三维开发包，从底层开发了试验原型系统。试验硬件环境为Pentium(R) Dual-core CPU T4200 @2.00 GHz，2.00 GB内存，采用SRTM 地形数据作为DEM数据源，天地图高分辨率影像为影像数据源，在给定淹没水位的情形下对海平面上升后研究区域进行动态淹没模拟。试验中首先构建了一定层次的QTM三角格网集作为研究区域的模拟图层，然后利用海陆分离思想提取试验区域的初始淹没单元，如图7所示。

图7　试验区域的初始淹没单元

进而利用本文提出的算法在给定淹没水位下获取所有被淹没的格网存入淹没格网集，并将淹没格网集中所有格网进行着色和渲染，然后与无源淹没算法的淹没效果进行了对比。图8为淹没水位为5 m的试验区域淹没范围的模拟结果，图9为利用无缘淹没算法得到的试验区域淹没范围。通过淹没模拟对比图可以看到，利用本文算法得到的淹没区域具有连通性，更加符合淹没的实际情况。

为了进一步比较本文算法与无源淹没算法的区别，并验证本文算法在大区域范围进行淹没模拟的有效性，在相同的试验硬件环境下，应用不同的淹没水位分别在不同面积的大范围试验区域进行有源淹没和无源淹没模拟，并比较了不同淹没算法的淹没面积，见表1。试验结果表明，随着淹没水位的升高，淹没范围逐渐增大；在利用相同淹没水位对同一的试验区域进行模拟时，本文算法得到的淹没范围相对于无源淹没算法较小，这是由于本文算法考虑到了淹没区域的连通性，只有高程低于给定水位且与初始淹没单元相连通的区域才被视为最终的淹没区域。此外，与传统的有源淹没算法相比，本文算法所应用的试验区域达到了数万平方千米的级别，可以有效地应用于更大范围的淹没分析中， 并且能将淹没结果动态地在三维球面上进行模拟和展示。

图8　本文算法得到的淹没模拟图

图9　无源淹没算法得到的淹没模拟图

km2

五、结论与展望

本文针对海平面上升后的沿海区域淹没模拟问题，提出了基于球面QTM 12邻近搜索的大范围有源淹没算法，并设计开发了相应的原型试验系统。试验证明该方法可以用于球面大范围海平面上升后的水淹动态模拟，能够为防灾减灾提供决策依据。但该算法也有一定的局限性，海平面上升是个渐变的过程，其年变化范围在毫米量级，而本文假定的海平面上升水位实际是叠加了平均潮位和极值潮位的最终上升水位。如何根据真实的海平面上升高度模拟淹没范围，并对海平面上升后给淹没区域带来的社会、经济影响等进行分析将是下一步的研究重点。

参考文献：

[1]李响，刘克修， 董军兴，等. 中国沿海海平面上升风险区划方法研究[J].海洋预报， 2014，31(2)：41-49.

[2]周子鑫.我国海平面上升研究进展与前瞻[J]. 海洋地质动态，2008， 24(10)：14-18.

[3]王慧，刘克修，范文静，等. 2012年中国沿海海平面上升显著成因分析[J].海洋学报， 2014，36(5)：8-17.

[4]刘仁义，刘南. 基于GIS的复杂地形洪水淹没区计算方法[J].地理学报，2001，56(1)：1-6.

[5]刘小生，陈英俊，黄玉生.基于GIS技术的洪水淹没区确定[J]. 测绘科学，2007，32(5)：136-137.

[6]王建军，吴北平，朱万里，等.基于ArcGIS的风暴潮灾害疏散系统的研究与开发[J].测绘通报，2013(9)：91-94.[15]国家测绘局．GB/T13990—921∶5000、1∶10000地形图航空摄影测量内业规范[S].北京：中国标准出版社，1992.

[7]韩丽蓉，张涛.基于RS/GIS的大南川水库溃坝灾害模拟方法研究[J].测绘通报，2014(1)：97-99.

[8]张东华，刘荣，张咏新，等.一种基于DEM的洪水有源淹没算法的设计与实现[J].东华理工大学学报:自然科学版，2009，32(3)：181-184.

[9]丁志雄，李纪人，李琳.基于GIS 格网模型的洪水淹没分析方法[J].水利学报，2004(6)：56-60.

[10]聂汉江，刘彬，何立新，等. 基于库群联合调度和DEM的区域洪水淹没范围模拟[J].水电能源科学，2014，32(5)：43-46.

[11]DUTTON G. Encoding and Handling Geospatial Data with Hierarchical Triangular Meshes [C]∥Proceeding of 7th International Symposium on Spatial Data Handling. Netherlands: Delft University of Technology, 1996. 34-43.

[12]SAHR K, WHITE D，KIMERLING A. Geodesic Discrete Global Grid Systems[J]. Cartography & Geographical Information Science, 2003,30(2): 121-134.

[13]LEE M，SAMET H. Navigating through Triangle Meshes Implemented as Linear Quadtree[J]. ACM Transactions on Graphics, 2000,19(2): 79-121.

[14]陈军，赵仁亮.GIS 空间关系的基本问题与研究进展[J]. 测绘学报，1999，28(2)：96-102.

[15]侯妙乐，邢华桥，赵学胜，等.球面四元三角网的复杂拓扑关系计算[J].武汉大学学报:信息科学版，2012，37(4)：468-471.

[16]赵学胜，侯妙乐，白建军.全球离散格网的空间数字建模[M].北京：测绘出版社，2007：40-46.引文格式： 颜金彪，胡最，禹信. CORS与静态相对定位技术下的坐标转换研究[J].测绘通报，2015(12)：54-56.DOI:10.13474/j.cnki.11-2246.2015.377

作者简介：邢华桥(1988—)，男，博士生，研究方向为全球离散格网、服务计算。E-mail：xinghuaqiao@126.com

基金项目：国家自然科学基金(41171304；40701152；41171306)

收稿日期：2015-02-26； 修回日期： 2015-10-27

中图分类号：P208

文献标识码：B

文章编号：0494-0911(2015)12-0046-04