留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

采用点质量模型方法反演中国大陆及周边地区陆地水储量变化

郭飞霄 苗岳旺 肖云 汪菲菲

郭飞霄, 苗岳旺, 肖云, 汪菲菲. 采用点质量模型方法反演中国大陆及周边地区陆地水储量变化[J]. 武汉大学学报 ● 信息科学版, 2017, 42(7): 1002-1007. doi: 10.13203/j.whugis20150031
引用本文: 郭飞霄, 苗岳旺, 肖云, 汪菲菲. 采用点质量模型方法反演中国大陆及周边地区陆地水储量变化[J]. 武汉大学学报 ● 信息科学版, 2017, 42(7): 1002-1007. doi: 10.13203/j.whugis20150031
GUO Feixiao, MIAO Yuewang, XIAO Yun, WANG Feifei. Recovery Water Storage Variation in China and Its Adjacent Area by Method of Point-Mass Model[J]. Geomatics and Information Science of Wuhan University, 2017, 42(7): 1002-1007. doi: 10.13203/j.whugis20150031
Citation: GUO Feixiao, MIAO Yuewang, XIAO Yun, WANG Feifei. Recovery Water Storage Variation in China and Its Adjacent Area by Method of Point-Mass Model[J]. Geomatics and Information Science of Wuhan University, 2017, 42(7): 1002-1007. doi: 10.13203/j.whugis20150031

采用点质量模型方法反演中国大陆及周边地区陆地水储量变化

doi: 10.13203/j.whugis20150031
基金项目: 

国家973计划 2013CB733303

国家自然科学基金 41374083

大地测量与地球动力学国家重点实验室开放研究基金 SKLGED2013-3-3-E

详细信息
    作者简介:

    郭飞霄, 硕士, 主要从事卫星重力数据处理研究。gravity_gfx@126.com

  • 中图分类号: P223;P228.41

Recovery Water Storage Variation in China and Its Adjacent Area by Method of Point-Mass Model

Funds: 

The National Basic Research Program(973 Programs) of China 2013CB733303

the National Natural Science Foundation of China 41374083

the State Key Laboratory of Geodesy and Earth's Dynamics SKLGED2013-3-3-E

More Information
    Author Bio:

    Guo Feixiao, master, specializes in satellite gravimetry data processing.E-mail:gravity_gfx@126.com

图(4) / 表(1)
计量
  • 文章访问数:  1419
  • HTML全文浏览量:  71
  • PDF下载量:  461
  • 被引次数: 0
出版历程
  • 收稿日期:  2015-05-21
  • 刊出日期:  2017-07-05

采用点质量模型方法反演中国大陆及周边地区陆地水储量变化

doi: 10.13203/j.whugis20150031
    基金项目:

    国家973计划 2013CB733303

    国家自然科学基金 41374083

    大地测量与地球动力学国家重点实验室开放研究基金 SKLGED2013-3-3-E

    作者简介:

    郭飞霄, 硕士, 主要从事卫星重力数据处理研究。gravity_gfx@126.com

  • 中图分类号: P223;P228.41

摘要: 研究了反演区域陆地水储量变化的点质量模型方法,采用Tikhonov正则化方法解决了反演过程中参数估计病态问题。利用GRACE(gravity recovery and climate experiment)时变重力场模型数据,用点质量模型方法反演了中国大陆及其周边地区陆地水储量变化,将点质量模型反演结果与球谐系数法反演结果、GLDAS(global land data assimilation)水文模型数据进行了验证分析,并选取了4个特征点计算了陆地水储量变化时间序列。实验结果表明,由于点质量模型方法将研究区域内不同网格质量变化对地球重力场的影响分离开来,所得区域陆地水储量变化局部信号更明显,并且点质量模型方法反演结果与GLDAS水文模型数据相关性更强。

English Abstract

