基于多重分形理论的海州湾海底地形构建

吉玮，朱文菊，李超，毕建涛

(1．国家海洋局数字海洋科学技术重点实验室，天津300171;2．中国科学院 遥感与数字地球研究所 数字地球重点实验室，北京100094;3．北京中科数遥信息技术有限公司，北京100101;4．山东省建设发展研究院，山东济南250001)

海上交通、海洋调查、海洋资源开发、海洋工程建设、海洋疆界勘定和海洋环境保护等，都需要了解海洋水下地形这一重要要素。但由于海上多变的自然环境造成的较高作业难度以及多波束测深仪造价昂贵等原因，中国海洋地形测量目前仍以单波束测深仪获取点源数据为主。这就需要对点源数据进行插值或估计而形成连续变化的曲面来对水下地形的变化进行表达。

地形高程插值方法最常见且使用较多的是克里格 (Kriging)插值法［1-2］，但是一般的Kriging方法为了体现某种趋势而采用加权平均估计，将实际高程值较大的地方消减，实际高程值较小的地方增大，对数值造成一定的平滑，结果使高程局部结构特征的表示与实际地形有较大差距。Mandelbrot在20世纪60年代提出的分形理论作为处理地形非线性和复杂性研究的有力工具，能够对系统演变过程中的复杂性、不规则性和不均匀性进行度量，且已在海岸地貌、流水地貌和喀斯特地貌等方面的研究中得到了广泛应用［3-6］，其不仅能够用来描述分形的复杂特征，还能特征描述分形本身的几何支撑度量。对于高程、水流等诸多非均匀的分形现象，可以使用多重分形测度或维数的连续谱来表示［7］。目前，多重分形方法已应用到地球科学的诸多领域［8-13］，但应用于海底地形插值方面的研究还不多见。本研究中，以海州湾作为研究区域，通过数字化海图水深数据，将多重分形理论引入到Kriging这一常见的方法中构建海底地形，并将插值结果与普通Kriging法进行比较，旨在对海洋水深这一非均匀分形现象进行更精确地预测。

1 研究区域概况与分析方法

1.1 研究区域概况

海州湾位于苏鲁隆起与苏北之南黄海拗陷的过渡地带，北起山东日照市岚山区的佛手咀 (35°05'55″N，119°21'53″E)，南至江苏连云港市高公岛(34°45'25″N，119°29'45″E)，是一个濒临黄海的大喇叭口形开敞海湾，湾口宽42 km，岸线长86．81 km，海湾面积为 876．39 km2［14-15］。由于供沙条件、水动力条件和岸坡形态的不同，海州湾地貌特征和冲淤动态各异。自20世纪30年代建港以来，该地区各种海洋建设工程持续进行，各类近海岸海洋活动频繁，特别是1994年拦海西大堤的建成进一步影响了海底的沉积环境［16-17］。海州湾作为江苏沿海开发的龙头型增长极以及连云港实施海湾型经济的重要载体，对这一区域海底地形的构建能对相关基础设施建设、旅游资源的利用与开发，以及渔业养殖等海洋经济产业的规划提供较为准确的科学依据。

1.2 研究数据

本研究中，海底地形构建使用的水深数据来源于2000年12月出版的日照港至灌河口海图，比例尺为1∶120 000，采用墨卡托投影，坐标系为1954北京坐标系，高程基准为1985国家高程基准，基本等高距为40 m。海图数据的采集单位为中国人民解放军海军司令部航海保证部，图中水深数据由海军海洋测量船使用测深仪测得。将海图扫描之后利用ArcGIS软件完成投影转换与配准的预处理工作，再在ArcGIS中采用人工矢量化的方式采集水深点数据，在约4265 km2研究区范围内共采集909个水深点，水深为0．05～37．00 m，平均水深为11．85 m。总体样本水深有一定差异，符合本研究需求，从中随机抽取709个作为插值样点用于插值研究，其余200个水深点用作检验 (图1)。

1.3 分析方法

