2. 海军大连舰艇学院 海洋测绘系, 辽宁 大连 116018
2. Department of Hydrography and Cartography, Dalian Naval Academy, Dalian 116018, China
构建高精度、高分辨率的海域重力场模型,对于大地测量、海洋地球物理勘探具有深远的意义[1-2].自20世纪80年代以来,卫星测高技术为海洋重力场、平均海面高的构建做出了巨大贡献[3-6].近年来,Cryosat-2卫星以其高精度海面高测量和稠密的空间分布,为构建海洋重力场提供了新的数据[7].
基于Cryosat-2数据研究海域重力场的工作在我国正在兴起.张胜军等[8-9]联合Cryosat-2与Jason-1GM数据研究了低、中、高纬度区域的海洋垂线偏差特征,讨论了该卫星数据的分辨率.万剑华等[10]利用Cryosat-2数据使用加权最小范数最小二乘解计算网格剩余垂线偏差分量的方法计算了南海的重力异常,精度达到4.5 mgal.本文基于Cryosat-2最近三年半的卫星测高GDR(geophysical data records)数据,使用海面高梯度,依据最小二乘配置方法求得了海域垂线偏差格网模型,结合移去恢复技术采用逆VeningMeinesz的球面一维傅里叶变换算法得到了我国海域重力异常.将所得重力异常格网经过船载重力测量的融合后,得到改善的我国的近海海域重力异常模型,优于国际上最近发布的DTU13模型与Sandwell V23.1模型.
1 数据与方法 1.1 Cryosat-2测高数据Cryosat-2卫星于2010年4月发射升空.该卫星的主要载荷是合成孔径干涉雷达测高仪SIRAL(synthetic aperture interferometric radar altimeter),提供的测高数据覆盖范围纬度可达到±88.0°,测量精度与之前测高卫星搭载的传统星载测高仪相比有明显提高.此外,Cryosat-2的运行轨道间距较小,约为8 km.其测高数据在我国近海海域的分布如图 1所示.文中选用了SIR_GDR_2A最近三年半的数据.海面高数据提取的编辑准则参考Cryosat-2的数据手册[11].
本文收集到的船载重力测量数据来源于我国近海海洋综合调查与评价专项(简称908专项),由不同部门在多个时间段测量获取,分布于不同区域,空间分辨率约10′,其分布如图 2所示.图中R1~R9为船载重力数据的测量区域.考虑到船载重力测量数据处理不是本文的重点,精细处理方法主要参考文献[12-16]等提出的方法进行.本文根据船测数据的分布,确定研究区域经纬度范围为东经110°~126°,北纬10°~42°.其精度指标可参考文献[12]中的相关内容.
文中选用DTU13重力场模型以及Sandwell V23.1模型作为外部检核的数据.DTU13模型由DTU08逐渐发展而来,在世界范围内得到了广泛应用,该模型仍在持续更新中,相关信息可参见文献[17].Sandwell V23.1模型于2014年发布,其中采用部分Cryosat-2和Jason1的大地测量周期内的数据,精度达到了1~2 mgal[7].
1.3 计算模型与方法 1.3.1 海面高梯度的计算先将观测到的海面高减去模型海面地形, 得到大地水准面起伏值Nres,将Nres沿卫星轨迹作距离微分,轨道误差将减小或削弱.由式(1)、式(2) 近似求得沿迹测高大地水准面梯度和方位角如下:
$ \varepsilon = - \partial N/\partial S \approx \left( {{N_{{\rm{res2}}}} - {N_{{\rm{res1}}}}} \right)/d $ | (1) |
$ \alpha = \arctan \left( {\Delta \lambda \left( {\cos {\varphi _0}} \right)/\Delta \varphi } \right) $ | (2) |
式(1)、(2) 中:ε为大地水准面梯度;Nres1、Nres2分别为计算大地水准面梯度前后两星下点处的大地水准面高;φ0、d、Δλ、Δφ分别为前后两星下点起点的纬度值、两点之间的距离、两点的经度之差、两点的纬度之差;α为沿迹大地水准面梯度方位角.
1.3.2 剩余梯度的计算重力场模型计算垂线偏差表达式[18]如下:
子午分量
$ \begin{array}{l} \xi = \frac{G}{{\gamma {\rho ^2}}}\sum\limits_{n = 2}^N {\sum\limits_{m = 0}^n {\left( {{{\bar C}_{nm}}\cos \left( {m\lambda } \right) + {{\bar S}_{nm}}\sin \left( {m\lambda } \right)} \right) \cdot } } \\ \;\;\;\;\;\;\frac{{{\rm{d}}{{{\rm{\bar P}}}_{nm}}\left( {\cos \theta } \right)}}{{{\rm{d}}\theta }} \end{array} $ | (3) |
卯酉分量
$ \begin{array}{l} \eta = \frac{G}{{\gamma {\rho ^2}}}\sum\limits_{n = 2}^N {{{\left( {\frac{R}{\rho }} \right)}^n}\sum\limits_{m = 0}^n {m\left( {{{\bar C}_{nm}}\cos \left( {m\lambda } \right) - } \right.} } \\ \;\;\;\;\;\;\left. {{{\bar S}_{nm}}\sin \left( {m\lambda } \right)} \right)\frac{{{{{\rm{\bar P}}}_{nm}}\left( {\cos \theta } \right)}}{{\cos \theta }} \end{array} $ | (4) |
式(3)、(4) 中:ξ、η分别为由重力场模型计算的垂线偏差子午分量与卯酉分量;G为万有引力常数;γ为计算点的正常重力;ρ为计算点的地心距离;R为地球平均半径;N为地球重力场模型的最高阶次;Pnm为缔合勒让德函数;θ为计算点的余纬;λ为计算点的经度;Cnm与Snm分别为重力场模型位系数.
由式(2)~(4) 计算模型垂线偏差梯度为
$ {\varepsilon _{\rm{M}}} = \xi \cos \alpha + \eta \sin \alpha $ | (5) |
式中:α含义与式(2) 相同;εM为沿迹星下点的模型垂线偏差梯度,在沿迹大地水准面梯度ε中扣除模型垂线偏差的εM,得到沿迹剩余梯度.
1.3.3 剩余垂线偏差的计算采用最小二乘配置法时,假设扰动位协方差函数K(·)各向同性.则K(·)是P、Q两点的球面角距的函数.把全球协方差函数模型应用于局部区域时,顾及参考位系数误差,则协方差函数模型为[19]
$ \begin{array}{l} K\left( {P,Q} \right) = {\left( {\frac{G}{{r\rho }}} \right)^2}{\sum\limits_{i = 2}^{\max } {{\varepsilon _i}\left( {T,T} \right)\left( {\frac{{{R^2}}}{{r\rho }}} \right)} ^i}{{\rm{P}}_i}\left( {\cos \varphi } \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {\frac{G}{{r\rho }}} \right)^2}\sum\limits_{i = \max + 1}^\infty {{\sigma _i}\left( {T,T} \right){{\left( {\frac{{R_B^2}}{{r\rho }}} \right)}^i}{{\rm{P}}_i}\left( {\cos \varphi } \right)} \end{array} $ | (6) |
式中:εi(T, T)为位系数模型的误差阶方差;σi(T, T)为位系数的阶方差;r、ρ分别为P、Q点的地心距;R为平均地球半径;RB为本杰哈马球半径.Tscherning/Rapp给出的阶方差模型为
$ {\sigma _i}\left( {T,T} \right) = \frac{A}{{\left( {i - 1} \right)\left( {i - 2} \right)\left( {i + 24} \right)}}\frac{{R_B^2}}{{r\rho }} $ | (7) |
式中:A取值为425.282 mgal2.基于协方差函数模型,建立沿迹剩余梯度与剩余垂线偏差之间的协方差矩阵以及沿迹剩余梯度的方差阵.
1.3.4 重力异常的计算利用逆VeningMeinesz公式,采用一维球面傅里叶变换的算法实现快速积分运算[18],计算的格网分辨率为2′.
$ \begin{gathered} \Delta {g_p} = \frac{{{\gamma _0}}}{{4\pi {\text{ }}}}\iint {H'\left( {{\zeta _q}\cos {\alpha _{qp}} + {\eta _q}\sin {\alpha _{qp}}} \right){\text{d}}{\sigma _q}} = \hfill \\ \;\;\;\;\;\;\;\;\;\frac{{{\gamma _0}}}{{4\pi {\text{ }}}}\iint_\sigma {H'{\varepsilon _{qp}}{\text{d}}{\sigma _q}} \hfill \\ \end{gathered} $ | (8) |
式中:p为计算点;γ0为计算点的正常重力;ζq与ηq分别为流动点垂线偏差北分量与东分量;αqp是从q点到p点的方位角; 核函数H′为
$ H' = - \frac{{\cos \frac{{{\psi _{pq}}}}{2}}}{{2\sin \frac{{{\psi _{pq}}}}{2}}} + \frac{{\cos \frac{{{\psi _{pq}}}}{2}\left( {3 + 2\sin \frac{{{\psi _{pq}}}}{2}} \right)}}{{2\sin \frac{{{\psi _{pq}}}}{2}\left( {1 + \sin \frac{{{\psi _{pq}}}}{2}} \right)}} $ | (9) |
式中:ψpq是球面距离.图 3是计算ψpq示意图.
将剩余重力异常与模型重力异常相加,即可得到基于卫星测高计算的海洋重力异常.
$ \Delta g = \Delta {g_{{\text{res}}}} + \Delta {g_{\text{G}}} $ | (10) |
式中:Δgres为测高数据计算的剩余重力异常; ΔgG为由重力场模型计算得到的重力异常.
2 计算流程(1) 从SIR_GDR_2A的三年半数据中,提取海面高,并进行高度计测量设备误差、信号传播路径误差、地球物理参数等的改正.
(2) 对海面高数据进行粗差探测,并将其剔除,得到编辑后的海面高.
(3) 从编辑后的海面高数据中,扣除海面地形模型,并对结果进行滤波,得到沿迹大地水准面高.
(4) 利用重力场模型,计算沿迹的垂线偏差分量,并计算沿迹方向的模型垂线偏差梯度.
(5) 利用沿迹大地水准面高计算海面梯度,并扣除模型大地水准面梯度,得到沿迹剩余大地水准面梯度.
(6) 利用最小二乘配置法,将沿迹剩余大地水准面梯度转换为垂线偏差格网,得到剩余垂线偏差.
(7) 利用剩余垂线偏差,结合逆VeningMeinesz公式计算剩余重力异常.
(8) 利用重力场模型计算模型重力异常.
(9) 将模型重力异常与剩余重力异常相加,得到基于Cryosat-2测高数据的海域重力异常.
3 分析与讨论为了保证文章内容表述简洁,选择南海部分海域(船载重力数据覆盖的R7、R8、R9 3个区域,见图 2)的重力异常模型构建试验区,阐述其中的关键技术.
3.1 试验区沿迹海面高梯度分析从测高的沿迹垂线偏差中减去模型的垂线偏差,沿迹剩余大地水准面梯度.Cryosat-2的星下点轨迹较为密集,在南海部分区域有160条升降弧段数据,弧段之间在经度方向的距离约为8 km.故从升弧段与降弧段中挑选4个弧段,以讨论其沿迹垂线偏差的空间分布特征,见图 4,图 5.考虑到Cryosat-2的GDR数据中不区分上升弧段与下降弧段,故先对GDR数据中的海面高数据进行编号,以区分上升弧段与下降弧段.
上升弧段77、83、145在研究区域内的起始星下点均在南海岛礁附近,其沿迹垂线偏差明显受到岛礁的影响.弧段83与弧段145上升轨迹过程中,仍然受到不能忽略的岛礁影响.在数据处理的过程中对受到陆地影响的垂线偏差,将跳跃值直接去掉,以保持数据的连续一致性,避免出现峰值.从EGM2008模型与观测的沿迹梯度来看,在上升弧段期间两者符合得很好,标准差为1.79 s,偏差平均值为0.14 s,见图 6.图 6中EGM2008表示由重力场模型EGM2008计算的梯度,OBS为实际观测数据计算的梯度.下降弧段8、66、152、162同样受到南海岛礁的影响.在近岸区域,弧段66受到的影响明显,其垂线偏差值发生剧烈变化,标准差为1.32 s,平均偏差为-0.01 s.该精度相对于上升弧段略差,见图 7.
不同阶次的参考重力场计算的重力异常与船测数据的比较见表 1~3.
从表 1~3可以看出,对于同一重力场模型,选用高阶次时,对计算结果比较有利.相对EGM2008而言,同一阶次选用EIGEN6C4重力场模型作为参考场时,其结果计算精度优于EGM2008模型.为了反映Cryosat-2反演海域重力场的潜力以及测试本文技术方法与流程的合理性,文中后续研究均选用EGM2008 2 160阶作为参考场.
3.3 试验区不同阶次参考重力模型的构建协方差函数阵对重力异常结果的影响不同重力场模型的阶次的阶方差对重力异常结果的影响见表 4.表中A表示采用EGM2008模型的最高阶2 160阶;B表示采用EGM2008模型的900阶为最高阶;C表示采用EIGEN6C4模型的900阶为最高阶.
从表 4可以看出,选择不同阶次构建的方差阵对计算结果精度的影响不大.下文中EGM2008选择2 160阶次的阶方差构建方差与协方差函数.
3.4 试验区剩余海面高低通滤波半径的确定在计算过程,需要对剩余海面高梯度进行粗差剔除以及减小高频误差的影响.分析了不同的高斯低通滤波半径对最终结果的影响.结果分别见表 5~7.
不同的滤波半径影响分为两个方面.在平均值方面,随着滤波半径的增大,平均值在增大.在精度方面,随着滤波半径的增大,标准差在减小.因为当滤波半径太大时,将海面高的高频信息都过滤了,难免对重力场的认识产生失真,故平均值有增大的趋势.为此,需要从众多可采用的滤波半径中确定最优的滤波半径.文中对该问题的解决办法是以所得结果的标准差第1次优于参考场模型时所用的半径作为最优滤波半径,进而对剩余海面高梯度进行低通滤波.经多次试验,建议45 km作为最优滤波半径.
3.5 试验区海域重力异常及其检核 3.5.1 试验区海域重力异常文中利用移去恢复技术,对逆Vening-Meinesz公式采用快速傅里叶变换方法计算得到了南海部分海域的重力异常,见图 8.
将本文仅利用Cryosat-2数据得到的重力场模型定义为CASM_GRA_ALT.为了检验本文技术方法的可靠性,将CASM_GRA_ALT、DTU13重力场模型、Sandwell V23.1模型, 分别内插船载重力测量点处的重力异常,然后将其与实际的船测重力值相减,差值的统计结果如表 8~10所示.
从表 8~10可以看出,本文采用的技术方法所得结果精度略低于DTU13重力场模型和Sandwell V23.1模型,但在平均偏差方面略优于两者.
3.5.3 试验区船测重力与测高重力数据融合首先将船载重力测量数据进行重新抽样,提取空间间距为20 km的数据点,剩余点作为检核点.然后采用最小二乘配置法进行20 km间距船测重力与测高重力CASM_GRA_ALT模型融合,得到近海海域重力异常,将该模型定义为CASM_GRA.利用剩余点船测数据对CASM_GRA精度进行检核,统计结果如表 11所示.
从表 8~11可以看出,本文所采用的技术方法与技术路线得到的南海部分海域重力场模型,其精度要好于DTU13、Sandwell V23.1模型.
3.6 中国近海海域重力异常的比较与检核借鉴试验区内构建的重力异常模型的技术方法,构建了中国近海海域重力异常模型CASM_GRA.利用DTU13全球重力场模型以及Sandwell V23.1模型内插各区域的船载重力测量点上的重力异常,并与实际的观测结果求差.各区域差值的统计情况如表 12、13所示.
从表 12、13可以看出,在渤海湾,舟山群岛及浙江沿海区域,我国的船载重力测量数据与DTU13模型、Sandwell V23.1模型有明显的系统偏差,可能是因为渤海湾是一个封闭性的海域,地理环境对测高信号的影响明显,而对于舟山群岛及浙江沿海区域,受近岸的海洋动力学环境影响比较明显.
本文CASM_GRA模型结果与未参与优化的船载重力数据的比较如表 14所示.
从表 14中可以看出,经过20 km分辨率的数据处理优化后,基于Cryosat-2数据构建的我国海域重力异常模型,相比DTU13、Sandwell V23.1模型而言,其精度有显著的提高,如表 15所示.一方面说明,Cryosat-2的海面高观测数据精度较高,尤其是在地理环境复杂的区域.随着时间的推移,Cryosat-2对地观测数据的逐渐增多,基于Cryosat-2数据进行海洋大地测量相关的研究具有极大的潜力.另一方面说明,本文设计的技术思路、技术方法在大范围内使用是合理可行的.
由表 15可见,本文构建的我国近海海域重力异常模型,相比DTU13重力异常模型,平均值改进区间在-11.31~12.98 mgal之间,标准差改进区间在-3.39~0.07 mgal之间.相比Sandwell V23.1重力异常模型,平均值改进区间在-10.11~13.58 mgal之间,标准差改进区间在-1.09~0.71 mgal之间.
4 结论与建议从Cryosat-2的GDR数据中,提取海面高数据确定海面高梯度及其剩余梯度,并将剩余梯度采用配置法格网化为垂线偏差剩余分量,进而依据逆VeningMeinesz公式采用快速傅里叶变换的方法计算了剩余重力异常,恢复模型重力异常后,将其与船载重力测量数据进行融合,得到了我国近海海域的重力异常.文中以南海部分海域为例,分析了参考重力场模型的选择、对海面高梯度进行低通滤波时的最优滤波半径的选择、方差与协方差函数构建时的重力场模型阶次选择等问题.所得重力异常结果优于DTU13与Sandwell V23.1模型结果.采用迭代的方法确定最优滤波半径,滤波半径建议选择45 km.高精度的高阶重力场模型作为参考场有助于提高海域重力异常的精度.新构建的重力异常模型,可以部分满足海洋地球物理勘探、海洋测量大地等方面科学研究与工程应用需求.
[1] |
VIGNUDELLI S, KOSTIANOY A G, CIPOLLINI P, et al. Coastal altimetry[M]. Berlin: Springer, 2011
|
[2] |
郭金运, 常晓涛, 孙佳龙, 等. 卫星雷达测高波形重定及应用[M]. 北京: 测绘出版社, 2013 GUO Jinyun, CHANG Xiaotao, SUN Jialong, et al. Waveform retracking of satellite radar altimeter and applications[M]. Beijing: Surveying and Mapping Press, 2013 |
[3] |
金涛勇, 李建成, 姜卫平, 等. 基于多源卫星测高数据的新一代全球平均海面高模型[J]. 测绘学报, 2011, 40(6): 723 JIN Taoyong, LI Jiancheng, JIANG Weiping, et al. The new generation of global mean sea surface height model based on multi-altimetric data[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(6): 723 |
[4] |
HWANG C, PARSONS B. Gravity anomalies derived from Seasat, Geosat, ERS-1 and Topex/Poseidon altimeter and ship gravity: a case study over the Reykjanes Ridge[J]. Geophysical Journal International, 2010, 122(2): 551 |
[5] |
HWANG C. Inverse Vening Meinesz formula and deflection-geoid formula: applications to the predictions of gravity and geoid over the South China Sea[J]. Journal of Geodesy, 1998, 72(5): 304 DOI:10.1007/s001900050169 |
[6] |
HWANG C, KAO E C, PARSONS B. Global derivation of marine gravity anomalies from Seasat, Geosat, ERS-1 and TOPEX/POSEIDON altimeter data[J]. Geophysical Journal International, 2015, 134(2): 449 |
[7] |
SANDWELL D T, MVLLER R D, SMITH W H F, et al. New global marine gravity model from CryoSat-2 and Jason-1 reveals buried tectonic structure[J]. Science, 2014, 346(6205): 65 DOI:10.1126/science.1258213 |
[8] |
张胜军, 李建成, 褚永海, 等. 基于Cryosat和Jason1GM数据的垂线偏差计算与分析[J]. 武汉大学学报(信息科学版), 2015, 40(8): 1012 ZHANG Shengjun, LI Jiancheng, CHU Yonghai, et al. Calculation and analysis of the deflection of vertical derived from cryosat and jason1 gm data[J]. Geomatics and Information Science of Wuhan University, 2015, 40(8): 1012 |
[9] |
张胜军, 金涛勇, 褚永海, 等. Cryosat-2数据的大地水准面分辨能力研究[J]. 武汉大学学报(信息科学版), 2016, 41(6): 759 ZHANG Shengjun, JIN Taoyong, CHU Yonghai, et al. Estimation of the resolution capability of the cryosat-2 altimeter[J]. Geomatics and Information Science of Wuhan University, 2016, 41(6): 759 |
[10] |
万剑华, 李家军, 刘善伟, 等. 基于Cryosat -2数据的南海2' × 2'重力异常计算与分析[J]. 中国石油大学学报(自然科学版), 2015, 39(3): 70 WAN Jianhua, LI Jiajun, LIU Shanwei, et al. Calculation and analysis of 2' × 2' gravity anomalies over the South China Sea based on Cryosat-2 satellite altimeter data[J]. Journal of China University of Petroleum (Natural Science), 2015, 39(3): 70 |
[11] |
SRIN, ESA and Mullard Space Science Laborator. Cryosat product handbook [R].London:University College London, 2012.
|
[12] |
柯宝贵, 章传银, 郭春喜, 等. 船载重力测量数据不同测区系统偏差纠正方法研究[J]. 武汉大学学报(信息科学版), 2015, 40(3): 417 KE Baogui, ZHANG Chuanyin, GUO Chunxi, et al. Research on the system error correction for shipborne gravimetric data of different region in china offshore[J]. Geomatics and Information Science of Wuhan University, 2015, 40(3): 417 |
[13] |
黄谟涛. 海洋重力测量半系统差检验、调整及精度计算[J]. 海洋通报, 1990, 9(4): 81 HUANG Motao. Examniation, adjustment and precision estimation of half-systematic error in marine gravity surveying[J]. Marine Science Bulletin, 1990, 9(4): 81 |
[14] |
黄谟涛. 海洋重力测线网平差[J]. 测绘学报, 1993, 22(2): 103 HUANG Motao. Marine gravity survey line network adjustment[J]. Acta Geodaetca et Cartographica Sinica, 1993, 22(2): 103 |
[15] |
刘雁春, 李明叁, 黄谟涛. 海洋测线网系统误差调整的秩亏网平差模型[J]. 武汉大学学报(信息科学版), 2012, 26(6): 533 LIU Yanchun, LI Mingsan, HUANG Motao. The rank_defect adjustment model for survey_line systematic errors in marine survey net[J]. Geomatics and Information Science of Wuhan University, 2012, 26(6): 533 |
[16] |
李明叁, 刘雁春, 黄谟涛, 等. 海洋测线网系统误差确定的3种模型[J]. 测绘学院学报, 2002, 19(3): 157 LI Mingsan, LIU Yanchun, HUANG Motao, et al. Three models for determination of survey-line systematic errors in marine survey net[J]. Journal of Institute of Surveying and Mapping, 2002, 19(3): 157 |
[17] |
ANDERSEN O B, KNUDSEN P, BERRY P. The DNSC08GRA global marine gravity field from double retracked satellite altimetry[J]. Journal of Geodesy, 2010, 84(3): 191 DOI:10.1007/s00190-009-0355-9 |
[18] |
管铮, 黄谟涛, 翟国君, 等. 局部重力场逼近理论和方法[M]. 北京: 测绘出版社, 1997 GUAN Zheng, HUANG Motao, ZHAI Guojun, et al. Local gravity field approximation theory and method[M]. Beijing: Surveying and Mapping Press, 1997 |
[19] |
TSCHERNING C C, RAPP R H. Closed covariance expressions for gravity anomalies, geoid undulations and deflections of the vertical implied by anomaly degree variance models[R]. Columbus: The Ohio State University, 1974. https://link.springer.com/chapter/10.1007/978-3-642-72245-5_24
|