CN113267822B - 基于地形约束因子权重优化提高海底地形反演精度的方法 - Google Patents

基于地形约束因子权重优化提高海底地形反演精度的方法 Download PDF

Info

Publication number
CN113267822B
CN113267822B CN202110484809.6A CN202110484809A CN113267822B CN 113267822 B CN113267822 B CN 113267822B CN 202110484809 A CN202110484809 A CN 202110484809A CN 113267822 B CN113267822 B CN 113267822B
Authority
CN
China
Prior art keywords
point
inversion
topography
gravity anomaly
terrain
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110484809.6A
Other languages
English (en)
Other versions
CN113267822A (zh
Inventor
郑伟
李钊伟
孙永进
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Academy of Space Technology CAST
Original Assignee
China Academy of Space Technology CAST
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Academy of Space Technology CAST filed Critical China Academy of Space Technology CAST
Priority to CN202110484809.6A priority Critical patent/CN113267822B/zh
Publication of CN113267822A publication Critical patent/CN113267822A/zh
Application granted granted Critical
Publication of CN113267822B publication Critical patent/CN113267822B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00Measuring gravitational fields or waves; Gravimetric prospecting or detecting
    • G01V7/02Details
    • G01V7/06Analysis or interpretation of gravimetric records

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种基于地形约束因子权重优化提高海底地形反演精度的方法,包括:确定Pi点对应位置处的自由空气重力异常和短波重力异常根据解算得到Pi点对应位置处的长波重力异常根据和分配给每个Pi点的权重值λi,通过地形约束因子权重优化法,求解得到P′点对应位置处的预测长波重力异常根据解算得到预测短波重力异常根据进行海底地形反演,得到海底地形反演结果Epp(x0,y0)。本发明提高了低频重力场模型构建的精度,为高精度海底地形反演提供较好的基础数据。

Description

基于地形约束因子权重优化提高海底地形反演精度的方法
技术领域
本发明属于海洋测绘学、地球重力学等交叉技术领域,尤其涉及一种基于地形约束因子权重优化提高海底地形反演精度的方法。
背景技术
海洋面积约占全球总面积的71%,并且蕴藏着丰富的自然资源。熟知海底地形将为合理开发和利用海洋资源提供基础信息;同时,丰富的海底地形信息结合重力场信息为航行安全将提供基础保障。传统水深测量方法以测船为载体搭载声纳设备开展,测深精度高,但是对于大范围海底地形测绘,存在耗时耗力的缺陷。而借助卫星测高技术或者重力卫星技术能够高效地恢复全球海洋重力场信息。基于重力信息与海底地形之间的相关性可以实现海底地形的快速构建,具有全天时、全天候且覆盖范围广、经济实用的优势,从而在一定程度上弥补了传统测深方式的不足。
自1978年,第一颗海洋测高卫星Seasat成功发射以来,众多学者试图利用获取的海面测高数据反演海底地形。1983年,Dixon等利用线性响应函数技术对大地水准面高度与海底地形的相关性进行了分析,并利用卫星测高反演的海洋大地水准面对夏威夷北部海山区域的海底地形进行了反演。随后部分学者对该研究开展了相应分析,并在此技术基础上进行了相应改进与应用。随着全球测高卫星数量的激增,例如:Geosat、ERS-1/2、TOPEX/Poseidon、Jason-1/2/3、ENVISAT、Cryosat-1/2、HY-2等,利用测高卫星反演高精度海洋重力场模型研究逐步受到重视,从而利用卫星测高技术开展海底地形反演应用也得以快速发展。现阶段基于卫星测高反演海底地形方法主要包括:导纳函数法、重力地质法(gravitygeologic method,GGM)、最小二乘配置法等。
重力地质法在船载测深控制点的约束条件下,将重力场信息划分为高频与低频两部分,基于高频信息与海底地形较强的相关性反演海底地形,公式简单且精度较好,因而受到各国学者青睐,并利用GGM法开展了大量海底地形反演研究工作。2011年,Kim等基于调优密度差GGM法得到德雷克海峡区域2'×2'空间分辨率的海底地形模型,相比传统GGM法精度提升了29m;2013年,Hu等基于GGM法在低频重力场计算过程中直接引入地形因子对算法进行了改进,并将其应用于印度洋Rodriguez三联点海域与南中国海区域,反演模型精度优于ETOPO1模型;2014年,Ouyang等利用GGM法获取了中国南海1'×1'空间分辨率海底地形,通过分析船载测深控制点分布情况和海底地形起伏情况对反演精度的影响,得出GGM法反演精度受海底地形变化影响的结论,最后综合ETOPO1模型与船载测深数据构建了该区域内海底地形模型;2017年,Xiang等针对GGM法低频重力场模型构建展开研究,提出非均匀控制点的自适应三角网有限元逼近法对低频重力场进行网格化处理,利用改进的GGM法在中国东海与南海区域进行实验验证,说明了该算法具备良好的稳定性;2018年,Kim和Yun尝试将GGM法应用于浅海区域海底地形反演研究,预测了朝鲜半岛塞慢吉姆堤海域的海底地形,与船载测深数据之间的均方根误差达到0.6m,验证了GGM法在浅海区域的有效性;2020年,Annan和Wan利用改进的GGM法反演构建了几内亚湾区域海底地形模型,同NGDC、ETOPO1和SIO模型对比该算法反演精度良好,研究表明在海山区域需要较多的船载测深控制点,且在整个区域不能采用单一密度差进行计算。
受限于船载测深技术固有的弊端,全球范围内船载测深控制点分布相对稀疏,而GGM法受控制点的约束,特别对低频重力场的格网化构建存在较大影响,从而影响海底地形构建的精度。然而,已有研究中的格网化通常利用GMT软件自带的张力样条算法处理,鲜有对低频重力模型构建的研究。而现阶段常用的格网化构建方法众多,主要包括张力样条插值、反距离加权插值法、多项式插值法、径向基函数插值法及各自的改进算法等,上述各方法直接通过周围控制点的值,或者利用特定数学公式进行,较少考虑已知点属性的结构特点及整体空间分布情况,易出现过拟合或拟合不足问题。
发明内容
本发明的技术解决问题:克服现有技术的不足,提供一种基于地形约束因子权重优化提高海底地形反演精度的方法,旨在提高低频重力场模型构建的精度,为高精度海底地形反演提供较好的基础数据。
为了解决上述技术问题,本发明公开了一种基于地形约束因子权重优化提高海底地形反演精度的方法,包括:
确定Pi点对应位置处的自由空气重力异常和短波重力异常其中,Pi点为已知控制点,Pi=(xi,yi),xi和yi分别表示Pi点的纬度和经度,i=1,2,...n,n为测区范围内已知控制点的总数;
根据和/>解算得到Pi点对应位置处的长波重力异常
根据和分配给每个Pi点的权重值λi,通过地形约束因子权重优化法,求解得到P′点对应位置处的预测长波重力异常/>其中,P′点为预测控制点,P′=(x0,y0),x0和y0分别表示P′点的纬度和经度;
根据解算得到预测短波重力异常/>
根据进行海底地形反演,得到海底地形反演结果Epp(x0,y0)。
在上述基于地形约束因子权重优化提高海底地形反演精度的方法中,确定Pi点对应位置处的自由空气重力异常包括:通过双线性内插,从测高卫星反演的重力异常模型中求取得到Pi点对应位置处的自由空气重力异常/>
在上述基于地形约束因子权重优化提高海底地形反演精度的方法中,确定Pi点对应位置处的短波重力异常包括:
获取Pi点对应位置处的水深Ecp(xi,yi);
通过如下式(1),解算得到Pi点对应位置处的短波重力异常
其中,G表示万有引力常数,Δρ表示海水与海底基岩的密度差,D表示测区范围内控制点的最大深度值。
在上述基于地形约束因子权重优化提高海底地形反演精度的方法中,根据和/>解算得到Pi点对应位置处的长波重力异常/>包括:
根据和/>通过如下式(2),解算得到Pi点对应位置处的长波重力异常/>
在上述基于地形约束因子权重优化提高海底地形反演精度的方法中,地形约束因子权重优化法求解的表达式如下:
其中,表示优化后的权重值,Pj点为另一已知控制点,Pj=(xj,yj),xj和yj分别表示Pj点的纬度和经度,j=1,2,...n,i≠j,/>表示优化后的Pi点与Pj点之间的变异函数γ(Lij,Zij)的最佳拟合值,/>表示优化后的Pi点与P′点之间的变异函数γ(Li0,Zi0)的最佳拟合值,Lij表示Pi点和Pj点在水平方向上的滞后距离,Zij表示Pi点和Pj点在垂直方向上的滞后距离,φ表示拉格朗日系数,Li0表示Pi点和P′点在水平方向上的滞后距离,Zi0表示Pi点和P′点在垂直方向上的滞后距离。
在上述基于地形约束因子权重优化提高海底地形反演精度的方法中,Lij、Zij、Li0和Zi0的解算公式如下:
Zij=|zi-zj|···(5)
Zi0=|zi-z0|···(7)
其中,zi、zj和z0分别为Pi点对应的地形因子、Pj点地形因子和P′点对应的地形因子。
在上述基于地形约束因子权重优化提高海底地形反演精度的方法中,γ(Lij,Zij)的解算流程如下:
确定Pi点和Pj点在水平方向上的实验变异函数γ(Lij):
其中,E[·]为数学期望计算符号;
确定Pi点和Pj点在垂直方向上的实验变异函数γ(Zij):
其中,和/>分别为zi和zj对应位置处的长波重力异常;
确定γ(Lij,Zij):
在上述基于地形约束因子权重优化提高海底地形反演精度的方法中,γ(Li0,Zi0)的解算流程如下:
确定Pi点和P′点在水平方向上的实验变异函数γ(Li0):
其中,表示P′点的长波重力异常实际值;
确定Pi点和P′点在垂直方向上的实验变异函数γ(Zi0):
其中,表示z0对应位置处的长波重力异常实际值;
确定γ(Li0,Zi0):
在上述基于地形约束因子权重优化提高海底地形反演精度的方法中,根据解算得到预测短波重力异常/>包括:
通过双线性内插,从测高卫星反演的重力异常模型中求取得到Pi点对应位置处的自由空气重力异常
通过如下式(14),解算得到预测短波重力异常
在上述基于地形约束因子权重优化提高海底地形反演精度的方法中,根据进行海底地形反演,得到海底地形反演结果Epp(x0,y0),包括:
根据通过如下式(15),进行海底地形反演,得到海底地形反演结果Epp(x0,y0):
本发明具有以下优点:
(1)本发明公开了一种基于地形约束因子权重优化提高海底地形反演精度的方法,基于新型TCFWO法可有效提高海底地形反演精度,且在船载测深控制点分布稀疏区域更能体现其优势,计算速度快。
(2)本发明公开了一种基于地形约束因子权重优化提高海底地形反演精度的方法,充分考虑海底地形信息对普通Kriging法变异函数的影响,分别从水平与深度方向构建变异函数模型,并顾及到水平方向变异函数各向异性影响,通过增加外部相关信息重新优化Kriging法格网化的权重,以提高低频重力场模型构建的精度,为高精度海底地形反演提供较好的基础数据。
附图说明
图1是本发明实施例中一种基于地形约束因子权重优化提高海底地形反演精度的方法的步骤流程图;
图2是本发明实施例中一种研究区域数据的示意图;2(a)为船载测深控制点及检核点分布,2(b)为卫星测高重力异常;
图3是本发明实施例中一种低频重力异常数据分位数图;
图4是本发明实施例中一种低频重力数据的空间结构分析结果示意图;4(a)为不同方向实验变异函数曲线,4(b)为各向异性变异函数拟合;
图5是本发明实施例中一种深度方向变异函数拟合图;
图6是本发明实施例中一种低频重力场模型示意图;6(a)为KR模型,6(b)为TCF模型;
图7是本发明实施例中一种交叉验证拟合回归曲线;7(a)为KR模型交叉验证结果,7(b)为TCF模型交叉验证结果;
图8是本发明实施例中一种密度差确定曲线;8(a)为相关系数曲线,8(b)为标准差曲线;
图9是本发明实施例中一种海底地形模型示意图;9(a)为TCFWO模型,9(b)为KR_GGM模型,9(c)为ETOPO1模型,9(d)为V19.1模型;
图10是本发明实施例中一种海底地形模型局部海底地形模型示意图;图10从左向右依次为ETOPO1模型、V19.1模型、KR_GGM模型及TCFWO模型;
图11是本发明实施例中一种船测检核点差值绝对值分布直方图;11(a)为KR_GGM模型,11(b)为TCFWO模型。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚,下面将结合附图对本发明公开的实施方式作进一步详细描述。
本发明涉及一种基于地形约束因子权重优化提高海底地形反演精度的方法。采用重力地质法GGM反演海底地形时,受船载测深控制点分布的约束,需要对低频重力异常数据进行格网化构建,构建模型的准确性直接影响海底地形反演精度。第一,本发明在普通Kriging法基础上,通过引入地形约束因子,提出了用于低频重力场模型构建的新型地形约束因子权重优化法(terrain constraint factor weight optimization,TCFWO)。第二,采用新型TCFWO法与普通Kriging法分别构建低频重力场模型(TCF模型与KR模型),利用交叉验证分析,TCF模型较KR模型的精度提升近40%,从而验证了新型TCFWO法的可靠性。第三,将构建的低频重力场模型分别应用到海底地形反演中,与船载测深检核点验证结果表明:研究区域内基于新型TCFWO法反演的TCFWO模型的精度优于ETOPO1模型和利用Kriging法构建的KR_GGM模型;在研究区内的控制点分布稀疏部分,TCFWO模型精度均优于国际常用模型V19.1和ETOPO1以及KR_GGM模型;相较KR_GGM模型,TCFWO模型精度提升约26%。结果表明,基于新型TCFWO法可有效提高海底地形反演精度,且在船载测深控制点分布稀疏区域更能体现其优势。新型地形约束因子权重优化方法的优点为海底地形反演精度高,计算速度快。
如图1,在本实施例中,该基于地形约束因子权重优化提高海底地形反演精度的方法,包括:
步骤101,确定Pi点对应位置处的自由空气重力异常和短波重力异常
在本实施例中,利用GGM反演海底地形时,可基于稀疏船载测深数据将卫星测高获得的自由空气重力异常Δgf(x,y),划分为长波重力异常Δgl(x,y)和短波重力异常Δgs(x,y),即本发明实施例所使用的通用基本公式如下:
Δgf(x,y)=Δgl(x,y)+Δgs(x,y)···(0-1)
其中,x和y分别表示纬度和经度。
进一步的,短波重力异常Δgs(x,y)与地形的关系的通用基本公式如下:
其中,E(x,y)表示(x,y)点对应位置处的水深;G表示万有引力常数,G=6.672×10-11N·m2/kg2;Δρ表示海水与海底基岩的密度差,用于调整地形和重力之间的近似线性关系;D表示测区范围内控制点的最大深度值。
可见,如果已知Pi点对应位置处的水深Ecp(xi,yi),则可基于公式0-2,解算得到Pi点对应位置处的短波重力异常
其中,Pi点为已知控制点(船载测深控制点),Pi=(xi,yi),xi和yi分别表示Pi点的纬度和经度,i=1,2,...n,n为测区范围内已知控制点的总数。
进一步的,可以通过双线性内插,从测高卫星反演的重力异常模型中求取得到Pi点对应位置处的自由空气重力异常
步骤102,根据和/>解算得到Pi点对应位置处的长波重力异常/>
在本实施例中,基于公式0-1,可以确定Pi点对应位置处的长波重力异常的解算公式如下:
步骤103,根据和分配给每个Pi点的权重值λi,通过地形约束因子权重优化法,求解得到P′点对应位置处的预测长波重力异常/>
在本实施例中,显然,采用船载测深控制点计算获取的长波重力异常是离散分布的。因此,有必要对其进行网格化。相比其他格网化方法,Kriging算法可以基于变异函数模型得到预测点位置处的无偏最优估值。因此,在Kriging算法基础上,本文提出了地形约束因子权重优化法构建长波重力异常模型。具体的:首先,基于根据式(2)解算得到然后,基于Kriging算法,利用下述公式对P′点位置处的长波重力异常进行估计:
其中,P′点为预测控制点,P′=(x0,y0),x0和y0分别表示P′点的纬度和经度。
由此可见,如何确定权重值λi是至关重要的,也是TCFWO法的关键。
TCFWO法同Kriging方法一样,其解算的估计值具有无偏与最优的特点。因此,应满足以下方程:
其中,表示P′点对应位置处的长波重力异常实际值,Var[·]为方差解算符号,E[·]为数学期望计算符号。
此外,TCFWO法是在假设变量满足二阶平稳的基础上的,这个假设需要满足两个基本条件:
基本条件1:任意控制点对应位置处的长波重力异常的平均值存在且为常数。表达式如下:
基本条件2:任意两个控制点对应位置处的长波重力异常之间的协方差只与两个控制点之间的距离和方向有关。表达式如下:
其中,K表示任意常数,K∈D;表示点(x,y)对应位置处的长波重力异常;Cov[·]是协方差计算符号;/>表示点(xi,yi)对应位置处的长波重力异常,表示点(xj,yj)对应位置处的长波重力异常,h表示/>之间的距离,f(h)是只与h有关的函数。
根据TCFWO法估计值的无偏与误差最优特点,联合公式0-3和公式(0-4),通过构造代价函数F将解算权重值λi问题转化为计算带有约束条件的方程最优化问题,构造的代价函数F如下:
其中,φ表示拉格朗日系数,γi0(d)表示Pi点和P′点之间的变异函数值,γij(d)表示Pi点和Pj点之间的变异函数值。变异函数值是描述数据空间相关性的一个统计量,它被定义为两点之间差异的方差。γi0(d)和γij(d)可以通过如下实验变异函数方程来确定:
其中,和/>分别表示m点和m+d点位置处的长波重力异常,d为m点和m+d点之间的欧氏距离,即滞后距离。
随后,对代价函数F分别求λi和φ的偏导数,并令求导后的结果为零,最终得到确定权重值的方程组如下:
根据公式0-9可知,变异函数值决定了插值的权重。而普通Kriging算法在解算变异函数时,只考虑了二维变量因素对该值的影响。为了优化公式0-9中的权重值,TCFWO法通过引入地形因子zi作为基于经纬度坐标的第三变量因素。因此,可以分别建立水平方向和垂直方向上的重力异常的变异函数模型。
基于以上分析可知,所述步骤103中所述的地形约束因子权重优化法求解的表达式如下:
其中,表示优化后的权重值,Pj点为另一已知控制点,Pj=(xj,yj),xj和yj分别表示Pj点的纬度和经度,j=1,2,...n,i≠j,/>表示优化后的Pi点与Pj点之间的变异函数γ(Lij,Zij)的最佳拟合值,/>表示优化后的Pi点与P′点之间的变异函数γ(Li0,Zi0)的最佳拟合值,Lij表示Pi点和Pj点在水平方向上的滞后距离,Zij表示Pi点和Pj点在垂直方向上的滞后距离,Li0表示Pi点和P′点在水平方向上的滞后距离,Zi0表示Pi点和P′点在垂直方向上的滞后距离。
优选的,Lij、Zij、Li0和Zi0的解算公式如下:
Zij=|zi-zj|···(5)
Zi0=|zi-z0|···(7)
其中,zi、zj和z0分别为Pi点对应的地形因子、Pj点地形因子和P′点对应的地形因子。
优选的,γ(Lij,Zij)的解算流程如下:
首先,确定Pi点和Pj点在水平方向上的实验变异函数γ(Lij):
然后,确定Pi点和Pj点在垂直方向上的实验变异函数γ(Zij):
其中,和/>分别为zi和zj对应位置处的长波重力异常。
最后,确定γ(Lij,Zij):
优选的,γ(Li0,Zi0)的解算流程如下:
首先,确定Pi点和P′点在水平方向上的实验变异函数γ(Li0):
其中,表示P′点的长波重力异常实际值。
然后,确定Pi点和P′点在垂直方向上的实验变异函数γ(Zi0):
其中,表示z0对应位置处的长波重力异常实际值。
最后,确定γ(Li0,Zi0):
步骤104,根据解算得到预测短波重力异常/>
在本实施例中,可通过双线性内插,从测高卫星反演的重力异常模型中求取得到Pi点对应位置处的自由空气重力异常然后,基于通用基本公式0-1可知,可以通过如下式(14),解算得到预测短波重力异常/>
步骤105,根据进行海底地形反演,得到海底地形反演结果Epp(x0,y0)。
在本实施例中,基于通用基本公式0-2可知,可以通过如下式(15),进行海底地形反演,得到海底地形反演结果Epp(x0,y0):
综上所述,本发明重点对低频重力场网格化构建方法展开探究,提出了新型地形约束因子权重优化法,TCFWO法充分考虑海底地形信息对普通Kriging法变异函数的影响,分别从水平与深度方向构建变异函数模型,并顾及到水平方向变异函数各向异性影响,通过增加外部相关信息重新优化Kriging法格网化的权重,以提高低频重力场模型构建的精度,为高精度海底地形反演提供较好的基础数据。
在上述实施例的基础上,下面对新型地形约束因子权重优化法进行验证说明。
1、新型地形约束因子权重优化法的验证
1.1、实验区域概述及数据准备
高精度的海底地形模型对自然资源的开发利用及海洋地质学的研究意义重大,本发明选取西太平洋马库斯-韦克海山群区域(156.00°E—164.47°E,17.88°N—26.26°N)作为研究对象,该地区地形变化复杂,地壳运动活跃,火山活动强烈,生物多样性繁多,自然资源富足遍布钴、锰、磷等矿产资源。
实验数据来源主要包括船载测深数据集、卫星测高反演的海洋重力数据及结果对比的其他海底地形模型。(1)由美国国家地球物理数据中心(NGDC,http://www.ngdc.noaa.gov/mgg/geodas/trackline.html)提供的船载测深数据(https://www.noaa.gov),区域内共计船载测深点53863个,最大深度为-6190.00m,最浅为-1038.00m,平均水深-4782.68m。船载数据点分布如图2(a)所示,每隔4个点选取1个检核点,构成水深反演结果及密度差确定的检核数据集,如图中圆点所示,共计10777个检核点。图中星点作为控制点共计43086个,以此对海底地形反演过程中涉及的参数及高频重力场进行解算;(2)由美国斯克里普斯海洋研究所与加州大学(SIO,Scripps Institution ofOceanography;USGD,University of California San Diego)提供的卫星测高反演的海洋重力异常数据(https://topex.ucsd.edu)如图2(b)所示,版本为V28.1,空间分辨率为1′×1′,区域内最大重力值为254.80mGal,最小为-63.00mGal,平均值为-4.69mGal。此外,还包括用于结果对比分析的海底地形模型数据ETOPO1(https://www.ngdc.noaa.gov)及V19.1(https://topex.ucsd.edu)。
1.2、低频重力信息地统计学分析
采用新型TCFWO法对离散低频重力异常数据进行网格化重建过程中,首先应保证区域变量满足二阶平稳假设。其次,需要通过变异函数了解变量的空间变异结构,从而保证构建的结果有效性与真实性。因此,需要对数据进行相应的地统计学分析,以全面掌握变量的结构特征。一般分别采用分位数图(QQ图,Quantile-Quantile Plot)和变异函数模型对数据的平稳性假设及数据的结构信息进行定量描述。
1.2.1、二阶平稳假设分析
二阶平稳性假设是TCFWO法的前提条件,该假设需要满足两个基本条件:(1)随机函数的均值存在且为常数;(2)任意两个空间变量之间的协方差仅与他们之间的距离和方向有关。数学公式表述如前所述,在此不再赘述。
对低频重力异常数据进行二阶平稳假设统计时,若近似满足正态分布,即认为满足上述两个基本条件。数据的正态分布检验常采用QQ图进行分析,基于分位图的思想绘制的分位数曲线如图3所示。图3中实线表示正态分布情况,散点表示低频重力样本数据,若满足正态分布,则散点数据应与直线基本重合。据图3可知,数据基本分布在直线附近,说明低频重力数据近似满足正态分布。因此,可以认为低频重力数据满足二阶平稳假设。
1.2.2、空间变异结构分析
变量在空间分布的结构特征,在构建格网化模型中发挥重要作用。地统计学中可通过变异函数模型进行分析,并且模型决定了TCFWO法插值的精度与准确性。因此,格网化之前需要首先对变量的变异函数结构进行分析,主要涉及变异函数的计算拟合及空间各向异同性分析。
(1)水平方向变异函数模型构建
空间各向异性表示在空间中各个方向的空间变异性存在差异,对应于变异函数模型中,不同方向上的变程与基台值都不完全相同,若相同则称其满足各向同性。利用公式(4)和(5)计算各低频重力异常数据在水平方向上的距离作为滞后距;根据公式(8)~(10)分别计算在0°、45°、90°及135°四个方向的实验变异函数值。最终,绘制的不同方向上实验变异函数曲线如图4(a)所示。根据变化曲线可知,在不同方向上变异函数的变程及基台值均存在一定差异,可以认为低频重力场在水平方向为各向异性。
因此,在变异函数构建过程中采用各向异性进行处理。本发明将四个方向的实验变异函数取平均值,对水平实验变异函数进行构建。实验样本分布如图4(b)中星点所示,分别基于理论变异函数指数函数模型、球状函数模型及高斯函数模型,利用最小二乘算法对实验变异函数曲线进行拟合。三种理论变异函数拟合的曲线结果分别见图4(b)中三角线、方格线及实线。
根据图4(b)中水平方向变异函数拟合图及参数可知,采用球状函数模型及高斯函数模型拟合的结果基本一致,拟合度分别为0.45与0.44;而采用指数函数模型拟合的变异函数曲线的拟合度为0.59,相对另外两种函数模型,此模型拟合精度更优。因此,本发明采用指数函数模型进行变异函数值解算。
(2)深度方向变异函数模型构建
首先,分别基于公式(4)~(10)解算滞后距与深度方向的实验变异函数值。假设数据在深度方向满足各向同性,然后按照水平方向变异函数拟合方法,低频重力数据在垂直方向上的实验样本数据及拟合曲线如图5所示。
据图5可知,采用球状函数模型及高斯函数模型的拟合度分别为0.94与0.93,其拟合结果近似;相较而言,采用指数函数模型拟合度为0.97,其拟合精度更优。故在深度方向上变异函数模型构建中采用指数函数模型。最后,基于公式(8)~(10)对两方向的变异函数进行处理。
1.3、低频重力场模型精度评定
为了验证新型TCFWO法的适用性,通过对比TCFWO法与Kriging法在低频重力场格网化中的效果,采用两算法对低频重力场数据进行格网化处理。首先,基于两方法分别构建的空间分辨率1′×1′的低频重力场模型为TCF模型(图6(a))与KR模型(图6(b))。然后,对格网化结果采用交叉验证法评定精度,本发明选取数据集中1000个数据点作为真实值,分别与预测值进行对比分析。以预测值为横轴,真实值为纵轴,两方法构建的散点图及其回归拟合曲线如图7所示。
据图6可知,KR模型最大值、最小值、平均值及标准差分别为:-97.18mGal、49.70mGal、-43.94mGal与16.15mGal,而TCF模型的最值均大于KR模型,标准差略低于KR模型,可以看出TCFWO法虽然提高了反演低频重力场的变化范围,但是模型稳定性有所提升。通过对比两交叉验证的回归曲线图7可知,KR模型交叉验证结果分布相比TCF模型的结果更加离散,可见TCFWO法的格网化结果与实际数据更加接近。
进一步分析交叉验证残差统计,可知KR模型的交叉验证的残差最小值为-48.42mGal,最大值为39.08mGal,平均值为-0.13mGal,标准差为6.08mGal,线性回归的拟合度为0.92;而采用TCF模型残差的最大值为20.23mGal,最小值为-31.51mGal,平均值0.07mGal,标准差为3.62mGal,线性回归的拟合度为0.98。相较而言,TCF模型相较KR模型的精度提高了约40%。由此可知,TCFWO法对低频重力场格网化的效果更优。下面将基于两种算法构建的低频重力场数据进行海底地形反演研究,进一步阐明TCFWO法的优越性。
2、海底地形反演实验应用
利用不同算法构建低频重力场模型,其最终应用即为海底地形反演提供基础数据源。据GGM原理可知,低频重力场模型的优劣直接关系到海底地形反演精度,因此判定TCFWO法的优劣还需要继续对地形反演结果进行更深入的分析。
2.1、密度差确定
密度差即海水与基岩的密度差别,基于GGM法的海底地形反演时,密度差Δρ的确定至关重要,且反演精度受该值影响较大。因此,需要选择较为合理的密度差。现阶段常用方法包括迭代法和向下延拓法,其中迭代法因其简单方便受到研究者的广泛应用。
本发明基于Kim J.B.提出的迭代密度差常数确定方法,选取43086个测深控制点和10777个检核点,利用控制点计算不同密度差常数下的海底地形模型,采用双线性内插得到检核点位置的预测水深,通过比较预测值与实际值之间相关系数与差值,分别绘制相关系数变化曲线图8(a)和差值标准差变化曲线图8(b)。
据图8可知,密度差常数在一定范围内,随着数值的增大相关系数随之增大,差值标准差逐渐减小;当密度差大于某一常数时,相关系数随密度差增大而降低,标准差逐渐增大。通过上述解算拟合的两曲线,最终确定当密度差常数取值约为1.11g/cm3时,拟合曲线的相关系数取最大值,同时标准差达到最小值。因此,本发明中研究区域的最优密度差常数约为1.11g/cm3,并以此开展海底地形反演研究。
2.2、海底地形模型精度分析
为进一步验证TCFWO法在低频重力场构建过程中的优势,下面重点从反演的海底地形结果展开分析。首先,利用卫星测高重力异常模型V28.1分别减去构建的TCF与KR低频重力场模型,得到高频重力异常模型;然后,基于GGM法反演得到1′×1′空间分辨率的海底地形模型TCFWO(图9(a))及KR_GGM(图9(b))。最后,同现今国际上通用的海底地形模型ETOPO1(图9(c))与V19.1模型(图9(d))对比,并利用10777个船载测深检核点分析四种模型的精度。
对比上述四种海底地形模型可知,基于GGM法反演的海底地形与现有模型的地形变化趋势基本一致。因此,GGM方法在海底地形反演方面具备较好优势。为进一步分析不同海底地形模型差异,分别计算各模型之间差值。TCFWO模型与KR_GGM模型之间的差异标准差为42.50m,然而TCFWO模型同V19.1模型和ETOPO1模型的标准差分别为213.85m和208.83m。KR_GGM模型与V19.1和ETOPO1模型的差值均值分别为-4.90m和-14.30m,标准差为221.75m和210.75m,其结果均大于TCFWO模型的统计值。由此可知,TCFWO模型与KR_GGM模型相比,更接近国际通用模型。
随后,采用船载测深检核点分析海底地形模型精度,即将四种不同海底地形模型利用双线性内插算法插值到10777个测船检核点上,通过比较预测值与实际值之间的差异,TCFWO模型和KR_GGM模型与检核点差异的最小值分别为-1677.20m和-1739.77m,最大值分别为2335.16m和2569.20m,两模型差异最值的结果明显大于V19.1模型,而小于ETOPO1模型;对比上述四种模型的标准差可以发现,研究区域内模型精度由高到低依次为V19.1模型、TCFWO模型、KR_GGM模型和ETOPO1模型,标准差分别为117.86m、129.14m、173.51m和191.71m,通过模型与测船检核数据的相关系数也可看出其精度的高低。相较而言,TCFWO模型较KR_GGM模型反演精度提升约26%。通过上述分析,反映出本发明提出的TCFWO法在海底地形反演过程中,相较普通Kriging算法具备明显优势,有效提高了海底地形反演精度。
利用GGM方法进行海底地形反演时,海底地形精度受控制点分布密集情况的影响较大。为进一步分析上述四种模型受控制点密集程度影响的情况,根据研究区域控制点分布的特点,截取船载测深点分布相对稀疏的范围(160.00°E—164.47°E,17.88°N—26.26°N)展开研究。并截取范围内的海底地形结果如图10所示,利用双线性插值计算测船检核点(3572个)位置处水深值。深入分析KR_GGM模型与TCFWO模型的差异情况,分别对两模型的检核差值的绝对值进行统计,绘制频率分布直方图分别如下图11(a)和11(b)所示。TCFWO模型的极值小于其余三种模型,并且TCFWO模型与检核点差值的标准差为152.88m,可见TCFWO模型的精度均优于其余3种模型,同KR_GGM模型相比精度提升约26.35%。尽管ETOPO1模型的精度最低,但是他与整体模型之间的标准差变化率最低,说明此模型的整体精度相近,表现出较好的稳定性。TCFWO模型的标准差变化率为18.35%,略低于ETOPO1模型。模型精度波动最大的是V19.1模型,标准差变化率为35.47%,反映出海底地形反演的精度受到船测控制分布的疏密情况影响较大。相较而言,四类模型精度均低于各自整体的精度,因此若要构建更高精度的海底地形加密船载测深点更为重要。
图11中横坐标表示差值绝对值,数值表示该差值不超过下标标注的范围,纵轴表示差异的频数。据图10可知,检核点差值绝对值不大于50m的点数,KR_GGM模型占比约44.57%,TCFWO模型占比约57.33%。TCFWO模型中约70%的点,其差值绝对值不大于100m;KR_GGM模型的差值绝对值小于500m占比约96.40%,TCFWO模型占比近98.16%,可见TCFWO模型更接近测船检核点深度。因此,TCFWO法在控制点分布稀疏区域反演的海底地形在精度及稳定性上优于Kriging法。
综合以上海底地形反演的结果可知,在整个研究区域,TCFWO模型精度优于KR_GGM模型与ETOPO1模型,精度与V19.1模型相近;在控制点分布稀疏区域TCFWO模型精度高于国际通用的海底地形模型(ETOPO1及V19.1)和KR_GGM模型;此外,TCFWO模型表现出较好的稳定性。因此,TCFWO方法可以有效地提高海底地形反演精度与模型稳定性,并且在船载测量区域稀疏的情况下优势更明显。
本发明虽然已以较佳实施例公开如上,但其并不是用来限定本发明,任何本领域技术人员在不脱离本发明的精神和范围内,都可以利用上述揭示的方法和技术内容对本发明技术方案做出可能的变动和修改,因此,凡是未脱离本发明技术方案的内容,依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化及修饰,均属于本发明技术方案的保护范围。
本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。

Claims (9)

1.一种基于地形约束因子权重优化提高海底地形反演精度的方法,其特征在于,包括:
确定Pi点对应位置处的自由空气重力异常和短波重力异常/>其中,Pi点为已知控制点,Pi=(xi,yi),xi和yi分别表示Pi点的纬度和经度,i=1,2,...n,n为测区范围内已知控制点的总数;
根据和/>解算得到Pi点对应位置处的长波重力异常
根据和分配给每个Pi点的权重值λi,通过地形约束因子权重优化法,求解得到P′点对应位置处的预测长波重力异常/>其中,P′点为预测控制点,P′=(x0,y0),x0和y0分别表示P′点的纬度和经度;
根据解算得到预测短波重力异常/>
根据进行海底地形反演,得到海底地形反演结果Epp(x0,y0);
地形约束因子权重优化法求解的表达式如下:
其中,表示优化后的权重值,Pj点为另一已知控制点,Pj=(xj,yj),xj和yj分别表示Pj点的纬度和经度,j=1,2,...n,i≠j,/>表示优化后的Pi点与Pj点之间的变异函数γ(Lij,Zij)的最佳拟合值,/>表示优化后的Pi点与P′点之间的变异函数γ(Li0,Zi0)的最佳拟合值,Lij表示Pi点和Pj点在水平方向上的滞后距离,Zij表示Pi点和Pj点在垂直方向上的滞后距离,φ表示拉格朗日系数,Li0表示Pi点和P′点在水平方向上的滞后距离,Zi0表示Pi点和P′点在垂直方向上的滞后距离。
2.根据权利要求1所述的基于地形约束因子权重优化提高海底地形反演精度的方法,其特征在于,确定Pi点对应位置处的自由空气重力异常包括:通过双线性内插,从测高卫星反演的重力异常模型中求取得到Pi点对应位置处的自由空气重力异常
3.根据权利要求2所述的基于地形约束因子权重优化提高海底地形反演精度的方法,其特征在于,确定Pi点对应位置处的短波重力异常包括:
获取Pi点对应位置处的水深Ecp(xi,yi);
通过如下式(1),解算得到Pi点对应位置处的短波重力异常
其中,G表示万有引力常数,Δρ表示海水与海底基岩的密度差,D表示测区范围内控制点的最大深度值。
4.根据权利要求3所述的基于地形约束因子权重优化提高海底地形反演精度的方法,其特征在于,根据和/>解算得到Pi点对应位置处的长波重力异常包括:
根据和/>通过如下式(2),解算得到Pi点对应位置处的长波重力异常/>
5.根据权利要求4所述的基于地形约束因子权重优化提高海底地形反演精度的方法,其特征在于,Lij、Zij、Li0和Zi0的解算公式如下:
Zij=|zi-zj|···(5)
Zi0=|zi-z0|···(7)
其中,zi、zj和z0分别为Pi点对应的地形因子、Pj点地形因子和P′点对应的地形因子。
6.根据权利要求5所述的基于地形约束因子权重优化提高海底地形反演精度的方法,其特征在于,γ(Lij,Zij)的解算流程如下:
确定Pi点和Pj点在水平方向上的实验变异函数γ(Lij):
其中,E[·]为数学期望计算符号;
确定Pi点和Pj点在垂直方向上的实验变异函数γ(Zij):
其中,和/>分别为zi和zj对应位置处的长波重力异常;
确定γ(Lij,Zij):
7.根据权利要求6所述的基于地形约束因子权重优化提高海底地形反演精度的方法,其特征在于,γ(Li0,Zi0)的解算流程如下:
确定Pi点和P′点在水平方向上的实验变异函数γ(Li0):
其中,表示P′点的长波重力异常实际值;
确定Pi点和P′点在垂直方向上的实验变异函数γ(Zi0):
其中,表示z0对应位置处的长波重力异常实际值;
确定γ(Li0,Zi0):
8.根据权利要求7所述的基于地形约束因子权重优化提高海底地形反演精度的方法,其特征在于,根据解算得到预测短波重力异常/>包括:
通过双线性内插,从测高卫星反演的重力异常模型中求取得到Pi点对应位置处的自由空气重力异常
通过如下式(14),解算得到预测短波重力异常
9.根据权利要求8所述的基于地形约束因子权重优化提高海底地形反演精度的方法,其特征在于,根据进行海底地形反演,得到海底地形反演结果Epp(x0,y0),包括:
根据通过如下式(15),进行海底地形反演,得到海底地形反演结果Epp(x0,y0):
CN202110484809.6A 2021-04-30 2021-04-30 基于地形约束因子权重优化提高海底地形反演精度的方法 Active CN113267822B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110484809.6A CN113267822B (zh) 2021-04-30 2021-04-30 基于地形约束因子权重优化提高海底地形反演精度的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110484809.6A CN113267822B (zh) 2021-04-30 2021-04-30 基于地形约束因子权重优化提高海底地形反演精度的方法

Publications (2)

Publication Number Publication Date
CN113267822A CN113267822A (zh) 2021-08-17
CN113267822B true CN113267822B (zh) 2024-05-31

Family

ID=77229886

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110484809.6A Active CN113267822B (zh) 2021-04-30 2021-04-30 基于地形约束因子权重优化提高海底地形反演精度的方法

Country Status (1)

Country Link
CN (1) CN113267822B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113960690B (zh) * 2021-09-03 2023-05-05 中国人民解放军战略支援部队信息工程大学 一种海面重力数据测量精度对海底地形反演结果影响计算方法及装置
CN114137624B (zh) * 2021-10-27 2024-02-27 中国海洋大学 一种基于卫星高度计反演海底地形的方法和系统
CN114663640B (zh) * 2022-05-20 2022-08-26 自然资源部第二海洋研究所 基于地形地貌和构造特征的海底地理实体划定与分类方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104459795A (zh) * 2014-12-08 2015-03-25 中国科学院南海海洋研究所 一种深度变密度的地壳伸展系数热校正重力异常反演方法
CN110765686A (zh) * 2019-10-22 2020-02-07 中国人民解放军战略支援部队信息工程大学 利用有限波段海底地形进行船载声呐测深测线设计的方法
CN111142169A (zh) * 2020-02-25 2020-05-12 中国地质大学(北京) 一种基于重力梯度数据的海底地形反演方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9229121B2 (en) * 2012-03-13 2016-01-05 Seoul National University R&Db Foundation Seismic imaging system for acoustic-elastic coupled media using accumulated Laplace gradient direction

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104459795A (zh) * 2014-12-08 2015-03-25 中国科学院南海海洋研究所 一种深度变密度的地壳伸展系数热校正重力异常反演方法
CN110765686A (zh) * 2019-10-22 2020-02-07 中国人民解放军战略支援部队信息工程大学 利用有限波段海底地形进行船载声呐测深测线设计的方法
CN111142169A (zh) * 2020-02-25 2020-05-12 中国地质大学(北京) 一种基于重力梯度数据的海底地形反演方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
利用卫星重力场模型反演区域地壳密度;刘国仕;《中国优秀硕士学位论文全文数据库》;全文 *
基于重力地质法反演南海海底地形;彭聪;王鹏;周兴华;杨磊;;海洋测绘(第02期);全文 *
马里亚纳海沟Challenger Deep的岩石圈流变结构与动力学分析;高玲举;张健;吴时国;;中国科学院大学学报(第03期);全文 *
高精度测高重力场反演南海海底地形;李倩倩;鲍李峰;;海洋测绘(第02期);全文 *

Also Published As

Publication number Publication date
CN113267822A (zh) 2021-08-17

Similar Documents

Publication Publication Date Title
CN113267822B (zh) 基于地形约束因子权重优化提高海底地形反演精度的方法
Jin et al. Surface wave phase-velocity tomography based on multichannel cross-correlation
Shen et al. Crustal and uppermost mantle structure beneath the United States
Grana Bayesian linearized rock-physics inversion
Linde et al. Geological realism in hydrogeological and geophysical inverse modeling: A review
Berdichevsky et al. On two-dimensional interpretation of magnetotelluric soundings
CN111538075B (zh) 干热岩勘探方法、装置、电子设备及存储介质
Darcel et al. Stereological analysis of fractal fracture networks
US20110153285A1 (en) Method of developing a petroleum reservoir from a facies map construction
CN111506861B (zh) 一种目的层有利区域裂缝强度计算方法
Ritzwoller et al. Ability of a global three‐dimensional model to locate regional events
US20110155389A1 (en) System and Method For Providing A Physical Property Model
CN108957554B (zh) 一种地球物理勘探中的地震反演方法
US20110125469A1 (en) Method of developing a petroleum reservoir by reservoir model reconstruction
CN113296166A (zh) 裂缝模型的构建方法
CN110989021B (zh) 水深反演方法、装置和计算机可读存储介质
Xing et al. Bathymetry inversion using the modified gravity-geologic method: application of the rectangular prism model and Tikhonov regularization
Zhou et al. New Framework of Combining Observations with Topographic Slope to Estimate VS 30 and Its Application on Building a VS 30 Map for Mainland China
Koulakov et al. Three‐dimensional seismic anisotropy in the crust and uppermost mantle beneath the Taiwan area revealed by passive source tomography
Gilder et al. Geostatistical Framework for Estimation of VS 30 in Data‐Scarce Regions
He et al. Bayesian frequency-dependent AVO inversion using an improved Markov chain Monte Carlo method for quantitative gas saturation prediction in a thin layer
Ren et al. Bayesian inversion of seismic and electromagnetic data for marine gas reservoir characterization using multi-chain Markov chain Monte Carlo sampling
Spichak Modern methods for joint analysis and inversion of geophysical data
Ng et al. Reconstructing ice‐flow fields from streamlined subglacial bedforms: A kriging approach
Kamer et al. KaKiOS‐16: a probabilistic, nonlinear, absolute location catalog of the 1981–2011 Southern California Seismicity

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant