王 琳,鐘鶴翔,謝志仁
(1.中國科學院南京地理與湖泊研究所,江蘇南京 210008;2.中國科學院研究生院,北京 100039; 3.江蘇省信息中心,江蘇南京 210013;4.南京師范大學地理科學學院,江蘇南京 210097)
海面-地面系統(tǒng)中三維古地面反演計算的剝蝕算法研究
王 琳1,2,鐘鶴翔3,謝志仁4
(1.中國科學院南京地理與湖泊研究所,江蘇南京 210008;2.中國科學院研究生院,北京 100039; 3.江蘇省信息中心,江蘇南京 210013;4.南京師范大學地理科學學院,江蘇南京 210097)
三維古地面反演計算是海面-地面系統(tǒng)研究中的一個基礎研究。該文在前人研究的基礎上提出了三維古地面反演計算的剝蝕算法,并與傳統(tǒng)的加積算法做了比較分析,得出剝蝕算法比加積算法有對時間序列離散判斷簡單、方便計算等優(yōu)點。同時,在GIS平臺下實現(xiàn)了該算法,并以長江三角洲地區(qū)為試驗區(qū),進行了古地面虛擬反演計算,取得了較好的效果。
地理信息系統(tǒng);海面-地面系統(tǒng);古地面反演計算
地理環(huán)境的表達和認知、古地理環(huán)境的恢復和再現(xiàn),始終是國內(nèi)外地學、地圖學、認知科學領域普遍關注的熱點;研究并建立基于 GIS技術(shù)和虛擬環(huán)境技術(shù)的沿海地區(qū)古地貌演變模型是海面-地面系統(tǒng)變化研究的一個重要方向[1-3]。三維古地面反演計算是古地理環(huán)境恢復和再現(xiàn)的重要研究內(nèi)容之一,也是海面-地面系統(tǒng)研究中的基礎內(nèi)容。海面-地面系統(tǒng)研究就是將海面升降變化和地面升降變化作為一個統(tǒng)一的整體加以考察,通過對保存在地面外部形態(tài)特征和內(nèi)部地層構(gòu)成中有關信息的提取,用來演繹海面和地面變化的歷史過程、當前動態(tài)和未來趨勢[4]。在海面-地面系統(tǒng)中,如何進行古海面-地面變化的反演以及未來海面-地面變化的預演是整個研究的核心。對于海面-地面相對運動的模擬計算、古地面高程的反演計算、古地形模擬已有一些研究成果[5-13]。本文在前人研究成果的基礎上,設計出海面-地面系統(tǒng)中三維古地面的剝蝕算法,并在GIS平臺下實現(xiàn)了該算法;然后選擇長江三角洲地區(qū)進行三維古地面反演模擬,取得了較好的效果。
傳統(tǒng)的三維古地面反演計算模型是以沉積物隨時間的推移在沉積基底之上逐漸加積從而導致地面升高為主線,同時考慮到海面的升降、地殼升降運動與沉積物壓實沉降等地球物理作用而設計的概念模型,在此基礎上通過試驗設計出三維古海面計算模型。該算法在時間尺度離散上存在著時間斷點判斷復雜、計算困難等缺點。鑒于此,通過逆向思維方式提出了三維古地面計算的剝蝕算法。算法的基本思路如下:現(xiàn)代地形中空間上離散點的高程值,都是由地殼運動、沉積作用及沉積壓實作用等在一定的地質(zhì)時間內(nèi)共同作用的結(jié)果。因此從現(xiàn)代地形中層層剝離掉這些作用造成的地面增高,便可以反演古地形的高度,同時與當時的古海面高程進行基準面轉(zhuǎn)換,從而得到各離散點的古高程值。
對于單一離散空間點的古高程計算式如下:
式中:O ld ET為歷史反演目標時間 T時期的古高程值;L an E為現(xiàn)代的地形高度;U pperSedDT為歷史反演目標時間T到現(xiàn)在這一連續(xù)時段內(nèi)的沉積物厚度;U pperConSed T為歷史反演目標時間 T到現(xiàn)在這一連續(xù)時間長度;TecV、EusV、Com V分別為歷史反演目標時間T到現(xiàn)在這一連續(xù)時段的平均地殼升降運動速率、地殼均衡運動速率和沉積壓實作用速率;O ldSLT為歷史反演目標時間 T時期的古海面高程,用來進行古海面基準轉(zhuǎn)換計算。其計算流程如圖1所示。
計算式(1)中用到的各數(shù)據(jù),假設將10 000 aB.P.至今分為 N個時段,并給出每個時間分段的時間點。具體計算流程如下。
(1)計算各時段的沉積持續(xù)時間:ConSed T[i]。每個地面離散點的沉積持續(xù)時間首先取決于它開始接受沉積的第一時段的位置。在精細的古地貌模擬中,不能將地層界面視為等時面。例如,沿海地帶不同地點的全新世沉積起始時間應該是一個與古地面高程和遭受全新世海侵相關的變量,不能簡單地一律設定為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個時段內(nèi)的沉積延續(xù)時間;Sub T[i]為時間段的分段時間點;StartT為該離散點上開始接受沉積的起始時間。3個公式分別表示3種情況下ConSed T[i]的值。計算時先判斷StartT位于哪個時間分段內(nèi),假設其位于第i個分段內(nèi),則第i個分段以前的各段沒有沉積,沉積持續(xù)時間為0,即:ConSed T[i-1]=0;第i個分段內(nèi)的沉積持續(xù)時間為第i個時段的分段時間點與沉積開始時間之差,即:Start T-Sub T[i];第i個分段后的所有分段的沉積持續(xù)時間為相鄰兩個分段時間點的差,即:Sub T[i]-Sub[i+1]。
(2)計算不同時段內(nèi)的沉積速率:SedV[i]。由于整個沉積時期的不同時段內(nèi)的沉積速度不同,要想達到最佳的歷史時期古地貌的擬合效果,就必須把整個沉積厚度離散分段,根據(jù)專家知識對不同時期賦以不同的比值,這樣在計算中不同的時期就會有沉積速率的差異。SedV[i]的計算公式為:
式中:SedV[i]為第i個時間分段內(nèi)的沉積速率; Sed T為整個模擬時段內(nèi)形成的沉積厚度;ConSed T [1]…ConSed T[N]為每個分段的沉積延續(xù)時間;k [i]為沉積速率分量控制參數(shù)。
(3)計算沉積厚度(即每個時段的沉積速度與沉積持續(xù)時間積之和):U pperSedD T。公式如下:
式中:SedV[i]T為所求目標時間T(或沉積開始時間StartT)所在時間分段的沉積速度;ConSed T[i]T為所求目標時間T(或沉積開始時間Start T)所在時間分段的沉積持續(xù)時間。
(4)計算升降校正:Rec H T。公式為:
式中:TecV、EusV和Com V分別是該離散點的垂直方向上構(gòu)造運動速率、均衡運動速率和壓實作用速率;U pperConSed T則為所求目標時間到目前的沉積持續(xù)時間。
(5)計算古地面高程?;竟綖?
式中:O ld ET為歷史反演目標時間 T以現(xiàn)代高程系為基準的古地面高程;Land E為該地面離散點的現(xiàn)代地面高程;U pperSedDT為歷史反演目標時間 T到現(xiàn)在這一連續(xù)時段內(nèi)的沉積物厚度;R ec HT為升降校正高度。
(6)古海面基準轉(zhuǎn)換計算?;竟綖?
式中:OldSLT為基于現(xiàn)代高程系的古海面高程,當古海面值低于現(xiàn)代海面時取古海面高程為負值,否則為正值;SLO ld ET為基于古海面基準高程的古地面高程值,該值體現(xiàn)的是反演目標時間的海陸關系。
通過以上計算可以得到三維古地面的高程,然后通過檢驗,對驗證結(jié)果進行判斷,基本符合驗證數(shù)據(jù)的結(jié)果便是所求結(jié)果,對于不符合的則重新進行計算。
從算法模型的古高程計算方程可以看出,該算法所需的原始數(shù)據(jù)主要有4類:沉積啟動時間數(shù)據(jù)、沉積厚度數(shù)據(jù)、沉積校正數(shù)據(jù)和古海面高程數(shù)據(jù)。其中沉積校正數(shù)據(jù)又包含:地殼構(gòu)造運動、地殼均衡運動和沉積壓實作用產(chǎn)生的沉積校正數(shù)據(jù)。在整個反演算法計算過程中,通過對原始數(shù)據(jù)進行自動半自動的時空離散來達到算法的自適應計算模擬,并通過對模擬結(jié)果的驗證逐步向真實的古海面-地面地形逐步逼近,以得到最佳反演結(jié)果。因此,算法中除原始數(shù)據(jù)外,還有時間離散控制參數(shù)和空間離散控制參數(shù)。
(1)沉積啟動時間數(shù)據(jù)。沉積啟動時間數(shù)據(jù)是控制反演區(qū)域沉積開始時間的數(shù)據(jù),它是通過空間離散控制參數(shù)把整個反演區(qū)域離散為柵格點沉積起始時間。在晚更新世末期的低海面時期,現(xiàn)在的三角洲平原區(qū)是一個大面積處于水面以上的受侵蝕的丘陵崗地地貌,起伏的地形決定了遭受海水淹沒并開始接受沉積的時間會有很大差別。因此每個柵格點開始接受沉積的起始時間是不同的,該值可以通過滑動底界模型求得[4]。
(2)沉積厚度數(shù)據(jù)。對于三角洲近海及平原區(qū)域,第四紀以來的沉積物加積作用是造成古地面變化的最為重要的因素,因此沉積厚度數(shù)據(jù)決定了整個反演計算結(jié)果的精度,也決定了時空離散參數(shù)的設定。因為整個算法都是建立在對沉積厚度數(shù)據(jù)數(shù)理分析的基礎上,只有通過對沉積厚度的分析,才能獲取相對精確的時間分段并計算出各個分段的平均沉積速率。
(3)沉積校正數(shù)據(jù)。沉積校正數(shù)據(jù)是用來校正地殼構(gòu)造運動、地殼均衡運動和沉積壓實作用對海面-地面高程的影響。但3種運動計算比較困難,在實際的試驗中可以由簡單模式確定,如平行狀等速率線所代表的向海掀斜運動模式;地殼均衡運動、沉積壓實作用與沉積厚度呈線性關系。算法在實現(xiàn)時留有數(shù)據(jù)接口和線性關系式參數(shù)接口。
(4)古海面高程控制。古海面高程是反演古海面-地面系統(tǒng)運動演變的重要控制參數(shù),它是相對古地面的一個重要的高度參考數(shù)據(jù)。在實際試驗中可以對已知離散時間上的古海面高度進行內(nèi)插從而求得目標時刻的古海面高度,以此作為相對古地面的基本海面標高。
(5)空間離散控制參數(shù)。通過對原始數(shù)據(jù)(包括現(xiàn)代地形數(shù)據(jù)、沉積厚度數(shù)據(jù)等)進行柵格化來對整個研究區(qū)進行空間離散。在整個研究區(qū),通過數(shù)據(jù)分析,將沉積過程基本一致、數(shù)據(jù)采集良好的區(qū)域分割成一個個單一試驗區(qū)進行試驗;然后對每個試驗區(qū)進行柵格化,柵格的大小可以根據(jù)試驗的要求和原始數(shù)據(jù)的精度來確定。
(6)時間離散控制參數(shù)。在整個第四紀地質(zhì)時期的沉積過程中不同時段的沉積速度有差異,但不可能把每個時段的沉積速度都求出??梢愿鶕?jù)專家經(jīng)驗,在對原始的沉積厚度、鉆井資料分析的基礎上,把試驗區(qū)的整個沉積過程大致分成幾個相應的沉積間段,在每個沉積間段內(nèi)假設有著同樣的沉積速率。不同時段間的沉積速率有個比值,在模型計算過程中可以調(diào)整此比值來進行不同沉積間段的沉積速度控制,在算法中可以通過時間分段技術(shù)實現(xiàn)整個計算過程的自動、半自動反演計算。
從數(shù)據(jù)分析可以看出,沉積厚度數(shù)據(jù)是影響反演計算結(jié)果的最重要的數(shù)據(jù)。而對研究區(qū)進行時空離散,通過時間分段技術(shù)進行的時間維度的離散,計算不同時段內(nèi)的平均沉積速率,并可以通過模型對每個時段內(nèi)的沉積速率自動調(diào)節(jié),達到模型自我不斷反演逼近真實地形的算法思想與數(shù)據(jù)無關。
在傳統(tǒng)的算法中以沉積基底上的沉積物加積作用為基本計算依據(jù),同時考慮到海面的升降、地殼運動與沉積壓實等地球物理作用等來計算古地面,姑且把該算法稱為“加積算法”,其基本方程如下:
式中:Lan E-Sed T項的地學意義是現(xiàn)代的地形(Lan E)剝除模擬時段內(nèi)形成的沉積總厚度(Sed T)之后,以現(xiàn)代高程基面計算的古地形高程,即模擬的起始地貌面;Low erConSed T為沉積持續(xù)時間:TStartT;Low erSedDT為從沉積起始時間到特定反演目標時間T時段內(nèi)的沉積厚度,其計算公式為:
式中:Low erConSed T[i]T為目標時間T所在的第i個時間分段到第i-1個時間分段的沉積延續(xù)時間。
從式(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基本相同,即每個時段的沉積速度與沉積持續(xù)時間積的和。SedV[1]s為開始沉積的所在時間分段的沉積速度;Low erConSed T[1]s為沉積開始時間StartT的所在時間分段沉積持續(xù)時間;SedV[i]T為所求目標時間T所在的時間分段的沉積速度;Low er-ConSed T[i]T為所求目標時間T所在的時間分段的沉積持續(xù)時間。
比較式(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
在加積算法實現(xiàn)時,可以先判斷Start T與T的先后位置。如果StartT在T之后,則沒有沉積,整個Low erSedD T為零;如果StartT在T之前,需分別判斷StartT和T所在的時間分段,這樣將N×(N+1)/2種可能性的判斷減少為2N次判斷。
在剝蝕算法實現(xiàn)時,在判斷StartT與T的先后位置后,如果StartT早于T,則從所求目標時間T開始計算沉積厚度;如果StartT晚于T,則從沉積起始時間StartT開始計算沉積厚度。這樣只需判斷StartT或T的時間分段的位置,最多判斷N次。與加積算法相比對于一個離散的柵格單元就減少一半的判斷次數(shù),當對整個區(qū)域進行計算時可以大大減少循環(huán)判斷的次數(shù),提高運算速度。
綜上所述,剝蝕算法和加積算法雖然在計算步驟上基本相同,而且兩個計算公式可以相互轉(zhuǎn)換,因此它們在原始數(shù)據(jù)相同,控制參數(shù)數(shù)據(jù)項目、反演目標時間相同的情況下,反演計算的結(jié)果是相同的。但是它們在時間分段技術(shù)上采用算法思想不同,在計算對于沉積持續(xù)時間位置判斷的次數(shù)上有很大差異。從上面的比較可以看出,如果把整個沉積時間分為N段,則剝蝕算法只需進行2N次判斷,而加積算法則需進行N×(N+1)/2次判斷,剝蝕算法能夠大大提高計算速度。
三角洲是陸海相互作用的產(chǎn)物,海面-地面系統(tǒng)變化是三角洲發(fā)育演變的控制性因素[1]。研究三角洲地區(qū)的海面-地面系統(tǒng)變化可以為三角洲平原古環(huán)境反演研究提供依據(jù)。本文選取長江三角洲地區(qū)作為試驗區(qū)域進行三維古地面反演計算的剝蝕算法試驗。試驗區(qū)范圍為:北緯32°10′~30°30′,東經(jīng)119°20′~122°。根據(jù)剝蝕算法,在A rcGIS平臺下,通過VBA+Desktop的開發(fā)方式對該算法進行開發(fā),在開發(fā)的試驗平臺上進行了試驗區(qū)域的剝蝕算法試驗研究。在試驗數(shù)據(jù)收集準備上,采用現(xiàn)代地面高程、全新世沉積厚度古海面標高、全新世分時段沉積速率比值、構(gòu)造升降速率、均衡調(diào)整速率、壓實作用速率等數(shù)據(jù)作為古地面反演計算模擬的原始數(shù)據(jù),將考古遺址點數(shù)據(jù)、古海濱砂堤數(shù)據(jù)作為試驗結(jié)果驗證數(shù)據(jù)。在試驗中空間離散參數(shù)選取了50 m ×50 m的柵格;在時間離散上對1萬年前的地質(zhì)歷史時期劃分了早、中、晚3個時期,分別為1 000~7 500 aB.P.、7 500~2 500 aB.P.、2 500 aB.P.至今;壓實作用速率和地殼均衡調(diào)整速率采用與沉積厚度呈線性關系的控制關系式參數(shù)數(shù)據(jù)。在不同沉積速率比等控制參數(shù)條件下對試驗區(qū)域進行三維古地面反演計算,并通過驗證數(shù)據(jù)對模擬結(jié)果進行自動、半自動的驗證分析,最終得出較好的模擬結(jié)果。其模擬結(jié)果的部分古地面時間切片如圖4。試驗結(jié)果顯示了該研究區(qū)的現(xiàn)代海陸格局大體上在3 000 aB.P.奠定基礎。在全新世海侵的大背景下,盡管這次模擬所采用的海面控制參數(shù)在全新世中期曾經(jīng)幾度高于現(xiàn)代海面,但從模擬結(jié)果看,研究區(qū)從未成為汪洋大海,有大片的陸地可供古人棲息。但是也不排除研究區(qū)內(nèi)歷史上可能有不少地帶曾受海洋的影響,幾度成為0~-2 m的海濱濕地。
圖4 部分模擬結(jié)果的二維渲染Fig.4 Some simulation results of two-dimensional renderingmap
(1)本文在前人研究的基礎上提出了三維古地面反演計算的剝蝕算法,并與傳統(tǒng)的加積算法相比,發(fā)現(xiàn)剝蝕算法在對時間序列進行離散、分時段求沉積速率和沉積持續(xù)時間,從而求得古地面的計算模擬上,有著沉積時間位置判斷簡單、方便計算的優(yōu)點。
(2)本算法的核心思想是對研究區(qū)進行時空離散逐步逼真的反演模擬,其中最重要的是通過時間分段技術(shù)所進行的時間維度的離散。該算法思想與原始數(shù)據(jù)無關,可以廣泛應用于主要由沉積作用造成的古地形變化的反演計算過程中,目前已有學者把該算法應用于黃土高原的古地形反演模擬中。已實現(xiàn)的算法模型可以廣泛應用到大河三角洲沖積地區(qū)的古地形反演計算。該算法的計算結(jié)果精度取決于研究區(qū)數(shù)據(jù)的精度,尤其是沉積厚度數(shù)據(jù)的精度。
(3)利用GIS平臺,實現(xiàn)了三維古地面反演計算的剝蝕算法,并以長江三角洲地區(qū)作為試驗區(qū),進行了古地面虛擬反演計算,取得較好的模擬效果。
(4)在本研究的算法實現(xiàn)和地區(qū)試驗中,仍需進一步探討和完善。例如,本文僅對全新世以來的沉積延續(xù)時間根據(jù)早、中、晚三個時期離散成3段,并分別求出這3段的沉積速率,在今后的研究中應對時段進一步離散后再進行計算,并進行比較分析,得出在計算效率和結(jié)果上最佳的分段組合。
[1] 謝志仁,王琳.三角洲地區(qū)海面—地面系統(tǒng)整合研究[J].南京師范大學學報(自然科學版),2006,29(2):110-114.
[2] 鐘鶴翔,謝志仁,崔樹紅,等.海面—地面升降綜合影響模擬和評估系統(tǒng)研究[J].海洋測繪,2004,24(4):32-39.
[3] 鐘鶴翔,謝志仁.基于 GIS的海面—地面系統(tǒng)變化影響評價模型研究[J].海洋測繪,2007,27(1):30-33.
[4] 謝志仁.海面變化與環(huán)境變遷:海面—地面系統(tǒng)和?!獨狻到y(tǒng)初探[M].貴陽:貴州科技出版社,1995.
[5] 信忠保,謝志仁,蔣衛(wèi)國,等.長江三角洲地區(qū)海面—地面系統(tǒng)模型構(gòu)建與試驗[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-),女,博士研究生,主要從事資源利用與生態(tài)效應、區(qū)域環(huán)境建模及其應用等研究。E-mail:wanglinnt@126.com