柴永杰

〔深圳市勘察測繪院(集團)有限公司,廣東 深圳 518028〕

0 引 言

地球外部重力場屬于地球的空間物理場范疇,其表征了地球內(nèi)部物質(zhì)的密度、體積等物理和幾何量的整體空間分布情況。因此,重力方法理論上來說可以作為一種手段來研究地球的內(nèi)部物質(zhì)分布。隨著重力衛(wèi)星GRACE[1]和GOCE[2]的成功發(fā)射,人們獲取了史無前例的高分辨率、高精度和(幾乎)全球覆蓋的靜態(tài)和時變重力場信息,從而為采用重力方法研究地球內(nèi)部物質(zhì)分布注入了新的活力和手段。Moho面為地殼與地幔之間的分界面。該面一個顯著的特點是密度跳躍(也稱不連續(xù)),因而采用重力方法研究Moho面的基本原理即是根據(jù)該界面密度跳躍特性所產(chǎn)生的重力異常建立起重力異常與Moho面參數(shù)之間的函數(shù)關系。

確定Moho面深度具有很廣泛的地學應用。例如可以更好地理解和認識地震觸發(fā)的機制、地熱的流向、熱平衡及其演化和板塊構造等地學過程。在定義與擴展大陸架方面,Moho面的深度也是重要的依據(jù)[3]。在大地測量方面,Moho深度也有很多重要的應用,例如Moho面深度可以用來計算均衡重力異常。均衡重力異常相對于自由空氣重力異常和布格重力異常更加平滑,因此為了減小重力異常內(nèi)插和外推的誤差通常采用均衡重力異常用于內(nèi)插和外推[4]。另一方面,為了提高GOCE重力梯度觀測值向下延拓反演平均海平面上重力異常值的穩(wěn)定性和精度,可以利用地形-均衡信息對GOCE沿軌觀測的重力梯度值進行歸算,而該歸算需要用到Moho面的深度信息。地形-均衡信息還可以用來構建高階重力場模型例如EGM2008[5]的高階部分的球諧系數(shù)主要來源于地形-均衡的位系數(shù)。在地球物理領域,Moho面結(jié)合地形信息可以用來反演巖石圈的有效彈性厚度[6];皮嬌龍等[7]研究了青藏高原區(qū)域Moho面起伏對通道流模型的動力學響應。在地質(zhì)方面Moho面深度與成礦之間也存在關聯(lián),如陳安國等的研究所述[8]。

1 研究進展

現(xiàn)階段主要有地震和重力兩種方法來確定Moho面的深度。當然也可采用重力和地震聯(lián)合的方法。地震方法的基本原理是基于地震波在Moho面處的反射和折射方程來建立數(shù)學模型,而重力方法則是基于Moho面的密度不連續(xù)性所引起的重力異常大小來建立觀測方程。地震和重力方法是基于不同的原理和假設,因此這兩種方法得出的結(jié)果必然存在差異[9]。

隨著新一代重力衛(wèi)星GRACE[1]和GOCE[2]的成功發(fā)射,人類跨入了衛(wèi)星重力遙感時代,從而可以快速、高效地獲取高分辨率、高精度和(幾乎)全球覆蓋的均勻重力場信息。在此衛(wèi)星重力遙感時代背景下,采用重力方法研究地殼結(jié)構能夠有效的克服地震數(shù)據(jù)的不均勻性和分辨率低等缺點,尤其在地震數(shù)據(jù)稀疏甚至完全缺失的情況下。

采用重力模型反演Moho面深度,國內(nèi)外學者做出了大量的工作。本文認為這些方法可以概括為4種:① 以地殼均衡假設為基礎導出的算法,例如Vening Meinesz-Moritz(VMM)方法[9-10];② 以FFT算法為基礎的Parker-Oldenburg 區(qū)域算法[11-12];③ 以球諧分析和球諧綜合為基礎的Parker-Oldenburg擴展算法[13],且該算法也能利用FFT技術;④ Moho引力信號提取的直接方法其包括一階泰勒線展開性化法、凝聚法。需要指出的是,以上各種算法之間也是相互聯(lián)系和一脈相承的,并非完全獨立。本文對以上各種算法進行了回顧和分析,并指出其優(yōu)缺點。

