将整个海域进行分区,分为3*3共9个区域,并分别对每个区域进行最小二乘法拟合,完成图像输出以及相关数据输出(平面方程,拟合质量评估
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.spatial import KDTree
from tqdm import tqdm
from concurrent.futures import ProcessPoolExecutor
import matplotlib
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
#读取数据
def read_data(file_path):
df = pd.read_excel(file_path, header=None)
data = df.values
x_coords_nautical = data[0, 1:].astype(float)
x_coords_meter = x_coords_nautical * 1852
y_coords_meter = []
depth_data_meter = []
for i in range(1, len(data)):
y_coords_meter.append(float(data[i, 0]) * 1852)
depth_data_meter.append(data[i, 1:].astype(float))
return np.array(depth_data_meter), np.array(x_coords_meter), np.array(y_coords_meter)
#最小二乘法拟合
def least_squares_fit(x_coords, y_coords, depth_data):
X, Y = np.meshgrid(x_coords, y_coords)
A = np.column_stack((X.flatten()[:, np.newaxis],
Y.flatten()[:, np.newaxis],
np.ones((X.size, 1))))
B = depth_data.flatten()
coefficients, residuals, rank, s = np.linalg.lstsq(A, B, rcond=None)
a, b, c = coefficients
Z_fit = a * X + b * Y + c
residuals = depth_data - Z_fit
ssr = np.sum(residuals**2)
y_mean = np.mean(depth_data)
sst = np.sum((depth_data - y_mean)**2)
ssr_model = np.sum((Z_fit - y_mean)**2)
r_squared = 1 - ssr / sst
mse = ssr / (depth_data.size - 3)
return a, b, c, ssr, r_squared, mse
#评估
def evaluate_fit(depth_data, Z_fit):
residuals = depth_data - Z_fit
#残差平方和(SSR)
ssr = np.sum(residuals**2)
#总平方和(SST)
y_mean = np.mean(depth_data)
sst = np.sum((depth_data - y_mean)**2)
#回归平方和(SSR_model)
ssr_model = np.sum((Z_fit - y_mean)**2)
#决定系数R²
r_squared = 1 - ssr / sst
#调整后的决定系数R²_adj(考虑变量数量)
n = depth_data.size # 样本数量
p = 2 # 自变量数量
r_squared_adj = 1 - (1 - r_squared) * (n - 1) / (n - p - 1)
#均方误差(MSE)
mse = ssr / (n - p - 1) # 自由度调整
#平均绝对误差(MAE)
mae = np.mean(np.abs(residuals))
#平均绝对百分比误差(MAPE)
mape = np.mean(np.abs(residuals / depth_data)) * 100
return {
'SSR': ssr,
'SST': sst,
'SSR_model': ssr_model,
'R_squared': r_squared,
'R_squared_adj': r_squared_adj,
'MSE': mse,
'MAE': mae,
'MAPE': mape
}
#3d
def plot_seafloor(x_coords, y_coords, depth_data):
X, Y = np.meshgrid(x_coords, y_coords)
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Y, depth_data, cmap='terrain',
linewidth=0, antialiased=True)
fig.colorbar(surf, shrink=0.5, aspect=5, label='水深 (米)')
ax.set_xlabel('东西方向位置 (米)')
ax.set_ylabel('南北方向位置 (米)')
ax.set_zlabel('水深 (米)')
ax.set_title('三维海底地形图')
ax.view_init(elev=35, azim=-45)
plt.tight_layout()
plt.savefig("seafloor_terrain.png", dpi=300, bbox_inches='tight')
plt.show()
#2d
def plot_depth_contour(x_coords, y_coords, depth_data):
X, Y = np.meshgrid(x_coords, y_coords)
plt.figure(figsize=(12, 8))
contour = plt.contourf(X, Y, depth_data, levels=100, cmap='terrain')
plt.colorbar(contour, label='水深 (米)')
cline = plt.contour(X, Y, depth_data, levels=100, colors='black', linewidths=0.5)
plt.clabel(cline, inline=True, fontsize=8)
plt.xlabel('东西方向位置 (米)')
plt.ylabel('南北方向位置 (米)')
plt.title('二维水深等高线图')
plt.tight_layout()
plt.savefig("depth_contour.png", dpi=300, bbox_inches='tight')
plt.show()
#最小二乘法拟合3d
def plot_fitted_plane(x_coords, y_coords, depth_data, a, b, c):
X, Y = np.meshgrid(x_coords, y_coords)
Z = a * X + b * Y + c
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
surf_real = ax.plot_surface(X, Y, depth_data, cmap='terrain', linewidth=0, antialiased=True, alpha=0.7)
surf_fit = ax.plot_surface(X, Y, Z, cmap='viridis',
linewidth=0, antialiased=True, alpha=0.7)
fig.colorbar(surf_real, shrink=0.5, aspect=5, label='水深 (米)')
ax.set_xlabel('东西方向位置 (米)')
ax.set_ylabel('南北方向位置 (米)')
ax.set_zlabel('水深 (米)')
ax.set_title('海底地形与最小二乘法拟合平面')
ax.view_init(elev=35, azim=-45)
plt.tight_layout()
plt.savefig("fitted_plane.png", dpi=300, bbox_inches='tight')
plt.show()
#最小二乘法拟合2d
def plot_fitted_contour(x_coords, y_coords, depth_data, a, b, c):
X, Y = np.meshgrid(x_coords, y_coords)
Z = a * X + b * Y + c
plt.figure(figsize=(12, 8))
contour_real = plt.contourf(X, Y, depth_data, levels=50, cmap='terrain', alpha=0.7)
contour_fit = plt.contour(X, Y, Z, levels=20, cmap='viridis', linewidths=2)
plt.colorbar(contour_real, label='水深 (米)')
plt.clabel(contour_fit, inline=True, fontsize=8)
plt.xlabel('东西方向位置 (米)')
plt.ylabel('南北方向位置 (米)')
plt.title('原始海底地形与最小二乘法拟合等高线')
plt.tight_layout()
plt.savefig("fitted_contour.png", dpi=300, bbox_inches='tight')
plt.show()
#主程序
if __name__ == "__main__":
file_path = "C:/Users/Yeah/Desktop/数模/第七题/B题/附件(数据).xlsx"
depth_data, x_coords, y_coords = read_data(file_path)
#实际地形图
plot_seafloor(x_coords, y_coords, depth_data)
plot_depth_contour(x_coords, y_coords, depth_data)
#最小二乘法拟合
a, b, c, ssr, r_squared, mse = least_squares_fit(x_coords, y_coords, depth_data)
print(f"\n最小二乘法拟合结果:")
print(f"平面方程: z = {a:.4f}x + {b:.4f}y + {c:.4f}")
print(f"其中x为东西方向,y为南北方向,z为水深")
#评估指标
print(f"\n拟合质量评估:")
print(f"残差平方和 (SSR): {ssr:.4f}")
print(f"决定系数 (R²): {r_squared:.4f}")
print(f"均方误差 (MSE): {mse:.4f}")
#拟合图像
plot_fitted_plane(x_coords, y_coords, depth_data, a, b, c)
plot_fitted_contour(x_coords, y_coords, depth_data, a, b, c)