海面之下：HKUST 近岸波高的間接量度

封存文章

本文源自 IMMC 2025 Problem B。舊版本用合成 ADCP signal 產生波高，再把均勻 10 分鐘觀測錯配到 Chebyshev nodes，不能證明現場精度或一分鐘處理能力。本版本改為可部署的量度鏈，所有精度仍須以現場 reference gauge 驗證。

題目要求估計 HKUST Sai Kung 現有岸線外 10 米位置的 wave height，每 10 分鐘輸出一次，連續 24 小時；影像可能在颱風、暴雨和低照度下取得。這不是單純 interpolation 問題，而是從原始壓力／流速／影像訊號，經過 transfer function、品質控制和現場校準，得到統計波參數。

主要輸出應先定義：

surface elevation time series $\eta(t)$ ；
significant wave height $H_{m0}=4\sqrt{m_0}$ ；
peak／mean period；
wave direction（若設備支持）；
每個 10 分鐘 burst 的 quality flag 和 uncertainty。

一、量度架構

較穩健的設計是多層感測：

海床壓力感測器／wave-capable ADCP： 主要水下訊號；
岸上相機： 在可見度足夠時提供 wave run-up／surface tracking；
水位／潮位 reference： 分離 tide 與 wave；
氣象站： 風、氣壓、雨量；
短期 reference deployment： 同位置或鄰近位置的 calibrated wave gauge，用於驗證而非永久依賴。

如果目標是「間接方法」，可以把壓力或影像當 proxy，再以短期直接量度建立 transfer model。把遠處 Kau Yi Chau 的平均值直接乘 seasonal／cyclone factors，不能還原 HKUST 近岸 10 米處的即時波高。

二、壓力到自由水面

海床 pressure sensor 測得：

p(t)=p_{\rm atm}(t)+\rho g[d+\eta_p(t)]+\epsilon_p(t).

先減 atmospheric pressure 和低頻 tide，得到 wave-band pressure fluctuation。在線性波理論下，某頻率 $f$ 的 pressure response factor 為

K_p(f) = \frac{\cosh k(f)(z+d)} {\cosh k(f)d},

其中 $z$ 是 sensor 相對 mean water level 的垂向位置， $d$ 是水深，wave number $k$ 解 dispersion relation：

\omega^2=gk\tanh(kd), \qquad \omega=2\pi f.

自由水面頻譜近似：

S_{\eta\eta}(f) = \frac{S_{pp}(f)} {(\rho gK_p(f))^2}.

在 $K_p$ 很小的高頻，除法會放大 noise，所以要根據 depth、sensor noise 和 coherence 設 cut-off，而不是無限制反演。

import numpy as np
from scipy.optimize import brentq
from scipy.signal import welch

G = 9.81

def wave_number(frequency_hz, depth_m):
    omega = 2 * np.pi * frequency_hz
    if omega == 0:
        return 0.0
    fn = lambda k: G * k * np.tanh(k * depth_m) - omega**2
    return brentq(fn, 1e-8, max(omega**2 / G * 100, 10))

def pressure_response(frequency_hz, depth_m, sensor_z_m):
    k = wave_number(frequency_hz, depth_m)
    return np.cosh(k * (sensor_z_m + depth_m)) / np.cosh(k * depth_m)

舊程式用

4 * np.std(velocities)

估 wave height。 $4\sigma$ 可由 surface-elevation variance 估 $H_{m0}$ ，但 orbital velocity 必須先經 frequency-dependent transfer function 轉成 $\eta$ ；不能直接把任意 Doppler-like array 的標準差稱為波高。

三、每 10 分鐘的 spectral estimate

每 10 分鐘收集一個 high-rate burst，例如 2–4 Hz：

timestamp synchronization；
移除 missing／spike；
atmospheric 與 tide correction；
detrend；
Welch spectrum；
pressure／velocity transfer；
在可信 frequency band 積分。

頻譜矩：

m_n=\int_{f_{\min}}^{f_{\max}}f^nS_{\eta\eta}(f)\,df.

因此：

H_{m0}=4\sqrt{m_0}, \qquad T_{m02}=\sqrt{\frac{m_0}{m_2}}.

def significant_wave_height(eta, fs_hz):
    f, spectrum = welch(
        eta,
        fs=fs_hz,
        window="hann",
        nperseg=min(len(eta), 256),
        noverlap=128,
        detrend="linear",
    )
    band = (f >= 0.04) & (f <= 0.5)
    m0 = np.trapezoid(spectrum[band], f[band])
    return 4 * np.sqrt(max(m0, 0))