郭飞霄, 苗岳旺, 肖云, 汪菲菲. 采用点质量模型方法反演中国大陆及周边地区陆地水储量变化[J]. 武汉大学学报 ● 信息科学版, 2017, 42(7): 1002-1007. doi: 10.13203/j.whugis20150031
引用本文: 郭飞霄, 苗岳旺, 肖云, 汪菲菲. 采用点质量模型方法反演中国大陆及周边地区陆地水储量变化[J]. 武汉大学学报 ● 信息科学版, 2017, 42(7): 1002-1007. doi: 10.13203/j.whugis20150031
GUO Feixiao, MIAO Yuewang, XIAO Yun, WANG Feifei. Recovery Water Storage Variation in China and Its Adjacent Area by Method of Point-Mass Model[J]. Geomatics and Information Science of Wuhan University, 2017, 42(7): 1002-1007. doi: 10.13203/j.whugis20150031
Citation: GUO Feixiao, MIAO Yuewang, XIAO Yun, WANG Feifei. Recovery Water Storage Variation in China and Its Adjacent Area by Method of Point-Mass Model[J]. Geomatics and Information Science of Wuhan University, 2017, 42(7): 1002-1007. doi: 10.13203/j.whugis20150031
  • 地球系统的质量及其分布是随时间不断变化的,质量的重新分布将导致地球重力场的变化。在一定的时间尺度上,地球重力场的变化主要来自大气压力、海底压力和陆地水储量的变化[1-2]。GRACE卫星的发射开创了高精度全球重力场观测与气候变化试验的新纪元,为获取全球尺度的物质质量迁移提供了一种有效途径。Wahr等建立了利用时变重力场模型估计地球表层质量变化的球谐系数算法模型[3],目前已广泛应用于陆地水储量变化、地震同震、海洋环流和冰川融化等领域研究[4-7]。由于信号泄露、混淆效应以及时变重力场模型误差等因素影响,采用球谐系数法只能反演较大空间尺度上地球表层质量随时间的变化。因此,提高地球表层质量变化反演结果的空间分辨率成为当前GRACE时变重力场研究的热点之一。

    地球重力场可看成是由正常重力场和扰动重力场两部分组成。假设地球质量异常都是聚集在地球表面上的一个个小质量异常块,那么扰动重力场就是这些质量异常块共同作用所产生的。可将扰动重力进行向下延拓,视质量异常块为点质量,进而求得区域质量变化,以期望提高反演结果的空间分辨率。Forsberg等建立了地球表层点质量变化引起空间点扰动重力变化的数学模型[8]。Baur利用点质量模型方法反演了格林兰岛地区冰盖变化[9],反演结果分辨率有所提高,反演结果优于球谐系数法结果,验证了方法的有效性。

    本文推导了利用点质量模型方法反演地球表层质量变化的数学模型,并将该方法应用于反演区域陆地水储量变化,采用GRACE卫星Level-2时变重力场模型数据计算了2010年中国大陆及其周边地区的陆地水储量变化,并与球谐系数法反演结果进行了对比分析。

    • 地球表层质量变化会引起空间点的扰动重力径向分量值的改变。而扰动重力径向分量值的改变可由时变重力场模型球谐系数的改变量来表示[10]:

      $$ \begin{array}{l} \delta {g_\rho } = - \frac{{{\rm{GM}}}}{{{r^2}}}\sum\limits_{l = 2}^{l\max } {\left( {l + 1} \right){{\left( {\frac{a}{r}} \right)}^l}} \sum\limits_{m = 0}^l {{{\bar P}_{lm}}\left( {\cos \varphi } \right) \times } \\ \left[ {\Delta {C_{lm}}\cos m\lambda + \Delta {S_{lm}}\sin m\lambda } \right] \end{array} $$ (1)

      式中,G为引力常数;M为地球质量;a为地球平均半径;lm分别为重力场模型的阶数和次数;(λ, φ)分别为地心经度和地心余纬;Plm(cosφ)为完全规格化的勒让德函数;r=a+h为空间点到地心的向径;h为空间点到地面上的高度;lmax是重力场球谐模型截断阶数,{ΔClm, ΔSlm}={ClmClmSlmSlm}是时变重力场模型{Clm, Slm}相对于参考重力场模型{Clm, Slm}的球谐系数改变量。

      扰动重力变化由两部分组成,包括地球表层物质变化直接引起的扰动重力变化和固体地球受到物质作用形变而引起的扰动重力变化。因此,式(1) 体现的是这两种变化总的影响之和,故地球重力场模型球谐系数改变量可写为两部分之和:

      $$ \left\{ \begin{array}{l} \Delta {C_{lm}} = \Delta C_{lm}^{{\rm{mass}}} + {k_1}\Delta C_{lm}^{{\rm{mass}}}\\ \Delta {S_{lm}} = \Delta S_{lm}^{{\rm{mass}}} + {k_l}\Delta S_{lm}^{{\rm{mass}}} \end{array} \right. $$ (2)

      式中,ΔClmmass是地球表层质量变化引起的重力场模型球谐系数改变量;kll阶负荷勒夫数。因此,仅由地球表层质量变化引起的扰动重力径向分量变化为:

      $$ \begin{array}{*{20}{c}} {\delta g_\rho ^{{\rm{mass}}} = - \frac{{{\rm{GM}}}}{{{r^2}}}\sum\limits_{l = 2}^{l\max } {\frac{{\left( {l + 1} \right)}}{{\left( {1 + {k_1}} \right)}}{{\left( {\frac{a}{r}} \right)}^l}} \sum\limits_{m = 0}^l {{{\bar P}_{lm}} \times } }\\ {\cos \varphi \times \left[ {\Delta {C_{lm}}\cos m\lambda + \Delta {S_{lm}}\sin m\lambda } \right]} \end{array} $$ (3)

      另一方面,如图 1所示,地球表面一点P的点质量变化δmj引起的空间一点S的扰动重力径向分量改变值δgρij为:

      图  1  点质量模型几何示意

      Figure 1.  Point-mass Modeling Geometry

      $$ {{\delta g}}_\rho ^{ij} = \frac{{G\delta {m_j}}}{{d_{ij}^2}}\cos \alpha $$ (4)

      式中,dij是空间点S和地球表面点P的距离。由图 1中几何关系, 可知:

      $$ \cos \alpha = \frac{{{r_i} - a\cos {\psi _{ij}}}}{{{d_{ij}}}} $$ (5)
      $$ \begin{array}{l} \cos {\psi _{ij}} = \sin {\varphi _i}\sin {\varphi _j} + \cos {\varphi _i} \times \\ \cos {\varphi _j}\cos \left( {{\lambda _i} - {\lambda _j}} \right) \end{array} $$ (6)
      $$ {d_{ij}} = \sqrt {{a^2} + r_i^2 - 2a \times {r_i}\cos {\psi _{ij}}} $$ (7)

      于是,式(4) 可写为:

      $$ { {\delta g}}_\rho ^{ij} = G\delta {m_j}\frac{{{r_i} - a\cos {\psi _{ij}}}}{{{{\left( {{a^2} + r_i^2 - 2a \times {r_i}\cos {\psi _{ij}}} \right)}^{3/2}}}} $$ (8)

      因此,地球表层质量变化引起的空间点S的扰动重力径向分量改变值为:

      $$ \begin{array}{l} {{\delta g}}_\rho ^i = G\sum\limits_{j = 1}^n {{\rm{\delta }}{m_j}} \times \\ \frac{{{r_i} - a\cos {\psi _{ij}}}}{{{{\left( {{a^2} + r_i^2 - 2a \times {r_i}\cos {\psi _{ij}}} \right)}^{3/2}}}} \end{array} $$ (9)

      式中,n是研究区域质量异常的点的个数;ri为空间点S到地心的向径。进一步,可将待求质量变化划为等效水柱高形式。设δmj=10×sj×Δhjλj1λj2为研究区域经度范围,φj1φj2为研究区域纬度范围,sj为待求研究区域面积,具体形式为:

      $$ {s_j} = \left( {a\Delta \varphi } \right) \times \left( {a\Delta \lambda } \right)\sin {\varphi _{j1}} $$ (10)

      于是,可将式(9) 所求区域质量变化换算成等效水柱高度Δhj形式:

      $$ \begin{array}{l} {{\delta g}}_\rho ^i = 10 \times G\sum\limits_{j = 1}^n {{s_j}} \times \\ \frac{{{r_i} - a\cos {\psi _{ij}}}}{{{{\left( {{a^2} + r_i^2 - 2a \times {r_i}\cos {\psi _{ij}}} \right)}^{3/2}}}} \times \Delta {h_j} \end{array} $$ (11)

      联立式(3) 和式(11) 即得采用点质模型反演区域质量变化的数学模型,i=1, 2…, k为空间离散点的个数。

    • 进一步,将式(11) 进行线性化,得:

      $$ \mathit{\boldsymbol{y = Ax + e}} $$ (12)

      式中,e为误差向量;y=[δgρi]k×1为观测向量可按式(3) 计算求得;x=[Δhj]n×1是待求量;A(k×n)为线性化后系数矩阵,可按式(11) 右端计算得到。

      研究表明, 采用最小二乘准则求解式(12) 是一个病态问题,其设计矩阵A近似秩亏,法方程组系数矩阵ATA将会奇异[9]。因此,为了求得式(12) 最优解,采用Tikhonov正则化原理[11-12],构造如下估计准则:

      $$ \mathit{\Omega = }{\left\| {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{Ax}}} \right\|^2} + \alpha {\left\| \mathit{\boldsymbol{x}} \right\|^2} = \min $$ (13)

      式中,α≥0是正则化参数;‖·‖表示欧氏2-范数。将式(13) 对待求参数x求偏导数,并令该偏导数等于0,即可得:

      $$ \mathit{\boldsymbol{\hat x = }}{\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A + }}\alpha \mathit{\boldsymbol{I}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{y}} $$ (14)

      与最小二乘解相比,式(14) 右端求逆部分增加了αI一项。Tikonov正则化方法解决这类病态问题的思想就是在法方程组系数矩阵的对角线增加这一阻尼项,抑制法方程组的病态性,以得到稳定可靠的参数解。从式(14) 可以看出,解算结果与正则化参数α相关,正则化参数的取值不同,所得估值可能不同。因此求解式(14) 的关键在于如何选择合适的正则化参数α。常见的方法有GCV方法、L-曲线法、岭迹法等。

      Hansen针对不适定问题提出的L-曲线法比较成熟和有影响的一种方法,其核心是定位L曲线上曲率最大的那个点[12]。在式(14) 中,‖yAx‖和‖x‖都是正则化参数的函数α的函数,以正则化参数α为变量,以μ=lg(‖yAx‖)作横轴,ν=lg(‖x‖)作纵轴,经过曲线拟合得到一条“L”形状的曲线,该曲线上曲率最大的那点所对应的α值就是L-曲线法所确定的正则化参数,该方法确定的正则化参数不是最优的,只是近似最优的。从式(14) 可以看出,应用L-曲线法选择正则化参数的合理性在于该方法强调数据拟合部分‖yAx‖和解部分‖x‖之间的平衡,这种平衡是通过正则化参数来实现的。

    • k个时变重力场模型改变量ΔClmk为:

      $$ \begin{array}{l} \Delta C_{lm}^k = C_{lm}^k - {{\bar C}_{lm}} = \\ C_{lm}^k - \left( {C_{lm}^l + \cdots + C_{lm}^k + \cdots + C_{lm}^n} \right)/n = \\ - \left[ {C_{lm}^l + \cdots + \left( {1 - n} \right)C_{lm}^k + \cdots + C_{lm}^n} \right]/n \end{array} $$ (15)

      根据误差传播规律,已知时变重力场模型球谐系数中误差,由式(15) 可得球谐系数改变量ΔClmk的中误差为:

      $$ \begin{array}{*{20}{c}} {{m_{\Delta C_{lm}^k}}^2 = }\\ {\frac{{m{{_{C_{lm}^1}^{}}^2} + \cdots + {{\left( {1 - n} \right)}^2}m{{_{C_{lm}^k}^{}}^2} + \cdots + m{{_{C_{lm}^n}^{}}^2}}}{{{n^2}}}} \end{array} $$ (16)

      同理可得:

      $$ \begin{array}{*{20}{c}} {m{{_{\Delta S_{lm}^k}^{}}^2} = }\\ {\frac{{m{{_{S_{lm}^1}^{}}^2} + \cdots + {{\left( {1 - n} \right)}^2}m{{_{S_{lm}^k}^{}}^2} + \cdots + m{{_{S_{lm}^n}^{}}^2}}}{{{n^2}}}} \end{array} $$ (17)

      按照误差传播定律,可由式(3) 可以得出点质量模型的扰动重力径向分量中误差:

      $$ \begin{array}{*{20}{c}} {m_{{g_\rho }}^{}2 = }\\ {{{\left( {\frac{{GM}}{{{r^2}}}} \right)}^2}\left[ {\sum\limits_{l = 2}^{l\max } {{f_l}} \sum\limits_{m = 0}^l {\left( {f_{lm}^cm{{_{\Delta C_{lm}^k}^{}}^2} + f_{lm}^sm{{_{\Delta S_{lm}^k}^{}}^2}} \right)} } \right]} \end{array} $$ (18)

      其中,

      $$ {f_l} = {\left( {\frac{{l + 1}}{{1 + {k_l}}}{{\left( {\frac{a}{r}} \right)}^l}} \right)^2} $$ (19)
      $$ f_{lm}^c = {\left[ {{{\bar P}_{lm}}\left( {\cos \varphi } \right) \times \cos m\lambda } \right]^2} $$ (20)
      $$ f_{lm}^s = {\left[ {{{\bar P}_{lm}}\left( {\cos \varphi } \right) \times \sin m\lambda } \right]^2} $$ (21)

      根据式(18) 可得到每个空间点的扰动重力误差。因此,可得观测向量的协方差阵Qy。然后根据式(14),按照误差传播定律可得估值的协方差阵为:

      $$ {\mathit{\boldsymbol{Q}}_{\hat x}} = \mathit{\boldsymbol{K}}{\mathit{\boldsymbol{Q}}_y}{\mathit{\boldsymbol{K}}^{\rm{T}}} $$ (22)

      式中, K=(ATA+αI)-1AT

      综上所述,可得单位权中误差μ

      $$ \mu = \pm \sqrt {\frac{{{\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{V}}}}{{k - n}}} $$ (23)

      式中,$ \mathit{\boldsymbol{V = A}}\hat x-\mathit{\boldsymbol{y}} $。

      需要指出的是,完整的误差估计还要包括泄露误差和GIA误差,后续工作将深入研究。

    • 实验中GRACE时变重力场模型采用德国地学中心发布的Level-2 Realease05版本数据,时间从2010-01~2010-12共12个时变重力场模型。考虑到GRACE时变重力场模型阶方差随着阶数增加而迅速增大,解算中对GRACE时变重力场模型截断至60阶。计算时首先对时变重力场求平均值作为稳态背景场,然后计算每个月时变重力场模型相对于稳态背景场的球谐系数改变量,再由每个月时变重力场模型球谐系数改变量根据§1所述点质量模型方法计算了中国大陆及其周边地区(10°~55°N,72°~136°E)的陆地水储量的月变化。

    • 根据点质量模型方法计算2010-01中国大陆及其周边地区陆地水储量变化。计算时,h取GRACE卫星轨道平均高度480 km,地面网格大小取1°×1°。由于线性化后的方程组系数矩阵呈病态,按照Tikhonov正则化原理,构造正则化准则,采用L-曲线方法求解正则化参数。图 2所示为采用点质量模型方法反演了2010-01中国大陆及周边地区陆地水储量变化。根据§1.3所述,得到反演结果单位权中误差为±0.01 cm。

      图  2  点质量模型反演2010-01中国大陆陆地水储量变化

      Figure 2.  2010-01 Water Storage Variation in China and Its Adjacent by Means of Point-Mass Model

      为验证反演结果,分别采用球谐系数法和GLDAS水文模型数据计算了相同区域范围2010-01陆地水储量变化,如图 3所示。其中,图 3(a)球谐系数法反演结果,计算所用数据及处理过程与点质量模型方法反演计算相同,计算时采用了滤波半径为400 km的扇形滤波[13]图 3(b)是GLDAS水文模型结果。对比图 2图 3(b)可以看出,采用点质量模型方法反演结果在整体上与GLDAS水文模型结果呈现出很好的相符性,GLDAS模型和点质量模型反演陆地水储量变化在滇藏交界地区、新疆西南部、中南以及辽宁地区呈减小趋势,在京津地区和内蒙古中部地区呈增长趋势。但是,利用点质量模型方法反演陆地水储量变化结果要大于GLDAS模型数据,主要可能原因在于GLDAS水文模型是以雷达和微波辐射计资料为主,主要反映了浅层地表土壤湿度及地表积雪厚度随时间的变化情况,体现了降雨和蒸发等季节信号特征[14]。而GRACE时变重力场则包括了总的陆地水储量变化和其他可能的地球物质迁移。

      图  3  2010-01中国大陆陆地水储量变化

      Figure 3.  2010-01 Water Storage Variation in China and Its Adjacent

      对比图 2图 3(a)可以看出:采用球谐系数法反演陆地水储量变化结果与点质量模型方程反演结果基本相一致,点质量模型方法得到的陆地水储量变化信号更明显。但两者同GLDAS水文模型结果相比,特别是在新疆西南部、内蒙古中部和辽宁地区,球谐系数法反演的陆地水储量变化信号不明显,而采用点质量模型方法则较好地反映出这些地区的水储量变化信号。这表明采用点质量模型方法,反演结果与GLDAS水文模型结果吻合程度更好。

    • 选取质量变化信号明显的4个点,其分布如图 2所述。分别采用点质量模型方法、球谐系数法和GLDAS模型数据计算了特征点从2010-01~2010-12的陆地水储量变化时间序列,如图 4所示。

      图  4  特征点质量变化时间序列

      Figure 4.  Water Storage Variation Time Series of Characteristic Points

      图 4为4个特征点2010年~12月的陆地水储量变化时间序列,其中红、蓝和黑色线条分别是点质量模型方法(p1)、球谐系数法(p2)和GLDAS水文模型数据(p3)结果,3种方法所得时间序列相关性如表 1所示。从图 4表 1可以看出:① 点质量模型方法反演结果和球谐系数法反演结果相一致,4个特征点的点质量模型方法反演结果和球谐系数法反演结果的相关系数均在0.90以上;② 点质量模型方法反演结果和GLDAS水文模型数据的相关性更强,4个特征点的点质量模型方法反演结果和GLDAS水文模型数据的相关系数分别为0.80、0.94、0.85和0.74,均大于球谐系数法反演结果和GLDAS水文模型数据的相关系数。

      表 1  不同反演结果相关系数统计

      Table 1.  Correlation Coefficient Statistics of Different Inversion Result

      特征点 cov(p1, p2) cov(p2, p3) cov(p1, p3)
      1 0.94 0.78 0.80
      2 0.95 0.85 0.94
      3 0.98 0.79 0.85
      4 0.90 0.70 0.74
    • 本文详细推导了空间点扰动重力径向分量改变量和点质量变化间的数学关系,并应用于中国大陆及其周边区域陆地水储量变化反演,将反演结果与球谐系数法和GLDAS水文模型数据进行比较分析。由于点质量模型方法将研究区域内不同网格质量变化对地球重力场的影响分离开来,理论上能够提高反演结果空间分辨率。由推导过程及反演结果可得出结论:点质量模型方法反演所得区域陆地水储量变化局部信号更明显,反演结果与GLDAS水文模型数据相关性更强,进一步表明了陆地表层质量变化主要是陆地水储量的变化。

参考文献 (14)

目录

    /

    返回文章
    返回