shipxy、shipxycom

《测绘学报》

构建与学术的桥梁 拉近与权威的距离

利用多元回归分析反演西南印度洋区域海底地形

范雕¹, 李姗姗¹ , 杨军军^2,3, 孟书宇⁴, 邢志斌⁵, 张驰¹, 冯进凯¹

1. 信息工程大学, 河南 郑州 450001;

2. 华中科技大学精密重力测量国家重大科技基础设施地球物理研究所, 湖北 武汉 430074;

3. 华中科技大学物理学院基本物理量测量教育部重点实验室, 湖北 武汉 430074;

4. 32022部队, 湖北 武汉 430074;

5. 航天工程大学, 北京 102206

收稿日期：2018-11-16；修回日期：2019-08-22

基金项目：国家自然科学基金（41774021；41774018；41504018；41674026；41674082；41574020）；国家重点研发计划（2016YFB0501702）；地理信息工程国家重点实验室开放基金（SKLGIE2016-M-3-2）

第一作者简介：范雕(1991-), 男, 博士生, 研究方向为物理大地测量和空间化海洋测绘工程。E-mail:fandiao2311@mails.jlu.edu.cn

通信作者：李姗姗, E-mail:zzy_lily@

摘要：针对海底地形与重力异常和重力异常垂直梯度在相应频段呈现强线性相关的特点，引入多元回归分析技术，提出并详细推导了联合多元重力数据的海底地形建模方法。然后，在西南印度洋SWIR（Southwest India Ridge）所在部分海域开展了海底地形反演试验及地形地貌分析研究。试验结果表明：6种海深模型中，基于多元回归分析技术构建的海深模型（BDVG模型）检核精度最高，相较于S&S V18.1模型和ETOPO1模型精度分别提高了11.51%和57.81%左右；2000 m以上水深海域，各个海深模型的检核精度较高，相对误差波动较小，反映了深海海域具有良好的反演效果；地形起伏剧烈海域或者浅海海域，BDVG海深模型，相较于以重力异常和重力异常垂直梯度作为单一输入源建立的BDG模型和BVGG模型相对误差及相对误差波动变化较小，反映了BDVG模型拥有更好的稳定性，从而体现了联合反演的必要性和优势。Indomed FZ-Gallieni FZ上唯一轴部缺失裂谷洋脊段（27洋脊段）目前属于岩浆供应充足阶段，构造作用的海底扩张对其影响较小；同时由于对称裂离方式影响，27洋脊段沿轴南北对称分布有地形隆起。

关键词：多元回归分析 海底地形 重力异常 重力异常垂直梯度 相干性

Predicting bathymetry by applying multiple regression analysis in the Southwest Indian Ocean Region

FAN Diao¹, LI Shanshan¹, YANG Junjun^2,3, MENG Shuyu⁴, XING Zhibin⁵, ZHANG Chi¹, FENG Jinkai¹

1. Information Engineering University, Zhengzhou 450001, China;

2. Institute of Geophysics and PGMF, Huazhong University of Science and Technology, Wuhan 430074, China;

3. MOE Key Laboratory of Fundamental Physical Quantities Measurement & Hubei Key Laboratory of Gravitation and Quantum Physics, PGMF and School of Physics, Huazhong University of Science and Technology, Wuhan 430074, China;

4. 32022 Troops, Wuhan 430074, China;

5. Space Engineering University, Beijing 102206, China

Foundation support: The National Natural Science Foundation of China(Nos. 41774021; 41774018; 41504018; 41674026; 41674082; 41574020); The National Key Research and Development Program of China(No. 2016YFB0501702); The Fund of State Key Laboratory of Geo-Information Engineering(No. SKLGIE2016-M-3-2)

First author: FAN Diao (1991—), male, PhD candidate, majors in physical geodesy and spatial marine surveying and mapping engineering.E-mail:fandiao2311@mails.jlu.edu.cn.

Corresponding author: LI Shanshan, E-mail:zzy_lily@.

Abstract: Considering the fact that the sea floor topography and gravity anomaly or vertical gravity gradient anomaly show strong linear correlation in the corresponding frequency bands, the method based on using multivariate regression analysis technique to combine multi-gravity data to construct the seafloor model was proposed. Then, the inversion test and analysis were carried out in the part of SWIR(Southwest India Ridge) in the Southwest Indian Ocean. The results showed that the bathymetry model (BDVG model) based on multiple regression analysis has the highest accuracy compared with other models, which is 11.51% and 57.81% higher than the S&S V18.1 model and ETOPO1 model respectively. The accuracy of each bathymetry model is higher, and the relative error fluctuation is small where the water depth is above 2000 m, reflecting the good inversion effect in the deep sea area. In places where the seafloor is fluctuated drastically or in shallow seaarea, BDVG model has less variation in relative error and relative error fluctuation than BDG model and BVGG model established by gravity anomaly and vertical gravity gradient anomaly as a single input source, reflecting the BDVG model has better stability and the necessity and advantage of joint inversion. The only shaft-deficient rift oceanic ridge section (27 oceanic ridge section) on the Indomed FZ-Gallieni FZ is currently in the stage of sufficient magma supply, and the seafloor expansion has less influence on it. At the same time, due to the influence of the symmetric splitting, several rises are symmetrically distributed along the north and south of the axis.

Key words: multiple regression analysis seafloor topography gravity anomaly vertical gravity gradient coherence

海洋测量按照内容和任务可分为海洋大地测量、海道测量、海岸地形测量、海底地形测量、海洋重力测量、海洋磁力测量、海洋工程测量和海洋遥感测量等。其中，海底地形测量是海洋测量最核心、最具特色的海洋信息获取方式^[1]。海底地形反映了海床的起伏变化，在导弹坐底发射、水下匹配导航、海底板块运动、沉船打捞、海洋油气勘探和海洋环境监测等方面发挥着重要作用^[2]。

