张 玲 杨晓平 魏占玉 黄伟亮
(中国地震局地质研究所,活动构造与火山重点实验室,北京 100029)
数字高程模型(Digital Elevation Model,DEM)作为地形定量描述的载体,是地貌学定量研究的主要对象,自20世纪50年代后期被首次提出以来,就受到了科学界和工程界极为广泛的关注。航空摄影测量一直是地形图测绘和更新的有效手段,其所获取的影像数据是高精度大范围的DEM生产最有价值的数据源(周启鸣等,2006)。航天科技、计算机技术等蓬勃发展以及遥感影像分辨率的不断提高,尤其是近年来商业航空LiDAR系统的最高分辨率已经达到了每个激光点0.3~0.5m(Tatsuro et al.,2008),这些软、硬件的改进都使衍生的DEM数据质量不断提高。怎样挖掘隐含于高分辨率三维数据中的信息并将它们应用于科研和生产实践中至关重要,这也是地形可视化所要解决的主要问题。在二维平面上,读图者并不能调整观看视角,因此研究人员通过不断地努力将描述地形的各种参数叠加在一起试图增强图形的可读性。传统的方法,如等高线图、分层设色图、地貌晕渲图和地貌参数叠加图等,已经得到了相当广泛的应用。下面首先概述它们的优点和局限性,然后介绍一种新的地形参数(openness),以及使用它得到的三维数据二维可视化方法(RRIM)在地貌研究中的应用。
根据实际用途的不同,研究人员发展了各种三维数据的二维可视化方法,最常用的是地形等高线图(图1a)。在一定比例尺的地形图上,用等高线的密度来表现地形坡度的大小,密度越大,地形坡度越大;反之,坡度越小。在地形图上面状地貌的判定特征为具有一定面积连续区域。其边缘等高线密集且近于闭合,指示边坡坡度陡倾;内部等高线稀疏,指示内部坡度较缓。而线状地貌可以被理解为其在一个方向上的长度远大于在其他方向上的长度。因此在等高线图上,若线状地貌的延伸方向与等高线垂直或斜交,等高线表现为顺地貌延伸方向上的连续突起。若延伸方向与等高线平行,则表现为等高线的密集分布。通过这种方法,可以大致确定诸如冲洪积扇、河流阶地等面状构造的边界,以及诸如山脊、冲沟等线状构造的延伸方向。地形的高程、坡度、坡向等地貌参数能够在地形图上通过量测计算得出。等高线图最大的缺点是立体感不强,因此有人提出明暗等高线图(图1b)(Tanaka,1950)与阶梯等高线图的概念(周启鸣等,2006)。另外,随着等高线密度的增加,已无法分辨出每一条等高线。因此,当三维数据分辨率不断提高,等高线图就无法满足二维显示的要求。
图1 Sheep Mountain背斜部分地段等高线图(a)与明暗等高线图(b)Fig.1 The contour map(a)and the illuminated contour map(b)of part region of Sheep Mountain anticline.DEM是1m分辨率的LiDAR数据,数据来源于NSF Open Topography Facility,
晕渲图是通过模拟太阳光对地面照射所产生的明暗程度,并用相应灰度色调或彩色得到随光度近似连续变化的色调。这种方法的最大优点就是符合人眼观察实际地形的习惯,最接近真实地模拟了地形面。由于地面的凹凸起伏和阴影取决于太阳光的入射角和太阳高度角,在常规制图中常常根据人眼的视觉习惯,太阳方位角选择为NW方向,太阳高度角为45°。图2中,太阳高度角均为45°,太阳光入射方向分别为黄色箭头方向:NW、NE、SW、SE 4个方向。不难看出,调整太阳光的入射方向可以得到完全不同的凹凸感。而且沿太阳光入射方向由于较大的地形高差会产生面积不同的阴影,这些阴影可能会遮盖相对较小的地貌现象。因此,在实际地貌研究中,需要不断地调整入射光的方向和太阳高度角来突显不同位置、不同大小的地貌现象。这个过程的实现就要求作图者具备一定的制图经验。
图2 Sheep Mountain背斜部分地段晕渲图Fig.2 Hill shading map of part region of Sheep Mountain anticline.
高程分层设色图(图3)是应用不同的颜色或者灰度级来表示不同的高度带,一般包括基于常规的高程分带设色和基于高程数据的灰度影像(半色调符号表示法)(周启鸣等,2006)。三维数据中的高程等信息隐含在颜色的渐变中。由于人眼对颜色过渡的分辨能力并不高,地形的轮廓往往不清晰。因此,这种方法不太适用于大比例尺图和地势平坦的地区。但是,这种方法与其他可视化方法叠加,则可改善图像的显示效果。如它与晕渲图(图4)、坡度图的叠加,都增强了影像的直观性和立体感。分层设色图与坡度图的叠加很好地显示了边界信息(图5),可以用来区分不同的地貌单元。但对于地形的细微构造表现得不够好,因此分层设色图更适用于表达大尺度的区域。
图3 Sheep Mountain背斜部分地段分层设色图Fig.3 Hypsometric tint map of part region of Sheep Mountain anticline.
图4 Sheep Mountain背斜部分地段彩色晕渲图Fig.4 Colored hill shading map of part region of Sheep Mountain anticline.
RRIM(Red Relief Image Map)是由日本学者Tatsuro(2008)首次提出的,旨在在二维平面中突显地貌上的微小构造特征,同时利用彩色和灰度的渐变来表现立体感。RRIM是由包含地形坡度、正openness和负openness 3个地形参数图叠加融合的平面图。
RRIM图中的核心参数为正、负openness,是由日本学者Yokoyama等(2002)提出的。这种方法是通过计算一定步长L内的最大天顶角(zenith)和天底角(nadir)来表现不规则表面的凹凸程度(图6)。
图5 Sheep Mountain背斜部分地段坡度图与分层设色图的融合Fig.5 Fusion of slope map and hypsometric tint map of part region of Sheep Mountain anticline.
openness的计算方法是基于格网的DEM。首先,在一个栅格点上,以L为半径,寻找最大高程角DβL与最小高程角DδL,以水平面为角度的零点,向上为正,向下为负,然后计算得到天顶角与天底角(图6)。其中,天顶角:DφL=90-DβL;天底角:DψL=90+DδL。
扩展到三维计算时,把目标栅格点与邻近的栅格点用直线连接,将分析区域分割为8条方向线和8个扇区(图7)。以半径L为步长,分别计算8个方向上目标点的最大天顶角与天底角。则正openness为8个方向上最大天顶角的平均值,即:φL=(0φL+45φL+…+315φL)/8;负openness为8个方向上最大天底角的平均值,即:ψL=(0ψL+45ψL+…+315ψL)/8。
图6 最大高程角、最小高程角、天顶角和天底角的定义(据Yokoyama et al.,2002)Fig.6 Definition of the maximum elevation angle,the minimum elevation angle,zenith angle and nadir angle(after Yokoyama et al.,2002).
图7 目标栅格点计算区域划分示意图(据 Yokoyama et al.,2002 改绘)Fig.7 The schematic drawing of regionalization of the grid point of interest(adapted after Yokoyama et al.,2002).
正、负openness与其他地貌参数最大的区别是分别突显了正、负地形。图8中,A1、A2和A3点分别位于山脊、平地和山谷。图8a,d,g分别为三者0°~180°方向的地形剖面,即对于A1点来说,为沿山脊线方向的地形剖面;对于A2点来说,可选为任意方向的地形剖面;对于A3点来说,为沿山谷线方向的地形剖面。图8b,e,h分别为三者90°~270°方向用来计算正openness的地形剖面,图8c,f,i分别为三者90°~270°方向用来计算负openness的地形剖面。由于90°~270°方向与0°~180°方向垂直,因此,90°~270°方向的地形剖面对于 A1点来说,为垂直于山脊线方向的地形剖面;对于A2点来说,为与0°~180°方向垂直的地形剖面;对于A3点来说,为垂直于山谷线方向的地形剖面。
在0°~180°方向的地形剖面中,3种地形的天顶角的平均值都在0°~90°之间,相差不大。但在90°~270°方向的地形剖面中,山脊两侧为向下的陡坡,因此天顶角DφL>90°,天底角DψL<90°;平地两侧地形平缓,因此,天顶角 0°<DφL<90°,天底角 0°<DψL<90°;山谷两侧为向上的陡坡,因此天顶角DφL<90°,天底角DψL>90°。将8个方向的天顶角与天底角进行求和平均计算求得正、负openness。在理想情况(山脊和山谷两侧地形足够陡峭)下,能够得出:0°<山谷 φL<平地 φL<90°<山脊 φL;0°<山脊 ψL<平地 ψL<90°<山谷 ψL。在实际地形中,地面凸起点在0°~360°方向上天顶角均>90°,天底角均<90°;地面下凹点在0°~360°方向上天顶角也均<90°,天底角均>90°,故能得出:0°<下凹点 φL<平地 φL<90°<凸起点 φL;0°<凸起点ψL<平地 ψL<90°<下凹点 ψL。
图8 山脊、平地和山谷3个地点不同方向剖面正负openness计算示意图(据Yokoyama et al.,2002改绘)Fig.8 The schematic drawing of calculation of positive and negative openness of hilltop,plain and depression along profiles in different directions(adapted after Yokoyama et al.,2002).
因此,正openness代表凸状表面,值越大,凸度越大,最高值代表最高点;而负openness代表凹状表面,值越大,凹度越大,最大值代表最低点。正、负openness的大小受步长半径L的影响,步长半径越大,愈突显大的构造;步长半径L越小,愈突显小的构造。这样就可以根据研究尺度的不同而设定不同的步长半径L。
openness本身具有与光源无关的特性,因此利用其制图就能够去除晕渲图等可视化方法中地形阴影的制约因素,而且计算结果受DEM噪音干扰的影响要比其他参数小得多(Yokoyama et al.,2002)。
可视化方法的进步在一定程度上依赖于地形参数的创新。openness不只是应用于地形可视化方面。Prima等(2006)通过计算日本本州北部地区DEM数据的地形坡度和openness等4个地形参数,实现了对该区域地形的监督分类,得出大区域的地形阶梯状构造与火山分布特征,指示了本州北部深部的岩浆管道系统是由上地幔上涌至上地壳。Peter(2009)将openness参数应用于可视域分析中,认为和地形坡度、断面曲率(cross sectional curvature)、最大曲率(maximum curvature)、最小曲率(minimum curvature)、平面凸度(plan convexity)、剖面凸度相比(profile convexity),正openness能给出最佳的可视域预测。由于破碎度很大程度上是由自然界的地表特征决定的(周启鸣等,2006),因此若将可视域分析扩展计算区域的破碎度,可判断可视区域的连通程度。
Tatsuro等(2008)规定了在RRIM中2个openness的计算遵从以下规则:
其中:Op是正openness,On是负openness。Tatsuro的openness计算规则进一步增强了单一图幅上对凹凸特性的双重显示。在RRIM图中openness参数的计算结果与地形坡度分别用灰度影像层、红色图层来表示(图9)。经验表明,红色对于人的眼睛来说有最丰富的色调,尤其是电脑下的色彩空间,因此能够最大限度地表现地貌信息。RRIM不仅仅适用于LiDAR数据,它还可以广泛地应用于其他各种三维地形数据,如SRTM、GTOPO30和ETOPO2(Tatsuro et al.,2008)。
虽然RRIM图能有效地在二维平面上显示三维数据的立体感,但它缺少一定的定量化的信息,如高程、坡度、坡向等。一种改进方式就是将RRIM图与其他地形参数图层进行叠加融合(Tatsuro et al.,2008),如将它与等高线图进行融合就增加了高程信息(图10),通过计算测量就能得到坡度和坡向。
图9 Sheep Mountain背斜部分地段RRIMFig.9 RRIM of part region of Sheep Mountain anticline.
图10 Sheep Mountain背斜部分地段RRIM与等高线融合图Fig.10 Fusion of RRIM and contour map of part region of Sheep Mountain anticline.
可以轻易地从RRIM图(图9)中直观地解译出背斜区总体的形态展布、水系分布情况以及背斜区由构造运动与水流共同作用造成的地表粗糙度。图11b,c,d,f,g,h,i表示了应用各种可视化方法对于相同小范围区域进行的地貌对比。通过对比可知,RRIM在细微构造的识别上有其他可视化方法无法比拟的优越性。图11h中显示的小冲沟和小路是在其他可视化方法图中无法分辨出来的。同时,RRIM图中还凸显了面状构造的边界信息,因此能够准确地勾勒出面状构造的边界,如图中的沟谷边界。对于1m分辨率的DEM实现了对m级构造信息的识别。在实际应用中,对于活动构造研究中需要重点分析的各种地貌识别标志,如线性山谷、断错水系、断层陡坎、断层三角面、河流阶地、冲洪积扇等,在RRIM图中可以很容易地进行判读。它们的解译是确定断层系统活动性质,测量断层滑移量,分析各构造单元相互关系,进一步计算地震复发周期等研究第四纪地表变形工作的基础,也是野外踏勘的前期准备工作。加之,三维数据本身具有可量测性,因此研究人员不必亲临地质构造现场就能够完成对目标现象的计算和分析。因此,对于那些由于自然条件恶劣、不易到达或者大尺度的研究区域,这种可视化方法具有很强的实用性,能在一定程度上填补一些未能开展工作的区域的研究空白。同时,也减少了野外工作量、提高了工作效率。Lin等(2012)在对日本中部山脉的Neodani断层开展机载LiDAR勘测中就应用了RRIM图来进行微小构造地貌的识别。目标研究区位于断层的北段,高程范围为156.2~1540.9m,平均坡度角35°左右。虽然主断层与周围伴生断层的活动模式一直受到特别的关注,但是由于断层展布区没有道路,加之浓密的植被覆盖和地形陡峭等因素的制约,因此一直都没有得到良好的航片解译和野外勘察结果。通过将RRIM图与分辨率为0.5m的LiDAR数据联合应用,发现了大量的微小构造现象,如植被覆盖下的1~6m高的断层陡坎、被整体错断的河道以及错断河道两侧基岩上的新鲜微小破裂。这些构造现象是无法在1/10000的航片、以前的野外调查和之前的2m分辨率LiDAR生成的DEM中解译出来的。这些微小构造现象的发现与解译对于该地区未来可能发生地震的发震时间和震级大小的评估都是至关重要的。如果有些目标区域因森林覆盖而使便携式GPS设备无法正常使用,而RRIM图却能发挥很好的定位功能。
图11 用不同可视化方法所得效果的对比Fig.11 The contrast of different visualization methods.
通过以上对各种可视化方法的简要介绍与对比,可知传统的可视化方法已经不能完全满足对分辨率不断提高的三维数据信息的判读。在二维平面上实现信息增强和要素提取一直是可视化方法不断改进的方向。等高线图、晕渲图和分层设色图都有各自的优点,可以相互补充,并通过叠加融合以达到信息增强的目的。新的可视化方法(RRIM)在很大程度上增强了三维数据二维显示时的立体感,去除了光源条件的制约,最重要的是可以根据研究尺度来突显各种细微构造,其直观性是传统方法无法比拟的。随着三维数据分辨率和精度的不断提高,良好的地貌解译工作不仅能够提高工作效率而且能够在一定程度上减少野外工作量。对于有研究价值而又植被覆盖严重、自然条件恶劣、不易到达的区域提供了一种更好的工作方法。如果未来三维数据信息的分辨率和精度足够高,能够和野外手持测量设备获得的测量精度相当,则有望实现用室内解译工作在一定程度上代替传统的野外工作。因此,RRIM图在活动构造研究、地质地貌填图以及地震危险性评价等工作中有着非常重要的实用价值。通过地形参数的创新可以从根本上改进可视化方法、拓宽研究区范围和研究尺度,期望在不久的将来openness参数和RRIM能够被应用到活动构造研究和其他更广泛的领域中。
周启鸣,刘学军.2006.数字地形分析[M].北京:科学出版社.
ZHOU Qi-ming,LIU Xue-jun.2006.Digital Terrain Analysis[M].Science Press,Beijing(in Chinese).
Lin Z,Heitaro K,et al.2013.Detection of subtle tectonic-geomorphic features in densely forested mountains by very high-resolution airborne LiDAR survey[J].Geomorphology,182:104—115.
Peter L G.2009.Using upward openness for viewshed prediction[A].In:Proceedings of ASPRS 2009 Annual Conference,Baltimore,Maryland.
Prima Oda,Echigo A,et al.2006.Supervised landform classification of Northeast Honshu from DEM-derived thematic maps[J].Geomorphology,78(3-4):373—386.
Tanaka K.1950.The relief contour method of representing topography on maps[J].Geographical Review,40:444.—456.
Tatsuro C,Shin-ichi K,et al.2008.Red relief image map:New visualization method for three dimensional data[A].In:The International Archives of the Photogrammetry,Remote Sensing and Spatial Information Sciences,Vol.XXXⅦ.Part B2.
Yokoyama R,Shirasawa M,Richard J P.2002.Visualizing topography by openness:A new application of image processing to digital elevation models[J].Photogrammetric Engineering and Remote Sensing,68(3):257—265.