刘书豪 个人
0%

SAGEA卫星重力数据处理：数据类型SHC和GRID介绍

发表于 分类于

GRACE(-FO)二级重力产品以球谐系数形式给出了地球重力场异常在不同频段的分布情况，其后续使用以及科学分析等仍然需要进一步处理（如物理量分离、滤波去噪等）。SAGEA (Satellite Gravity Error Assessment) 项目基于Python，提供了相应的数据下载、数据处理和误差分析等功能，项目开源地址为https://github.com/NCSGgroup/SaGEA 。本系列作为使用手册，希望可以帮助你更好地使用该工具，如果你在使用过程中有任何问题，欢迎与我讨论（liushuhao@hust.edu.cn）。本文将介绍SAGEA中的数据类型SHC（Spherical Harmonic Coefficient）和GRD (Gridded data)，并对它们的基本用法作简单介绍。

按照项目中README.md的指引安装后，路径pysrc/auxiliary/data_class/下的源文件SHC.py和GRD.py将提供对该工具包核心数据类型的定义。本篇文章主要介绍它们的结构和基本用法。

球谐系数SHC

属性和生成

SHC被实例化后，其核心属性SHC().value以一个numpy.ndarray二维数组给出了一组或多组球谐系数，其中第一个维度是数据的组数，第二个维度则是对该组球谐系数的有序排列，具体为

1
2
3
4
5
6
>>> shc = SHC().value
>>> shc

[[c1[0,0], s1[1,1], c1[1,0], c1[1,1], s1[2,2], s1[2,1], c1[2,0], c1[2,1], c1[2,2], s1[3,3] ...],
 [c2[0,0], s2[1,1], c2[1,0], c2[1,1], s2[2,2], s2[2,1], c2[2,0], c2[2,1], c2[2,2], s2[3,3] ...],
 [                                        ...                                         ]]

注意，即使被存储的数据仅含有一组时，其属性.value仍然为一个二维数组，此时的第一维长度为1。另外，如果存储的数据是多组的，那么这些组球谐系数的最大阶次一定是一致的。

一般来说，SHC不会直接拿来实例化，而是通过load_SHC等函数得到（具体请查看该系列相关文章）。不过仍然可以通过传入数据来得到其实例。具体操作为

1
shc = SHC(c, s=None)

有四种方式来传入参数：

  1. 参数c和s都是二维的数组，均以(l,m)的索引方式分别给出cos系数和sin系数的单组球谐系数，其中l是阶索引，m是次索引；
  2. 参数c和s都是二维的数组，均以(q,l,m)的索引方式分别给出cos系数和sin系数的多组球谐系数，其中q是组数索引，l是阶索引，m是次索引；
  3. 参数s为None, c是表示一组完整球谐系数的一维数组，以[c[0,0], s[1,1], c[1,0], c[1,1], s[2,2], s[2,1], ...]的顺序给出；
  4. 参数s为None, c是表示多组完整球谐系数的二维数组，以[[c1[0,0], s1[1,1], c1[1,0], ...], [c2[0,0], s2[1,1], c2[1,0], ...], ...]的顺序给出。

不论以哪种方式传入，shc.value的结构与按照前面所描述的一致。

基本信息获取

如前所述，SHC().value是该类型以二维数组的形式给出，其第一维代表组数，可以用方法SHC().get_length() -> int来得到组数个数，也可以用SHC().is_series() -> bool来判断其是否为多组数据。SHC().value的第二维代表该组的球谐系数，因此，其长度一定是\((\text{lmax}+1)^2\)，而不能是任意整数。其中\(\text{lmax}\)是整数，代表最大的阶次数，可以用方法SHC().get_lmax() -> int来得到最大的阶次数。

用法如下：

1
2
3
4
5
6
7
8
9
10
11
>>> sh_array = np.array([
    [c1[0,0], s1[1,1], c1[1,0], c1[1,1], s1[2,2], s1[2,1], c1[2,0], c1[2,1], c1[2,2], s1[3,3] ..., c1[96,96]],
	[c2[0,0], s2[1,1], c2[1,0], c2[1,1], s2[2,2], s2[2,1], c2[2,0], c2[2,1], c2[2,2], s2[3,3] ..., c2[96,96]],
])
>>> shc = SHC(sh_array)
>>> shc.get_length()
2
>>> shc.is_series()
True
>>>shc.get_lmax()
96

有时用户需要得到以二维矩阵形式给出的c, s系数，可以用方法SHC().get_CS2d()实现。该方法返回的是一个二维元组，分别是系数c构成的numpy.ndarray三维数组和s构成的numpy.ndarray三维数组。其中第一维度是组数，第二维度是阶索引l，第三维度是次索引m。注意即使SHC().is_series() is False，方法SHC().get_CS2d()返回的c和s仍然是三维的，此时第一维长度是1。具体如下：

