1 簡介

AIDW 主要是針對 IDW 的缺點進行了改進，考慮了樣本點與預測點的位置，即方向和距離，具體見下圖：

2 改進

IDW 公式：

從IDW算法可看出，插值點的估算值僅與插值樣本距插值點的遠近相關，并未考慮樣本點的方位性（即樣本點被表示為各向同性）。

IDW 插值的基本假設是樣點在插值區呈均勻分布。但眾多情況下，樣點在各向分布并非均勻，甚至會出現樣點集中于某一方向的現象，違背了基本假設，其插值合理性就難被保證。針對 IDW 這一插值局限，作者提出了調和反距離權重（AIDW）插值算法。

AIDW 增加了可反映插值點與樣本點方位關系的調和權重系數 K，其基本假設是：距插值點近的樣本點，對其后方的樣本點有遮蔽效應，當兩樣本點與插值點的連線夾角 $\alpha <360°/n$（n 為插值搜索鄰域內的樣點個數）時，遮蔽效應存在，當 $\alpha \ge 360°/n$ 時，遮蔽效應消失。在 AIDW 插值過程中，受遮蔽影響的樣本點，其插值權重將被削弱，削弱的程度取決于該樣點 K 值的大小。

按上述假設：

• 圖1(a) 所示的 5 個樣點在方向上均勻地分布在插值點（中心點）周圍，任意兩樣點與插值點的連線夾角均大于或等于 72°（即α $\alpha \ge 360°/5$），即認為該５個樣點間相互不存在遮蔽效應；
• 在圖1?中，任意兩樣點與插值點的連線夾角均小于72°，即認為距插值點的近樣點，對其后的樣點均具有遮蔽效應；在大多情況下，樣點在插值點周圍的分布應類似圖1(b)，既不像圖1(a)均勻分布，也不像圖1?集中分布。
• 圖1(b)中 $Z_1$、$Z_3$ 對任一樣點均無遮蔽，$Z_2$ 對 $Z_4$、$Z_5$ 有遮蔽，$Z_4$ 對 $Z_5$ 也有遮蔽。

將 IDW 傳統的算法思想與本文的基本假設結合，提出了 AIDW 算法:

${\sin }^p\theta$的理解：

• 從下圖(a)可以看出，$Z_i$ 逐漸移向 $Z_0$ 的過程種，$\theta$ 逐漸增大，當三者形成等腰三角形時，$\theta =90°$，此時最大，即 ${\sin }^p\theta =1$，$Z_j$ 不會對 $Z_i$ 產生遮蔽影響。
• 從下圖(b)可以看出，$Z_i$ 保持與 $Z_0$相同的距離向 $Z_j$ 移動，當兩者位于同一條線時，$Z_i$的影響完全被遮蔽，即$\theta =0°,{\sin }^p\theta =0$。

計算樣例：

按AIDW算法，在圖3(b)中因$Z_1$對$Z_6$、$Z_3$對$Z_7$和$Z_8$、$Z_4$對$Z_7$有遮蔽影響，這些受遮蔽樣點的插值權重被削減，$Z_{10}、Z_{11}、Z_{12}$分別被$Z_4$、$Z_3$、$Z_7$ 完全遮蔽，它們的插值權重降至為０。依照式(2)和式(3)，最終插值點估算值的計算式為：

• $Z$ 為插值點（中心點）估算值
• $Z_1?Z_9$為樣點觀測值
• $d_1?d_9$為樣點與插值點的歐氏距離
• p 是冪指數
• $\theta$ 角如圖3(b)所示

3 程序設計流程

AIDW 的插值程序可分為插值前準備和插值計算兩個過程：

• 插值前準備主要是用于搜索合適的插值樣點，并為下一步的插值計算提供 $d_i$ 和 $f_{ij}$ 值；
• 插值計算過程主要是求算反映樣點遮蔽程度的 ${\sin }^p{\theta }_{ij}$ 值，并結合 $d_i、z_i$ 值求算插值點的Z值。

• 插值搜索鄰域的大小以格點數(ｋ×ｋ)表示
• m 是搜索鄰域內的樣點數
• n 是插值所需的樣點數
• d、f 分別為樣點與插值點的歐氏距離和兩樣點間的歐氏距離
• i、j、u、v 均為插值樣點的序號

4 插值結果

• 對比 M 點(黑色框)，IDW 出現孤立圓現象(站點集中于一側)，AIDW 消除了孤立圓現象
• 對比 C 點(紅色框)，IDW 出現同心圓現象，AIDW 消除了同心圓現象
• 對比 K 點(黃色框)，兩者均出現孤立圓，通過分析，K 點周圍的站點分布均勻。

從上圖可以看出 AIDW 有效改進了 IDW 的缺點，同時又能保留 IDW 的插值思想，但 AIDW 需要計算 $\theta$ ，因此在插值時間上要比 IDW 慢。

5 python 實現

from sklearn.neighbors import NearestNeighbors"""類函數"""
class AIDW:def __init__(self, x, y, m,p=2):self.m = m # 搜索鄰域內的樣點數self.nbrs = NearestNeighbors(n_neighbors=m, algorithm='ball_tree').fit(x)self.thresh = 360/mself.p = pself.y = yself.x = xdef fit(self, x_new):indices = self.nbrs.kneighbors(x_new, return_distance=False)x_sample = self.x[indices[0]]y_sample = self.y[indices[0]]x_ = x_sample-x_newzi = []ki = 1for i in range(1, self.m-1):for j in range(i):cos_ = np.sum(x_[i]*x_[j])/((np.sum(x_[i]**2)**(1/2))*(np.sum(x_[j]**2)**(1/2)))radian = np.arccos(cos_)angle = radian*180/np.pi if angle>=self.thresh:continueelse:ki*=np.sin(radian)**self.pdi = np.sum(x_[i]**2)**(1/2)zi.append(ki/(di**self.p))z = np.sum(np.array(zi)*y_sample[1:-1])/np.sum(zi)return z"""demo"""
import numpy as np
import matplotlib.pyplot as plt# create sample points with structured scores
X1 = 10 * np.random.rand(1000, 2) -5def func(x, y):return np.sin(x**2 + y**2) / (x**2 + y**2 )z1 = func(X1[:,0], X1[:,1])# 'train'
aidw = AIDW(X1, z1, m=15)# 'test'
spacing = np.linspace(-5., 5., 100)
X2 = np.meshgrid(spacing, spacing)
grid_shape = X2[0].shape
X2 = np.reshape(X2, (2, -1)).T
z2 = []
for x2 in X2:z2.append(aidw.fit(x2.reshape(1,-1)))
z2 = np.array(z2)# plot
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(10,3))
ax1.contourf(spacing, spacing, func(*np.meshgrid(spacing, spacing)))
ax1.set_title('Ground truth')
ax2.scatter(X1[:,0], X1[:,1], c=z1, linewidths=0)
ax2.set_title('Samples')
ax3.contourf(spacing, spacing, z2.reshape(grid_shape))
ax3.set_title('Reconstruction')
plt.show()

參考：
顧及方向遮蔽性的反距離權重插值法.李正泉.
An Adjusted Inverse Distance Weighted Spatial Interpolation Method.