1.1 地殼均衡假設

地殼均衡假設是指地殼處于靜態(tài)平衡狀態(tài),即地殼的密度和厚度在不同區(qū)域內(nèi)達到平衡。Moho面是地球的地殼和地幔之間的界面,是地殼和地幔密度和速度差異的表現(xiàn)。地殼均衡假設和Moho面之間存在密切的關系。根據(jù)地殼均衡假設,地殼的密度和厚度是在不同的地質(zhì)時間尺度上達到平衡的,而Moho面是地殼和地幔之間的物理界面,隨著地質(zhì)時間的演化,其深度和形態(tài)也在不斷變化。因此,研究地殼均衡假設可以幫助我們更好地理解Moho面的形成和演化過程。此外,地殼均衡假設也可以用于解釋地球內(nèi)部的一些現(xiàn)象,例如地震波傳播的速度和路徑,地殼形變等。通過對地殼均衡假設和Moho面的研究,我們可以更深入地了解地球內(nèi)部的結(jié)構和演化過程。

在19世紀中期,Pratt在喜馬拉雅山脈地區(qū)通過測量發(fā)現(xiàn)實測的垂線偏差小于可見的地形質(zhì)量所計算的理論結(jié)果,因此他推斷地球內(nèi)部一定存在某種質(zhì)量虧損來解釋上述差異。另外,假定地形僅僅是簡單的疊加在勻質(zhì)的地殼之上,則布格重力歸算將會扣除地形起伏所引起的大部分不規(guī)則重力場信號,因此布格重力異常理論上將會非常的平滑且其大小應該在零值附近擺動,然而事實卻并不如此且布格重力異常具有系統(tǒng)性特征,即在山區(qū)表現(xiàn)為大的負值。對于以上理論與實際不符的情況,唯一的解釋就是在山區(qū)地形以下有質(zhì)量的虧損,換句話說,地形質(zhì)量存在某種補償。本文簡單回顧一下3種經(jīng)典的地殼均衡模型。

1.1.1 普拉特-海福特模型

Pratt首先提出了該模型的核心思想,后來Hayford推導了詳細的數(shù)學計算公式并且把該模型廣泛的應用于大地測量學中,因此通常把該模型稱之為Pratt-Hayford模型。Pratt-Hayford模型(圖1)的基本原理是位于補償面以下的物質(zhì)密度是均勻的,位于補償面以上每個截面的密度隨著地形高度的變化而變化。即位于補償面以上每個截面的質(zhì)量相等,因此高出海水面的地形越高,則該截面的密度越小;低于海水面的地形越深,則該截面的密度越大。

圖1 Pratt-Hayford均衡模型

1.1.2 艾黎-海斯卡涅模型

該模型是Airy首先提出,為了便于大地測量學的應用,Heiskanen給出了精確的數(shù)學表達式,并對其進行了廣泛的應用。因此通常把該模型稱之為Airy-Heiskanen模型(圖2)。該模型的基本思想是地形漂浮于密度較大的巖漿(地幔)上面,因此地形越高,則在重力的引導下其下沉得越深。在山區(qū),其截面密度為常數(shù)ρc=2.67 g/cm3,下面密度較大的地幔密度為ρm=3.27 g/cm3,且定義地殼和地幔密度差為Δρc/m=ρm-ρc。根據(jù)Airy模型,在山區(qū)下面會形成山根,海洋下面會形成反山根。

圖2 Airy-Heiskanen 均衡模型

地震研究表明,Airy-Heiskanen均衡模型比Pratt-Hayford均衡模型更加的接近實際地球狀態(tài),但是Pratt-Hayford均衡模型也存在一定的合理性,所以有些學者提出了聯(lián)合均衡模型,即在陸地上采用Airy-Heiskanen均衡模型,在海洋上采用Pratt-Hayford均衡模型。盡管如此,實際情況也比聯(lián)合均衡模型要更復雜,例如地球的曲率影響、地殼的彈性性質(zhì)、冰川均衡調(diào)整(GIA)等因素都會對均衡狀態(tài)有所影響。