1
2
3
4
5
6
7
8
>>> sh_array = np.array([
    [c1[0,0], s1[1,1], c1[1,0], c1[1,1], s1[2,2], s1[2,1], c1[2,0], c1[2,1], c1[2,2], s1[3,3] ..., c1[96,96]],
	[c2[0,0], s2[1,1], c2[1,0], c2[1,1], s2[2,2], s2[2,1], c2[2,0], c2[2,1], c2[2,2], s2[3,3] ..., c2[96,96]],
])
>>> shc = SHC(sh_array)
>>> c, s = SHC().get_CS2d()
>>> c.shape, s.shape
(1, 97, 97), (1, 97, 97)

球谐系数的加减

SHC添加了魔法方法__add__(self, other), __sub__(self, other)以实现两个SHC实例的加减运算，实现的本质是二者属性value的的加减。因此，它们的运算规则应该与相应的数组加减规则一致。如这句代码

1
>>> shc_1 + shc_2

返回的结果是一个新的SHC实例，其属性value应该等于shc_1.value + shc_2.value。并且计算过程中，应该满足下列条件之一：

  1. shc_2.value.shape == shc_1.value.shape，即二者拥有同样的组数和最大阶次。此时分别计算每组对应系数的相加；
  2. shc_2.is_series is False，即shc_2仅代表一组球谐系数。此时shc_1中的每一组系数分别与shc_2代表的这一组系数相加。

这些规则同样适用于相减的情况。

实际应用中，有时候可能只需要从某一阶开始或结束加减，此时可以使用方法SHC().add(other, lbegin, lend)和SHC().value(other, lbegin, lend)。其中参数other是SHC的实例，lbegin和lend均以int类型给出，分别代表从哪一阶开始（含）加减，从哪一阶结束（含）加减。返回的结果是自身按照相应规则加减被传入参数other后的SHC实例。具体如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
>>> shc1 = SHC(np.array([
    [c1[0,0], s1[1,1], c1[1,0], c1[1,1], s1[2,2], s1[2,1], c1[2,0], c1[2,1], c1[2,2], s1[3,3] ..., c1[61,61]],
]))

>>> shc2 = SHC(np.array([
    [c2[0,0], s2[1,1], c2[1,0], c2[1,1], s2[2,2], s2[2,1], c2[2,0], c2[2,1], c2[2,2], s2[3,3] ..., c2[61,61]],
]))

>>> shc1.add(shc2, lbegin=0, lend=1)
>>> shc1.value
np.array([
    [c1[0,0]+c2[0,0], s1[1,1] + s2[1,1], c1[1,0] + c2[1,0], c1[1,1] + c2[1,1], s1[2,2], s1[2,1], ...],
])

基于加减方法，SHC()还定义了.de_background(background: SHC=None)方法，即扣除背景场。参数background以SHC实例的形式传入，或者为None。如果。如果参数background是SHC实例，需要确保background.is_series() if False，即传入的背景场必须是单组球谐系数，此时返回的结果等价于自身每一组球谐系数减去背景场后构成的SHC实例。如果back_ground is None，则扣除自身所有组球谐系数的平均值。具体如下

1
2
3
4
5
>>> shc.de_background(shc_bg).value == (shc - shc_bg).value
True

>>> shc.de_background().value == shc.value - numpy.mean(shc.value, axis=0)
True

物理量转换和阶方差RSS、RMS的获取

在评估重力场时，我们可能需要通过将无量纲的球谐系数转换为大地水准面高或等效水高后，通过阶方差信息来判断。

一般来说，重力场球谐产品以无量纲球谐系数的形式描述地球重力场在不同频段上的信息，而SHC().convert_type(from_type, to_type)可以将其两个转换为我们所需要的物理量。其中，参数from_type和to_type的类型应该是物理量纲类型PhysicalDimensions。

SAGEA为一些固定的名称，如物理量纲、滤波器方法等等提供了对应的枚举类型，在../pysrc/auxiliary/preference/EnumClasses.py文件中。PhysicalDimensions定义了这些值：

1
2
3
4
5
6
7
8
9
10
class PhysicalDimensions(Enum):
    Dimensionless = 0  # 无量纲系数
    EWH = 1  # 等效水高 [m]
    Pressure = 2  # 压力 [bar]
    Density = 3  # 质量密度 [kg/m^3]
    Geoid = 4  # 大地水准面高 [m]
    Gravity = 5  # 重力 [mGal]
    HorizontalDisplacementEast = 6  # [m]
    HorizontalDisplacementNorth = 7  # [m]
    VerticalDisplacement = 8  # [m]

顾名思义，from_type定义自身目前代表的物理量类型，to_type定义了转换后的物理量类型。它们缺省值均为None，并且当传入None时，程序将为其赋值为PhysicalDimensions.Dimensionless。使用中，如果一个SHC实例在调用.convert_type()时会改变自身的值.value。