1．3．1 Kriging插值法 Kriging插值法又称空间局部插值法，是以变异函数理论和结构分析为基础，在一定的区域内对区域化变量进行无偏最优估计的一种方法，是地质学研究的主要方法之一。Kriging法是根据待插值点与临近实测点的空间位置，对待插值点的高程值进行线性无偏最优估计，通过生成一个关于高程的Kriging插值图来表达研究区域的原始地形。本研究中，采用普通Kriging方法进行插值，具体的原理与插值步骤参见文献 ［2］，具体计算过程在ArcGIS软件中完成。

1．3．2 多重分形插值法 多重分形 (Multifractal)是指用多个维数来描述非均匀的复杂几何体，以此来全面刻画其特征。在地形变化复杂的地区或者相对较大的区域，高程异常的变化不规则，用规则的曲面取代实际上不规则的似大地水准面，必然会导致大地水准面的不规则变化，导致拟合后高程异常的精度降低、误差较大。多重分形可通过高程值计算出高程局部异常的奇异性指数，通过奇异性指数得到多重分形Kriging插值的结果。

局部奇异性分析方法实际上是将样点高程值段的密度在分形空间中进行度量，以确定分形密度和分形维数，奇异性所度量的是场值随量度范围大小的变化规律。奇异性指数 (Z)表示基于对场值的某种滑动加权平均，计算公式如下:

其中:Ω(x0，ε)为围绕中心点x0半径ε的小滑动窗口;ω(‖x0-x‖)为对在Ω (x0，ε)中与中心点x0相隔‖x0-x‖距离的任意点x的加权函数，它往往与距离呈负相关。加权函数的选择不仅与距离有关，与场的空间性自相关，还与处理目的有关。

成秋明［8］给出的多重分形方法将滑动平均关系表达为

其中，α(x0)为x0点处的局部奇异性指数。可以看出，式 (2)中不仅包含了空间相关性的成分，还有度量奇异性指数。如果α(x0)=2时，通过该方法所计算的加权平均值与通常的加权平均无异;但当处于高程值差异较大地段而且局部地区具有奇异性，即α(x0)＜2时，通过该方法所得的结果将高于通常的加权平均结果;相反，当处于高程值差异较低地段，即α(x0)＞2时，通过该方法所得的结果将低于通常的加权平均结果。由此可见，该方法有利于加强峰谷值，对于具有奇异性的空间结构来说，传统的插值方法不能很好地估计出峰谷值。指数α对应分形空间维数，用分形空间维数与正常的欧氏空间维数的差Δα=2-α即可表示分形密度与正常密度的空间维数的差异。

2 结果与分析

2.1 普通Kriging法的插值结果

利用探索性数据分析工具得知插值样点集(图1中黑色实心点)偏态值为-0．157，可认为样点呈近正态分布。再对数据进行试验半方差函数的计算和球面模型拟合 (步长取3500 m，组数为12组)，得到研究区域的海底地形插值结果。从普通Kriging法的插值结果 (图2)可以看出，海州湾水下地形起伏不大，呈明显条带状分布。提取栅格中检查点位置的水深值，并与检查点的初始值进行线性拟合 (图3)。结果表明:水深为-15～0 m范围内的点比较集中，说明插值结果与实际值比较接近;-25～-15 m范围内的插值结果与实际值有一定偏离，但插值结果与实际值均分布在斜率k=1的直线附近，表明总体估计较为准确。

2.2 多重分形Kriging法的插值结果

图1 样本点分布示意图Fig.1 The distribution of sampling points

图2 普通Kriging法的插值结果Fig.2 Map of interpolation by Kriging method

图3 基于普通Kriging法插值的交叉验证图Fig.3 Cross－examination by Kriging method