1.1.3 維寧曼尼斯區(qū)域模型

通過考察Pratt-Hayford和Airy-Heiskanen模型,Vening Meinesz發(fā)現(xiàn)上述兩種模型是局部補償機制,也就說補償僅僅沿著界面的垂直方向,顯然這不符合實際情況,因此Vening Meinesz在Airy模型的基礎上進行了改進,提出了區(qū)域補償模型。局部補償和區(qū)域補償?shù)牟町惾鐖D3所示。根據(jù)Vening Meinesz假設,地形被視為一種負載施加于剛性的彈性地殼上面,也就說其考慮了地殼本身的彈性性質(zhì),因此部分地形負載不僅僅會通過垂直方向也會通過水平方向進行補償。通過以上分析,可知Airy-Heiskanen模型是Vening Meinesz模型的一種特殊情況。

圖3 區(qū)域補償與局部補償之間的差異

現(xiàn)實中地球均衡補償機制非常復雜,涉及的因素很多,例如地質(zhì)構造運動、冰川均衡調(diào)整等地球動力學過程。為了盡可能逼真的對實際地球均衡狀態(tài)進行描述,在Moho面反演建模的過程中充分吸收Pratt-Hayford、Airy-Heiskanen、Vening Meinesz和Moritz的思想,即地形的均衡補償是由于密度補償、山根和反山根在全球尺度上的綜合補償所致,該思想可以抽象成是一種“廣義均衡模型”。

1.2 Vening Meinesz-Moritz算法研究現(xiàn)狀

Pratt模型是以地殼密度的虧損或者剩余來補償?shù)匦蔚馁|(zhì)量。而Airy模型是以地殼厚度的虧損或者剩余去補償?shù)匦蔚馁|(zhì)量。Pratt和Airy均衡補償模型都是一種局部(Local)補償模式。因此,從物理角度分析,這兩種補償模型都忽略了地殼張力等其他因素對地殼均衡的影響,為此Vening Meinesz(1931)引入巖石圈的撓曲模型對Airy均衡模型進行了進一步的拓展,提出了區(qū)域(Regional)補償模型。后來,Moritz[10]進一步將Vening Meinesz區(qū)域均衡模型擴展成全球(Global)均衡模型。Sj?berg[9]對Moritz[10]全球均衡模型進行了數(shù)學上的重新表述,并將求解Moho面深度問題歸納為求解第一類非線性Fredholm積分方程,并進一步推導了顧及非線性二次項的直接反演公式。Moritz[10]和Sj?berg[9]方法都是以全球均衡思想為基礎,但不同的是前者給出的是迭代解,而后者給出的是基于球近似導出的直接解。

Vening Meinesz-Moritz(VMM)算法雖然成功的反演了不同區(qū)域的Moho面深度,但是本文認為該算法存在3個缺點:① 當包括到非線性項的二次項時該反演方程包括一個二維空間積分項,即任意一點的積分理論上都需要進行全球積分導致算法的計算效率較低;② 當計算點和積分點之間的球面距離為零的時候上述空間積分會出現(xiàn)奇異值,因此還需要進一步近似處理才能消除奇異值的影響(例如,計算點附近采用平面近似);③ 當顧及非線性項的三次項及更高次項的時候?qū)С龅乃惴〞拥膹碗s,因此限制了該算法的反演效率和精度。

1.3 Parker-Oldenburg算法研究現(xiàn)狀