而阶方差的形式在不同情景下会使用RMS (Root Mean Squre)或RSS (Root Sum Squre)，即 \[ \text{RMS}(l) = \left(\cfrac{1}{2l+1} \left( C_{l,0}^2 + \sum_{m=1}^l \left( C_{l,m}^2 + S_{l,m}^2 \right) \right) \right) ^ {\cfrac{1}{2}} , \] 或 \[ \text{RSS}(l) = \left( C_{l,0}^2 + \sum_{m=1}^l \left( C_{l,m}^2 + S_{l,m}^2 \right) \right) ^ {\cfrac{1}{2}} . \] 对应的SHC()方法分别是.get_degree_rms()和.get_degree_rss()。对于.value.shape == (n, (lmax+1)^2)情况，这两种方法返回的结果形状为(n, lmax + 1)，其中n是组数, lmax是最大阶数。具体使用如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
import pysrc.auxiliary.preference.EnumClasses as Enums
import ...


shc = load_SHC(...)  # 读取文件内容详见本系列对应文章

shc.de_background()  # 扣除平均场

shc.convert_type(from_type=Enums.PhysicalDimensions.Dimensionless， to_type=Enums.PhysicalDimensions.Geoid)  # 将无量纲系数转换为大地水准面高

degree_array = np.arange(0, shc.get_lmax() + 1, 1)  # 得到阶序列，以便后续画图

rms = shc.get_degree_rms()  # 获取RMS
rss = shc.get_degree_rss()  # 获取RSS

plot_index = 10  # 随便定义一组来画图，不过当然不能超过组数

fig = plt.figure(figsize=(5, 6))
ax = fig.add_axes([0.2, 0.1, 0.7, 0.8])

ax.plot(degree_array[2:], rms[10][2:], label="RMS")  # GRACE二级重力场产品不提供0阶和1阶项信息，需要额外技术得到。所以从第2阶开始作图。
ax.plot(degree_array[2:], rss[10][2:], label="RSS")

ax.legend()

ax.set_xlabel('Degree')
ax.set_ylabel('Geoid [m]')

ax.set_yscale("log")

plt.show()

上述代码得到的图片为

格网系数GRD

作为对球面物理信号的表示，球谐系数（SHC）提供了信号在频谱中的信息，即不同波长下信号的强度。而格网产品作为以地理坐标为基底，则更加直观地表示信号在不同地理位置的强度。SAGEA则提供了另一种数据类型GRD类，来存储（多组的）这类数据。

属性和生成

GRD有三个核心属性：.value， .lat和.lon。它们分别是三维的、一维的和一维的np.ndarray类型，分别代表不同组数的格网值、地理纬度（\(^\circ\)）和地理经度（\(^\circ\)），并且它们的形状分别是(nset, nlat, nlon)，(nlat,)和(nlon,)。注意，GRD类型数据只能存储规则格网（即固定经纬度间隔）下的全球数据，即nlon始终为nlat的两倍。

本系列文章在未特别说明情况下，地理经度范围均在\(-180^\circ\)到\(180^\circ\)范围内定义。

一般来说，GRD很少被用户主动生成，不过生成时需要三个必要参数和一个缺省参数，即grid = GRID(value, lat, lon, option=1)。前三个分别对应其属性.value， .lat和.lon，而第四个则代表传入经纬度的类型。如果option为非0的任何值则代表传入的lat，lon是以\(^\circ\)为单位的地理经纬度，而如果option为0这代表传入的lat，lon是以\(\text{rad}\)为单位的地心经纬度。当然option只是用来指明传入参数的参考系信息，并不会改变GRD.lat和GRD.lat属性的参考系和单位。

基本信息获取

类似于SHC数据，GRD也提供一些函数来查看其基本信息。如GRD().get_length()获取组数，GRD().get_grid_space()获取格网密度等。在此不过多叙述。

球谐系数和格网系数的相互转换

SAGEA提供全球的球谐综合和球谐分析程序，以实现球谐系数和格网数据之间的相互转换。在包装为SHC方法和GRD方法后，使用方法为

1
2
3
4
5
6
7
8
9
import ...

shc = load_SHC(...)

grid_space = 1  # 定义格网密度
grid = shc.to_GRD(grid_space = grid_space, special_type: Enums.PhysicalDimensions = None)

lmax = 60  # 定义最大阶次数
shc_new = grid.to_SHC(lmax = lmax)

使用过程中，参数grid_space和lmax分别定义了要得到GRD实例的格网间隔和要得到SHC实例的最大阶次数，并且它们都有缺省值None。如果传入shc.to_GRD(grid_space = None)，则默认得到的格网间隔为int(180 / self.get_lmax())；而如果传入grid.to_SHC(lmax = None)，则默认得到的最大阶次数为lmax = int(180 / self.get_grid_space())。

另外，在通过球谐系数得到水平位移的格网数据Enums.PhysicalDimensions.HorizontalDisplacementEast或Enums.PhysicalDimensions.HorizontalDisplacementNorth，需要以其作为参数传入special_type。注意：这里的special_type仅提供指定球谐综合方式的作用，并不转换球谐系数本身，转换球谐系数仍需要使用SHC().convert_type()方法。