1.本发明涉及海底地形反演技术领域，尤其涉及一种恢复海陆交界区海底地形的改进重力地质方法。
背景技术：2.目前，重力反演海底地形业已成为全球海底地形模型构建的主要手段，国内外发布的全球海底地形数据资料，其中绝大部分海域海底地形数据均依赖重力数据恢复得到，如sto_ieu2020、bat_whu2020、sandwell系列模型、geoco系列模型等。应用重力反演海底地形方法研究主要集中在全海洋覆盖空间的规则矩形海域海底地形恢复，鲜有涉及海陆交界不规则海域海底地形恢复研究。原因可能是，1)导纳函数方法和sa(simulated annealing)方法一般在频率域处理重力数据和地形数据，涉及傅里叶正反变换操作，从而通常应用于规则矩形海域海底地形反演，对待海陆交界区不规则海域海底地形恢复效果不佳；2)海底地形反演方法实施过程中一般涉及数据格网化过程，而海陆交界区海洋区域数据格网化结果势必受陆地区域数据干扰，导致采用常规海底地形反演方法恢复陆海交界区海洋部分地形数据鲁棒性较弱。然而，一方面海陆交界区海底地形作为全球海底地形重要组成部分，是自主建立全球海底地形模型不可回避的重要内容；另一方面随着邻海区域多源重力数据不断累积丰富和数据处理技术不断进步，海陆交界区海洋部分相关可用数据质量将得到极大提升，可预见临海区域海洋地形相关数据的开发与利用势在必行。
技术实现要素：3.本发明针对重力反演海陆交界不规则形状海域海底地形算法实施困难和鲁棒性较弱等缺点，基于海底地形数据与重力数据间近似线性关系，提出一种恢复海陆交界区海底地形的改进重力地质方法，对于发挥重力数据在海陆交界区海底地形资料获取，以及全球海底地形模型精准构建方面优势大有裨益。
4.为了实现上述目的，本发明采用以下技术方案：
5.一种恢复海陆交界区海底地形的改进重力地质方法，包括：
6.步骤1：将声学水深测量数据按比例分为两类，一类为水深离散控制点dk，用于参与海底地形模型构建；另一类为水深离散检核点dj，用于评估海底地形模型质量以获取最佳密度参数；
7.步骤2：设置密度差异常数初值
△
ρ0，基于dk、
△
ρ0及参考水深求解离散控制点处短波重力异常
8.步骤3：基于海面重力异常数据及计算离散控制点处长波重力异常
9.步骤4：将格网化获得目标区域格网化长波重力异常
10.步骤5：基于获得海陆交界区海洋部分网格长波重力异常
11.步骤6：基于海陆交界区海面重力异常数据及获取海陆交界区海洋部分网格
短波重力异常
12.步骤7：利用反演得到海陆交界区海洋部分网格化海底地形结果d
gwo
，即海陆交界区海底地形模型；
13.步骤8：基于dj、采用均方根评估得到的海陆交界区海底地形模型精度；
14.步骤9：在设置的密度差异常数初值基础上调整密度差异，重复步骤2至步骤9，直至得到的海陆交界区海洋部分网格化海底地形模型精度最高为止；
15.步骤10：将海底地形模型检核均方根最小对应的密度差参量作为最佳密度差常数，基于步骤7相同方式得到对应的海陆交界区海洋部分网格化海底地形结果；
16.步骤11：使用最邻近内插方法将步骤10得到的海底地形结果中数值大于零的数替换，得到最终海陆交界区海洋部分网格化海底地形结果。
17.进一步地，所述步骤2包括：
18.利用式(1)获得离散控制点处短波重力异常：
[0019][0020]
式中，为控制点短波重力异常；g为万有引力常数；δρ为研究点地壳密度和海水密度之差；dk是控制点处水深；d是参考水深。
[0021]
进一步地，所述步骤3包括：
[0022]
步骤3.1：将海面重力异常
△
g内插到离散控制点处获得控制点处重力异常
△gk
，其插值公式为
[0023][0024]
式中
[0025][0026]
式中，(x0,y0)表示离散控制点位置；x1、x2、y1和y2表示重力异常模型值点位置；
[0027]
步骤3.2：利用和
△gk
，按照下式计算离散控制点处长波重力异常
[0028][0029]
其中，为控制点长波重力异常。
[0030]
进一步地，所述步骤6包括：
[0031]
步骤6.1：使用步骤3.1中描述的内插方法将海陆交界区海面重力异常数据内插得到海陆交界区海洋部分网格重力异常
△ggwo
；
[0032]
步骤6.2：将
△ggwo
与作差得到海陆交界区海洋部分网格短波重力异常
[0033]
进一步地，所述步骤7包括：
[0034]
按照下式反演得到海陆交界区海洋部分网格化海底地形结果d
gwo
[0035][0036]
其中，g为万有引力常数；d是参考水深。
[0037]
进一步地，所述步骤8中，按照下式评估得到的海陆交界区海底地形模型精度：
[0038][0039]
式中，d
′i表示第i个检核点处海底地形模型插值；表示第i个检核点海深值。
[0040]
与现有技术相比，本发明具有的有益效果：
[0041]
本发明提供了解决海陆交界区不规则形状海域海底地形建模策略。利用短波重力异常与海深间近似线性关系，以外部检核差值统计结果为参考，通过不断调整整个区域海水与地壳的密度差，致使海深与短波重力异常达到平衡，进而获得整个不规则形状海域最佳密度差常数；然后基于短波重力异常与海深间函数关系，反演得到海陆交界区不规则形状海域海底地形。
[0042]
本发明具有较强的鲁棒性，且最终反演得到的海陆交界区不规则形状海域海底地形相较于scs_dtu18和scs_etopo1精度分别提高了78％和71％左右。
附图说明
[0043]
图1为本发明实施例一种恢复海陆交界区海底地形的改进重力地质方法的基本流程图；
[0044]
图2为本发明实施例船测水深点分布图；
[0045]
图3为本发明实施例离散控制点处短波重力异常和长波重力异常分布图；
[0046]
图4为本发明实施例海陆交界区海洋部分网格长波重力异常分布图；
[0047]
图5为本发明实施例海陆交界区海洋部分重力异常和短波重力异常分布图；
[0048]
图6为本发明实施例得到的最终海陆交界区海洋部分网格化海底地形结果scs_iggm；
[0049]
图7为dtu18海底地形模型和etopo1海底地形模型示例图。
具体实施方式
[0050]
下面结合附图和具体的实施例对本发明做进一步的解释说明：
[0051]
如图1所示，一种恢复海陆交界区海底地形的改进重力地质方法，包括：
[0052]
步骤一：将声学水深测量数据按照一定比例分为两类；一类为水深离散控制点dk，用于参与海底地形模型构建；另一类为水深离散检核点dj，用于评估海底地形模型质量以获取最佳密度参数。
[0053]
步骤二：设置密度差异常数初值
△
ρ0，求解离散控制点处短波重力异常
[0054]
利用式(1)获得离散控制点处短波重力异常：
[0055][0056]
式中，为控制点短波重力异常；g为万有引力常数；δρ为研究点地壳密度和海水密度之差；dk是控制点处水深；d是参考水深，通常选择研究区域最大水深作为参考水深。
[0057]
步骤三：计算离散控制点处长波重力异常
[0058]
所述步骤三包括：
[0059]
步骤3.1：将海面重力异常
△
g内插到离散控制点处获得控制点处重力异常
△gk
，其插值公式为
[0060][0061]
式中
[0062][0063]
式中，(x0,y0)表示离散控制点位置；x1、x2、y1和y2表示重力异常模型值点位置。
[0064]
步骤3.2：利用步骤二得到控制点短波重力异常
△gsk
和步骤3.1获取的控制点处重力异常
△gk
，计算控制点处长波重力异常；公式如下
[0065][0066]
其中，
△gk
为海面控制点重力异常；为控制点长波重力异常；为控制点短波重力异常。
[0067]
步骤四：将步骤三控制点处长波重力异常格网化获得目标区域格网化长波重力异常具体地，使用的格网化方法如gmt(generic mapping tools)软件“surface”命令将步骤三控制点处长波重力异常格网化为目标区域格网化长波重力异常
[0068]
步骤五：使用步骤3.1中描述的内插方法，将目标区域格网化长波重力异常内插获得海陆交界区海洋部分网格长波重力异常
[0069]
步骤六：获取海陆交界区海洋部分网格短波重力异常
[0070]
所述步骤六包括：
[0071]
步骤6.1：使用步骤3.1中描述的内插方法将目标海域海面重力异常
△
g内插得到海陆交界区海洋部分网格重力异常
△ggwo
。
[0072]
步骤6.2：将步骤6.1获得的
△ggwo
与步骤五获得的作差得到海陆交界区海洋部分网格短波重力异常
[0073]
步骤七：利用步骤六获得的反演得到海陆交界区海洋部分网格化海底地形d
gwo
；计算式如下
[0074][0075]
式中，参数意义同上。
[0076]
步骤八：评估步骤七构建的海陆交界区海底地形模型精度，公式如下
[0077][0078]
式中，d
′i表示第i个检核点处海底地形模型插值；表示第i个检核点海深值。
[0079]
步骤九：在设置的密度差异常数初值基础上以一定步长调整密度差异，重复步骤二至步骤九，直至构建的海陆交界区海洋部分网格化海底地形模型精度最高为止。
[0080]
步骤十：将海底地形模型检核均方根最小对应的密度差参量作为最佳密度差常数，利用式(5)反演得到海陆交界区海洋部分网格化海底地形结果。
[0081]
步骤十一：使用最邻近内插/推估方法将步骤十得到的海底地形结果中数值大于零的数替换，得到最终海陆交界区海洋部分网格化海底地形结果。具体地，最邻近内插/推估方法是利用已知点坐标和待求点坐标计算出已知点和待求点间空间距离，选取5个距离待求点最近的已知点海深作为输入数据，通过距离加权方法内插/推估得到待求点海深值。
[0082]
作为一种具体可实施方式，选取中国南海陆海交界区域2
°×2°
(6
°
n～8
°
n,115
°
e～117
°
e)范围为研究对象，使用本发明方法，融合船载水深测量数据和grav_31.1重力异常数据反演恢复海陆交界区海底地形数据。
[0083]
步骤一，将声学水深测量数据按照一定比例分为两类。一类为水深离散控制点，另一类为水深离散检核点。实施例中，检核点和控制点分别为2467个和9871个，控制点和检核点的分布如图2所示，其中灰色圆点为控制点，黑色三角形为检核点。
[0084]
步骤二，利用步骤一中控制点数据得到试验实施海域参考水深为2971m。设置密度差异常数初值
△
ρ0为1.00g/cm3，获得离散控制点处短波重力异常和长波重力异常分别如图3(a)和图3(b)所示。
[0085]
步骤三，将步骤二控制点长波重力异常格网化得到目标区域格网化长波重力异常然后将目标区域格网化长波重力异常内插获得海陆交界区海洋部分网格长波重力异常海陆交界区海洋部分网格长波重力异常如图4所示。
[0086]
步骤四，将目标海域grav_31.1重力异常数据(称为scs_grav_31.1)与步骤三作差，获取海陆交界区海洋部分网格短波重力异常命名为scs_short_grav。scs_grav_31.1和scs_short_grav分别如图5(a)和图5(b)所示。
[0087]
步骤五，利用步骤四获得的scs_short_grav反演得到海陆交界区海洋部分网格化海底地形，相应海底地形模型检核rms值为827.43m。
[0088]
步骤六，以0.01g/cm3为步长调整密度差异重复步骤二至步骤五，直至构建的海陆交界区海洋部分网格化海底地形模型精度最高(对应的检核rms最小)为止。最后，最佳密度差常数(海底地形模型检核rms最小对应的密度差参量)结果为2.01g/cm3。
[0089]
步骤七，利用步骤六得到的最佳密度差常数，反演得到实施例海域海底地形模型。
[0090]
步骤八，使用最邻近内插/推估方法将步骤七得到的海底地形结果中数值大于零的数替换，得到最终海陆交界区海洋部分网格化海底地形结果，命名为scs_iggm。scs_iggm如图6所示。
[0091]
引入目前广泛使用的dtu18海底地形模型和etopo1海底地形模型参与模型横向比对，dtu18海底地形模型和etopo1海底地形模型如图7所示。直观感受而言，三款海底地形模型(海深模型)均能较好吻合，其中scs_dtu18模型西北方向相较于其他模型存在4处明显的地形隆起。
[0092]
scs_iggm、与scs_dtu18、scs_etopo1间差值统计结果见表1。
[0093]
表1不同海底地形模型间差值统计(单位：m)
[0094][0095]
依据表1可以看出，scs_iggm与scs_dtu18、scs_etopo1相关系数均超过了0.99，表明本发明方法恢复的海陆交界部分海底地形数据与其整体吻合较好。进一步比较分析scs_iggm、scs_dtu18、scs_etopo1模型效能，以未参与建模的3086个外部水深点评估海底地形模型外符合精度，将上述海深模型内插到实测水深点，内插结果与实测水深点处船测海深值作差，相应的差值统计结果见表2。
[0096]
表2海底地形模型外符合精度统计结果(单位：m)
[0097][0098]
从表2可以看出，scs_iggm、scs_dtu18、scs_etopo1与外部实测水深值的rms分别为40.35m、183.79m和138.15m。可见，scs_iggm相较于scs_dtu18和scs_etopo1精度分别提高了78％和71％左右，进一步验证了本发明方法的有效性。
[0099]
以上所示仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
技术特征：1.一种恢复海陆交界区海底地形的改进重力地质方法，其特征在于，包括：步骤1：将声学水深测量数据按比例分为两类，一类为水深离散控制点d
k
，用于参与海底地形模型构建；另一类为水深离散检核点d
j
，用于评估海底地形模型质量以获取最佳密度参数；步骤2：设置密度差异常数初值
△
ρ0，基于d
k
、
△
ρ0及参考水深求解离散控制点处短波重力异常步骤3：基于海面重力异常数据及计算离散控制点处长波重力异常步骤4：将格网化获得目标区域格网化长波重力异常步骤5：基于获得海陆交界区海洋部分网格长波重力异常步骤6：基于海陆交界区海面重力异常数据及获取海陆交界区海洋部分网格短波重力异常步骤7：利用反演得到海陆交界区海洋部分网格化海底地形结果d
gwo
，即海陆交界区海底地形模型；步骤8：基于d
j
、采用均方根评估得到的海陆交界区海底地形模型精度；步骤9：在设置的密度差异常数初值基础上调整密度差异，重复步骤2至步骤9，直至得到的海陆交界区海洋部分网格化海底地形模型精度最高为止；步骤10：将海底地形模型检核均方根最小对应的密度差参量作为最佳密度差常数，基于步骤7相同方式得到对应的海陆交界区海洋部分网格化海底地形结果；步骤11：使用最邻近内插方法将步骤10得到的海底地形结果中数值大于零的数替换，得到最终海陆交界区海洋部分网格化海底地形结果。2.根据权利要求1所述的一种恢复海陆交界区海底地形的改进重力地质方法，其特征在于，所述步骤2包括：利用式(1)获得离散控制点处短波重力异常：式中，为控制点短波重力异常；g为万有引力常数；δρ为研究点地壳密度和海水密度之差；d
k
是控制点处水深；d是参考水深。3.根据权利要求1所述的一种恢复海陆交界区海底地形的改进重力地质方法，其特征在于，所述步骤3包括：步骤3.1：将海面重力异常
△
g内插到离散控制点处获得控制点处重力异常
△
g
k
，其插值公式为式中
式中，(x0,y0)表示离散控制点位置；x1、x2、y1和y2表示重力异常模型值点位置；步骤3.2：利用和
△
g
k
，按照下式计算离散控制点处长波重力异常其中，为控制点长波重力异常。4.根据权利要求3所述的一种恢复海陆交界区海底地形的改进重力地质方法，其特征在于，所述步骤6包括：步骤6.1：使用步骤3.1中描述的内插方法将海陆交界区海面重力异常数据内插得到海陆交界区海洋部分网格重力异常
△
g
gwo
；步骤6.2：将
△
g
gwo
与作差得到海陆交界区海洋部分网格短波重力异常5.根据权利要求1所述的一种恢复海陆交界区海底地形的改进重力地质方法，其特征在于，所述步骤7包括：按照下式反演得到海陆交界区海洋部分网格化海底地形结果d
gwo
其中，g为万有引力常数；d是参考水深。6.根据权利要求1所述的一种恢复海陆交界区海底地形的改进重力地质方法，其特征在于，所述步骤8中，按照下式评估得到的海陆交界区海底地形模型精度：式中，d
i
'表示第i个检核点处海底地形模型插值；表示第i个检核点海深值。
技术总结本发明涉及海底地形反演技术领域，公开一种恢复海陆交界区海底地形的改进重力地质方法，利用短波重力异常与海深间近似线性关系，以外部检核差值统计结果为参考，通过不断调整整个区域海水与地壳的密度差，致使海深与短波重力异常达到平衡，进而获得整个不规则形状海域最佳密度差常数；然后基于短波重力异常与海深间函数关系，反演得到海陆交界区不规则形状海域海底地形。本发明对于发挥重力数据在海陆交界区海底地形资料获取，以及全球海底地形模型精准构建方面优势大有裨益。型精准构建方面优势大有裨益。型精准构建方面优势大有裨益。
技术研发人员：范雕 李姗姗 赵东明 杜斌 胡兴猛 范昊鹏 冯进凯 单建晨 王傲明
受保护的技术使用者：中国人民解放军战略支援部队信息工程大学
技术研发日：2022.07.06
技术公布日：2022/11/1
