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

張 玲 梁詩明

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

引言

地形分析以地形、地質(zhì)圖、數(shù)字圖像資料及詳細野外地質(zhì)調(diào)查等為研究基礎(chǔ)，是地質(zhì)地貌研究的重要手段，在新構(gòu)造研究中占有舉足輕重的地位（張會平等，2004）。通過垂直于地表面的截面切割地面，以反映地面起伏曲線的圖形即為帶狀地形剖面圖（王春范，2018）。由于帶狀地形剖面圖可提供更廣泛的大地形態(tài)特征視圖，有助于識別活動構(gòu)造中指示地貌特征的異?；虻匦乌厔荩℅rohmann，2004），是道路、管線、勘探等工程設(shè)計的基礎(chǔ)文件之一，在工程領(lǐng)域得到廣泛應(yīng)用。通常情況下，帶狀地形剖面的實現(xiàn)是通過連續(xù)步驟將數(shù)字地形數(shù)據(jù)中的點投影到具有一定寬度的樣條中（Fielding 等，1994；Burbank 等，2013），投影至每個條帶中的數(shù)據(jù)可用于統(tǒng)計分析，以定義其屬性，通常包括投影點的最大、平均、最小海拔高度（Burbank 等，2013）。隨著高精度數(shù)字地形數(shù)據(jù)的日益普及，很多GIS 軟件已具備自動生成功能。然而，這種方法忽略了投影方式引入的非構(gòu)造變形，如圖1 所示。在圖1（a）球面橢球坐標(biāo)系統(tǒng)中繪制了2 個標(biāo)準(zhǔn)圓，它們的南部間隙約束了喜馬拉雅弧形造山帶的南北邊界。將圖1（a）中的2 個標(biāo)準(zhǔn)圓利用圓柱等積投影至平面坐標(biāo)系后，它們在橫向和縱向均發(fā)生了變形（圖1（b））；將圖1（a）中的2 個標(biāo)準(zhǔn)圓利用橫軸墨卡托投影至平面坐標(biāo)系后同樣發(fā)生了畸變（圖1（c）），2 個標(biāo)準(zhǔn)圓之間的間隙在高緯度地區(qū)變大，在低緯度地區(qū)變小，越靠近高緯度地區(qū)畸變越大。因此，高精度數(shù)字地形數(shù)據(jù)特別適用于構(gòu)造變形強度較小的小型地質(zhì)構(gòu)造。對于類似于喜馬拉雅弧形造山帶等大型構(gòu)造，在球面坐標(biāo)系下獲得的帶狀地形剖面相比于通過上述傳統(tǒng)方法得到的結(jié)果，更接近真實地形，具有更多的構(gòu)造地貌特征細節(jié)。

圖1 不同投影方式引起的非構(gòu)造變形示意（以穿過喜馬拉雅弧形造山帶的環(huán)形為例）Fig. 1 The unreal deformation introduced by different projections (we take the annulus covering the Himalayan orogenic belt as an example)

筆者在實際工作中，探索了基于球面坐標(biāo)系統(tǒng)制作帶狀地形剖面圖的方法。本文以典型的喜馬拉雅弧形造山帶為例，利用SRTM 90 m 分辨率的DEM 數(shù)據(jù)，通過Python 語言實現(xiàn)了提取過程，并對制作結(jié)果進行展示。

1 提取流程

喜馬拉雅弧形造山帶是地學(xué)界研究陸陸碰撞和板塊俯沖模型等重大構(gòu)造變形問題的經(jīng)典地區(qū)之一。喜馬拉雅造山帶的地震活動性、區(qū)域應(yīng)變、地形起伏和河流裂點等的空間分布特征均印證了其“完美”的弧形分布形態(tài)，如圖2（a）所示（Bendick 等，2001）。Bendick 等（2001）曾采用最小二乘法確定出最佳球圓函數(shù)以匹配上述數(shù)據(jù)，計算結(jié)果表明，各數(shù)據(jù)擬合的球面圓具有高度相似性，其平均圓心位置為（（42.4°±2.1°）N，（91.6°±1.6°）E），平均半徑為（1 696±55）km。

圖2 喜馬拉雅造山帶位置及條帶狀地形剖面提取結(jié)果Fig. 2 Location and the topographic swath profile of the Himalayan orogenic belt

