引言

地形分析以地形、地质图、数字图像资料及详细野外地质调查等为研究基础，是地质地貌研究的重要手段，在新构造研究中占有举足轻重的地位（张会平等，2004）。通过垂直于地表面的截面切割地面，以反映地面起伏曲线的图形即为带状地形剖面图（王春范，2018）。由于带状地形剖面图可提供更广泛的大地形态特征视图，有助于识别活动构造中指示地貌特征的异常或地形趋势（Grohmann，2004），是道路、管线、勘探等工程设计的基础文件之一，在工程领域得到广泛应用。通常情况下，带状地形剖面的实现是通过连续步骤将数字地形数据中的点投影到具有一定宽度的样条中（Fielding等，1994；Burbank等，2013），投影至每个条带中的数据可用于统计分析，以定义其属性，通常包括投影点的最大、平均、最小海拔高度（Burbank等，2013）。随着高精度数字地形数据的日益普及，很多GIS软件已具备自动生成功能。然而，这种方法忽略了投影方式引入的非构造变形，如图1所示。在图1（a）球面椭球坐标系统中绘制了2个标准圆，它们的南部间隙约束了喜马拉雅弧形造山带的南北边界。将图1（a）中的2个标准圆利用圆柱等积投影至平面坐标系后，它们在横向和纵向均发生了变形（图1（b））；将图1（a）中的2个标准圆利用横轴墨卡托投影至平面坐标系后同样发生了畸变（图1（c）），2个标准圆之间的间隙在高纬度地区变大，在低纬度地区变小，越靠近高纬度地区畸变越大。因此，高精度数字地形数据特别适用于构造变形强度较小的小型地质构造。对于类似于喜马拉雅弧形造山带等大型构造，在球面坐标系下获得的带状地形剖面相比于通过上述传统方法得到的结果，更接近真实地形，具有更多的构造地貌特征细节。

图 1  不同投影方式引起的非构造变形示意（以穿过喜马拉雅弧形造山带的环形为例）

Figure 1.  The unreal deformation introduced by different projections (we take the annulus covering the Himalayan orogenic belt as an example)

下载: 全尺寸图片 幻灯片

笔者在实际工作中，探索了基于球面坐标系统制作带状地形剖面图的方法。本文以典型的喜马拉雅弧形造山带为例，利用SRTM 90 m分辨率的DEM数据，通过Python语言实现了提取过程，并对制作结果进行展示。

1.   提取流程

喜马拉雅弧形造山带是地学界研究陆陆碰撞和板块俯冲模型等重大构造变形问题的经典地区之一。喜马拉雅造山带的地震活动性、区域应变、地形起伏和河流裂点等的空间分布特征均印证了其“完美”的弧形分布形态，如图2（a）所示（Bendick等，2001）。Bendick等（2001）曾采用最小二乘法确定出最佳球圆函数以匹配上述数据，计算结果表明，各数据拟合的球面圆具有高度相似性，其平均圆心位置为（（42.4°±2.1°）N，（91.6°±1.6°）E），平均半径为（1 696±55）km。

图 2  喜马拉雅造山带位置及条带状地形剖面提取结果

Figure 2.  Location and the topographic swath profile of the Himalayan orogenic belt

下载: 全尺寸图片 幻灯片

本研究将不规则地形——喜马拉雅东构造结引入“完美”弧线，以突出地形变化。使用最小二乘法，以喜马拉雅造山带内最高地形点作为约束，确定了最小化误差平方和情况下的最佳球面圆参数，圆心为（40.89°N，88.8°E），半径为1 500 km，这接近于喜马拉雅弧4 000 m等高线代表的区域应力发生转换的边界（Molnar等，1989）。

以已有学者估计的弧形山脉小圆半径1 200 km为内径，1 540 km为外径，以弧形山脉西端点（33.40°N，75.75°E）为计算起点，构建覆盖喜马拉雅弧形造山带的球面圆弧带状地形剖面（图2（a）、（b））。其中，带状地形剖面外径和计算正方向分别定义为1 540 km、北方向。每隔0.2°作为1个阶跃划分地形计算单元，分别进行最大、平均、最小高程统计（图3）。

