<h2>TL;博士</h2>
<p>在这个问题上,您需要解决多个挑战,主要是:</p>
<ul>
<li>从梯度(矢量场)重建电势(标量场)</li>
</ul>
<p>而且:</p>
<ul>
<li>非矩形网格凹面船体的观测</李>
<li>数值二维线积分和数值精度</李>
</ul>
<p>似乎可以通过选择一种特殊的插入剂和一种智能的集成方式来解决这个问题(正如<code>@Aguy</code>所指出的)</p>
<h2>MCVE</h2>
<p>首先,让我们构建一个MCVE来突出上述关键点</p>
<h3>数据集</h3>
<p>我们重新创建一个标量场及其梯度</p>
<pre><code>import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
def f(x, y):
return x**2 + x*y + 2*y + 1
Nx, Ny = 21, 17
xl = np.linspace(-3, 3, Nx)
yl = np.linspace(-2, 2, Ny)
X, Y = np.meshgrid(xl, yl)
Z = f(X, Y)
zl = np.arange(np.floor(Z.min()), np.ceil(Z.max())+1, 2)
dZdy, dZdx = np.gradient(Z, yl, xl, edge_order=1)
V = np.hypot(dZdx, dZdy)
</code></pre>
<p>标量字段如下所示:</p>
<pre><code>axe = plt.axes(projection='3d')
axe.plot_surface(X, Y, Z, cmap='jet', alpha=0.5)
axe.view_init(elev=25, azim=-45)
</code></pre>
<p><a href="https://i.stack.imgur.com/JBM2N.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/JBM2N.png" alt="enter image description here"/></a></p>
<p>向量场如下所示:</p>
<pre><code>axe = plt.contour(X, Y, Z, zl, cmap='jet')
axe.axes.quiver(X, Y, dZdx, dZdy, V, units='x', pivot='tip', cmap='jet')
axe.axes.set_aspect('equal')
axe.axes.grid()
</code></pre>
<p><a href="https://i.stack.imgur.com/MNqbf.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/MNqbf.png" alt="enter image description here"/></a></p>
<p>事实上,梯度与势能水平是垂直的。我们还绘制了梯度大小:</p>
<pre><code>axe = plt.contour(X, Y, V, 10, cmap='jet')
axe.axes.set_aspect('equal')
axe.axes.grid()
</code></pre>
<p><a href="https://i.stack.imgur.com/y9nDQ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/y9nDQ.png" alt="enter image description here"/></a></p>
<h3>原始场重建</h3>
<p>如果我们天真地从梯度重建标量场:</p>
<pre><code>SdZx = np.cumsum(dZdx, axis=1)*np.diff(xl)[0]
SdZy = np.cumsum(dZdy, axis=0)*np.diff(yl)[0]
Zhat = np.zeros(SdZx.shape)
for i in range(Zhat.shape[0]):
for j in range(Zhat.shape[1]):
Zhat[i,j] += np.sum([SdZy[i,0], -SdZy[0,0], SdZx[i,j], -SdZx[i,0]])
Zhat += Z[0,0] - Zhat[0,0]
</code></pre>
<p>我们可以看到,总体结果大致正确,但在梯度幅值较低的情况下,水平精度较低:</p>
<p><a href="https://i.stack.imgur.com/DsJYz.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/DsJYz.png" alt="enter image description here"/></a></p>
<h2>插值场重建</h2>
<p>如果我们提高栅格分辨率并选择特定的插值(通常在处理网格栅格时),我们可以获得更精细的场重建:</p>
<pre><code>r = np.stack([X.ravel(), Y.ravel()]).T
Sx = interpolate.CloughTocher2DInterpolator(r, dZdx.ravel())
Sy = interpolate.CloughTocher2DInterpolator(r, dZdy.ravel())
Nx, Ny = 200, 200
xli = np.linspace(xl.min(), xl.max(), Nx)
yli = np.linspace(yl.min(), yl.max(), Nx)
Xi, Yi = np.meshgrid(xli, yli)
ri = np.stack([Xi.ravel(), Yi.ravel()]).T
dZdxi = Sx(ri).reshape(Xi.shape)
dZdyi = Sy(ri).reshape(Xi.shape)
SdZxi = np.cumsum(dZdxi, axis=1)*np.diff(xli)[0]
SdZyi = np.cumsum(dZdyi, axis=0)*np.diff(yli)[0]
Zhati = np.zeros(SdZxi.shape)
for i in range(Zhati.shape[0]):
for j in range(Zhati.shape[1]):
Zhati[i,j] += np.sum([SdZyi[i,0], -SdZyi[0,0], SdZxi[i,j], -SdZxi[i,0]])
Zhati += Z[0,0] - Zhati[0,0]
</code></pre>
<p>哪一种性能更好:</p>
<p><a href="https://i.stack.imgur.com/UIj0E.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/UIj0E.png" alt="enter image description here"/></a></p>
<p>因此,基本上,使用临时插值增加网格分辨率可以帮助您获得更准确的结果。插值还解决了从三角形网格获取规则矩形网格以执行积分的需要</p>
<h2>凹凸壳</h2>
<p>您还指出了边缘的不精确性。这些都是插值选择和集成方法相结合的结果。当标量场到达内插点较少的凹区域时,积分方法无法正确计算标量场。选择能够外推的无网格插值时,问题就消失了</p>
<p>为了说明这一点,让我们从MCVE中删除一些数据:</p>
<pre><code>q = np.full(dZdx.shape, False)
q[0:6,5:11] = True
q[-6:,-6:] = True
dZdx[q] = np.nan
dZdy[q] = np.nan
</code></pre>
<p>然后,插入剂的构造如下所示:</p>
<pre><code>q2 = ~np.isnan(dZdx.ravel())
r = np.stack([X.ravel(), Y.ravel()]).T[q2,:]
Sx = interpolate.CloughTocher2DInterpolator(r, dZdx.ravel()[q2])
Sy = interpolate.CloughTocher2DInterpolator(r, dZdy.ravel()[q2])
</code></pre>
<p>执行积分时,我们发现,除了经典的边缘效应外,我们在凹面区域(外壳为凹面的旋转点虚线)中的准确值较低,并且我们在凸面外壳外没有数据,因为Clough-Tocher是基于网格的插值:</p>
<pre><code>Vl = np.arange(0, 11, 1)
axe = plt.contour(X, Y, np.hypot(dZdx, dZdy), Vl, cmap='jet')
axe.axes.contour(Xi, Yi, np.hypot(dZdxi, dZdyi), Vl, cmap='jet', linestyles='-.')
axe.axes.set_aspect('equal')
axe.axes.grid()
</code></pre>
<p><a href="https://i.stack.imgur.com/J6xFO.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/J6xFO.png" alt="enter image description here"/></a></p>
<p>基本上,我们在拐角处看到的误差,很可能是由于积分问题,再加上限于凸包的插值</p>
<p>为了克服这一问题,我们可以选择不同的插值函数,如RBF(径向基函数核),它能够在凸包之外创建数据:</p>
<pre><code>Sx = interpolate.Rbf(r[:,0], r[:,1], dZdx.ravel()[q2], function='thin_plate')
Sy = interpolate.Rbf(r[:,0], r[:,1], dZdy.ravel()[q2], function='thin_plate')
dZdxi = Sx(ri[:,0], ri[:,1]).reshape(Xi.shape)
dZdyi = Sy(ri[:,0], ri[:,1]).reshape(Xi.shape)
</code></pre>
<p>请注意这个插值器的界面稍有不同(注意parmaters是如何传递的)</p>
<p>结果如下:</p>
<p><a href="https://i.stack.imgur.com/71SGA.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/71SGA.png" alt="enter image description here"/></a></p>
<p>我们可以看到凸包外的区域可以外推(RBF是无网格的)。因此,选择adhoc Interplant绝对是解决问题的关键。但我们仍然需要意识到,外推可能表现良好,但在某种程度上是毫无意义和危险的</p>
<h2>解决你的问题</h2>
<p>由<code>@Aguy</code>提供的答案是非常好的,因为它设置了一种巧妙的集成方式,而不受凸包外缺失点的干扰。但正如你所提到的,在凸面外壳内部的凹面区域存在不精确性</p>
<p>如果您希望删除检测到的边缘效果,则必须求助于能够进行外推的插值工具,或者找到另一种集成方法</p>
<h3>中间剂变化</h3>
<p>使用RBF插值似乎可以解决您的问题。以下是完整的代码:</p>
<pre><code>df = pd.read_excel('./Trial-Wireup 2.xlsx')
x = df['X'].to_numpy()
y = df['Y'].to_numpy()
z = df['Delay'].to_numpy()
r = np.stack([x, y]).T
#S = interpolate.CloughTocher2DInterpolator(r, z)
#S = interpolate.LinearNDInterpolator(r, z)
S = interpolate.Rbf(x, y, z, epsilon=0.1, function='thin_plate')
N = 200
xl = np.linspace(x.min(), x.max(), N)
yl = np.linspace(y.min(), y.max(), N)
X, Y = np.meshgrid(xl, yl)
#Zp = S(np.stack([X.ravel(), Y.ravel()]).T)
Zp = S(X.ravel(), Y.ravel())
Z = Zp.reshape(X.shape)
dZdy, dZdx = np.gradient(Z, yl, xl, edge_order=1)
SdZx = np.nancumsum(dZdx, axis=1)*np.diff(xl)[0]
SdZy = np.nancumsum(dZdy, axis=0)*np.diff(yl)[0]
Zhat = np.zeros(SdZx.shape)
for i in range(Zhat.shape[0]):
for j in range(Zhat.shape[1]):
#Zhat[i,j] += np.nansum([SdZy[i,0], -SdZy[0,0], SdZx[i,j], -SdZx[i,0]])
Zhat[i,j] += np.nansum([SdZx[0,N//2], SdZy[i,N//2], SdZx[i,j], -SdZx[i,N//2]])
Zhat += Z[100,100] - Zhat[100,100]
lz = np.linspace(0, 5000, 20)
axe = plt.contour(X, Y, Z, lz, cmap='jet')
axe = plt.contour(X, Y, Zhat, lz, cmap='jet', linestyles=':')
axe.axes.plot(x, y, '.', markersize=1)
axe.axes.set_aspect('equal')
axe.axes.grid()
</code></pre>
<p>以图形方式呈现如下所示:</p>
<p><a href="https://i.stack.imgur.com/Yck6p.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Yck6p.png" alt="enter image description here"/></a>
<a href="https://i.stack.imgur.com/kktBR.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/kktBR.png" alt="enter image description here"/></a></p>
<p>边缘效应消失了,因为RBF插值可以在整个gr上进行外推id。您可以通过比较基于网格的插值的结果进行确认</p>
<h3>线性</h3>
<p><a href="https://i.stack.imgur.com/jyZT9.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/jyZT9.png" alt="enter image description here"/></a>
<a href="https://i.stack.imgur.com/g86m6.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/g86m6.png" alt="enter image description here"/></a></p>
<h3>克拉夫·托彻</h3>
<p><a href="https://i.stack.imgur.com/smGCD.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/smGCD.png" alt="enter image description here"/></a>
<a href="https://i.stack.imgur.com/bhYl8.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/bhYl8.png" alt="enter image description here"/></a></p>
<h3>积分变量顺序变化</h3>
<p>我们还可以尝试找到更好的方法来整合和缓解边缘效应,例如,让我们更改整合变量顺序:</p>
<pre><code>Zhat[i,j] += np.nansum([SdZy[N//2,0], SdZx[N//2,j], SdZy[i,j], -SdZy[N//2,j]])
</code></pre>
<p>使用经典的线性插值。结果非常正确,但左下角仍有边缘效果:</p>
<p><a href="https://i.stack.imgur.com/YNNkh.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/YNNkh.png" alt="enter image description here"/></a></p>
<p>正如您所注意到的,问题发生在积分开始的区域中的轴的中间,并且缺少参考点</p>