目前海底地形测量技术主要包括船基声呐测深技术、机载激光雷达测量技术(airborne LiDAR bathymetry，ALB)、海岸带一体化地形测量技术、潜基海底地形测量技术等。其中船基海底地形测量主要借助船舶搭载多波束/单波束回声测深仪开展水深测量^[3]，是目前海洋地理空间信息获取的主要方式；机载激光雷达测量技术主要作业方式在50 m的水域，是船基测量的重要补充；潜基海底地形测量以AUV(autonomous underwater vehicle，AUV)/ROV(remote operated vehicle，ROV)/深拖系统为平台，携载多波束测深系统、侧扫声呐系统、压力传感器等对海域进行测量，以提高海底地形地貌的分辨率和精度，从而满足科学研究和工程应用需求。

然而上述海底地形测量作业模式的测量区域范围相对于全球海域而言，数据覆盖范围十分有限，文献[4]认为深海测绘将花费200多船年(ship-year)时间，耗费大量人力物力才能绘制完成。同时以船舶为载体进行测量的测量结果分布很不均匀，特别在纬度高于30°S的海域，船测数据十分稀疏，存在大面积的数据空白。随着卫星测高技术的发展，特别执行大地测量任务(geodetic mission，GM)的Geosat等测高卫星的数据积累，卫星测高数据现已能够完全覆盖这些海区^[5-6]。应用卫星测高技术获得的覆盖全球且均匀分布的海面高数据，依据大地测量方法可恢复全球海洋重力场信息，并且随着测高和定轨精度的大大提高，测高数据不断更新并大量积累，海洋重力场模型也不断优化改进^[7]。另外，海面重力信息主要由海水质量亏损、地球质量异常和海深结构、地壳以下质量均衡补偿引起，其中海水质量亏损影响的贡献占主要，其余因素对海面重力信息的贡献由于向上延拓而被平滑。基于此，利用海深与海面重力信息的相关性特点^[8-9]，采用相应反演方法可获得深海大洋大尺度海底地形，进而可形成航空、航天、地面、水面、水下五位一体海底地形测量手段，构成多元立体测量平台，达到立体、高效、高精度观测海洋目的，为海洋基础地理信息获取提供多种可选择的技术手段，完善海底地形测量全要素多平台测量体系。

凭借海面重力数据与海底地形在某些波段高度相关的特点开展海底地形预测最早由文献[10]提出。文献[11]证明了该方法的可行性。随后国内外学者进行了大量的试验研究，文献[12]研究了大地水准面高反演海底地形的技术途径，并详细分析了不同参数对转换函数的影响。文献[13]研究了应用向下延拓方法和最小二乘配置方法反演海底地形的技术手段。文献[14]回顾了主要的海底地形反演方法，包括线性滤波方法和导纳函数方法等。文献[15]在南极德雷克海峡(Drake Passage)采用自由空间重力异常(free-air gravity anomalies，FAGA)，应用重力地质方法(gravity-geologic method，GGM)开展了反演试验。试验结果表明：应用GGM方法单纯将船测数据格网化检核精度更高。文献[16]推导了利用重力异常垂直梯度的海底地形反演方法。文献[17]基于重力异常垂直梯度反演了中国南海2′的海底地形模型。文献[18]采用线性方法分别反演了以重力异常和重力异常梯度为输入源的海底地形模型，取得了良好的反演效果。文献[19]联合Geosat、ERS-1/2、T/P、Jason-1等测高卫星，反演了中国南海海域的重力异常，并以此为基础采用解析法反演了中国南海海底地形。文献[20—21]利用重力异常在皇帝海山链展开了试验分析，获得了良好的反演效果。经以上分析可知，国内外学者的研究主要集中在单独利用重力异常、重力异常垂直梯度等重力数据构建海底地形模型。然而一个现象往往与多个因素相联系，海底地形变化在重力异常和重力异常垂直梯度中均会有所反映。因此，建立重力异常和重力异常垂直梯度的最优组合共同预测海底地形，比单一利用重力异常或者重力异常垂直梯度构建海底地形模型似乎更加有效，也更符合实际。例如，文献[22]分别以重力异常和重力异常垂直梯度采用线性方法反演不同波段海底地形，进而将不同波段反演结果求和获得最终的海底地形模型，模型精度相较于V15.1模型提高了29.5%左右；文献[23]基于自适应赋权技术赋予重力异常和重力异常垂直梯度反演海深模型最佳权值组合，进而构建最终的海底地形模型，检核精度结果相较于单一输入源有较大提高，该方法的缺点是需要外部控制点作为约束方能获得不同输入源的有效权比。

然而上述融合方法本质上是以重力异常或者重力异常垂直梯度作为单一输入源，从而获得相应的不同输入源的反演结果，然后将不同输入源的反演结果进行某种组合，所以实质上仍是在技术层面达到联合反演海底地形目的。有别于此，本文基于重力异常或者重力异常垂直梯度经滤波和延拓处理后结果与海底地形表现明显线性相关特点，提出基于多元回归分析技术，以重力异常和重力异常垂直梯度为数据源，联合重力异常和重力异常垂直梯度的海底地形建模方法(该方法不需外部控制点约束即可获得不同输入数据的最佳权重组合)。并将联合重力数据反演结果与目前国际应用广泛的Smith and Sandwell V18.1海深模型(简称S & S V18.1模型)、ETOPO1海深模型进行精度比较、评估模型效能。最后以本文建立的海底地形模型，基于岩浆-构造旋回模型研究分析相关海域的地形地貌及地质构造特点。

