在直角坐标系中旋转检测版

最近遇到一个需求，需要在检测版测试中将检测版旋转一定角度，来避免射线方向的Smearing。

定义研究区域和网格点#

首先定义研究区域范围为x方向 $( 0 - 1000m )$ ，y方向 $( 0 - 800m )$ , 间隔为 $12m$ ，并生成网格点坐标为xgrid和ygrid

1
import numpy as np
2
import matplotlib.pyplot as plt
3
from scipy.interpolate import RegularGridInterpolator
4

5
xmin, xmax = 0, 1000
6
ymin, ymax = 0, 800
7
x0 = (xmin+xmax)/2
8
y0 = (ymin+ymax)/2
9
dx = 12
10
dy = 12
11
x = np.arange(xmin, xmax, dx)
12
y = np.arange(ymin, ymax, dy)
13
ygrid, xgrid = np.meshgrid(y, x, indexing='ij')

网格点坐标旋转#

我们以在原坐标中定义了规则网格 $\left( x, y \right)$ 即xgrid和ygrid，选择中心位置 $( x_0, y_0 )$ （也可以是任意位置）和角度 $\theta$ ，进行坐标旋转，获得新的坐标系中对应点的坐标 $( x_1, y_1 )$ ：

\begin{pmatrix} x_1 \\\\ y_1 \end{pmatrix} = \begin{pmatrix} \cos\theta & \sin\theta \\\\ -\sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} x - x_0 \\\\ y - y_0 \end{pmatrix}

代码如下

1
def rotate(x, y, x0, y0, theta):
2
    x1grid = np.zeros([y.size, x.size])
3
    y1grid = np.zeros([y.size, x.size])
4
    rad = np.deg2rad(theta)
5
    rotmat = np.array([[np.cos(rad), np.sin(rad)],
6
                       [-np.sin(rad), np.cos(rad)]])
7
    for i, xx in enumerate(x):
8
        for j, yy in enumerate(y):
9
            x1grid[j, i], y1grid[j, i] = rotmat @ np.array([xx-x0, yy-y0])
10
    return x1grid, y1grid

这里定义 $( \theta=60^\circ )$ ，并绘制出坐标旋转后新坐标系下 $( x, y )$ 对应的坐标 $( x_1, y_1 )$

1
theta=60
2
x1grid, y1grid = rotate(x, y, x0, y0, theta)
3
plt.scatter(x1grid, y1grid, s=2, color='C1')

在新坐标系中设置检测版#

我们在新坐标系下x方向 $( -500 \sim 500m )$ 和y方向 $( -600 \sim 600m )$ 范围内设置检测版，其中x方向设置6个异常体，y方向设置7个异常体，异常大小为12%

1
x_lim_cb = [-500, 500]
2
y_lim_cb = [-600, 600]
3
hx = 2
4
hz = 2
5
x_cb = np.arange(*x_lim_cb, hx)
6
x_pert = np.zeros_like(x_cb)
7
y_cb = np.arange(*y_lim_cb, hx)
8
y_pert = np.zeros_like(y_cb)
9
pert=0.12
10
num_per_x = 6
11
num_per_y = 7
12
x_pert = np.sin(num_per_x*np.pi*np.arange(x_cb.size)/x_cb.size)
13
y_pert = np.sin(num_per_y*np.pi*np.arange(y_cb.size)/y_cb.size)
14
yy, xx = np.meshgrid(y_pert, x_pert, indexing='ij')
15
xy_pert = xx*yy*pert
16

17
plt.pcolor(x_cb, y_cb, xy_pert)
18
plt.gca().set_aspect('equal')

二维插值获得旋转后的检测版#

我们用 $( x_1, y_1 )$ 在xy_pert中插值，即可获得原坐标 $( x, y )$ 中每个点对应的异常体大小。

1
interp = RegularGridInterpolator((y_cb, x_cb), xy_pert, bounds_error=False, fill_value=0.)
2
xy_pert_inter = interp(( y1grid, x1grid))
3
plt.pcolor(xgrid, ygrid, xy_pert_inter)
4
plt.gca().set_aspect('equal')