# tbathy-CircularArrayCross-Correlation

**Repository Path**: dragonlalala/tbathy-circular-array-cross-correlation

## Basic Information

- **Project Name**: tbathy-CircularArrayCross-Correlation
- **Description**: 基于圆形阵列互相关算法的近海岸水深反演算法
- **Primary Language**: Unknown
- **License**: Not specified
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 0
- **Created**: 2026-01-09
- **Last Updated**: 2026-04-12

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# 基于圆形阵列互相关法的水深反演算法实现

> 根据《Unmanned Aerial Vehicle Depth Inversion to Monitor River-Mouth Bar Dynamics》论文复现，论文见docs文件夹。



## 第一步：**算法总结与公式推导**

该论文提出了一种基于 **UAV 视频影像** 的 **二维水深反演方法**，核心是 **利用涌浪（swell）在视频中的传播特性**，通过 **时域互相关方法** 提取波速向量场与特征频率，再利用 **线性波色散关系** 反演水深。

### 1. 核心假设
- 海面像素亮度时间序列与水面波面高度具有高度相关性；
- 视频已通过 GCP（地面控制点）正射校正，消除 UAV 摆动影响；
- 波场近似为窄带涌浪（频带限定为 0.05–0.2 Hz）；
- 线性波理论适用（即忽略非线性效应，如破碎波区除外）。

---

### 2. 水深反演公式（基于线性色散关系）

论文采用经典的线性波色散关系：

\[
c = \frac{gT}{2\pi} \tanh\left(\frac{2\pi h}{gT^2}\right)
\]

等价地，用频率 \( f = 1/T \) 和波速 \( c \) 表示水深 \( h \)：

\[
\boxed{
h = \frac{c}{2\pi f} \tanh^{-1}\left( \frac{2\pi f c}{g} \right)
}
\quad \text{(Eq. 1)}
\]

其中：
- \( c \)：局部波速（celerity，单位 m/s）；
- \( f \)：局部主导波频率（Hz）；
- \( g = 9.81 \, \text{m/s}^2 \)；
- \( h \)：水深（m）。

> 注意：该公式隐式依赖 \( c \) 和 \( f \)，需先估计这两个量。

---

### 3. 波速与传播方向估计（圆形阵列互相关法）

#### 3.1 方法概述
对每个分析点 \( A \)，在其周围半径为 \( r \) 的圆周上设 \( N \) 个参考点 \( B_n \)（\( n = 1,...,N \)）。通过互相关计算波从前端 \( A \) 传播到 \( B_n \) 的时间延迟 \( t_n \)。

#### 3.2 几何关系
若波以方向角 \( \theta \)（从正东逆时针测量，或按论文定义）传播，则投影距离为：

\[
\text{投影距离} = r \cos(\theta - \theta_n)
\]

因此时间延迟为：

\[
\boxed{
t_n = \frac{r \cos(\theta - \theta_n)}{c}
}
\quad \text{(Eq. 3)}
\]

其中 \( \theta_n \) 是 \( B_n \) 相对于 \( A \) 的方位角。

#### 3.3 非线性最小二乘求解
将 Eq. (3) 线性化（围绕初值 \( c_0, \theta_0 \)）：

\[
t_n = F_n(c_0, \theta_0) + \left( \frac{\partial F_n}{\partial c} \right)_0 dc + \left( \frac{\partial F_n}{\partial \theta} \right)_0 d\theta + v_n
\]

写成矩阵形式：

\[
\mathbf{Y} \, d\mathbf{X} = \mathbf{L} - \mathbf{V}
\]

其中：
- \( d\mathbf{X} = [d\theta, dc]^T \)
- \( \mathbf{Y}_{n,:} = [ (\partial F_n / \partial \theta)_0,\ (\partial F_n / \partial c)_0 ] \)
- \( \mathbf{L}_n = t_n - F_n(c_0, \theta_0) \)

#### 3.4 加权最小二乘（考虑信号质量）
引入权重矩阵 \( \mathbf{P} \)，基于互相关系数 \( R_n \)：

\[
P_{nn} = \frac{R_n}{\bar{R}}, \quad \bar{R} = \frac{1}{N} \sum_{n=1}^N R_n
\]

正规方程为：

\[
\boxed{
(\mathbf{Y}^T \mathbf{P} \mathbf{Y}) d\mathbf{X} = \mathbf{Y}^T \mathbf{P} \mathbf{L}
}
\quad \text{(Eq. 7)}
\]

迭代求解直至 \( |dc| < 0.01 \, \text{m/s} \) 且 \( |d\theta| < 0.01^\circ \)。

---

### 4. 主导频率估计（加权平均交叉谱法）

对每个点 \( A \)，利用其与周围 \( B_n \) 的交叉谱 \( C_{AB,n}(f) \) 计算加权平均频率：

\[
\boxed{
f_r = \frac{1}{N} \sum_{n=1}^N \frac{ \int_{f_1}^{f_2} |C_{AB,n}(f)| f \, df }{ \int_{f_1}^{f_2} |C_{AB,n}(f)| \, df }
}
\quad \text{(Eq. 8)}
\]

其中 \( f_1 = 0.05 \, \text{Hz}, f_2 = 0.2 \, \text{Hz} \)。

---

### 5. 后处理
- 仅保留满足 \( t_n \) 沿圆周呈正弦变化的点（否则视为噪声/破碎区/多向波）；
- 反演水深 \( h \) 后，减去潮位（TWL）得海床高程；
- 最终结果用 \( 5 \times 5 \) 高斯滤波平滑。

---

### 6. 推荐参数（来自论文验证）
- 圆形半径 \( r = 20 \, \text{m} \)
- 参考点数 \( N = 8 \)
- 视频长度 ≥ 5 分钟（采样率 6 Hz）
- 网格分辨率：2 m
- 有效水深范围：0 – 8 m（受波周期限制）

---


### 7. 实验结果
phantom4RTK_202205183数据集：
![](Results/phantom4RTK_202205183/depth_map.png)

## 待办清单
1. 高斯平滑结果图
2. 研究波速反演不出来的区域(✔)
3. 高程图无法显示的问题(✔)
4. 流程优化，减少重复的计算
5. 了解频率计算原理(✔)
6. 引入参数文件，将参数集中管理(✔)
7. 保存相关的mat，用于多方法对比(✔)