1 原理与方法

图 1为质量亏损与重力异常的关系。如图 1所示，以区域平均海深面作为O-xy面，以指向海面为z轴建立笛卡尔坐标系。其中海面一点P的坐标为(x,y,z₀)，区域平均水深为d，海底地形高于平均海深面的高度为b(向上为正)，对应的海洋深度为h(向下为负)，r为质量体指向P点的向径，海水密度和洋壳密度分别为ρ_w和ρ_c。根据海底地形和海面重力数据的关系，可以将其关系描述为关于地形波长的函数^[24-25]。频率域内重力异常和海底地形的关系利用Parker公式可表示为^[26]

图 1 质量亏损与重力异常的关系Fig. 1 Relationship between mass loss and gravity anomaly

图选项

(1)

式中，F(Δg(x,y,z₀))表示海面重力异常的傅里叶变换；G为地球引力常数，通常取6.672×10^-8cm³/(g·s²)；Δρ为洋壳与海水的密度差异；Z₀为研究海域参考海深；f为径向频率(

文献[27]研究发现，式(1)在平坦海域高次项的影响较小，通常情况下, 主要考虑式(1)中一次项的贡献。当仅考虑一次项时，式(1)可整理为

(2)

令

(3)

式中，(2πGΔρ)为布格常数(Bouguer constant)。则式(2)可理解为若以海底地形作为输入，经过函数Z(f_x, f_y)的转换可以输出重力异常。定义函数Z(f_x, f_y)为表征将海底地形转换为重力异常能力的重力导纳函数^[28-29]，式(2)可表示为

(4)

导纳函数Z(f_x, f_y)中由于向下延拓因子exp(-2πfz₀)的衰减效应，随着波数的增大，导纳函数值将减小，从而延拓因子exp(-2πfz₀)具有平滑正演结果的作用；同时重力异常结果的振幅也将比海底地形振幅的(2πGΔρ)倍小，导纳函数Z(f_x, f_y)曲线如图 2黑色实线所示。转换函数如图 3所示。

图 2 导纳函数Fig. 2 Admittance function

图选项

图 3 转换函数Fig. 3 Transform function

图选项

根据挠曲均衡理论，岩石圈作为弹性薄板漂浮在地幔上达到静力平衡状态，那么海底地形作为岩石圈弹性板块上质量载荷将会在莫霍面得到反映，导致莫霍面发生相应变化，因此海面重力异常将受海底地形和地壳挠曲均衡的综合影响。当地形波长远小于挠曲波长时，海面重力异常主要受海底地形影响；然而当地形波长接近甚至大于挠曲波长时，需要考虑地壳均衡对重力异常的贡献。假设由于海底地形b(x,y)的影响导致莫霍面的变化为m(x,y)，那么莫霍面变化与海底地形在频率域的关系^[30]可表示为

(5)

式中，M(f_x, f_y)和B(f_x, f_y)分别表示m(x,y)和b(x,y)的傅里叶变换；ρ_m为地幔密度；Φ_e(f_x, f_y)的表达式如下

(6)

式中，λ为挠曲波长，其大小等于Φ_e(f_x, f_y)=0.5时对应的波长。

当挠曲均衡波长^[31]，从而满足上述情况的导纳函数可表示为

(7)

式中，T_c为地壳厚度。不同挠曲波长λ对应的导纳函数如图 2所示。

基于以上分析，定义转换函数Q(f_x, f_y)具有线性和各向同性的特点，从而可将重力异常作为输入，输出海底地形结果

(8)

对比式(8)和式(4)可知，转换函数和导纳函数互为倒数关系^[32]，不同挠曲波长对应的转换函数示意图如图 3所示。由图 3可知，在短波或者长波波段范围，转换函数数值较大，将会导致该波段范围内重力场的微小噪声被显著放大，进而严重影响海底地形的反演结果；同时在短波或者长波波段范围内，转换函数将不存在Hankel变换对；图 3结果显示不同的挠曲波长将显著影响长波段海底地形反演结果。基于以上分析可知，由于转换函数的固有特点，短波段或者长波段的海底地形解将很不稳定。从而选择使用带通滤波器限制中短波段导纳函数，使得海底地形解算结果保持稳定，同时保证转换函数Q(f_x, f_y)存在Hankel变换对

(9)

式中，W(f_x, f_y)为引入的带通滤波器。通过式(9)获得的转换函数可解算特定波段的海底地形。因而利用重力异常反解某一波段海底地形b_s(x,y)可表示为

(10)

式中，exp(2πfz₀)可视为向下延拓因子。由此可见，频率域内重力异常经滤波并向下延拓后与海深呈现线性关系。由于海底地质条件复杂，一般不使用式(10)的理论参数获取线性系数。实际计算中，通常利用对应波段重力数据和海底地形数据进行线性拟合，利用线性回归分析技术可将式(10)表示为

(11)

式中，a和c分别表示线性回归的比例因子和常数项，经滤波和延拓后的重力异常与海底地形关系如图 4所示(图 4中数据为2.1节试验海区数据)，由图 4可以看出二者存在明显的线性关系(1 Gal=10^-2m/s²)。B_s(f_x, f_y)和ΔG(f_x, f_y)如下

图 4 重力异常和海底地形关系Fig. 4 Gravity anomaly and seafloor topography

图选项

(12)

重力异常垂直梯度与重力异常关系为

(13)

式中，负号表示与z轴指向地球外部，与垂线方向相反。依据傅里叶变换求导法则对式(4)求导得

(14)

因此利用重力异常垂直梯度解算特定波段海底地形可表示为

(15)

式(10)和式(15)显示，重力异常或者重力异常垂直梯度经滤波并向下延拓等处理后与海深呈现线性关系。类似式(10)，由于海底地质条件复杂，一般也不使用式(15)计算理论线性系数。实际操作中，通常利用对应波段重力数据和海底地形数据进行线性拟合，即

(16)

式中，a_z和c_z表示线性回归的比例因子和常数项，经滤波和延拓等处理后的重力异常垂直梯度与海底地形关系如图 5所示(图 5中数据为2.1节试验海区数据)，由图 5可以看出二者存在明显的线性关系。B_s(f_x, f_y)和ΔG_z(f_x, f_y)如下

图 5 重力异常垂直梯度与海底地形关系Fig. 5 Vertical gravity gradient anomaly and seafloor topography

图选项

(17)

综合式(11)和式(16)，将两式相加得

(18)

式中，e为线性拟合常数项。式(18)即为通过多元回归分析技术联合重力异常和重力异常垂直梯度反演海底地形的函数模型。即以经滤波和向下延拓后的重力异常Δg′和经滤波和向下延拓等处理后的重力异常垂直梯度Δg′_z为自变量，以海底地形b为因变量的多元线性回归方程可表示为

(19)

式中，β₀为回归方程的常数项；β₁、β₂为偏回归系数，表示在其他自变量不变的情况下，重力异常或者重力异常垂直梯度每改变一个单位时，单独引起b的平均变化量。至此，经过以上详细的推导流程，建立了联合重力异常和重力异常垂直梯度反演海底地形的可行性方案。实践过程中，基于该理论的海底地形模型构建方法将融合重力异常信号和重力异常垂直梯度信号中蕴含的各自不同频段海底地形特征。

假设利用经滤波和向下延拓过的重力异常或者重力异常垂直梯度反演某一波段的海底地形为b_s(x,y)，则反演波段海底地形与其他波段海底地形(b_r(x,y))求和可得到最终的海底地形(b(x,y))为

(20)

2 数值分析试验2.1 试验区概况和数据准备

印度洋是地球第3大洋，总面积约为75 Mkm²，平均海洋深度3800 m左右。印度洋海底中最为独特的地形地貌特征为大洋底存在巨大的倒Y形洋中脊分布，分别为北部的中印度洋海岭(Central India Ridge，CIR)、东南印度洋海岭(Southeast India Ridge，SEIR)和西南印度洋海岭(Southwest India Ridge，SWIR)。GIR、SEIR和SWIR 3个海岭交汇处称为罗德里格斯三联点(Rodriguez Triple Juction，RTJ)。RTJ北部的CIR向西北方向延伸与红海和东非大裂谷相连，其上分布有中央裂谷。RTJ东南方向的SEIR向东南伸展到阿姆斯特丹岛(Amsterdam Island)，并继续向东延伸，形成澳大利亚-南极海隆(Australia-Antarctic Rise)。位于RTJ西南方向的SWIR向西南方向延伸至Bouvet三联点(Bouvet Triple Juction，BTJ)，并绕过非洲大陆南端与大西洋中脊相连接，全长大约8000 km，总体为斜向扩张洋脊，走向与扩张方向成约60°角，是非洲板块和南极洲板块的边界。SWIR上发育有大量的断裂带(fraction zone，FZ)，文献[33]根据SWIR的扩张历史和几何特征，SWIR自西向东分为7个段落，即沿RTJ-BTJ方向发育有Melville FZ、Atlantis FZ、Gallieni FZ、Indomed FZ、Discovery FZ和Eric Simpson FZ等。如此丰富的地貌分布导致SWIR成为目前海底构造地貌学研究热点，然而南半球船测海深测量数据稀疏且稀少，严重阻碍科学研究与实践。从而本文利用卫星测高重力数据反演海底地形的试验区域选择SWIR部分海域(34°S—41°S，46°E—53°E海域)作为试验海域，融合少量船测海深数据构建海底地形模型；并依据建立的海底地形模型，基于岩浆-构造旋回模型研究分析相关海域的地形地貌及地质构造特点。研究海域数据来源如下：

(1) 卫星测高重力异常数据：来源于丹麦科技大学(Technical University of Denmark)空间实验室(DTU Space)发布的1′×1′ DTU10模型，研究海域重力异常如图 6(a)所示。

图 6 研究海域试验数据Fig. 6 Test data in study area

图选项

(2) 卫星测高重力异常垂直梯度数据：来源于SIO，UCSD(Scripps Institution of Oceanography, University of California, San Diego)，V24.1版本，分辨率为1′。研究海域重力异常垂直梯度如图 6(b)所示。

(3) 船测海深数据：来源于NGDC(The National Geophysical Data Center)发布的研究海域实测数据。研究海域原始船测水深数据共11 892个，通过3σ准则对实测数据进行粗差剔除，最终得到11 392个实测水深数据。选取其中2/3数量的数据共7594个作为控制点，以期对构建的海深模型进行控制；剩余3978个实测水深数据作为外部检核点对最终的海深模型开展精度评估使用。控制点和检核点的分布如图 6(c)所示，其中白色六边形为控制点，黑色三角形为检核点，图中背景为S & S V18.1海深模型。

2.2 海底地形模型构建

根据利用卫星测高重力数据反演海底地形的原理和方法可知，为了达到以卫星测高重力数据作为输入源，输出海底地形解算结果稳定的目的，需要对转换函数进行相应的滤波处理，即使用特定波段重力数据反演相应波段海底地形。其中，重力数据输入波段范围通常选择与海底地形表现较强相关性波段。因此采用多元回归分析技术，利用卫星测高重力数据构建海底地形模型的步骤可分为：

(1) 选取反演波段范围。基于频率域相干性分析技术，分别对重力异常/重力异常垂直梯度与海底地形的相干性进行研究分析，获得重力异常/重力异常垂直梯度与海底地形呈现强相关性的波段范围。

(2) 卫星测高重力数据的滤波和延拓处理。以步骤(1)获得的波段范围作为参考，设计合适的带通滤波器对重力输入数据进行滤波处理，得到相应波段的重力数据；同时将反演波段重力数据经向下延拓等处理，获得经滤波和延拓等处理后的重力数据。

(3) 获取偏回归系数和常数项。以步骤(2)获得的经滤波和延拓后的重力数据与对应波段先验海底地形为基础，依据最小二乘原理对反演波段的重力异常、重力异常垂直梯度和海底地形开展拟合分析，获得多元回归的偏回归系数和常数项。

(4) 建立海底地形模型。以步骤(3)获得的多元回归的偏回归系数和常数项为基础，通过多元回归分析技术获得反演波段海底地形。再综合各个波段的海底地形结果，最终建立联合多元重力数据(重力异常和重力异常垂直梯度)构建的海底地形模型。

由于研究海域船舶测量点分布过于稀疏，单纯依靠控制点难以得到研究海域格网化海深，而反演实施过程中需利用先验海深信息补偿非反演波段海底地形信息。因此本文以S & S V18.1为先验海深模型填补研究海域数据空白区，融合海深测量控制点得到格网化先验海深模型。依据文献[34]中的相干性分析技术，分别对研究海域重力异常和重力异常垂直梯度与海底地形进行频率域相干性处理以获得与海底地形呈现强相关性波段范围，相干性分析结果如图 7所示，其中黑色点线和蓝色点虚线分别表示重力异常、重力异常垂直梯度与海底地形相干性结果。

图 7 重力数据与海底地形的相干性Fig. 7 Coherence between gravity data and seafloor topography

图选项

若认为重力数据与海底地形相干性在0.5以上为强相关，那么从重力异常、重力异常垂直梯度与海底地形相干性结果可以看出：重力异常垂直梯度与海底地形的强相干性波段范围为18~ 160 km，重力异常与海底地形在20~200 km表现出强相干性，结果表明重力异常垂直梯度对高频地形信息更加敏感。结合图 3转换函数在长波段和短波段数值分布特点，最终依据卫星测高重力数据反演海底地形的波段范围选择为45~160 km。

因利用卫星测高重力数据反演海底地形主要依据傅里叶变换在频率域上完成，为了消除边缘效应的影响，傅里叶逆变换结果的边界分别内缩10′。然后基于最小二乘原则，以反演波段重力异常/重力异常垂直梯度和先验海底地形数据为输入，以半径为10′的移动窗口对该窗口范围内的重力和海底地形数据进行单元回归和多元回归分析，获得窗口中央点的偏回归系数和常数项。再以研究点重力异常/重力异常垂直梯度数据作为输入，结合该点的偏回归系数和常数项获得反演波段海底地形。分别以重力异常和重力异常垂直梯度作为单独输入数据，经过上述解算流程得到的偏回归系数和常数项分别如图 8(a)—(d)所示。

图 8 重力异常和重力异常垂直梯度的回归系数和常数项Fig. 8 Regression coefficient and constant term of vertical gravity gradient anomaly and gravity anomaly

图选项

对比图 8(a)和图 8(c)重力数据与海底地形回归系数的分布情况，二者具有空间分布的一致性特点：在南北方向呈现线型分布。图 8(b)和图 8(d)显示重力异常和重力异常垂直梯度的常数项也具有空间分布一致性特征。联合重力异常和重力异常垂直梯度数据，基于多元回归分析技术获得的重力异常和重力异常垂直梯度的偏回归系数和常数项如图 9(a)—(c)所示。

图 9 多元回归偏回归系数和常数项Fig. 9 Multiple regression partial regression coefficient and constant term

图选项

对比图 9(a)和图 9(b)可知，多元回归分析中，重力异常和海底地形、重力异常垂直梯度和海底地形的偏回归系数在南北方向依然呈现出近线型分布，且二者的空间分布特征保持一致。对比基于多元回归分析获得的重力异常和重力异常垂直梯度与海底地形的偏回归系数(图 9(a)、(b))和基于单元回归技术获得的仅仅以重力异常/重力异常垂直梯度数据作为唯一输入源的回归系数(图 8(a)、(c))，重力异常和重力异常垂直梯度的回归系数数值大小差别较大。同时回归分析的常数项(图 8(b)、(d)和图 9(c))显示：仅仅以重力异常和重力异常垂直梯度作为数据输入的常数项二者差别较小；而以重力异常和重力异常垂直梯度数据作为共同输入源获得的多元回归分析的常数项与单元(重力异常/重力异常垂直梯度)输入常数项虽然空间分布特征一致，但是其数值大小差异明显。基于以上分析可以发现，联合多元数据(重力异常和重力异常垂直梯度)并非简单地将单元数据相加，进而验证了联合多元数据反演海底地形的有效性。最终基于多元回归分析技术，联合多元重力数据构建的海深模型(图 10)。为描述方便，将重力异常和重力异常垂直梯度构建的海深模型分别称为BDG模型和BVGG模型，联合多元重力数据反演得到的海深模型称为BDVG模型。

图 10 BDVG模型Fig. 10 BDVG model

图选项

2.3 试验结果精度分析

为进一步说明基于多元回归分析技术构建的海底地形模型与重力异常和重力异常垂直梯度构建的海底地形模型结果取算术平均的本质差异，将重力异常和重力异常垂直梯度构建的海底地形模型相加并取算术平均获得的海底地形模型称为BMDVG模型。比对评估BDG模型、BVGG模型、BDVG模型和BMDVG模型精度，以Smith and Sandwell发布的最新海深模型V18.1 (简称S & S V18.1模型)和国际通用的ETOPO1模型作为比对参考，通过模型内插方法将模型海深值内插到外部检核点处与检核点测量海深进行差异比较，进而对各个模型展开精度的检验评估，各个模型的检核统计结果见表 1。

表 1 海深模型检核统计结果Tab. 1 Bathymetry model statistical results

m
海深模型	|最大值|	|最小值|	平均值	均方差	相关系数
BDG模型	785.28	0.03	-1.83	100.65	0.993 1
BVGG模型	910.31	0.02	7.11	111.46	0.992 0
BDVG模型	959.02	0.04	-3.30	97.76	0.993 4
BMDVG模型	847.80	0.01	2.64	103.93	0.992 8
S & S V18.1	1 363.74	0.02	4.37	110.48	0.991 5
ETOPO1	1 249.39	0.13	-25.38	231.72	0.962 4

表选项

表 1中各模型检核统计结果显示：与检核点差值最大的海深模型是S& S V18.1海深模型，差值超过1300 m；ETOPO1海深模型值与检核点的最大互差结果也较大，也超过了1000 m；BDG海深模型检核差值为780 m左右。究其原因，笔者认为可能是因为差值结果最大值仅仅是单点比较结果，若测量点自身存在粗差，而在数据预处理粗差剔除阶段未被发现，势必会对海深模型的检核统计结果产生影响。如分析检核数据结果发现，各个模型差值较大值所共有的检核点共3个，剔除这3个认为存在粗差的检核点后的BDG模型、BVGG模型、BDVG模型和BMDVG模型检核差值绝对值最大值分别为509.48、708、571.53和581.36 m。对比分析6种海深模型的检核均方差表明：研究海域范围内，检核精度最优的是利用多元回归分析，联合重力异常和重力异常垂直梯度数据构建的BDVG模型，其余依次为BDG模型、BMDVG模型、BVGG模型、S & S V18.1模型、ETOPO1模型。其中本文构建的BDVG海深模型、BDG海深模型、BVGG海深模型和BMDVG模型检核均方差结果优于S & S V18.1和ETOPO1模型。相较于ETOPO1模型，本文构建的4种海深模型检核精度提高了近一倍有余。S & S V18.1模型检核精度与BVGG模型相当，低于BDG模型、BDVG模型和BMDVG模型。比对分析本文构建的4种海深模型(BDG模型、BVGG模型、BDVG模型和BMDVG模型)，联合多元重力数据构建的BDVG模型相较于BDG模型和BVGG模型检核精度分别提高了约2.87%和12.29%，BDVG模型相较BDG模型精度提升有限(不超过3%)，而较BVGG提升效果明显。笔者认为产生该结果的原因：因试验开始阶段考虑重力异常和重力异常垂直梯度与海底地形频谱相干性及转换函数等因素而选择的反演波段范围所致。理论上讲，重力异常垂直梯度主要反映高频地形信息，即重力异常垂直梯度对高频地形信息相比重力异常更加敏感，而本文反演波段范围在45~160 km(高频部分重力信息未充分利用)，该波段重力异常与地形符合较好，从而导致最终联合多元重力数据构建的BDVG模型偏向于重力异常为输入源构建的海底地形模型，重力异常垂直梯度在模型构建中对结果进行补充和修正，这也从侧面体现了采用多元回归分析技术，联合多元重力数据构建海底地形模型的优势和意义。

不同重力数据反演结果取算术平均构建的BMDVG模型精度低于采用多元回归分析构建的BDVG模型，说明了虽然多元回归模型推导过程中将不同重力数据构建的海底地形模型相加并取平均，但是实际操作并非简单地将反演模型权重各半，多元回归分析技术可通过自动调节偏回归系数以获得联合多元重力数据构建海底地形模型的最佳系数组合。6种海深模型插值与检核点的相关系数结果与检核精度结果一致，6种海深模型插值与检核点均表现出很强的相关性，ETOPO1海深模型插值与检核点相关系数略低于另外5种海深模型。

进一步分析比对本文构建的4种海深模型，统计BDG模型、BVGG模型、BMDVG模型和BDVG模型不同差值范围内检核点个数，结果如图 11所示。

图 11 模型检核差值统计Fig. 11 Each model checking statistics

图选项

图 11中4种海深模型不同检核差值范围检核点数量统计结果显示：随着检核差值增大，检核点数量不断减少，检核差值小于50 m范围内的检核点数量较多。比对4种海深模型，检核差值小于50 m左右范围内，采用多元回归分析构建的BDVG模型检核点数量明显多于其他3种海深模型；检核差值大于80 m左右，检核点数量小于其他3种模型，反映了联合多元重力数据，采用多元回归技术构建的海底地形模型优越性。比对BMDVG模型和BDVG模型检核差值曲线可以发现，BDVG模型检核差值小的检核点数量多于BMDVG模型，进一步说明了本文提出的多元回归分析技术并非简单将重力异常和重力异常垂直梯度反演海底地形模型权重各半。分析评价各个海深模型的模型效能，定义海深模型内插到检核点并与检核点海深的差值与检核点海深之比为相对误差，各个模型的相对误差统计结果见表 2。

表 2 各模型相对误差统计结果Tab. 2 Relative error of each model statistical results

(%)
海深模型	|最大值|	|最小值|	平均值	标准差
BDG模型	90.85	6.35e-04	-0.58	6.43
BVGG模型	180.97	6.20e-04	-1.51	11.23
BDVG模型	49.91	1.28e-03	-0.08	4.69
BMDVG模型	135.91	2.15e-04	-1.04	8.67
S & S V18.1	157.26	6.12e-04	-0.40	6.38
ETOPO1	159.90	6.83e-03	0.17	10.97

表选项

表 2中各模型的相对误差统计结果显示：6种海深模型的相对误差标准差最小为基于多元回归分析技术构建的BDVG模型，其标准差为4.69%；BVGG模型相对误差标准差与ETOPO1模型相当。各模型相对误差的标准差统计结果与检核统计精度结果大致一致。其中BVGG模型、ETOPO1模型的相对误差的绝对值最大分别为180%和159%左右，BDVG模型相对误差的绝对值最大值在5种海深模型中相对较小，但也为49.91%。造成该结果的原因笔者认为与检核差值统计结果类似，可能因为相对误差绝对值最大值统计结果仅仅是反映检核点中一个检核点的统计结果，若因为原始测量或者其他因素造成该测量点本身存在测量粗差，又由于粗差剔除不彻底未发现该点导致反映在最终统计结果中，这样将会对结果产生误导。

另外，BVGG模型相对误差平均值结果达到-1.51%，明显比其他模型结果大，说明该模型存在明显的系统性误差。笔者认为原因在于：在本文反演波段范围内重力异常相比重力异常垂直梯度而言与海底地形信息符合较好。重力异常垂直梯度主要反映高频地形信息，从而导致以重力异常垂直梯度为输入数据构建的BVGG模型存在较为明显的系统性误差。而本文提出采用多元回归分析技术，联合重力异常和重力异常垂直梯度数据构建海底地形模型中，重力异常垂直梯度信息主要对模型进行补充和修正。如将重力异常和重力异常垂直梯度反演结果权重各半构建的海底地形模型(BMDVG模型)，检核精度反而低于BDG模型，相对误差平均值也达到1.04%，说明了该方法重力异常垂直梯度信息对模型的补偿和修正不尽合理，从而进一步表明本文提出的采用多元回归分析技术与采用将重力异常和重力异常垂直梯度反演结果权重各半方法构建的海底地形存在本质上的差异。为更直观地观察不同海深模型的相对误差状况，将各模型相对误差的空间分布情况绘图如图 12所示。

图 12 相对误差空间分布Fig. 12 Relative error distribution

图选项

图 12相对误差的空间分布图中，各个模型相对误差的空间分布具有类似的特征。如在(-37°S, 52°E)和°S, 47°E)附近区域相对误差值明显大于其他区域，对比图 10可知该区域有海沟存在，地形变化较为剧烈。通过多元回归分析建立的BDVG模型在这两个区域的表现较好，相对误差值较小，表明基于多元回归技术建立的BDVG模型稳定性能较单独使用重力异常或者重力异常垂直梯度数据构建的海底地形模型较强。进一步评估分析5种海深模型在不同海深条件下的模型效能，比对不同水深环境下5种海深模型的相对误差如图 13所示。

