球面圓弧數字地形剖面提取方法研究
——以喜馬拉雅弧形造山帶為例1

張 玲 梁詩明

（中國地震局地質研究所, 地震動力學國家重點實驗室, 北京 100029）

引言

地形分析以地形、地質圖、數字圖像資料及詳細野外地質調查等為研究基礎，是地質地貌研究的重要手段，在新構造研究中占有舉足輕重的地位（張會平等，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 不同投影方式引起的非構造變形示意（以穿過喜馬拉雅弧形造山帶的環形為例）Fig. 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 喜馬拉雅造山帶位置及條帶狀地形剖面提取結果Fig. 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）。

計算過程中，確定每個計算單元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′地理坐標。

圖3 球面圓弧帶狀地形剖面計算原理示意Fig. 3 The schematic diagram of the mathematical procedure of arc swath topographic based on ellipsoid model of globe

以A′地理坐標計算為例，說明不斷向前推進地地形計算單元端點地理坐標的計算方法。球面上點A′的地理坐標是由相應的任意點A 在x，y，z軸上的分解運動旋轉變換得到的，僅考慮靜止的慣性參考系，點A的歐拉角（ φM,λM,Ω） 是已知的，其中， φM、 λM分 別為點A 經、緯度坐標， Ω為計算步長，即為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 結語

進行地形分析研究時，帶狀地形剖面分析可直觀地反映地球表面形態。通過與構造地質、沉積學和年代學等研究方法相結合，可深入分析造山帶變形、盆山耦合關系和動力學機制等問題。其中，可反映地形剖面平面位置信息的剖面投影平面是重要的地理要素。本文提出的球面圓弧數字地形剖面提取方法計算過程中，無須分塊計算平面位置信息的剖面投影平面，可避免提取地殼尺度大型構造帶地形剖面時的投影問題，剔除由此引入的非構造變形。與傳統方法相比，提取結果更接近真實地形，體現更具細節的地貌特征，以便更加深入地開展地學研究。同時，本文提出的方法也適用于小型構造。由于計算過程簡單，可高效地以不同尺度對地貌特征展開研究。