图 3  球面圆弧带状地形剖面计算原理示意

Figure 3.  The schematic diagram of the mathematical procedure of arc swath topographic based on ellipsoid model of globe

下载: 全尺寸图片 幻灯片

计算过程中，确定每个计算单元4个端点的地理坐标至关重要（图3）。由于喜马拉雅弧内、外径已知，以点A（34.02°N，76.68°E）和点B（32.77°N，74.85°E）为起点，以球面上的点M（91.6°E，42.4°N）为圆点，进行角度旋转，计算得到四边形端点A′、B′地理坐标。

以A′地理坐标计算为例，说明不断向前推进地地形计算单元端点地理坐标的计算方法。球面上点A′的地理坐标是由相应的任意点A在x，y，z轴上的分解运动旋转变换得到的，仅考虑静止的惯性参考系，点A的欧拉角（$ {{\varphi }_{\text{M}},\lambda }_{\text{M}},\varOmega $）是已知的，其中，$ {{\varphi }}_{\text{M}} $、$ {\lambda }_{\text{M}} $分别为点A经、纬度坐标，$ \varOmega $为计算步长，即为0.2°。欧拉角可通过下式转化为空间笛卡尔坐标系（Oxyz）中的等效旋转：

计算得到欧拉角的矩阵：

将点A的对应值围绕圆心点M以一定角度进行旋转，即可得到点A′的空间直角坐标：

根据同样的方法，求取点B′的坐标。至此，即可获得地形计算单元点A、A′、B、B′ 坐标和区域范围。

当将点A′、B′坐标作为已知时，利用相同的方法，以0.2°的步长可计算得到点A″、B″坐标，进而逐步计算得到弧形地形条带上所有计算单元边界的端点。在每个地形计算单元中，搜索提取最大、平均、最小高程。以计算步长0.2°为横轴，以统计结果（最大、平均、最小高程）为纵轴，投影绘图即可获得带状地形剖面。

2.   结果分析

由喜马拉雅弧形造山带地形剖面提取结果可得地形变化特征（图2（b）、（c））：①在30°～80°弧度区间内，印度板块与欧亚板块为正向碰撞，而向两侧斜向碰撞角度逐渐增大。同样地，地形（最大）在前文提到的中部“完美”弧形地段基本保持稳定，而在两侧的东、西构造结处显示出明显的下降。②在弧度30°位置处，从平均地形中可见明显的下降。该位置对应Leo Pargil裂谷所在位置，是青藏高原内部东向迁移和撕裂的起始位置（Bian等，2020）。根据相同的原理，可从地形剖面中平均地形突然发生下降的点位识别出Thakkola裂谷、孔错裂谷和亚东裂谷。③弧度80°所在区域是地形特征十分特殊的区域，既包含最大高程超过7 000 m的南迦巴瓦峰和加拉白垒峰，又是平均地形相对于两侧急速下降的区域。这一地区属于喜马拉雅东构造结的前缘地带，是整个喜马拉雅造山带中构造应力最强烈、隆升和剥蚀速率最快、新生代岩浆活动、深熔作用和变质作用最强的地区（Ding等，2001；丁林等，2013）。独特的地形地貌特征为该地区变形和动力学机制提供依据。

3.   结语

进行地形分析研究时，带状地形剖面分析可直观地反映地球表面形态。通过与构造地质、沉积学和年代学等研究方法相结合，可深入分析造山带变形、盆山耦合关系和动力学机制等问题。其中，可反映地形剖面平面位置信息的剖面投影平面是重要的地理要素。本文提出的球面圆弧数字地形剖面提取方法计算过程中，无须分块计算平面位置信息的剖面投影平面，可避免提取地壳尺度大型构造带地形剖面时的投影问题，剔除由此引入的非构造变形。与传统方法相比，提取结果更接近真实地形，体现更具细节的地貌特征，以便更加深入地开展地学研究。同时，本文提出的方法也适用于小型构造。由于计算过程简单，可高效地以不同尺度对地貌特征展开研究。