图 13 相对误差随海深分布Fig. 13 Relative error change with sea depth

图选项

从图 13不同海深条件下相对误差分布可以看出：ETOPO1模型相较于其他5种模型相对误差波动较大；海深低于1500 m的情况下，5种海深模型的相对误差均开始发生较大波动。基于多元回归技术构建的BDVG模型在该水深范围内，相比于BDG模型、BVGG模型和BMDVG模型相对误差波动较小，且均匀分布在零值附近，而BDG模型、BVGG模型和BMDVG模型均存在不同程度过低估计海深的情况，进一步验证了利用多元回归分析技术，联合重力异常和重力异常垂直梯度数据建立海底地形模型的优越性和科学性。BDVG模型在海深小于1500 m范围内，相对误差的波动幅度也较S & S V18.1低。基于以上分析可知，利用多元回归技术，联合多元卫星测高重力数据构建的海底地形模型具有比仅仅以一种重力数据作为输入建立的海底地形具有更高的精度和更好的稳定性，进一步证明了联合多元重力数据构建海底地形模型的重要性和可行性。

2.4 地形地貌分析

以基于多元回归分析，联合重力异常和重力异常垂直梯度数据构建的BDVG海深模型(图 14)为基础，分析研究试验海域地形地貌特征。图 14清晰显示了研究海域SWIR洋脊E-W斜向展布，同时能清晰观察到试验海区西部和东部(图 14中Ⅰ和Ⅱ所示)的Indomed断裂带和Gallieni断裂带。进一步分析SWIR洋脊可以发现，洋脊段发育着大量近E-W断裂，然而在°S，50.5°E)附近的洋脊段相比于毗邻的洋脊段出现裂谷缺失现象，而南北两侧对称发育着巨型隆起；同时文献[35]也发现该段洋脊在50.5°E附近存在明显的低磁异常带，由此推测该洋脊段存在着激烈的岩浆活动，该洋脊段通常称为27洋脊段，由东到西依次为28洋脊段和29洋脊段^[36]。

