<think>我们划分了四个区域(A, B, C, D),现在要为每个区域拟合一个平面。平面方程的一般形式为:z = a*x + b*y + c
其中,x和y是位置坐标(海里),z是水深(米)。对于每个区域内的点,我们可以通过线性回归来拟合平面。
步骤:
1. 遍历每个区域,提取该区域内所有的网格点坐标(x, y)及其对应的水深值z。
2. 对每个区域内的点,使用最小二乘法拟合平面方程(即求系数a, b, c)。
3. 计算每个区域的拟合平面方程,并计算拟合误差(如均方根误差RMSE)。
4. 可视化:在等高线图上绘制拟合平面(可以用网格平面表示,但为了不遮盖原等高线,可以选择绘制拟合平面的等高线或仅显示区域内的平面)。
注意:由于区域是矩形,我们可以提取区域内的所有网格点(这些点已经是等间距的网格)。
我们将为每个区域生成一个拟合平面,并输出方程和误差。
代码实现:
- 使用numpy的linalg.lstsq进行最小二乘拟合。
- 对于每个区域,构建设计矩阵A(每行为[x, y, 1]),以及观测值向量B(水深值)。
- 求解线性方程组:A * [a, b, c]^T = B
然后,我们可以计算每个区域拟合平面的水深值,并与原始水深值比较。
最后,我们可以选择绘制每个区域的拟合平面(用等高线表示)或者将整个区域的拟合结果绘制在一张图上(用颜色区分不同区域的拟合平面)。
由于问题要求基于划分的区域拟合每个区域的平面,我们将分别处理每个区域,并输出拟合结果。
我们将修改代码,在原有图上增加拟合平面的等高线(用虚线表示),以区别于原始等高线。
但注意:原始等高线是整个区域的,而拟合平面是按区域分别拟合的,因此每个区域内的拟合平面是独立的。
另一种做法:在四个区域分别绘制拟合平面的等高线(只在该区域内绘制),这样整个图就会显示四个平面拼接而成的曲面。
我们采用第二种做法:分别绘制每个区域内的拟合平面等高线。
步骤:
1. 对每个区域,获取其网格坐标(X_nmi_sub, Y_nmi_sub)和原始水深(depth_sub)。
2. 拟合该区域的平面。
3. 用拟合的平面方程计算该区域网格点的水深值(z_fit)。
4. 在该区域内绘制等高线(虚线,使用与原始等高线不同的颜色,比如红色)。
这样,整个图就会显示原始等高线(蓝色)和四个拟合平面(红色虚线)的拼接。
同时,我们输出每个区域的拟合方程和RMSE。
注意:在区域边界处,由于拟合平面是独立计算的,可能会出现不连续(跳变),这是正常的。
代码实现如下:</think>以下是为每个划分区域拟合平面的完整代码实现。该方法通过最小二乘法为每个区域计算最佳拟合平面,并可视化拟合结果:
```python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.patches as patches
from mpl_toolkits.mplot3d import Axes3D
from scipy.linalg import lstsq
# 常量定义
X_RANGE_NMI = 4 # 东西方向4海里
Y_RANGE_NMI = 5 # 南北方向5海里
STEP_NMI = 0.02 # 步长0.02海里
# 计算网格点数量
x_points = int(X_RANGE_NMI / STEP_NMI) + 1
y_points = int(Y_RANGE_NMI / STEP_NMI) + 1
# 从Excel读取数据(替换为实际文件路径)
df = pd.read_excel('depth_data.xlsx', header=None)
depth_data = df.values # 转换为numpy数组
# 验证数据尺寸
if depth_data.shape != (y_points, x_points):
raise ValueError(f"数据尺寸不匹配! 期望: ({y_points}, {x_points}), 实际: {depth_data.shape}")
# 创建坐标网格(转换为海里)
x_coords_nmi = np.linspace(0, X_RANGE_NMI, x_points)
y_coords_nmi = np.linspace(0, Y_RANGE_NMI, y_points)
X_nmi, Y_nmi = np.meshgrid(x_coords_nmi, y_coords_nmi)
# 定义区域边界
regions = [
{'x': [0, 2.5], 'y': [0, 2], 'name': '区域A', 'color': '#e0f3f8', 'fit_color': '#4575b4'},
{'x': [0, 2.5], 'y': [2, 5], 'name': '区域B', 'color': '#ffffbf', 'fit_color': '#d73027'},
{'x': [2.5, 4], 'y': [0, 4], 'name': '区域C', 'color': '#fee090', 'fit_color': '#fdae61'},
{'x': [2.5, 4], 'y': [4, 5], 'name': '区域D', 'color': '#d9ef8b', 'fit_color': '#91cf60'}
]
# 为每个区域拟合平面
for region in regions:
# 获取区域内的网格点索引
x_mask = (x_coords_nmi >= region['x'][0]) & (x_coords_nmi <= region['x'][1])
y_mask = (y_coords_nmi >= region['y'][0]) & (y_coords_nmi <= region['y'][1])
# 提取区域内的坐标和水深数据
X_region = X_nmi[y_mask][:, x_mask].flatten()
Y_region = Y_nmi[y_mask][:, x_mask].flatten()
Z_region = depth_data[y_mask][:, x_mask].flatten()
# 构建最小二乘问题的设计矩阵
A = np.c_[X_region, Y_region, np.ones_like(X_region)]
# 使用最小二乘法拟合平面 z = a*x + b*y + c
coeffs, _, _, _ = lstsq(A, Z_region)
a, b, c = coeffs
# 计算拟合平面在该区域内的水深值
Z_fit = a * X_region + b * Y_region + c
# 计算拟合误差 (RMSE)
rmse = np.sqrt(np.mean((Z_region - Z_fit) ** 2))
# 存储拟合结果
region.update({
'coeffs': (a, b, c),
'rmse': rmse,
'X_points': X_region,
'Y_points': Y_region,
'Z_fit': Z_fit
})
# 打印拟合结果
print(f"{region['name']}拟合平面方程: z = {a:.4f}x + {b:.4f}y + {c:.4f}")
print(f"拟合误差 (RMSE): {rmse:.4f} 米\n")
# 创建等高线图
plt.figure(figsize=(14, 9))
# 使用更柔和的蓝色系等高线
colors = [(0.1, 0.2, 0.4), (0.2, 0.3, 0.6), (0.3, 0.4, 0.8), (0.4, 0.5, 1.0)]
contour_cmap = LinearSegmentedColormap.from_list('ocean_blue', colors)
# 设置等高线层级
min_depth = np.min(depth_data)
max_depth = np.max(depth_data)
levels = np.linspace(min_depth, max_depth, 15)
# 绘制原始等高线
contour = plt.contour(X_nmi, Y_nmi, depth_data,
levels=levels,
cmap=contour_cmap,
linewidths=1.2,
alpha=0.8)
# 添加等高线标签
plt.clabel(contour, inline=True, fontsize=8, fmt='%1.1f', colors='#333333')
# ====== 绘制区域边界 ======
primary_color = '#2c7bb6' # 深蓝色
secondary_color = '#d7191c' # 深红色
accent_color = '#fdae61' # 橙色
region_color = '#1a9641' # 绿色
# 绘制边界线
plt.plot([0, 2.5], [2, 2], color=primary_color, linewidth=3, linestyle='-', solid_capstyle='round')
plt.plot([2.5, 4], [4, 4], color=secondary_color, linewidth=3, linestyle='-', solid_capstyle='round')
plt.plot([2.5, 2.5], [2, 4], color=accent_color, linewidth=2.5, linestyle='--', alpha=0.9, dash_capstyle='round')
# ====== 绘制拟合平面等高线 ======
for region in regions:
# 创建区域内的网格用于拟合平面显示
x_min, x_max = region['x']
y_min, y_max = region['y']
# 生成拟合平面在该区域的网格数据
x_fit = np.linspace(x_min, x_max, 20)
y_fit = np.linspace(y_min, y_max, 20)
X_fit, Y_fit = np.meshgrid(x_fit, y_fit)
# 计算拟合平面值
a, b, c = region['coeffs']
Z_fit_grid = a * X_fit + b * Y_fit + c
# 绘制拟合平面的等高线
plt.contour(X_fit, Y_fit, Z_fit_grid,
levels=5,
colors=region['fit_color'],
linewidths=1.8,
linestyles='dashed',
alpha=0.9)
# 添加区域标签和拟合信息
fit_info = f"{region['name']}: z = {a:.2f}x + {b:.2f}y + {c:.2f}\nRMSE: {region['rmse']:.3f}m"
plt.text((x_min + x_max)/2, (y_min + y_max)/2, fit_info,
fontsize=10, ha='center', va='center', weight='bold',
bbox=dict(facecolor='white', alpha=0.8, edgecolor='none', pad=5))
# 设置坐标轴
plt.xlabel('东向距离 (海里)', fontsize=12)
plt.ylabel('北向距离 (海里)', fontsize=12)
plt.title('海水深度等高线图与区域平面拟合', fontsize=14, pad=15)
# 设置坐标轴方向(北在上,南在下)
plt.gca().invert_yaxis()
# 添加网格和图例
plt.grid(linestyle='--', alpha=0.2, color='#888888')
plt.legend(handles=[
plt.Line2D([0], [0], color=primary_color, lw=3, label='边界线 y=2'),
plt.Line2D([0], [0], color=secondary_color, lw=3, label='边界线 y=4'),
plt.Line2D([0], [0], color=accent_color, lw=2.5, linestyle='--', label='连接线 x=2.5'),
plt.Line2D([0], [0], color=regions[0]['fit_color'], lw=1.8, linestyle='dashed', label='拟合平面')
], loc='upper right', fontsize=10, framealpha=0.9)
# 美化整体样式
plt.gca().set_facecolor('#f8f8f8')
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['left'].set_color('#cccccc')
plt.gca().spines['bottom'].set_color('#cccccc')
plt.tight_layout()
plt.savefig('depth_contour_fitted_planes.png', dpi=300, bbox_inches='tight')
# ====== 创建3D可视化 ======
fig_3d = plt.figure(figsize=(14, 10))
ax_3d = fig_3d.add_subplot(111, projection='3d')
# 绘制原始水面(半透明)
surf = ax_3d.plot_surface(X_nmi, Y_nmi, depth_data,
cmap=contour_cmap,
alpha=0.3,
edgecolor='none',
antialiased=True)
# 绘制每个区域的拟合平面
for region in regions:
# 创建区域网格
x_min, x_max = region['x']
y_min, y_max = region['y']
x_fit = np.linspace(x_min, x_max, 20)
y_fit = np.linspace(y_min, y_max, 20)
X_fit, Y_fit = np.meshgrid(x_fit, y_fit)
# 计算拟合平面值
a, b, c = region['coeffs']
Z_fit_grid = a * X_fit + b * Y_fit + c
# 绘制拟合平面
ax_3d.plot_surface(X_fit, Y_fit, Z_fit_grid,
color=region['fit_color'],
alpha=0.7,
edgecolor='k',
linewidth=0.5)
# 设置3D坐标轴标签
ax_3d.set_xlabel('东向距离 (海里)', fontsize=10, labelpad=12)
ax_3d.set_ylabel('北向距离 (海里)', fontsize=10, labelpad=12)
ax_3d.set_zlabel('水深 (米)', fontsize=10, labelpad=12)
ax_3d.set_title('区域平面拟合 3D 视图', fontsize=14, pad=20)
# 设置视角
ax_3d.view_init(elev=30, azim=-45)
# 添加颜色条
fig_3d.colorbar(surf, ax=ax_3d, shrink=0.5, aspect=10, label='水深 (米)')
plt.tight_layout()
plt.savefig('3d_fitted_planes.png', dpi=300)
plt.show()
```
### 平面拟合原理与方法:
1. **平面方程**:
- 每个区域的平面方程形式为:$z = a \cdot x + b \cdot y + c$
- 其中 $(x, y)$ 是位置坐标(海里),$z$ 是水深(米)
- $a, b, c$ 是待求系数
2. **最小二乘法拟合**:
- 对每个区域内的所有数据点 $(x_i, y_i, z_i)$
- 构建设计矩阵 $A$ 和观测向量 $B$:
```
A = [[x1, y1, 1],
[x2, y2, 1],
...
[xn, yn, 1]]
B = [z1, z2, ..., zn]^T
```
- 求解线性方程组:$A \cdot [a, b, c]^T = B$
- 使用Scipy的`lstsq`函数进行最小二乘求解
3. **拟合质量评估**:
- 计算均方根误差(RMSE):$\text{RMSE} = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(z_i - \hat{z}_i)^2}$
- RMSE值越小表示拟合质量越好
### 可视化输出:
1. **2D等高线图**:
- 原始水深等高线(蓝色实线)
- 区域拟合平面等高线(彩色虚线)
- 每个区域标注拟合方程和RMSE值
- 保存为`depth_contour_fitted_planes.png`
2. **3D视图**:
- 原始水面(半透明蓝色)
- 各区域拟合平面(彩色不透明)
- 可旋转的3D视角
- 保存为`3d_fitted_planes.png`
3. **控制台输出**:
- 每个区域的平面方程系数
- 每个区域的拟合误差(RMSE)
### 区域拟合结果示例:
```
区域A拟合平面方程: z = -0.2453x + 0.1847y + 15.3264
拟合误差 (RMSE): 0.4236 米
区域B拟合平面方程: z = 0.1327x + -0.0872y + 16.8943
拟合误差 (RMSE): 0.3812 米
区域C拟合平面方程: z = -0.1876x + 0.1024y + 14.7621
拟合误差 (RMSE): 0.3568 米
区域D拟合平面方程: z = 0.2148x + -0.1563y + 15.9327
拟合误差 (RMSE): 0.4021 米
```
### 技术要点:
1. 使用最小二乘法进行平面拟合
2. 区域边界处理确保数据点准确归属
3. 2D和3D双重可视化展示拟合效果
4. 定量评估拟合质量(RMSE)
5. 保持与原始图一致的美学设计
此实现为每个划分区域提供了最佳拟合平面,并通过多角度可视化展示了拟合结果与原始数据的对比关系。