在普通Kriging插值方法的基础上，为准确反映海州湾海底地形，度量水深分布的局部奇异性，引入多重分形理论。首先采用不同格网尺寸等级，滑动计算多重分形测度范围内高程分布所包围形体的体积，计算格网虚拟局部奇异性指数;然后回归计算双对数曲线的曲率，得到奇异性指数［即分形维数α(x0)］。具体计算过程在Matlab软件中完成。当曲率大于指定的阈值 (在此设为0．9)时，双对数关系成立;否则，取α(x0)=2。结合海州湾海岸线，得到了基于多重分形理论的Kriging插值结果。该插值结果 (图4)与普通Kriging插值结果趋势相同，都呈明显的条带状分布。将检查点位置的水深值与检查点的初始值进行线性拟合得到交叉验证图 (图5)。

图4 多重分形Kriging法的插值结果Fig.4 Map of interpolation by multifractal Kriging method

图5 多重分形Kriging法插值的交叉验证图Fig.5 Cross－examination of multifractal Kriging method

2.3 两种插值方法的比较

将两种插值方法的结果与验证数据集进行对比。从表1可以看出:实际值中最大值为-0．10 m，最小值为-29．50 m;利用普通Kriging法插值后的最大值为-0．09 m，最小值为-31．16 m，可见，其对数据集的拉伸较为明显;基于多重分形理论的Kriging法插值后的最大值为-0．10 m，最小值为 -3 0．0 2 m，该方法对数据集的拉伸较小。从标准差与均方根误差上看，基于多重分形的Kriging方法均比普通Kriging方法小，说明其得到的插值结果更为准确。

表1 两种插值方法的结果比较Tab.1 Comparison of two methods

另外，对比图2与图4的插值结果，发现图2较图4平滑，没有保留异常点的水深信息。由于后者弥补了半变异函数局部平滑的缺陷，其能够有效地勾勒局部细节，度量地表高程局部隆起和下陷的局部奇异特征，更好地反映了海州湾的局部地形。

3 海州湾水下地形的三维可视化与分析

为了更形象地刻画海州湾海底的地形，利用Phong简单光照模型，基于多重分形理论的Kriging方法的插值结果，通过调整光照 (设定光照方向为水平向135°，垂直向为45°)、颜色与选择观察方向，参考自然地物可视化的相关研究［18-22］，利用IDL语言开发绘制了海州湾水下三维地形图(图6)，从宏观上较好地展示了海州湾的整体水下地形;除此之外，海湾里的几个暗礁 (图6中紫色矩形区)也较好地呈现了出来。

由于古黄河的南徙，其携带的大量泥沙沉积不断使陆地向前推进，最终造成原来的孤岛云台山(图6中A点)与大陆相连［14，23］，从而形成了从A到B的“大喇叭口”［24］形开敞海湾。海岸地貌类型主要有海蚀地貌和海积地貌。如图6所示，海域宽阔但水深较浅，低潮线以下多为水下暗坡，地势向东北方向倾斜，平均比降在0．37‰。在湾口附近和口外浅海，水深为10～27 m，是一个起伏和缓的冲刷面，地势向东偏北方向倾斜。利用2012年获取的研究区域THEOS卫星融合影像 (图7)，对比图6中暗礁的坐标，确定其对应着湾口东北部的基岩型岛屿:达山岛、车牛山岛和牛背岛，周围还存在着达东礁和牛嘴礁等。虽然采集的水深数据无法模拟出海平面以上的岛屿，但通过多重分形Kriging方法准确地拟合出了海底暗礁的位置。

4 结论

利用研究区域的水深数据点，分别使用普通Kriging插值法和基于多重分形的Kriging插值方法进行了海州湾的水下地形构建与三维可视化。研究表明:(1)普通Kriging插值法通过半变异函数从预测点周围的观测值中生成权系数进行预测，可以较好地反映整个区域内的数据空间特征，但插值结果较平滑，丢失了局部信息;(2)多重分形理论由于采用空间自相似性过滤，改进了普通Kriging方法中半变异函数局部平滑的缺点，考虑了空间数据分布的局部奇异性，可以获得更好的空间插值效果。对于研究区域海州湾，采用多重分形Kriging插值方法在宏观上较细腻地反映了平坦的海底地形，在微观上也展示了达山岛和车牛山岛等周围的暗礁。证明该方法可广泛用于模拟复杂的几何形体，尤其适用于复杂粗糙地形的重建和对数据异常的敏感度量场的重建。

