内容摘要：书接上文，在了解了网格数据的基本格式后，今天我们来聊聊位场数据处理的基本方法。主要内容包括：如何计算正常重力，离散点插值成网格，去趋势和平滑，异常体正演，生成模拟数据等。本节课开始，将逐步围绕重磁位场数据的处理和解释展开，因此，需要一定的地球物理基础知识，如果哪里讲的不清楚，欢迎下方留言。
1、计算重力异常
应用重力数据研究地壳结构，通常需要进行一系列校正，才能得到重力异常。这些校正一般分四种：正常场校正、高度校正、中间层校正、地形校正。前两种称之为自由空气校正、全算上叫布格校正。
可能有人会问，为什么要做这么多校正呢？因为，重力随着场源与观测点之间的不同差别很大。地球实际上是一个椭球体，随着纬度不同，重力值变化很快（图1）。还有就是地球表面有山、有谷，高度不同测量得到的重力值也不一样。
图1 地球的正常重力值
要想看到由于地壳结构差异引起的重力异常，首先必须把这些与纬度、高度相关的重力变化都校正掉。这种定义在标准椭球体上的理论重力场，我们称之为正常重力场。对所有的数据，采用统一个标准进行校正，就可以看到重力异常啦。
图2 重力异常的定义
一般来说，计算好的自由空气异常样子看上去与地形高度相关；而布格异常与地形不相关，在山区有时候会因为存在低密度的山根，而出现与地形反相关的特征。
图3 剩余密度与重力异常
简单说了计算重力异常的原理后，我们开始动手吧，看看Geoist中，有那些函数可以用于重力数据校正，请看下面代码：
from geoist.pfm import normgra
lat = 45.0
height = 1000.0
grav = 976783.05562266 #mGal
gamma = normgra.gamma_closed_form(lat, height)
gamma1 = normgra.gamma_closed_form(lat, height, ellipsoid = GRS80)
somi = normgra.gamma_somigliana_free_air(lat, height)
somi1 = normgra.gamma_somigliana_free_air(lat, height, ellipsoid = GRS80)
bug = normgra.bouguer_plate(height)
print(gamma, gamma1)
print(somi, somi1)
print(bug)
运行后，三个print的结果如下：
980311.2896926827 980311.4329622437
980311.1769377294 980311.3202522187
111.94694713134103
计算正常场的代码在normgra模块中实现，支持两种计算正常重力值的公式，接口分别为gamma_closed_form和gamma_somigliana_free_air，需要输入纬度和高度两个参数，默认是WGS84椭球，还支持GRS80椭球参数。
另外，还有函数bouguer_plate可以计算布格板（中间层）重力异常（2pirouGh），输入参数为观测高度。
输入重力值单位为mGal，高度单位为m，经纬度为dd.ddd格式。布格异常计算公式如下：
bouguer = grav -gamma - bug
-3640.181017154048
2、插值到网格
实际陆地重力数据处理过程，首先是对每个测量点进行重力异常计算。计算好每个点上的重力异常后，要看重力异常场的特征或者进一步进行空间导数和延拓等变换。都需要基于网格化的数据来进行，这时候最常用的就是网格化。
如果您使用的是用球谐模型来描述的重力异常场，网格化就可以省略了，直接用球谐函数解算规则格点上的重力异常就可以形成网格数据。
网格化方法有很多，最小曲率法、克里金法、最小二乘配置法、反距离加权平均法等等，不同的方法优缺点我们这里不想展开讨论。但无论什么方法，核心目的其实都一样，就是尽量避免没有观测数据的地方出现虚假异常。因为一旦出现具有特征的虚假异常，无论是求导还是延拓都将针对不真实的数据进行，也会对将来的解释带来影响。
一般而言，如果评价一个插值方法是否合适，首先应该看插值结果是否产生出高频特征的畸变类异常。如果对于未覆盖观测的区域插值结果呈现出的是平缓的、低频的异常特征，往往是比较合理的。在整个频率成份上看，不同空间位置上的异常不会产生与测点疏密相关的高低频特征，这样的插值算法就比较好。
下面我们看看怎么实现插值吧！
2.1 投影变换
一般在插值前，需要将经纬度坐标，转换为投影坐标，那就需要先做个投影变换，结合前面介绍的网格投影知识，我们还是使用pyproj库。离散点投影变换相比网格投影变换要简单，没有重采样过程，点和点之间都是一一对应关系，记住Proj和transform两个函数即可。代码如下：
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pathlib import Path
from geoist.pfm import normgra
from geoist import DATA_PATH
# 1.读取数据
datapath = Path(Path(normgra.__file__).parent, 'data')
filename = Path(datapath, 'ynyx_grav.csv') #'ynyx_grav.csv') #ynp1_grav.csv
gradata = pd.read_csv(filename)
print('1. 重力数据已经读入，格式为：{}'.format(gradata.keys()))
# 2. 计算FGA
gradata['freeair'] = normgra.gamma_closed_form(gradata['lat'], gradata['elev'])
gradata['buglayer'] = normgra.bouguer_plate(gradata['elev'])
gradata['FGA'] = gradata['grav'] - gradata['freeair']
gradata['BGA_s'] = gradata['grav'] - gradata['freeair'] - gradata['buglayer']
# 2.1 正常场计算方法不同
#gradata['BGA_s1'] = gradata['grav'] - normgra.gamma_somigliana_free_air(gradata['lat'], gradata['elev']) - gradata['buglayer']
# 2.2 输出结果
gradata.to_csv(Path(DATA_PATH, 'ynp1_grav_anomaly.csv'), index = False, sep = ',')
print('2. 重力异常计算完成，已保存到：{}'.format(Path(DATA_PATH, 'ynp1_grav_anomaly.csv')))
import pyproj
p_jw = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
p_lcc = "+proj=lcc +lon_0=102.5 +lat_0=24.38 +lat_1=45 +ellps=WGS84 +datum=WGS84 +no_defs"
proj_xy = pyproj.Proj(p_lcc) #projection = pyproj.Proj(proj="merc", lat_ts=gradata['lat'].mean())
proj_coords = proj_xy(gradata['lon'].values, gradata['lat'].values)
gradata['x'] = proj_coords[0]
gradata['y'] = proj_coords[1]
proj_jw = pyproj.Proj(p_jw) # 目标坐标系统
origin_lon, origin_lat = pyproj.transform(proj_xy, proj_jw, gradata['x'].values, gradata['y'].values)
fig, (ax0, ax1) = plt.subplots(nrows=1, ncols=2,figsize=(12,6))
ax0.set_title("Locations of gravity anomlay")
ax0.plot(gradata['lon'], gradata['lat'], "ok", markersize=1.5)
ax0.set_xlabel("Longitude")
ax0.set_ylabel("Latitude")
ax0.set_aspect("equal")
ax1.set_title("Projected coordinates of gravity anomlay")
ax1.plot(gradata['x'], gradata['y'], "ok", markersize=1.5)
ax1.set_xlabel("Easting (m)")
ax1.set_ylabel("Northing (m)")
ax1.set_aspect("equal")
plt.tight_layout()
plt.show()
前面的代码是一段导入数据代码，直接从第3部分开始看即可，可视化效果如下：
图4 离散点投影变换
2.2 网格化
离散点网格化，需要使用到geoist的gridder模块，里面有两个类需要导入，分别为：spline，mask。直接看代码吧！
from geoist.gridder import spline, mask
BGA = spline.Spline().fit(proj_coords, gradata['BGA_s'].values)
res_BGA = gradata['BGA_s'].values - BGA.predict(proj_coords)
grid = BGA.grid(spacing=2e2, data_names=["BGAs"])
print('4. 网格化信息如下：',grid)
type(grid)
grid1 = mask.distance_mask(proj_coords, maxdist=5e2, grid=grid)
grid1.BGAs.plot.pcolormesh()
上面代码中，BGA是spline的类实例，在初始化的时候可以将坐标和需要网格化的点值输入。实例化后BGA支持predict函数，可以对每个点上的插值结果进行校验。grid函数可以生成网格对象，类型是xarray的dataset，如下：
4. 网格化信息如下： <xarray.Dataset>
想看一下数据长什么样，直接可以xarray实例的plot.pcolormesh()函数。细心的朋友们注意到还有一个mask，这个其实就是让距离大于500m的点不显示，避免数据太稀疏产生的插值错误信息。如图5所示：
图5 网格化后的布格重力异常
3、去趋势
在不同比例尺的重力异常图上，经常可以看到背景性的趋势变化特征。这种趋势性特征，根据研究目的的不同，通常要去除掉。
一般在重力位场解释中认为，深部场源体的异常特征以低频为主，而浅部的多呈现高频特征。这些特征都是定性的，没有严格上的定量关系。比如：一张百km尺度的重力异常图/曲线，整个界面上的趋势异常特征，可以解释为与莫霍面的起伏相关，因为地壳和地幔之间的密度差异最大。
如果您要通过重力资料研究浅部的构造，比如断裂的倾角、走向，这种趋势性的异常往往首先要考虑去掉。这就是去趋势，或场源分离，反之，如果要分离浅部的高频干扰，就是平滑。
下面我们就看看在怎么实现这些处理步骤，在上面网格化的基础上，我们再导入trend类，使用方法与上面的网格化十分类似，都是在类的初始化时设置好参数，开始拟合，代码如下：
from geoist.gridder import trend
trend = trend.Trend(degree = 1).fit(proj_coords, gradata['BGA_s'].values)
print('4. 拟合的趋势系数：'.format(trend.coef_))
trend_values = trend.predict(proj_coords)
residuals = gradata['BGA_s'].values - trend_values
resBGA = spline.Spline().fit(proj_coords, residuals)
grid2 = resBGA.grid(spacing=2e2, data_names=["resBGA"])
grid3 = mask.distance_mask(proj_coords, maxdist=5e2, grid=grid2)
plt.figure()
grid3.resBGA.plot.pcolormesh()
大家注意residuals = gradata['BGA_s'].values - trend_values，这段代码，这个计算得到的residuals数据就是剩余布格重力异常，这个异常通常可以直接用于资料解释。
小结：今天讲的内容都是位场数据处理是基础，下一节我们将继续讨论如何进行位场的空间解析延拓，导数计算等。掌握这些用法，对于基本的重力异常分析和定性解释就可以轻松实现，完成一般工程项目报告中的数据处理和图件绘制都是小case啦！
网友评论