本研究將不規(guī)則地形-喜馬拉雅東構(gòu)造結(jié)引入“完美”弧線，以突出地形變化。使用最小二乘法，以喜馬拉雅造山帶內(nèi)最高地形點作為約束，確定了最小化誤差平方和情況下的最佳球面圓參數(shù)，圓心為（40.89°N，88.8°E），半徑為1 500 km，這接近于喜馬拉雅弧4 000 m 等高線代表的區(qū)域應(yīng)力發(fā)生轉(zhuǎn)換的邊界（Molnar 等，1989）。

以已有學(xué)者估計的弧形山脈小圓半徑1 200 km 為內(nèi)徑，1 540 km 為外徑，以弧形山脈西端點（33.40°N，75.75°E）為計算起點，構(gòu)建覆蓋喜馬拉雅弧形造山帶的球面圓弧帶狀地形剖面（圖2（a）、（b））。其中，帶狀地形剖面外徑和計算正方向分別定義為1 540 km、北方向。每隔0.2°作為1 個階躍劃分地形計算單元，分別進行最大、平均、最小高程統(tǒng)計（圖3）。

計算過程中，確定每個計算單元4 個端點的地理坐標(biāo)至關(guān)重要（圖3）。由于喜馬拉雅弧內(nèi)、外徑已知，以點A（34.02°N，76.68°E）和點B（32.77°N，74.85°E）為起點，以球面上的點M（91.6°E，42.4°N）為圓點，進行角度旋轉(zhuǎn)，計算得到四邊形端點A′、B′地理坐標(biāo)。

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

以A′地理坐標(biāo)計算為例，說明不斷向前推進地地形計算單元端點地理坐標(biāo)的計算方法。球面上點A′的地理坐標(biāo)是由相應(yīng)的任意點A 在x，y，z軸上的分解運動旋轉(zhuǎn)變換得到的，僅考慮靜止的慣性參考系，點A的歐拉角（ φM,λM,Ω） 是已知的，其中， φM、 λM分 別為點A 經(jīng)、緯度坐標(biāo)， Ω為計算步長，即為0.2°。歐拉角可通過下式轉(zhuǎn)化為空間笛卡爾坐標(biāo)系（Oxyz）中的等效旋轉(zhuǎn)：

計算得到歐拉角的矩陣：

將點A 的對應(yīng)值圍繞圓心點M 以一定角度進行旋轉(zhuǎn)，即可得到點A′的空間直角坐標(biāo)：

根據(jù)同樣的方法，求取點B′的坐標(biāo)。至此，即可獲得地形計算單元點A、A′、B、B′ 坐標(biāo)和區(qū)域范圍。

當(dāng)將點A′、B′坐標(biāo)作為已知時，利用相同的方法，以0.2°的步長可計算得到點A″、B″坐標(biāo)，進而逐步計算得到弧形地形條帶上所有計算單元邊界的端點。在每個地形計算單元中，搜索提取最大、平均、最小高程。以計算步長0.2°為橫軸，以統(tǒng)計結(jié)果（最大、平均、最小高程）為縱軸，投影繪圖即可獲得帶狀地形剖面。

2 結(jié)果分析

由喜馬拉雅弧形造山帶地形剖面提取結(jié)果可得地形變化特征（圖2（b）、（c））：①在30°～80°弧度區(qū)間內(nèi)，印度板塊與歐亞板塊為正向碰撞，而向兩側(cè)斜向碰撞角度逐漸增大。同樣地，地形（最大）在前文提到的中部“完美”弧形地段基本保持穩(wěn)定，而在兩側(cè)的東、西構(gòu)造結(jié)處顯示出明顯的下降。②在弧度30°位置處，從平均地形中可見明顯的下降。該位置對應(yīng)Leo Pargil 裂谷所在位置，是青藏高原內(nèi)部東向遷移和撕裂的起始位置（Bian 等，2020）。根據(jù)相同的原理，可從地形剖面中平均地形突然發(fā)生下降的點位識別出Thakkola 裂谷、孔錯裂谷和亞東裂谷。③弧度80°所在區(qū)域是地形特征十分特殊的區(qū)域，既包含最大高程超過7 000 m 的南迦巴瓦峰和加拉白壘峰，又是平均地形相對于兩側(cè)急速下降的區(qū)域。這一地區(qū)屬于喜馬拉雅東構(gòu)造結(jié)的前緣地帶，是整個喜馬拉雅造山帶中構(gòu)造應(yīng)力最強烈、隆升和剝蝕速率最快、新生代巖漿活動、深熔作用和變質(zhì)作用最強的地區(qū)（Ding 等，2001；丁林等，2013）。獨特的地形地貌特征為該地區(qū)變形和動力學(xué)機制提供依據(jù)。

3 結(jié)語

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