Parker[12]提出了一種重力場正演的快速算法,該算法的基本思想是在平面直角坐標系下(平面近似)將任意界面所引起的重力異常表達為傅里葉級數(shù)。接著,Oldenburg[11]對Parker的正演算法進行了簡單的數(shù)學變換并提出了采用重力異常迭代求解界面深度的反演算法。同時,為了提高迭代收斂的速度和穩(wěn)定性,提出了一個濾波器,后來人們通常稱該算法為Parker-Oldenburg(P-O)算法。本文認為該算法的一個顯著優(yōu)點是采用了FFT技術,從而大大提高了計算效率,同時也提高了對分辨率較高的格網(wǎng)數(shù)據(jù)的處理能力。同時,該算法的一個明顯缺點是其基于平面近似,因此只能用來進行小范圍的局部研究,即范圍不能過大,否則地球的曲率對反演結(jié)果影響不容忽視。最初,Parker-Oldenburg算法只能用于二維重力數(shù)據(jù)的反演。后來,Gomez-Ortiz and Agarwal[14]對Parker-Oldenburg算法進行了擴展,將其由二維擴展到三維,因此大大提高了其應用范圍。從均衡的角度看,本文認為Parker-Oldenburg算法也可以看成在Vening Meinesz區(qū)域均衡模型基礎上導出的算法,只不過算法適用于平面近似的區(qū)域,但是其補償范圍比Airy均衡模型要廣泛。

1.4 Parker-Oldenburg擴展算法研究現(xiàn)狀

正如前面所論述,雖然Parker-Oldenburg計算效率高,但是由于其是基于平面近似導出的算法,因此該算法只能用于平面近似區(qū)域。為了能夠?qū)FT技術也能應用于全球Moho面的快速反演,Chen and Tenzer[13]借鑒Parker-Oldenburg的思想,在球近似的情況下提出了一種新的算法,暫且稱之為“Parker-Oldenburg擴展算法”,從而有效提高了反演全球Moho面深度的計算效率。和Vening Meinesz-Moritz方法比較,“Parker-Oldenburg擴展算法”表現(xiàn)出了明顯優(yōu)勢,主要體現(xiàn)在4點:① 避免了任意一點的Moho面求解都需要費時的全球積分計算;② 不存在求解方程的奇異值問題;③ 該算法的二次以及更高次項都可以利用FFT技術進行計算且一次、二次及更高次項的公式具有一致性,即不存在階次越高算法越復雜的情況;④ 該算法可以顧及Moho面的橫向變化(二維變化)從而有效的提高反演結(jié)果的精度。但是本文認為“Parker-Oldenburg擴展算法”仍然存在以下缺點:① 平均Moho面深度和密度差選擇不夠精確會導致迭代不收斂;② 該算法不能夠根據(jù)陸地和海洋的特征分別選擇先驗參數(shù),因此導致反演的結(jié)果在海洋和陸地上均具有一定的系統(tǒng)性誤差。

從均衡的角度去考察,本文認為“Parker-Oldenburg擴展算法”也可以理解為以“廣義均衡模型”為基礎導出的算法。因為該算法在數(shù)學建模的時候,既考慮了Moho面深度的變化,同時也考慮了Moho面密度差的變化,換句話說,高出平均海水面的地形,既不是單一的Pratt模型提出的密度虧損來補償,也不是單一的Airy模型提出的地殼厚度來補償,而是密度和地殼聯(lián)合補償,同時該算法也包含了Moritz的全球均衡思想。即可以認為該算法綜合了Pratt、Airy以及Moritz 3種均衡模型的思想,認為“廣義均衡模型”是“Parker-Oldenburg擴展算法”的物理基礎。眾所周知,地殼的均衡狀態(tài)是一個復雜的物理過程,其主要取決于巖石圈的負載、巖石圈有效的彈性厚度、硬度及其流變學特性、軟流圈的黏度等因素。同時,冰川均衡調(diào)整、冰川融化、板塊運動、地幔對流等地球動力學過程也對地殼均衡產(chǎn)生一定的影響。因此,本文認為基于地殼均衡模型的重力模型存在一定的物理缺陷,其主要體現(xiàn)在目前還難以模型化的地球動力學效應所引起的重力異常信號如何扣除,例如采用重力方法所反演的Moho面在俯沖帶區(qū)域具有比較大的缺陷。

1.5 直接法研究現(xiàn)狀