图 14 BDVG模型Fig. 14 BDVG model

图选项

进一步分析27洋脊段、28洋脊段和29洋脊段的地形走势，沿SWIR洋脊方向(图 14实线1)作3个洋脊段地形剖面，如图 15所示。图 15中从东到西分别为27洋脊段、28洋脊段和29洋脊段，大概对应图中黑色、绿色和蓝色方框位置。从图 15可以看出，27洋脊段沿轴方向大约40 km没有裂谷，然后沿东西方向海深均逐渐变深而重新出现裂谷。28洋脊段两侧水深依次增加，东侧水深又迅速减小，地形抬升超千米过渡到27洋脊段；西侧水深增加3500 m左右过渡到29洋脊段范畴。29洋脊段东西两侧水深均不断增加，近似对称分布。28洋脊段和29洋脊段地形抬升部分为发育的轴部火山脊(axial volcanic ridges，AVRs)剖面。从图 15可以看出，28洋脊段和29洋脊段在E-W方向呈脊槽现象，然而在27洋脊段脊槽分布现象消失，往东又再次出现。在29洋脊段沿南北方向做地形剖面(图 14中2号线位置)，结果如图 16所示。由图 16可知，29洋脊段南北两侧地形不对称，呈现南高北低的地貌形态，而由于反演结果分辨率的限制，对于南北谷壁细致的地质构造呈现不清晰。

