02_地球物理测点数据预处理

内容摘要： 本节课针对地球物理数据预处理过程中的插值、去趋势、坐标转换等常用的数据处理步骤进行讲解。常见的地球物理数据一般有矢量和标量两种，地形、重力这些都是标量，而GPS和地磁一般是矢量格式。由于观测环境的限制，一般采样都是非规则的，可能沿着多条测线进行，也可能是空间上不规则的离散点。在进行计算或画图之前，通常需要进行网格化或插值处理。

1、Verde工具包

在Python领域，用于插值的函数和工具包很多，考虑到地球物理数据的一些特殊性，我们选择Verde这套开源工具。

该工具包的插值实现基于经典的格林函数理论，类似于机器学习的线性回归，在scipy等工具包上扩展了一些地球物理领域的功能。

2、数据生成类—测试棋盘格CheckBoard

用途：在地球物理正反演过程中，经常需要验证数据的分辨率，这时候需要生成模拟数据，反演参数和算法的解释能力，棋盘格方法是常用的手段之一。棋盘格异常是指一组规则相间的异常形态，如下图所示：

图1 生成检测板异常

测试代码：

import verde as vd
synth = vd.datasets.CheckerBoard()
# Generating a grid results in a xarray.Dataset
grid = synth.grid(shape=(150, 100))
print("\nData grid:\n", grid)
synth.predict((2000, -2500))
grid = synth.grid(shape=(150, 100)) #生成一个150×150的网格数据对象
table = synth.scatter(size=100) #生成100个随机点，数据类型为pandas.dataframe对象。
grid.scalars.plot.pcolormesh()
table.head()
table.plot.scatter(x='northing', y='easting', c= 'scalars')

3、数据重采样（Data Decimation）

通常地球物理数据观测沿着测线进行，在一条侧线上的采样点十分密集，而要插值的网格远大于采样点间距，这时候需要按照一定规则抽取多条测线上的数据，可以使用verde下的BlockReduce类实现。

实例化时需要指定reduction方法，网格大小spacing；输出坐标类型center_coordinates为可选参数，如果设置为True，重采样坐标为网格中心点。实例化后filter函数实现数据从采样，coordinates参数要求是两列坐标tuple，data是待采样分量。

图2 重采样后的数据对比

reducer = vd.BlockReduce(reduction=np.median, spacing=5 / 60)
print(reducer)
coordinates, bathymetry = reducer.filter(
    coordinates=(data.longitude, data.latitude), data=data.bathymetry_m)
plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.Mercator())
ax.set_title("Locations of decimated data")
# Plot the bathymetry data locations as black dots
plt.plot(*coordinates, ".k", markersize=1, transform=crs)
vd.datasets.setup_baja_bathymetry_map(ax)
plt.tight_layout()
plt.show()
//如果输出采样后的网格中心点坐标
reducer_center = vd.BlockReduce(
    reduction=np.median, spacing=5 / 60, center_coordinates=True)
coordinates_center, bathymetry = reducer_center.filter(
    coordinates=(data.longitude, data.latitude), data=data.bathymetry_m)

4、数据去趋势（Trend Estimation）

通常在地球物理异常数据中，需要对趋势性的线性变换信号进行场分离，以便获取与背景无关的异常信号，即剩余异常，再进行反演和相应的定量或定性解释。在Verde包中，Trend类可以实现，实例化后fit、predict和filter函数可以分别完成拟合、预测和合并计算。

图3 数据去趋势

trend = vd.Trend(degree=1).fit(coordinates, data.bouguer)
print(trend.coef_)
#通过求得的trend函数，可以在其它坐标点上进行预测，以便获得剩余异常信号
trend_values = trend.predict(coordinates)
residuals = data.bouguer - trend_values
#也可以使用filter方法，一次性完成拟合，预测和剩余计算
__, res_filter, __ = vd.Trend(degree=1).filter(coordinates, data.bouguer)
print(np.allclose(res_filter, residuals))

5、坐标投影转换和插值

在数据插值、反演和模型计算时，一般需要使用笛卡尔坐标系，经纬度坐标需要经过投影变为km或m的长度单位。这时候需要利用pyproj库实现投影变换。这时候使用Spline类提供的方法来实现，fit，grid的用法如下：

图4 投影变换

import pyproj
import matplotlib.pyplot as plt
import geoist.others.fetch_data as fetch_data
data = fetch_data.fetch_gra_data()
projection = pyproj.Proj(proj="merc", lat_ts=data.Latitude.mean())
proj_coords = projection(data.Longitude.values, data.Latitude.values)
print(proj_coords)
plt.figure()
plt.title("Projected coordinates of gravity anomlay")
plt.plot(proj_coords[0], proj_coords[1], ".k", markersize=0.5)
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.gca().set_aspect("equal")
plt.tight_layout()
plt.show()
spline = vd.Spline().fit(proj_coords, data.Bouguer)
grid = spline.grid(spacing=10e3, data_names=["bouguer"]) #注意这个bouguer作为名字
print(grid)
grid.bouguer.plot.pcolormesh()
#通过distance_mask可以去除距离测线过远的数据
grid1 = vd.distance_mask(proj_coords, maxdist=30e3, grid=grid)
grid1.bouguer.plot.pcolormesh()

但是，在笛卡尔坐标系下处理完数据后可能需要再以经纬度形式标现，这时候需要进行投影反变换。首先需要生成坐标区域region和网格配置grid_geo。

# Get the geographic bounding region of the data
region = vd.get_region((data.Longitude, data.Latitude))
spacing = 1/60
print("Data region in degrees:", region)
# Specify the region and spacing in degrees and a projection function
grid_geo = spline.grid(
    region=region,
    spacing=spacing,
    projection=projection,
    dims=["latitude", "longitude"],
    data_names=["bouguer"],)
print("Geographic grid:")
print(grid_geo)
grid_geo = vd.distance_mask(
    (data.Longitude, data.Latitude), maxdist=30e3, grid=grid_geo, projection=projection)
grid_geo.bouguer.plot.pcolormesh() #还可以使用contour/contour/imshow等可视化方法

以上这些操作，也可以一次性完成，通过Chain这个类就可以实现。

chain = vd.Chain([("reduce", vd.BlockReduce(np.median, spacing * 111e3)),
        ("trend", vd.Trend(degree=1)),
        ("spline", vd.Spline()),])
print(chain)
chain.fit(proj_coords, data.Bouguer)
residuals = data.Bouguer - chain.predict(proj_coords)
plt.figure()
plt.title("Histogram of fit residuals")
plt.hist(residuals, bins="auto", density=True)
plt.xlabel("residuals (m)")
plt.xlim(-1500, 1500)
plt.tight_layout()
plt.show()
#与spline的grid类似，chain实例中可以直接调用grid，结果包含趋势项和预测分量，去除噪声。
grid = chain.grid(region=region,
    spacing=spacing,
    projection=projection,
    dims=["latitude", "longitude"],
    data_names=["bouguer"],)
print(grid)
print(chain.named_steps["trend"]) #也可以通过named_steps来访问grid方法。#"spline"

另外，该工具包还可以对矢量进行插值，设置模型光滑参数，根据数据方差定权等功能。对性能评估方面，包含了分组测试vd.train_test_split、评分spline.score、交叉验证vd.cross_val_score等，请自行学习和体会用法。

动手实践：进一步熟悉Verde软件包功能，分别使用重力、地磁、地形和GPS示例数据进行实验，并检验不同插值参数和模型的效果，完成Markdown笔记，上传到Gitee的项目库中（使用Git推送方式完成）。