frequency band 必須按 sampling、depth 和 site spectrum 驗證。

四、由 offshore reference 轉到近岸

如果主資料來自較外海站，需 wave transformation，而不是 spatial polynomial interpolation。對單一 component，energy flux 在無 dissipation 情況近似守恆：

E C_g b=\text{constant}, \qquad E=\frac18\rho gH^2.

沿 wave ray 的 shoaling coefficient：

K_s=\sqrt{\frac{C_{g0}}{C_g}},

refraction coefficient $K_r$ 由 bathymetry 和 Snell law 得到：

\frac{\sin\theta}{C} =\text{constant}.

近岸還有 bottom friction、breaking、diffraction 和 sheltering：

H_{\rm near} \approx H_{\rm off}K_sK_rK_fK_b.

在 landing step 附近的複雜岸線，應使用 survey bathymetry 和經驗證的 spectral wave model（例如 SWAN 類工具），再用現場 sensor calibration。

Breaking limit

淺水波高受限於

H_b\lesssim\gamma_b d,

$\gamma_b$ 常在約 0.7–0.8 的量級，但依 slope 和 wave conditions 改變。任何模型在 10 米離岸點預測超過 local depth 的巨大波高，都應被 QC flag。

五、影像方法

低照度暴雨影像不適合只靠 RGB edge detection。可使用：

固定、剛性相機和防水罩；
known control points／fiducial markers；
相機內外參 calibration；
polarizing／infrared illumination（符合法規與安全）；
temporal denoising；
optical flow 或 waterline segmentation；
rain streak／glare quality flag。

由像素 $(u,v)$ 到世界坐標平面的 homography：

\lambda \begin{bmatrix}u\\v\\1\end{bmatrix} = \mathbf H \begin{bmatrix}X\\Y\\1\end{bmatrix}.

但 wave surface 是三維且動態，單一 homography 只對近似平面有效。若由岸上相機估水線／run-up，需用 reference targets 轉換 vertical scale，再與壓力 sensor 同步比較。

影像在 lens 被雨水遮蔽或夜間 saturation 時應降權，而不是硬輸出。

六、感測器融合

令 latent wave height 為 $H_t$ ，不同 sensor observation：

\begin{aligned} y_t^{(p)}&=H_t+b_p+\epsilon_t^{(p)},\\ y_t^{(v)}&=a_vH_t+b_v+\epsilon_t^{(v)},\\ y_t^{(c)}&=a_cH_t+b_c+\epsilon_t^{(c)}. \end{aligned}

可用 Kalman／Bayesian state-space model 估計，但 noise variance 應由 calibration data 取得，並隨 quality flags 改變：

\operatorname{Var}(\epsilon_t^{(c)}) \uparrow \quad\text{when rain/visibility flag is poor}.

「三個 sensor」不會自動帶來可信度；它們要有獨立 failure modes 和可追查 calibration。

七、為何 Chebyshev 插值不適用

舊版本把 144 個每 10 分鐘觀測值配到

x_j=\cos\frac{(2j+1)\pi}{2n}

再縮放至 0–24 小時。問題是 sensor observation 本來在均勻時間：

t_j=j/6\ {\rm hour}.

把數值重新貼到 Chebyshev timestamps 會改寫資料的發生時間。 Chebyshev nodes 是設計可控制採樣節點時的 polynomial approximation 工具；題目已固定每 10 分鐘觀測，不能事後搬動節點。

此外，degree-143 global polynomial 對 noisy wave data 完全不必要。若只有短缺口，可用 local linear／spline 或 state-space smoother，並把 imputed points 標示；長缺口不應偽造。

Barycentric 權重

對固定 nodes $x_j$ ：

w_j =\frac{1}{\prod_{k\ne j}(x_j-x_k)}.

舊文寫成 $1/(x-x_j)^2$ ，它依 evaluation point $x$ 改變，並不是 barycentric weight。

八、錯誤的「95% confidence interval」

舊圖對每一時間點畫：

H_t \pm t_{0.975,n-1}\frac{s_H}{\sqrt n}.

這是整批樣本 mean 的標準誤近似，不能代表每個瞬時 wave-height estimate 的 uncertainty，也忽略 serial correlation。

每個 10 分鐘估計的不確定性應來自：

finite spectral degrees of freedom；
pressure transfer uncertainty；
depth／tide error；
sensor calibration；
fusion model；
data quality。

