海面-地面系统中三维古地面反演计算的剥蚀算法研究

王 琳,钟鹤翔,谢志仁

(1.中国科学院南京地理与湖泊研究所,江苏南京 210008;2.中国科学院研究生院,北京 100039; 3.江苏省信息中心,江苏南京 210013;4.南京师范大学地理科学学院,江苏南京 210097）

海面-地面系统中三维古地面反演计算的剥蚀算法研究

王 琳1,2,钟鹤翔3,谢志仁4

(1.中国科学院南京地理与湖泊研究所,江苏南京 210008;2.中国科学院研究生院,北京 100039; 3.江苏省信息中心,江苏南京 210013;4.南京师范大学地理科学学院,江苏南京 210097）

三维古地面反演计算是海面-地面系统研究中的一个基础研究。该文在前人研究的基础上提出了三维古地面反演计算的剥蚀算法,并与传统的加积算法做了比较分析,得出剥蚀算法比加积算法有对时间序列离散判断简单、方便计算等优点。同时,在GIS平台下实现了该算法,并以长江三角洲地区为试验区,进行了古地面虚拟反演计算,取得了较好的效果。

地理信息系统;海面-地面系统;古地面反演计算

地理环境的表达和认知、古地理环境的恢复和再现,始终是国内外地学、地图学、认知科学领域普遍关注的热点;研究并建立基于 GIS技术和虚拟环境技术的沿海地区古地貌演变模型是海面-地面系统变化研究的一个重要方向[1-3]。三维古地面反演计算是古地理环境恢复和再现的重要研究内容之一,也是海面-地面系统研究中的基础内容。海面-地面系统研究就是将海面升降变化和地面升降变化作为一个统一的整体加以考察,通过对保存在地面外部形态特征和内部地层构成中有关信息的提取,用来演绎海面和地面变化的历史过程、当前动态和未来趋势[4]。在海面-地面系统中,如何进行古海面-地面变化的反演以及未来海面-地面变化的预演是整个研究的核心。对于海面-地面相对运动的模拟计算、古地面高程的反演计算、古地形模拟已有一些研究成果[5-13]。本文在前人研究成果的基础上,设计出海面-地面系统中三维古地面的剥蚀算法,并在GIS平台下实现了该算法;然后选择长江三角洲地区进行三维古地面反演模拟,取得了较好的效果。

1 三维古地面反演计算的剥蚀算法

1.1 剥蚀算法

传统的三维古地面反演计算模型是以沉积物随时间的推移在沉积基底之上逐渐加积从而导致地面升高为主线,同时考虑到海面的升降、地壳升降运动与沉积物压实沉降等地球物理作用而设计的概念模型,在此基础上通过试验设计出三维古海面计算模型。该算法在时间尺度离散上存在着时间断点判断复杂、计算困难等缺点。鉴于此,通过逆向思维方式提出了三维古地面计算的剥蚀算法。算法的基本思路如下:现代地形中空间上离散点的高程值,都是由地壳运动、沉积作用及沉积压实作用等在一定的地质时间内共同作用的结果。因此从现代地形中层层剥离掉这些作用造成的地面增高,便可以反演古地形的高度,同时与当时的古海面高程进行基准面转换,从而得到各离散点的古高程值。

对于单一离散空间点的古高程计算式如下:

式中:O ld ET为历史反演目标时间 T时期的古高程值;L an E为现代的地形高度;U pperSedDT为历史反演目标时间T到现在这一连续时段内的沉积物厚度;U pperConSed T为历史反演目标时间 T到现在这一连续时间长度;TecV、EusV、Com V分别为历史反演目标时间T到现在这一连续时段的平均地壳升降运动速率、地壳均衡运动速率和沉积压实作用速率;O ldSLT为历史反演目标时间 T时期的古海面高程,用来进行古海面基准转换计算。其计算流程如图1所示。

计算式(1）中用到的各数据,假设将10 000 aB.P.至今分为 N个时段,并给出每个时间分段的时间点。具体计算流程如下。

(1）计算各时段的沉积持续时间:ConSed T[i]。每个地面离散点的沉积持续时间首先取决于它开始接受沉积的第一时段的位置。在精细的古地貌模拟中,不能将地层界面视为等时面。例如,沿海地带不同地点的全新世沉积起始时间应该是一个与古地面高程和遭受全新世海侵相关的变量,不能简单地一律设定为10 000 aB.P.,而必须通过一定的判别计算才能确定。基本计算公式为:

图1 古海面-地面变化中三维古地面反演计算流程Fig.1 The flow chart of three-dimensional paleosurface inversion algorithm in sea level-land surface changes

式中:i为第i个时段;ConSed T[i]为在第i个时段内的沉积延续时间;Sub T[i]为时间段的分段时间点;StartT为该离散点上开始接受沉积的起始时间。3个公式分别表示3种情况下ConSed T[i]的值。计算时先判断StartT位于哪个时间分段内,假设其位于第i个分段内,则第i个分段以前的各段没有沉积,沉积持续时间为0,即:ConSed T[i-1]=0;第i个分段内的沉积持续时间为第i个时段的分段时间点与沉积开始时间之差,即:Start T-Sub T[i];第i个分段后的所有分段的沉积持续时间为相邻两个分段时间点的差,即:Sub T[i]-Sub[i+1]。

(2）计算不同时段内的沉积速率:SedV[i]。由于整个沉积时期的不同时段内的沉积速度不同,要想达到最佳的历史时期古地貌的拟合效果,就必须把整个沉积厚度离散分段,根据专家知识对不同时期赋以不同的比值,这样在计算中不同的时期就会有沉积速率的差异。SedV[i]的计算公式为:

式中:SedV[i]为第i个时间分段内的沉积速率; Sed T为整个模拟时段内形成的沉积厚度;ConSed T [1]…ConSed T[N]为每个分段的沉积延续时间;k [i]为沉积速率分量控制参数。

(3）计算沉积厚度(即每个时段的沉积速度与沉积持续时间积之和）:U pperSedD T。公式如下:

式中:SedV[i]T为所求目标时间T(或沉积开始时间StartT）所在时间分段的沉积速度;ConSed T[i]T为所求目标时间T(或沉积开始时间Start T）所在时间分段的沉积持续时间。

(4）计算升降校正:Rec H T。公式为:

式中:TecV、EusV和Com V分别是该离散点的垂直方向上构造运动速率、均衡运动速率和压实作用速率;U pperConSed T则为所求目标时间到目前的沉积持续时间。

(5）计算古地面高程。基本公式为:

式中:O ld ET为历史反演目标时间 T以现代高程系为基准的古地面高程;Land E为该地面离散点的现代地面高程;U pperSedDT为历史反演目标时间 T到现在这一连续时段内的沉积物厚度;R ec HT为升降校正高度。

(6）古海面基准转换计算。基本公式为:

式中:OldSLT为基于现代高程系的古海面高程,当古海面值低于现代海面时取古海面高程为负值,否则为正值;SLO ld ET为基于古海面基准高程的古地面高程值,该值体现的是反演目标时间的海陆关系。

通过以上计算可以得到三维古地面的高程,然后通过检验,对验证结果进行判断,基本符合验证数据的结果便是所求结果,对于不符合的则重新进行计算。

1.2 算法模型所需数据分析

从算法模型的古高程计算方程可以看出,该算法所需的原始数据主要有4类:沉积启动时间数据、沉积厚度数据、沉积校正数据和古海面高程数据。其中沉积校正数据又包含:地壳构造运动、地壳均衡运动和沉积压实作用产生的沉积校正数据。在整个反演算法计算过程中,通过对原始数据进行自动半自动的时空离散来达到算法的自适应计算模拟,并通过对模拟结果的验证逐步向真实的古海面-地面地形逐步逼近,以得到最佳反演结果。因此,算法中除原始数据外,还有时间离散控制参数和空间离散控制参数。

(1）沉积启动时间数据。沉积启动时间数据是控制反演区域沉积开始时间的数据,它是通过空间离散控制参数把整个反演区域离散为栅格点沉积起始时间。在晚更新世末期的低海面时期,现在的三角洲平原区是一个大面积处于水面以上的受侵蚀的丘陵岗地地貌,起伏的地形决定了遭受海水淹没并开始接受沉积的时间会有很大差别。因此每个栅格点开始接受沉积的起始时间是不同的,该值可以通过滑动底界模型求得[4]。

(2）沉积厚度数据。对于三角洲近海及平原区域,第四纪以来的沉积物加积作用是造成古地面变化的最为重要的因素,因此沉积厚度数据决定了整个反演计算结果的精度,也决定了时空离散参数的设定。因为整个算法都是建立在对沉积厚度数据数理分析的基础上,只有通过对沉积厚度的分析,才能获取相对精确的时间分段并计算出各个分段的平均沉积速率。

(3）沉积校正数据。沉积校正数据是用来校正地壳构造运动、地壳均衡运动和沉积压实作用对海面-地面高程的影响。但3种运动计算比较困难,在实际的试验中可以由简单模式确定,如平行状等速率线所代表的向海掀斜运动模式;地壳均衡运动、沉积压实作用与沉积厚度呈线性关系。算法在实现时留有数据接口和线性关系式参数接口。

(4）古海面高程控制。古海面高程是反演古海面-地面系统运动演变的重要控制参数,它是相对古地面的一个重要的高度参考数据。在实际试验中可以对已知离散时间上的古海面高度进行内插从而求得目标时刻的古海面高度,以此作为相对古地面的基本海面标高。

(5）空间离散控制参数。通过对原始数据(包括现代地形数据、沉积厚度数据等）进行栅格化来对整个研究区进行空间离散。在整个研究区,通过数据分析,将沉积过程基本一致、数据采集良好的区域分割成一个个单一试验区进行试验;然后对每个试验区进行栅格化,栅格的大小可以根据试验的要求和原始数据的精度来确定。

(6）时间离散控制参数。在整个第四纪地质时期的沉积过程中不同时段的沉积速度有差异,但不可能把每个时段的沉积速度都求出。可以根据专家经验,在对原始的沉积厚度、钻井资料分析的基础上,把试验区的整个沉积过程大致分成几个相应的沉积间段,在每个沉积间段内假设有着同样的沉积速率。不同时段间的沉积速率有个比值,在模型计算过程中可以调整此比值来进行不同沉积间段的沉积速度控制,在算法中可以通过时间分段技术实现整个计算过程的自动、半自动反演计算。

从数据分析可以看出,沉积厚度数据是影响反演计算结果的最重要的数据。而对研究区进行时空离散,通过时间分段技术进行的时间维度的离散,计算不同时段内的平均沉积速率,并可以通过模型对每个时段内的沉积速率自动调节,达到模型自我不断反演逼近真实地形的算法思想与数据无关。

1.3 与传统算法比较分析

在传统的算法中以沉积基底上的沉积物加积作用为基本计算依据,同时考虑到海面的升降、地壳运动与沉积压实等地球物理作用等来计算古地面,姑且把该算法称为“加积算法”,其基本方程如下:

式中:Lan E-Sed T项的地学意义是现代的地形(Lan E）剥除模拟时段内形成的沉积总厚度(Sed T）之后,以现代高程基面计算的古地形高程,即模拟的起始地貌面;Low erConSed T为沉积持续时间:TStartT;Low erSedDT为从沉积起始时间到特定反演目标时间T时段内的沉积厚度,其计算公式为:

式中:Low erConSed T[i]T为目标时间T所在的第i个时间分段到第i-1个时间分段的沉积延续时间。

从式(10）可以看出,剥蚀算法与加积算法的主要差异为所求的沉积物厚度位于目标时间T之前还是之后(图2）。

图2 两种算法所求沉积物厚度的位置差异Fig.2 Location differentia between the aggradation thickness as results of the two kinds of algorithm

Low erSedD T的求解思路与U pperSedD T基本相同,即每个时段的沉积速度与沉积持续时间积的和。SedV[1]s为开始沉积的所在时间分段的沉积速度;Low erConSed T[1]s为沉积开始时间StartT的所在时间分段沉积持续时间;SedV[i]T为所求目标时间T所在的时间分段的沉积速度;Low er-ConSed T[i]T为所求目标时间T所在的时间分段的沉积持续时间。

比较式(6）和式(11）可以看出:式(11）要多求解一个SedV[1]s和Low erConSed T[1]s;此外,由于StartT与目标时间T在时间分段上的位置差异,使得LowerConSed T[1]s和LowerConSed T[i]T在时间分段上定位的情况可能比较多。仅考虑所求目标时间在沉积起始时间之后的情况,StartT与目标时间T可能性分布如图3,共有N×(N+1）/2种可能性。

图3 StartT与目标时间 T可能性分布Fig.3 The possible distribution of Sta rtT and target T

在加积算法实现时,可以先判断Start T与T的先后位置。如果StartT在T之后,则没有沉积,整个Low erSedD T为零;如果StartT在T之前,需分别判断StartT和T所在的时间分段,这样将N×(N+1）/2种可能性的判断减少为2N次判断。

在剥蚀算法实现时,在判断StartT与T的先后位置后,如果StartT早于T,则从所求目标时间T开始计算沉积厚度;如果StartT晚于T,则从沉积起始时间StartT开始计算沉积厚度。这样只需判断StartT或T的时间分段的位置,最多判断N次。与加积算法相比对于一个离散的栅格单元就减少一半的判断次数,当对整个区域进行计算时可以大大减少循环判断的次数,提高运算速度。

综上所述,剥蚀算法和加积算法虽然在计算步骤上基本相同,而且两个计算公式可以相互转换,因此它们在原始数据相同,控制参数数据项目、反演目标时间相同的情况下,反演计算的结果是相同的。但是它们在时间分段技术上采用算法思想不同,在计算对于沉积持续时间位置判断的次数上有很大差异。从上面的比较可以看出,如果把整个沉积时间分为N段,则剥蚀算法只需进行2N次判断,而加积算法则需进行N×(N+1）/2次判断,剥蚀算法能够大大提高计算速度。

2 三维古地面计算的剥蚀算法试验

三角洲是陆海相互作用的产物,海面-地面系统变化是三角洲发育演变的控制性因素[1]。研究三角洲地区的海面-地面系统变化可以为三角洲平原古环境反演研究提供依据。本文选取长江三角洲地区作为试验区域进行三维古地面反演计算的剥蚀算法试验。试验区范围为:北纬32°10′～30°30′,东经119°20′～122°。根据剥蚀算法,在A rcGIS平台下,通过VBA+Desktop的开发方式对该算法进行开发,在开发的试验平台上进行了试验区域的剥蚀算法试验研究。在试验数据收集准备上,采用现代地面高程、全新世沉积厚度古海面标高、全新世分时段沉积速率比值、构造升降速率、均衡调整速率、压实作用速率等数据作为古地面反演计算模拟的原始数据,将考古遗址点数据、古海滨砂堤数据作为试验结果验证数据。在试验中空间离散参数选取了50 m ×50 m的栅格;在时间离散上对1万年前的地质历史时期划分了早、中、晚3个时期,分别为1 000～7 500 aB.P.、7 500～2 500 aB.P.、2 500 aB.P.至今;压实作用速率和地壳均衡调整速率采用与沉积厚度呈线性关系的控制关系式参数数据。在不同沉积速率比等控制参数条件下对试验区域进行三维古地面反演计算,并通过验证数据对模拟结果进行自动、半自动的验证分析,最终得出较好的模拟结果。其模拟结果的部分古地面时间切片如图4。试验结果显示了该研究区的现代海陆格局大体上在3 000 aB.P.奠定基础。在全新世海侵的大背景下,尽管这次模拟所采用的海面控制参数在全新世中期曾经几度高于现代海面,但从模拟结果看,研究区从未成为汪洋大海,有大片的陆地可供古人栖息。但是也不排除研究区内历史上可能有不少地带曾受海洋的影响,几度成为0～-2 m的海滨湿地。

图4 部分模拟结果的二维渲染Fig.4 Some simulation results of two-dimensional renderingmap

3 结论与讨论

(1）本文在前人研究的基础上提出了三维古地面反演计算的剥蚀算法,并与传统的加积算法相比,发现剥蚀算法在对时间序列进行离散、分时段求沉积速率和沉积持续时间,从而求得古地面的计算模拟上,有着沉积时间位置判断简单、方便计算的优点。

(2）本算法的核心思想是对研究区进行时空离散逐步逼真的反演模拟,其中最重要的是通过时间分段技术所进行的时间维度的离散。该算法思想与原始数据无关,可以广泛应用于主要由沉积作用造成的古地形变化的反演计算过程中,目前已有学者把该算法应用于黄土高原的古地形反演模拟中。已实现的算法模型可以广泛应用到大河三角洲冲积地区的古地形反演计算。该算法的计算结果精度取决于研究区数据的精度,尤其是沉积厚度数据的精度。

(3）利用GIS平台,实现了三维古地面反演计算的剥蚀算法,并以长江三角洲地区作为试验区,进行了古地面虚拟反演计算,取得较好的模拟效果。

(4）在本研究的算法实现和地区试验中,仍需进一步探讨和完善。例如,本文仅对全新世以来的沉积延续时间根据早、中、晚三个时期离散成3段,并分别求出这3段的沉积速率,在今后的研究中应对时段进一步离散后再进行计算,并进行比较分析,得出在计算效率和结果上最佳的分段组合。

[1] 谢志仁,王琳.三角洲地区海面—地面系统整合研究[J].南京师范大学学报(自然科学版）,2006,29(2）:110-114.

[2] 钟鹤翔,谢志仁,崔树红,等.海面—地面升降综合影响模拟和评估系统研究[J].海洋测绘,2004,24(4）:32-39.

[3] 钟鹤翔,谢志仁.基于 GIS的海面—地面系统变化影响评价模型研究[J].海洋测绘,2007,27(1）:30-33.

[4] 谢志仁.海面变化与环境变迁:海面—地面系统和海—气—冰系统初探[M].贵阳:贵州科技出版社,1995.

[5] 信忠保,谢志仁,蒋卫国,等.长江三角洲地区海面—地面系统模型构建与试验[J].地理研究,2005,24(4）:513-519.

[6] TANABE S,SA ITO Y,SA TO Y,et a1.Stratigraphy and Holocene evolution of themud-dominated Chao Phraya delta,Thailand[J].Quaternary Science Review s,2003,22:789-807.

[7] ALL ISON M A,KHAN S R,GOODBRED S L J,et a1.Stratigraphic evolution of the late Holocene Ganges Brahmaputra lower delta p lain[J].Sedimentary Geology,2003,155:317-342.

[8] DIRSZOWSKYA R W,DESLOGESJ R.Evolution of the Moose Lake Delta,British Columbia:Implications for Holocene environmental change in the Canadian Rocky Mountains[J].Geomorphology,2004,57:75-93.

[9] 王成善,刘志飞,王国芝,等.新生代青藏高原三维古地形再造[J].成都理工学院学报,2000,20(1）:2-7.

[10] 张凤英,刘尚仁.沉积物压实量分式及其在研究海面高度中的应用[J].中山大学学报(自然科学版）,1996(S1）:38-44.

[11] 王凯,李国胜,方国洪,等.大河三角洲河口海岸机理模型研究: (I）模式理论与进展[J].地理研究,2003,22(1）:21-29.

[12] 许红昆,林畅松.沉积物压实的一种校正方法[J].地学前缘, 2000,7(2）:366.

[13] 杨建明.海平面变化研究中沉积物压实量的估算[J].台湾海峡,1988,7(4）:382-388.

Research on Denudation Arithmetic of Three-Dimensional Paleosurface Inversion Algorithm in Sea Level-Land Surface System

WANG Lin1,2,ZHONG He-xiang3,XIE Zhi-ren4
(1.N anjing Institute of Geography and L im nology,Chinese Academ y of Sciences,N anjing 210008;
2.Graduate University of Chinese Academ y of Sciences,Beijing 100039;3.Jiangsu Information Center, N anjing 210013;4.College of Geographical Science,N anjing N orm al University,N anjing 210097,China）

Three-dimensional paleosurface inversion algo rithm is a basic study in sea level-land surface system research.Based on p revious studies,a denudation arithmetic of three-dimensional paleosurface inversion calculation is p roposed.The algorithm is to peel off the layersof increased heightw hich was caused by crustalmovement,sedimentation,and sedimentary compaction,etc.in a certain geological period of time,and inverse the height of ancient terrain.Then,a datum conversion of the inversed height is conducted using the ancient sea surface height at that time,and get the value of discrete pointsof the ancient elevation.The algorithm is independent w ith the raw data,and can be w idely used in the calculation of the inversion p rocess mainly caused by deposition of ancient terrain changes.The calculation accuracy of denudation arithmetic depends on the accuracy of the data in study areas,especially the accuracy of sedimentary thickness.Compared w ith traditional aggradation arithmetic,the denudation arithmetic ismuch more convenient in time sequence discrete judgments.In addition,the Yangtze River Delta region is selected as a case study area to perfo rm the denudation arithmetic,and better results areobtained w hen an inversion calculation to simulate 10 thousand years terrain changes in this region p rocessed.

geographic information system;sea level-land surface system;paleosurface inversion algo rithm

P208

A

1672-0504(2010）02-0049-05

2009-09-27;

2009-11-24

国家自然科学基金项目(40171008）

王琳(1981-）,女,博士研究生,主要从事资源利用与生态效应、区域环境建模及其应用等研究。E-mail:wanglinnt@126.com