水深測量資料問題

這是一道年代較早的數學建模競賽題，主題是水深測量資料。雖然題目很舊，有限測點下重建空間表面仍然是現代建模的重要問題，也能展示由資料、假設到決策區域的完整推理。

問題背景

下表給出低潮時，平面座標 $(X,Y)$ 上的水深 $Z$。$X,Y$ 以 yard 計，$Z$ 以 ft 計。一艘船的吃水深度為 5 ft；在矩形 $(75,200)\times(-50,150)$ 內，哪些區域應避免？

分析

直接在每個位置量度水深可能昂貴而耗時，因此需要由有限測點估計區域內的水深。任務不是只畫一張平滑圖，而是決定吃水 5 ft 的船應避開哪些位置。

這是一個空間插值問題。資料稀疏、不規則且有兩個空間座標，傳統一維插值不合適。本文採用徑向基底函數插值（Radial Basis Function，RBF）。

數學模型

徑向基底函數只取決於輸入與中心之間的距離：

其中 $x$ 是輸入向量，$c$ 是中心，$\|\cdot\|$ 是 Euclidean 距離，$\phi$ 是指定的徑向核。插值通常寫成

並選擇權重 $w_j$，使 $s(x_i)=z_i$。因此，權重不是單純按每次查詢位置直接分配，而是由整個插值線性系統共同求得。

基本步驟是：

1. 以測點作 RBF 中心；
2. 選擇核函數及形狀參數；
3. 解出令插值穿過測量值的權重；
4. 在規則網格評估插值面；
5. 以水深界線標示避免及警戒區。

Python 實作

先匯入函式庫及定義資料：

import numpy as np
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt

points = np.array([
    [129.0, 7.5], [140.0, 141.5], [108.5, 28.0], [88.0, 147.0],
    [185.5, 22.5], [195.0, 137.5], [105.5, 85.5], [157.5, -6.5],
    [107.5, -81.0], [77.0, 3.0], [162.0, -66.5], [162.0, 84.0],
    [117.5, -38.5]
])

depths = np.array([4, 8, 6, 8, 6, 8, 8, 9, 9, 8, 9, 4, 9])

若只回答題目指定矩形，網格應使用 $x\in[75,200]$、$y\in[-50,150]$：

x = np.linspace(75, 200, 200)
y = np.linspace(-50, 150, 200)
X, Y = np.meshgrid(x, y)

原實驗採用 multiquadric RBF：

rbf = Rbf(
    points[:, 0],
    points[:, 1],
    depths,
    function='multiquadric',
    epsilon=2
)
Z = rbf(X, Y)

epsilon=2 會影響曲面的平滑與條件數，不能只憑單一圖片決定。較嚴謹研究應以留一測點交叉驗證或敏感度分析比較多個核與形狀參數。

為吃水加入 0.5 ft 額外餘量，等價於把避免界線設在預測原始水深 $5.5$ ft：

safety_margin = 0.5
Z_safe = Z - safety_margin

mask_unsafe = Z_safe <= 5
mask_caution = (Z_safe > 5) & (Z_safe < 7)

原分析亦計算插值面梯度：

zy, zx = np.gradient(Z)
gradient_magnitude = np.sqrt(zx**2 + zy**2)

再以測點距離及梯度方向建立各向異性距離代理：

def anisotropic_distance(x, y, point, zx, zy):
    dx = x - point[0]
    dy = y - point[1]
    direction = np.arctan2(dy, dx)
    gradient_direction = np.arctan2(zy, zx)
    angle_diff = np.abs(direction - gradient_direction)
    return np.sqrt(dx**2 + dy**2) * (
        1 + gradient_magnitude * np.cos(angle_diff)
    )

distances = np.array([
    anisotropic_distance(X, Y, point, zx, zy)
    for point in points
])
distance_proxy = np.min(distances, axis=0)
uncertainty = np.clip(
    distance_proxy / np.max(distance_proxy), 0, 1
)

這個 uncertainty 只是一個依測點距離與局部梯度而定的視覺代理，不是由誤差模型推導的標準誤、置信區間或超越機率。原公式也可能因 $\cos(\text{angle difference})<0$ 而縮小甚至扭曲「距離」，因此不應直接用於航行安全決策。

最後繪製水深、避免界線、警戒界線、代理陰影及測點：

plt.figure(figsize=(12, 10))

contour = plt.contourf(
    X, Y, Z_safe,
    levels=np.linspace(0, 10, 21),
    cmap='viridis',
    alpha=0.7
)
plt.colorbar(contour, label='Depth (feet) with safety margin')

plt.contour(
    X, Y, mask_caution,
    levels=[0.5], colors='yellow', linewidths=2
)
plt.contour(
    X, Y, mask_unsafe,
    levels=[0.5], colors='red', linewidths=4
)

uncertainty_contour = plt.contourf(
    X, Y, uncertainty,
    levels=np.linspace(0, 1, 11),
    cmap='Greys',
    alpha=0.3
)
plt.colorbar(
    uncertainty_contour,
    label='Distance proxy (darker = farther from support)'
)

plt.scatter(
    points[:, 0], points[:, 1],
    c=depths, cmap='viridis',
    edgecolor='k', s=50, zorder=3
)

plt.title('Water Depth Map with Safety Margin')
plt.xlabel('X (yards)')
plt.ylabel('Y (yards)')
plt.xlim(75, 200)
plt.ylim(-50, 150)
plt.gca().set_aspect('equal', adjustable='box')
plt.grid(True)
plt.tight_layout()
plt.show()

原專案輸出如下：

如何解讀結果

紅色 $5$ ft 安全餘量後界線，對應原始預測水深不超過 $5.5$ ft 的位置；這是模型下應避免的候選區域。黃色範圍則表示更寬鬆的警戒帶。

不過，RBF 插值會在測點之外仍輸出平滑數值，並不代表該處已有證據。尤其靠近邊界或遠離測點時，不同核、epsilon 或少量量度誤差都可能大幅移動 $5.5$ ft 等深線。真正安全分析應至少：

1. 檢查測點座標是否都位於題目範圍及潮位基準一致；
2. 比較 RBF、kriging、薄板樣條等替代方法；
3. 用留一法評估預測誤差；
4. 對水深誤差建立機率模型；
5. 報告「水深低於安全界線的機率」，而不只是一條插值等高線。

結語

RBF 是處理稀疏、不規則多維資料的有用工具，能把有限水深測點轉成連續表面，再提出應避免區域。但這張表面是模型輸出，不是真實海床的完整量度；安全餘量、核函數與不確定性方法都會影響決策邊界。

這道題最有價值的建模啟示，不只是選擇一個插值函數，而是把「哪裏資料不足」與「哪裏預測水深不足」同時放入最終答案。