图 15 位置1剖面Fig. 15 Profile at position 1

图选项

图 16 29洋脊段剖面Fig. 16 29 ocean ridge profile

图选项

沿图 14中3号线做27洋脊段的地形剖面，结果如图 17所示。结合图 14和图 17可以看出，27洋脊段沿轴两侧发育有呈对称分布的两组共轭隆起，且向轴侧坡度陡峭，离轴侧坡度较缓，同时轴部火山活动明显。文献[37]研究了SWIR轴部的岩浆-构造旋回，提出了轴部演化主要由岩浆活跃期和构造活跃期交替进行，文献[36]也将其称之为岩浆-构造动力学模式。岩浆活跃期，岩浆供应量充足导致建造成脊；构造活跃期，岩浆供应量较少，海底地形地貌主要由构造扩张所消化；所以岩浆-构造动力学模式本质上对应岩浆由多到少的供应周期。因此由岩浆-构造动力学模式，结合27洋脊段火山活动明显特点，目前27洋脊段属于岩浆活跃期阶段，轴部两侧的对称隆起推测是随海底扩张迁移所致；同时27洋脊段轴部的对称裂离也从侧面验证了目前该洋脊段岩浆供应充足，正是由于27洋脊段岩浆供应充足，从而削弱了构造活动在海底扩张中的影响。对比29洋脊段剖面图(图 16)，29洋脊段表现出了不同的裂离方式——非对称裂离。究其原因，笔者认为可能是27洋脊段和29洋脊段分属于不同的岩浆控制系统，图 15中27洋脊段脊槽分布模式与另外洋脊段存在差异也可作为侧面验证，因此可能存在其他叠加因素对27洋脊段产生相应影响。

