文 | 唐哲
编辑 | 唐哲
大地测量方法推算板块碰撞起始时间的基本思想
为了推算印度板块和欧亚板块的初始碰撞时间，可以假设如下物理模型：欧亚板块地壳（青藏高原块体）是一个弹性块体，在印度板块持续冲撞下发生弹性变形，根据地壳均衡理论，地表面隆升而莫霍面下沉，青藏高原块体总体上处于动态均衡状态且不断增厚，上述隆升机制和地壳增厚模型与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万年前与板块碰撞的百万年尺度比起来量级相对较小，故对计算板块的初始碰撞时间影响不大。