文 | 唐哲
編輯 | 唐哲
大地測量方法推算闆塊碰撞起始時間的基本思想
為了推算印度闆塊和歐亞闆塊的初始碰撞時間，可以假設如下實體模型：歐亞闆塊地殼（青藏高原塊體）是一個彈性塊體，在印度闆塊持續沖撞下發生彈性變形，根據地殼均衡理論，地表面隆升而莫霍面下沉，青藏高原塊體總體上處于動态均衡狀态且不斷增厚，上述隆升機制和地殼增厚模型與Sun等［30］所采用的模型一緻［圖2（a）］，具有合理性，隻是需要說明如下幾點：
彈性地殼問題，因為在目前大多數的地球模型中均假設地殼是彈性塊體，地表面的冰川均衡調整（GlacialIsostaticAdjustment，GIA）現象或者非彈性變形均是上地幔黏彈性變形的表現形式而已，本文的主要研究對象基本上針對地殼塊體的變形問題，是以考慮彈性地殼是比較容易接受的。
動态均衡問題，印度闆塊碰撞歐亞闆塊緻使其地表上升和莫霍面下降，是假設整個青藏高原塊體處于一個理想的均衡狀态，如常見的ParrtHayford模型與Airy-Heiskanen模型［圖2（b）］。
雖然很多研究結果表明，高山區或高原區往往沒有完全處于地殼均衡狀态，但是地殼厚度模型以及高程資料總體揭示高山區接近于地殼均衡狀态，作為近似，假設青藏高原塊體處于動态均衡狀态是可以接受的。
線性隆升問題，關于青藏高原塊體隆升曆史，學術界普遍認為它是一個長期複雜且分階段的過程（圖1），目前大地測量資料結果揭示了青藏高原高原動态是具有線性趨勢特征，例如Sun等在利用重力和GPS資料研究青藏高原地殼增厚時就是假設現今高原隆升是線性的。
與地質學手段不同，大地測量已經無法取得任何曆史上的高原變形測量資料，是以隻能對高原地殼變形進行線性假設限制了，參照地質學研究結果得出非線性變形的假設是困難且難以得到青藏高原古高度，是以本研究暫不做考慮。
于是，本研究的基本思想是，在上述青藏高原塊體隆升機制的基礎上，利用現今地表面GNSS觀測的隆升速率，以均衡理論為基礎獲得莫霍面的下沉速率，再通過計算地表面位移與莫霍面變化速率之和獲得地殼增厚速率。
考慮初始闆塊厚度，最後借此推算出兩個闆塊的初始碰撞時間，為達到此目的，利用現代大地測量特别是GNSS位移場資料，計算确定地表面隆升速率和莫霍面變化速率是重要的資料處理任務。
大地測量觀測資料和地殼模型
本研究主要采用GNSS觀測位移場資料确定構造垂直位移場，并采用CRUST1.0地殼模型确定莫霍面深度及其變化速率。
中國大陸構造環境監測網絡（簡稱“陸态網”），是以GNSS觀測為主，以甚長基線幹涉測量、人衛雷射測距等空間技術為輔，結合精密重力和水準測量等多種技術手段，由260個連續觀測和2000個不定期觀測站點構成的、覆寫中國大陸的高精度、高時空分辨率和自主研發資料處理系統的觀測網絡。
陸态網可監測大陸大陸岩石圈，近海和近地空間的物質結構和四維構造形态的變化，反映現今地殼運動和動力學的總體态勢，為本研究提供了極其重要的地表位移場，Wang等發表了中國大陸特别是青藏高原地區的首個水準位移場結果，揭示了青藏高原地區複雜的構造運動特征，具有先導意義。
Liang等根據陸态網GNSS觀測資料給出了中國青藏高原及周邊地區的垂直位移場的第一個結果，展示了青藏高原地區以隆升為主的特征，Pan等［27］利用陸态網資料對青藏高原隆升率問題進行了精細研究，基于GRACE時變重力對GNSS垂直位移場進行了負荷位移改正，并給出了青藏高原地區的經過負荷改正的水準與垂直位移場。
最近，Rao等對如何準确計算負荷改正問題進行研究，發現利用格林函數與水文負荷資料計算方法更為合理，并給出了最新的經過合理負荷改正的青藏高原垂直位移場，本文采用其研究方法和思路以獲得青藏高原表面構造位移場。
為了推算青藏高原塊體的增厚速率，我們需要知道地殼的厚度及其增厚率；增厚率可以由地表GNSS垂直位移以及莫霍面的變化率獲得，而地殼厚度可以由青藏高原高程資料以及莫霍面深度資料得到，為此，我們采用數字高程模型（DigitalElevationModel）和地殼模型資料CRUST1.0分别獲得青藏高原地區地表面高程資料，以及莫霍面深度的資料。
地形資料是通過有限的地形高程資料對真實地面地形的數字化表達，本文使用的地形資料來自美國奮進号航天飛機的雷達地形測繪（ShuttleRadarTopographyMission，SRTM）得到的SRTM15+衍生資料産品，分辨率為1°×1°。
全球地殼模型CRUST1.0于2013年釋出，分辨率為1°×1°，是目前為止較為詳細和系統的全球地殼模型，本文所使用的靜态莫霍面深度資料來自于CRUST1.0模型，它是根據全球最新地震研究資料得到的地殼厚度資料進行視窗平均得到。
地表面與莫霍面變形速率的估算
值得注意的是，由于GNSS觀測站點在整個青藏高原的分布是不均勻的，甚至是稀疏的，觀測得到的垂直位移場也是離散的，為了獲得整個青藏高原表面的GNSS垂直位移場，需要對觀測結果進行插值。
根據地理學第一定律：任何事物都是與其他事物相關的，隻不過相近的事物關聯更緊密，衍生出地學中最常用的2種插值方法——反距離權重插值（InverseDistanceWeighted，IDW）和克裡金插值（Kriging）。
IDW是簡便和常用的空間插值方法，它以插值點與樣本點的距離作為權重進行權重平均，權系數與距離成反比，樣本點離插值點越遠對插值點影響越小，權重函數為距離平方的倒數，在鄰近網格點内搜尋12個距離較短的資料進行計算。
克裡金插值法是一種基于統計學的插值方法，原理是利用區域化變量為基礎，以線性、高斯、球面、指數等變異函數為基本工具，對未知樣點進行線性無偏、最優化估計。
對于克裡金插值，必不可少的2個步驟為建立變異函數和協方差函數來估算取決于自相關模型的統計相關性值，進而預測未知值，克裡金插值方法給出了線性、高斯、球面、指數半變異函數模型，由于高原西部的GNSS測站點較少，使用高斯模型時會出現明顯的奇異點。
綜合考慮到IDW與克裡金插值方法的優缺點，本研究主要采用變異函數為球面模型的克裡金插值方法進行位移場插值計算。
莫霍面深度的變化率是根據均衡理論以及CRUST1.0模型資料計算得到的，具體的計算方法及結果在4.2節中給出。
如3.1節所述，根據Pan等給出的GNSSGRACE改正的青藏高原地區的垂直位移場和Rao等給出的GNSS-水文（GNSS-Hydrology）改正的垂直位移場為基礎理論進行插值計算。
這兩個結果的主要差異是在地表品質遷移負荷改正上的不同，前者利用GRACE球諧系數産品計算了負荷改正，僅反映了低階水文信号；而後者采用了格林函數積分方法計算了地表品質（含水文和構造）負荷效應，我們認為采用GNSS-水文改正的垂直位移場更為合理，而采用這兩種位移場資料開展研究，其目的是為了進行結果對比。
本小節的目的是把3.1節獲得的青藏高原GNSS站點經負荷改正的垂直位移，使用兩種插值方法，計算整個青藏高原地區的垂直位移場，所得結果如圖3所示。
圖3顯示，插值方法對結果具有明顯的影響，克裡金插值方法似乎更能突出局部效應，而IDW方法具有更強的平滑作用，這個現象說明選擇合适的插值方法是重要的；另外，GRACE與格林函數負荷改正的差別不明顯，可能是由于負荷改正相對于構造位移場仍然較小的原因。
我們利用均衡理論和地表位移場擷取莫霍面的變化速率，地殼均衡理論是地球實體學研究地殼的重要理論之一，其可靠性也由衆多的地球實體資料所證明，近年來在大地測量的研究中發揮了重要作用。
Parrt-Hayford模型與Airy-Heiskanen模型是最常見的2個地殼均衡理論模型，AiryHeiskanen均衡模型由于更接近地球實體原型，而被廣泛使用，是以，本研究采用該模型推算莫霍面的變化率，假設地殼垂直上升速率為T（mm/a），根據Airy-Heiskanen均衡模型，相應的莫霍面下降速率可由下式計算。
式中：ΔvMoho為莫霍面下降速率（mm/a），ρc為地殼密度（g/cm3），ρm為地幔密度（g/cm3），由地殼模型CRUST1.0求出的青藏高原的平均地殼密度為2.77g/cm3和地幔密度為3.37g/cm3，于是，利用公式（1），對上小節的地表面GNSS-負荷改正的位移場進行插值，就可以得到青藏高原莫霍面的變化（圖4）。
在前兩小節我們已經獲得了青藏高原地表與莫霍面的變化速率，将兩者求和就可以得到整個青藏高原區域地殼厚度的變化速率（圖5），該變化速率展現了青藏高原塊體增厚率，可用于推算青藏高原的初始隆升時間。
在前小節中我們推算了青藏高原地殼厚度的增厚率，為了證明推算結果的合理性，我們再通過青藏高原塊體的體積變化進行讨論，我們首先把青藏高原的各個子單元的增厚率進行求和，就可以獲得青藏高原塊體的體積變化率，結果表明，青藏高原塊體的體積變化率為2.5km3/a。
另一方面，我們通過計算青藏高原塊體的品質流入與流出的方法計算出體積變化率，即，從地殼物質平衡的角度考慮青藏高原的地殼物質體積變化。
根據Westaway的印度歐亞闆塊碰撞的體積流量平衡假設，印度地殼厚度為35km且印度與歐亞闆塊碰撞範圍為72°~97°E，碰撞接觸區域為正交于闆塊運動方向且長度約為2500km的區域，印度闆塊的移動速率為5cm/ar。
于是，印度闆塊每年流入歐亞闆塊的品質體積為4.4km3/a，其中經剝蝕而流出的品質體積約為1.7km3/a，剩餘的2.7km3/a體積增量則貢獻于青藏高原的地殼增厚。
綜合上述結果表明，本文計算的2.5km3/a體積變化率與Westaway的2.7km3/a體積增量基本符合，是以，我們估算的地殼增厚結果是可信的。
設Z是獨立觀測量x1，x2，…，xn的函數，即
式中：x1，x2，…，xn為直接觀測量，它們的中誤差為m1，m2，…，mn，則觀測值函數的中誤差為
式中：∂f/∂xi是Z對各變量xi的偏導數，mz是觀測函數的中誤差。
本文主要誤差來源于資料的誤差（包括GNSS測量資料以及DEM資料）和插值過程所産生的誤差，文中所使用為線性模型，則中誤差為。
式中：a1，a2，…，an均為常數。
印度闆塊與歐亞闆塊的初始碰撞時間結果
值得指出的是，如果上節中地殼增厚速度是線性的，我們就可以推算印度闆塊與歐亞闆塊的初始碰撞時間，但從另一方面考慮，若考慮增厚率為非線性的，且其他地質學、地球實體學等的方法與資料能夠得出青藏高原增厚規律，假設隆升的分階段性或加速減速趨勢特征。
根據我們推算的地殼增厚結果，重新推算其碰撞的初始時間是易于實作的，在本研究中，我們僅從大地測量技術的角度研究青藏高原初始碰撞時間推算，不考慮到整個青藏高原塊體隆升過程的複雜性（圖1），這是一種簡單且具有說服力的思路，本研究屬于初步探讨，在碰撞時間計算模拟中我們假設青藏高原地殼增厚率是線性的。
闆塊的初始碰撞指的是兩個大陸之間的大洋區域的岩石圈沿着縫合帶某一點或某一段開始俯沖并俯沖殆盡，造成兩側相連的大陸岩石圈直接接觸形成一個整體的岩石闆塊大陸。
青藏高原内部構造變化造成的地表形變效應可由大地測量手段測得，這就為利用大地測量手段估算高原的初始碰撞時間奠定了基礎，構造運動造成的地表效應可以展現在高原的地殼垂直形變上，這主要與GNSS觀測的垂直速度以及利用一些手段改正得到（如地表品質遷移的負荷效應改正）的構造運動垂直速度有關。
根據第2節中的基本思想和對青藏高原的變形與隆起機制的假設，估算初始碰撞時間的步驟如下。
首先，将青藏高原及其周邊區域（65°~110°E，20°~45°N）按照1°×1°的空間分辨率劃分格網，同時将莫霍面資料和地形資料STRM（ShuttleRadarTopographyMission，）也按照1°×1°分辨率劃分格網。
把GNSS位移場（未經任何負荷改正）、經GRACE負荷改正所得的地殼厚度增厚率，以及經水文模型負荷改正所得的地殼厚度增厚率，也都按照1°×1°分辨率劃分格網，其次，依據每個網格點，根據地殼厚度以及增厚率，計算出3類資料的開始增厚時間點，結果如圖6所示。
圖6中每個網格點呈現出明顯的波動，一方面可能是由GNSS觀測精度、DEM資料精度和插值方法所緻，也可能由于地表剝蝕分布不均勻所緻，每個網格點的波動性背後包含有非常複雜的造山運動和構造運動的原因。
還需要進一步分析，青藏高原的隆升本身就是一個分區域分階段的複雜過程，對于高原的隆升的初始時間仍然是學術界一個極具争議的問題，即使是上百萬年尺度的地質學研究也無法給出準确結果。
圖1所展示的地質學結果也具有很強的不連續性，如喜馬拉雅和印度—亞洲縫合線以及拉薩—羌塘縫合谷也出現了隆起時間差異很大，是以網格之間出現波動是合理的，為了得到總體隆升時間起點，我們隻要把圖6中整個青藏高原塊體進行平均即可，以求得印度闆塊與歐亞闆塊的初始碰撞時間，結果如表1所列。
表1是本研究的最終結果，由此可見，印度闆塊與歐亞闆塊的初始碰撞時間範圍為36~100Ma，表1中利用GNSS直接測量位移場所對應的結果，對比其他兩組經負荷改正的結果總體上偏小（分别為40Ma和36Ma）。
但是，這似乎與Najman等所給出的兩闆塊碰撞發生在40Ma的結論相一緻，利用經過GRACE負荷改正GNSS位移場所得到的初始碰撞時間的範圍為51~65Ma，這與現在科學界所普遍認為的兩個闆塊的初始碰撞時間的範圍較為符合。
而對于利用水文模型負荷改正GNSS位移場所得到的闆塊初始碰撞時間為83~100Ma，較為符合目前最長碰撞時間的結果，如圖1所示的利用拉薩—羌塘縫合谷所得到的結果（95~100Ma），總之，表1中的結果基本上符合現有的由地質學所得的結果範圍。
表1顯示，使用克裡金插值方法與IDW進行插值所得結果具有一定差異，達4~17Ma，這個差異展現了插值方法所帶來的誤差或不确定性，一般而言，IDW方法計算速度快，但是，克裡金插值精度要略高于IDW，考慮到這個因素，我們認為使用克裡金插值方法更好。
另外，在上述3類（GNSS、GNSS-GRACE和GNSS-Hydrology）結果中，GNSS-Hydrology結果，即83Ma更合理一些，因為GNSS-Hydrology結果是對GNSS觀測位移場做了更為合理的負荷改正，更為真實地反映了固體地殼的變形場，而GNSS直接測量位移場包含了地表品質遷移的負荷變形資訊，很容易了解其結果的不合理性。
GNSS-GRACE結果雖然對應于GRACE負荷改正，但是由于其截斷效應沒能反映出負荷變形的全頻域資訊，其合理性也是不充分的，資料的誤差（包括GNSS測量資料以及DEM資料）和插值過程所産生的誤差，根據線性函數誤差傳播，估算的隆起時間的中誤差大緻在±16Ma。
是以，綜合上述讨論，我們認為，利用大地測量觀測技術确定的印度闆塊與歐亞闆塊的初始碰撞時間為（83±16Ma），這個結果大于圖1中的60~65Ma的碰撞時間，小于其利用LQV所得到的95~100Ma的結果，總體上介于地質學結果的合理範圍内。
值得注意的是，相對于上述地質學結果，大地測量結果（83±16Ma）似乎偏大，然而，該值偏大具有一定的合理性，最直接的原因是本研究假設青藏高原是線性隆升的，當然也有上述讨論的觀測或插值誤差的原因。
如果考慮青藏高原是非線性隆升的，其特性與圖1所顯示的一樣，早期隆升速率較快而晚期隆升速率較慢，則印度闆塊與歐亞闆塊的初始碰撞時間要比83Ma更晚，更加接近于地質學的公認結果。
在上述研究中我們沒考慮冰GIA效應的影響，現在對此加以讨論，GIA是一種跟負荷相關的地球實體現象，它是指黏彈地球對始于末次冰期地表冰和海水負荷的響應。
GIA效應影響着冰後地殼運動、全球海平面變化、地球重力場、地球旋轉運動和應力狀态等方面，為研究地幔流變、末次冰期冰蓋厚度、衛星測高、衛星重力、全球氣候變暖以及闆塊運動、地殼垂直運動等提供改正。
科學界提出了一系列的GIA模型，本研究選取了以三維有限元可壓縮地球模型為基礎的Geruo13模型進行改正。
利用與上述相同的方法，我們計算了考慮GIA效應後的印度闆塊與歐亞闆塊的初始碰撞時間，其結果表明，在考慮GIA效應後，無論利用GNSS-GRACE還是GNSS-Hydrology所計算得到的初始碰撞時間差異與不考慮GIA效應時所得到的結果差異不大，大約在1Ma。
這可能是因為最近的一次GIA效應開始于10萬年前與闆塊碰撞的百萬年尺度比起來量級相對較小，故對計算闆塊的初始碰撞時間影響不大。