图 17 27洋脊段剖面图Fig. 17 27 ocean ridge profile

图选项

3 总结

本文提出了基于多元回归分析技术联合多元重力数据的海底地形构建方法，选取南半球西南印度洋SWIR所在的部分海域开展研究试验。首先依据试验海域重力和海底地形的相干性分析结果，结合转换函数在不同挠曲波长不同频段的函数变化差异特点，获得了依据试验海域卫星测高重力数据反演海底地形的有效波段范围。然后采用多元回归分析技术，联合卫星测高重力异常和重力异常垂直梯度数据构建了相应的海底地形模型(BDVG模型)；并将BDVG模型与仅采用重力异常和重力异常垂直梯度作为输入源建立的海底地形模型(BDG模型和BVGG模型)和将重力异常和重力异常垂直梯度反演海底地形结果分别权重各半建立的BMDVG模型以及S & S V18.1海深模型、ETOPO1海深模型进行了比对分析。最后，基于BDVG海深模型简要分析了SWIR上27洋脊段、28洋脊段和29洋脊段的地形地貌以及地质构造特征。

试验结果表明，基于多元回归分析技术构建的海底地形模型检核精度最高，相较于S & S V18.1模型和ETOPO1海深模型精度分别提高了11.51%和57.81%左右；水深2000 m以上海域，海深模型的检核精度较高，相对误差波动较小；水深小于1500 m的海域，海深模型的检核精度下降较快，相对误差波动较剧烈，反映了深海海域良好的反演效果。地形起伏剧烈海域和浅海海域，基于多元回归分析技术构建的BDVG海底地形模型相较于以重力异常和重力异常垂直梯度作为单一输入源建立的BDG模型和BVGG模型相对误差及相对误差波动变化较小，反映了BDVG海底地形模型拥有更好的稳定性能，从而体现了联合反演的必要性和优势；另外基于多元回归分析构建的海底地形模型并非简单地将反演模型权重各半，其可通过自动调节偏回归系数以获得联合多元重力数据构建海底地形模型的最佳系数组合，而不需外部检核点作为约束。27洋脊段是Indomed FZ—Gallieni FZ上唯一轴部缺失裂谷洋脊段；27洋脊段目前属于岩浆供应充足阶段，构造作用的海底扩张对其影响较小；同时27洋脊段的裂离方式为对称裂离，沿轴分布有南北对称的地形隆起。

基于卫星测高重力数据构建的海底地形模型目前只能反映大尺度的地形地貌特征，经前文对研究海域洋脊段的地形地貌分析可知，依据构建的海底地形模型可在一定程度上对海底地质构造演化展开推测演绎，进而熟悉了解目标海域的地质构造。最终若需要细致了解该地区的地质构造特征，可结合其他地质和水文等资料对目标海区开展更加深入的研究。

【引文格式】范雕, 李姗姗, 杨军军, 等. 利用多元回归分析反演西南印度洋区域海底地形. 测绘学报，2020，49(2)：147-161. DOI: 10.11947

悼念丨永远怀念宁津生院士

论文推荐 | 宋迎春：加权混合估计中权值的确定方法

资讯 | 在家学习效率低？收下这 9 个实用工具，让你自学成学霸

书讯 | 王家耀院士：地图制图数据处理的模型与方法

论文推荐 | 赵泉华：全卷积网络和条件随机场相结合的全极化SAR土地覆盖分类