從信號恢復的角度看,不基于任何均衡假設,而是從重力觀測值中利用地震、地球物理等輔助數(shù)據(jù)恢復Moho面所引起的重力異常信號,然后采用扣除各種重力異常之后的精化重力數(shù)據(jù)進行Moho面反演,本項目稱之為“直接法”。直接法大體上可分為兩步進行,首先將計算點至Moho面之間以及Moho面到地心之間的所有重力異常信號從觀測值中進行扣除;然后把扣除各種重力異常信號之后的精化重力異常當作輸入數(shù)據(jù)反演Moho面深度。

按照數(shù)學處理的不同,本文將其大致分為一階泰勒線展開性化方法和凝聚法等。從數(shù)學形式上看,泰勒展開線性化方法就是對第一類非線性Fredholm積分方程利用泰勒展開定理進行線性化,但是其展開的位置選擇平均Moho面深度處,從而可大大簡化算法的復雜度和計算量。凝聚法的基本思想是將Moho面相對于平均Moho面的起伏壓縮到平均Moho面深度所對應的球面上,然后導出相應的算法。需要指出的是直接法的一個核心問題是如何干凈、有效的扣除非Moho面所引起的重力異常信號。

隨著計算機技術的不斷發(fā)展和重力數(shù)據(jù)獲取的進一步完善,重力反演在Moho面深度研究中的應用越來越廣泛,而且具有很大的潛力。首先,重力數(shù)據(jù)的精度和覆蓋范圍將會得到進一步提高。其次,多源數(shù)據(jù)融合將會成為一種趨勢。重力反演方法可以與其他地球物理學方法結(jié)合使用,例如地震勘探、熱流場分析、地表形貌等數(shù)據(jù),這些數(shù)據(jù)的融合可以增加反演結(jié)果的可靠性和精度。第三,反演算法的改進和發(fā)展也將會推動重力反演Moho深度的研究,未來隨著新的反演算法的發(fā)展,反演結(jié)果的準確性和可靠性將會進一步提高。

2 結(jié) 語

通過對重力反演Moho算法的研究現(xiàn)狀剖析,未來可以對重力反演Moho面算法從如下幾個方面進行改進。

(1)重力數(shù)據(jù)的優(yōu)點是分辨率高、精度高和覆蓋范圍廣;地震方法的優(yōu)點是反演結(jié)果精度和可靠性高,缺點是數(shù)據(jù)采集成本高,導致地震數(shù)據(jù)分辨率和精度不均勻、某些區(qū)域分布稀疏甚至缺失導致這些區(qū)域的反演結(jié)果具有較大的不確定性。因此,第一個切入點是推導聯(lián)合重力和地震優(yōu)點的新算法。

(2)當前計算非Moho面所引起重力異常的時候采用的地殼密度通常為常密度(一維)或者橫向變化的密度(二維),但是我們知道地殼真實密度是隨著深度和位置三維變化,因此采用一維或者二維的密度值去替代三維的密度值必然會產(chǎn)生誤差。因此,第二個切入點是聯(lián)合多源數(shù)據(jù)精化地殼模型,然后推導顧及地形密度三維變化的重力場正演算法。

(3)基于地殼均衡模型的重力算法存在一定的物理缺陷,主要是地球動力學效應(例如地殼形變)難以模型化;重力算法根據(jù)研究范圍的大小采用一定的近似:局部的采用平面近似,區(qū)域和全球的采用球近似。由于真實地球的形狀更接近橢球,平面和球近似必然會產(chǎn)生近似誤差。因此,第三個切入點是在數(shù)學上采用更加接近真實地球形狀的橢球作為數(shù)學推導的基準面;在物理上聯(lián)合GNSS、衛(wèi)星重力和水文數(shù)據(jù)量化地殼形變與Moho面變化之間的關系。

重力反演Moho深度在地球科學中的應用將會更加廣泛。隨著對地球內(nèi)部結(jié)構和演化的研究不斷深入,重力反演Moho深度不僅可以用于地震、地質(zhì)勘探等領域,還可以用于地球物理學、地球化學等領域的研究,有望在未來取得更多的研究成果。