图6 研究区地形模拟结果Fig.6 Terrain simulation of the observed area

图7 海州湾THEOS融合影像Fig.7 THEOS fusion image of Haizhou Gulf

［1］申静，苏天赟，王国宇，等．基于Kriging算法的海底地形插值设计与实现［J］．海洋科学，2012，36(5):24-28．

［2］包世泰，廖衍旋，胡月明，等．基于 Kriging的地形高程插值［J］．地理与地理信息科学，2007，23(2):28-32．

［3］张捷，包浩生．分形理论及其在地貌学中的应用——分形地貌学研究综述及展望［J］．地理研究，1994，13(3):104-112．

［4］高义，苏奋振，周成虎，等．基于分形的中国大陆海岸线尺度效应研究［J］．地理学报，2011，66(3):331-339．

［5］金德生，陈浩，郭庆伍．河道纵剖面分形——分线性形态特征［J］．地理学报，1997，52(2):154-162．

［6］熊波，王建力，张天文．渝东南岩溶山区耕地利用变化下土壤颗粒体积分形特征研究［J］．中国岩溶，2011，30(3):295-301．

［7］刘鹏举，赵仁亮，朱金兆，等．保持地貌特征的数字高程模型生成方法研究［J］．中国矿业大学学报，2006，35(4):523-526．

［8］成秋明．多重分形与地质统计学方法用于勘查地球化学异常空间结构和奇异性分析［J］．地球科学，2001，26(2):161-166．

［9］Lee C K．Multifractal characteristics in air pollutant concentration time series［J］．Water，Air and Soil Pollution，2002，135(1/4):389-409．

［10］李锰，朱令人，龙海英．不同类型地貌的各向异性分形与多重分形特征研究［J］．地球学报，2003，24(3):237-242．

［11］曹汉强，朱光喜，李旭涛，等．多重分形及其在地形特征分析中的应用［J］．北京航空航天大学学报，2004，30(12):1182-1185．

［12］崔灵周，李占斌，郭彦彪，等．基于分形信息维数的流域地貌形态与侵蚀产沙关系［J］．土壤学报，2007，3(2):197-203．

［13］Wang D，Fu B J，Lu K S，et al．Multifractal analysis of landuse pattern in space and time:a case study in the Loess Plateau of China［J］．Ecological Complexity，2010(7):487-493．

［14］中国海湾志编纂委员会．中国海湾志．第四分册［M］．北京:海洋出版社，1993:354-380．

［15］左书华，庞启秀，杨华，等．海州湾海域悬沙分布特征及运动规律分析［J］．山东科技大学学报:自然科学版，2013，32(1):10-17．

［16］王宝灿，虞志英，刘苍字，等．海州湾岸滩演变过程和泥沙流动向［J］．海洋学报，1980，2(1):79-96．

［17］张存勇，冯秀丽，陈斌林．海州湾南部近岸柱状样沉积物重金属污染［J］．海洋地质与第四纪地质，2008，28(5):37-43．

［18］薛安，马蔼乃，李天宏．基于OpenGL实现真实感地形表现的研究［J］．中国图象图形学报，2001，6(2):800-805．

［19］孙博文．分形算法与程序设计——Delphi实现［M］．北京:科学出版社，2005:257-331．

［20］王梦，金文标．基于函数迭代系统的3-D分形插值算法［J］．计算机应用，2006，26(11):2701-2703．

［21］黄天云，张传武．分形插值算法在分形自然景物模拟中的应用［J］．计算机工程与设计，2007，28(16):3394-3397．

［22］张涛，徐晓苏，王其，等．基于分形插值的三维海底地图生成算法［J］．中国惯性科学学报，2008，16(2):171-173．

［23］黄志强．全新世云台山的海陆变迁［J］．徐州师范学院学报:自然科学版，1992，10(2):28-33．

［24］耿秀山．黄渤海地貌特征及形成因素探讨［J］．地理学报，1981，36(4):423-434．