可用 segment bootstrap 或 spectrum-specific confidence interval，再經 Monte Carlo 傳播 transfer parameters。

九、舊圖像的合適定位

舊 time series 和 3D surface 由 simulate_adcp_signal() 生成，並非現場測量。3D surface 又用

Z=H(t)e^{-0.1x}

人為製造 spatial decay，沒有 bathymetry、energy flux 或 breaking model。因此：

0.15–6.80 m range 是歷史值／設定，不是 model validation；
Kau Yi Chau 平均 0.41 m 不能驗證 HKUST site；
圖形只能標成 synthetic demonstration；
不能稱系統在颱風中 robust。

十、設備清單

海上

wave-capable pressure sensor 或 ADCP；
anti-corrosion mount／concrete frame；
redundant pressure sensor；
tide／water-level reference；
waterproof logger；
synchronized clock；
protected power／battery；
telemetry（若可安全可靠）；
retrieval marker 和 maintenance plan。

岸上

fixed calibrated low-light camera；
weatherproof housing、lens wiper／hood；
control points；
local edge computer；
weather station；
UPS 和 local storage。

軟件

acquisition daemon；
time synchronization；
spectral processing；
QC rules；
sensor-fusion service；
database；
monitoring／alert dashboard；
audit log 和 versioned calibration。

十一、實時 pipeline 與 latency

每個 10 分鐘 window 完結後：

步驟	目標時間
data close／checksum	5–10 s
QC／detrend	5 s
spectrum／transfer	5–15 s
image feature（如可用）	10–30 s
fusion／uncertainty	5 s
database／dashboard	5–10 s

在本地 edge hardware 上，計算少於一分鐘在工程上可能，但要用 benchmark 證明。如果 telemetry 中斷，系統應在本地儲存，稍後補傳；不可將 missing data 當成平浪。

十二、驗證計劃

Laboratory

pressure calibration；
clock drift；
power endurance；
waterproofing；
synthetic sinusoid／known spectrum unit tests。

Field intercomparison

在正常和較高浪況下，共置 reference wave gauge：

\operatorname{bias} =\frac1N\sum_t(\hat H_t-H_t^{\rm ref}),

\operatorname{RMSE} = \sqrt{\frac1N\sum_t(\hat H_t-H_t^{\rm ref})^2}.

另報 correlation、spectral agreement、event peak error、availability 和 interval coverage。

Extreme conditions

颱風時不應安排人員臨時下海。設備必須預先部署，並用 historical storm、 wave flume／synthetic injection 和實際遠端資料逐步驗證。

十三、結構設計不能只用一天資料

24 小時量度可測系統與捕捉一次 event，但 marine structure design 的 extreme wave 需要長期 hindcast／observation、extreme-value analysis、return period 和 joint water-level–wave conditions。可用 peaks-over-threshold：

\Pr(H-u\le y\mid H>u) \approx 1-\left(1+\xi\frac y\sigma\right)^{-1/\xi},

但 threshold、independence、nonstationarity 和 uncertainty 都需專業處理。本文的量度系統只供建立／驗證資料鏈，不直接給出 landing step design wave。

結論

可靠的近岸波高系統由「量度鏈」而不是高階插值構成：

壓力／流速訊號經 dispersion-based transfer 變成 surface spectrum；
每 10 分鐘以 spectral moments 得到 $H_{m0}$ ；
遠處觀測要經 bathymetry、shoaling、refraction 和 breaking 轉換；
低照度影像只在 calibration 和 quality flags 下輔助；
sensor fusion 同時輸出 estimate、uncertainty 和 failure state；
現場共置 reference 才能證明精度。

Chebyshev nodes、合成訊號和歷史範圍不能代替驗證。修正後的方案不再聲稱已在颱風下達到 robust real-time accuracy，而是提出一個可部署、可校準和可失敗得透明的系統。

參考資料

International Mathematical Modeling Challenge. (2025). Problem B: Nearshore wave measurement for the HKUST landing step.
Holthuijsen, L. H. (2007). Waves in Oceanic and Coastal Waters. Cambridge University Press.
Dean, R. G., & Dalrymple, R. A. (1991). Water Wave Mechanics for Engineers and Scientists. World Scientific.
Young, I. R. (1995). The determination of confidence limits associated with estimates of the spectral peak frequency. Ocean Engineering, 22(7), 669–686.
Booij, N., Ris, R. C., & Holthuijsen, L. H. (1999). A third-generation wave model for coastal regions. Journal of Geophysical Research: Oceans, 104(C4), 7649–7666.