检查
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
from matplotlib.patches import Polygon
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 build_coord_tree(x_coords, y_coords):
X, Y = np.meshgrid(x_coords, y_coords)
coords = np.column_stack((X.flatten(), Y.flatten()))
tree = KDTree(coords)
return tree, X, Y
# 查找最近点深度
def find_depth(x, y, tree, depth_data, X, Y):
_, idx = tree.query([x, y], k=1)
return depth_data.flatten()[idx]
# 计算覆盖宽度
def calculate_coverage_width(beam_angle, depth, slope_angle=0):
"""计算考虑坡度的多波束覆盖宽度"""
theta = np.radians(beam_angle/2)
alpha = np.radians(slope_angle)
# 基本覆盖宽度(平坦海底)
basic_width = 2 * depth * np.tan(theta)
# 考虑坡度的修正
if np.cos(theta - alpha) == 0:
return basic_width
# 坡度影响下的覆盖宽度
slope_effect_width = (depth / np.cos(alpha)) * (np.sin(theta) / np.cos(theta - alpha) +
np.sin(theta) / np.cos(theta + alpha))
return slope_effect_width
# 自适应测线布设
def adaptive_survey(x_coords, y_coords, depth_data, beam_angle=120,
overlap_range=(0.1, 0.2), base_spacing=None):
tree, X, Y = build_coord_tree(x_coords, y_coords)
a, b, c, _, _, _ = least_squares_fit(x_coords, y_coords, depth_data)
# 总等深线方向（即坡度方向）
slope_direction = np.arctan2(b, a) # 计算总等深线方向
slope_angle = np.degrees(slope_direction)
# 旋转坐标系，使X轴与等深线方向一致
cos_rot = np.cos(-slope_direction)
sin_rot = np.sin(-slope_direction)
# 存储测线点
all_lines = []
# 确定测线方向（垂直于等深线）
line_direction = slope_direction + np.pi/2
cos_line = np.cos(line_direction)
sin_line = np.sin(line_direction)
# 如果未指定基础测线间距，则根据平均水深计算
if base_spacing is None:
avg_depth = np.mean(depth_data)
avg_width = calculate_coverage_width(beam_angle, avg_depth)
base_spacing = avg_width * (1 - overlap_range[1]) # 使用最大重叠率计算间距
# 计算起始位置
x_center = (x_coords[0] + x_coords[-1]) / 2
y_center = (y_coords[0] + y_coords[-1]) / 2
# 计算需要的测线数量
total_range = np.sqrt((x_coords[-1]-x_coords[0])**2 +
(y_coords[-1]-y_coords[0])**2)
num_lines = int(np.ceil(total_range / base_spacing))
# 生成测线
for i in range(-num_lines//2, num_lines//2 + 1):
# 计算当前测线中心点（在等深线方向上偏移）
offset_distance = i * base_spacing
offset_x = offset_distance * np.cos(slope_direction)
offset_y = offset_distance * np.sin(slope_direction)
line_x_start = x_center + offset_x - cos_line * 10000 # 超出范围确保覆盖
line_y_start = y_center + offset_y - sin_line * 10000
line_x_end = x_center + offset_x + cos_line * 10000
line_y_end = y_center + offset_y + sin_line * 10000
# 在测线方向上采样点
num_points = 100
x_line = np.linspace(line_x_start, line_x_end, num_points)
y_line = np.linspace(line_y_start, line_y_end, num_points)
# 存储该测线上的有效点
valid_points = []
coverage_widths = []
depths = []
for x, y in zip(x_line, y_line):
if (x < x_coords[0] or x > x_coords[-1] or
y < y_coords[0] or y > y_coords[-1]):
continue
depth = find_depth(x, y, tree, depth_data, X, Y)
if depth <= 0:
continue
# 计算当前位置的坡度方向
grad_x, grad_y = np.gradient(depth_data)
local_slope_dir = np.arctan2(grad_y[int(y/100)][int(x/100)],
grad_x[int(y/100)][int(x/100)])
# 计算覆盖宽度
coverage_width = calculate_coverage_width(beam_angle, depth,
np.degrees(local_slope_dir))
valid_points.append((x, y))
coverage_widths.append(coverage_width)
depths.append(depth)
if valid_points:
all_lines.append({
'points': valid_points,
'coverage_widths': coverage_widths,
'depths': depths,
'offset': offset_distance
})
# 计算覆盖率和重叠率
coverage_map = build_coverage_map(x_coords, y_coords, all_lines)
return all_lines, coverage_map, slope_angle
# 构建覆盖地图
def build_coverage_map(x_coords, y_coords, all_lines):
coverage_map = np.zeros((len(y_coords), len(x_coords)))
for line in all_lines:
points = line['points']
coverage_widths = line['coverage_widths']
for (x, y), width in zip(points, coverage_widths):
# 找到最近的网格点
x_idx = np.argmin(np.abs(x_coords - x))
y_idx = np.argmin(np.abs(y_coords - y))
# 计算影响范围
half_width = width / 2
x_range = np.where(np.abs(x_coords - x) <= half_width)[0]
y_range = np.where(np.abs(y_coords - y) <= half_width)[0]
for i in y_range:
for j in x_range:
dx = x_coords[j] - x
dy = y_coords[i] - y
if np.sqrt(dx**2 + dy**2) <= width/2:
coverage_map[i, j] += 1
return coverage_map
# 绘制测线覆盖图
def plot_coverage(x_coords, y_coords, depth_data, all_lines, coverage_map):
plt.figure(figsize=(14, 10))
# 绘制地形
X, Y = np.meshgrid(x_coords, y_coords)
contour = plt.contourf(X, Y, depth_data, levels=50, cmap='terrain', alpha=0.7)
plt.colorbar(contour, label='水深 (米)')
# 绘制测线
for i, line in enumerate(all_lines):
points = np.array(line['points'])
plt.plot(points[:, 0], points[:, 1], '-', linewidth=2,
label=f'测线 {i+1}', alpha=0.8)
# 绘制覆盖不足区域
no_coverage = np.where(coverage_map == 0)
plt.scatter(x_coords[no_coverage[1]], y_coords[no_coverage[0]],
c='red', s=10, label='漏测区域', alpha=0.5)
plt.xlabel('东西方向位置 (米)')
plt.ylabel('南北方向位置 (米)')
plt.title('测线布设与覆盖效果')
plt.legend()
plt.tight_layout()
plt.savefig("survey_coverage.png", dpi=300, bbox_inches='tight')
plt.show()
# 计算评估指标
def calculate_metrics(all_lines, coverage_map, x_coords, y_coords, depth_data):
total_length = 0
for line in all_lines:
points = np.array(line['points'])
for i in range(1, len(points)):
total_length += np.sqrt((points[i][0] - points[i-1][0])**2 +
(points[i][1] - points[i-1][1])**2)
total_area = (x_coords[-1] - x_coords[0]) * (y_coords[-1] - y_coords[0])
no_coverage_area = np.sum(coverage_map == 0) / coverage_map.size * total_area
# 计算重叠过多区域
overlap_excess = np.where(coverage_map > 2)
overlap_excess_area = len(overlap_excess[0]) / coverage_map.size * total_area
# 计算平均重叠率
avg_overlap = np.mean(coverage_map[coverage_map > 0] - 1)
return {
'total_length': total_length,
'no_coverage_percent': no_coverage_area / total_area * 100,
'overlap_excess_area': overlap_excess_area,
'avg_overlap_rate': avg_overlap,
'num_lines': len(all_lines)
}
# 主程序
if __name__ == "__main__":
file_path = "C:/Users/Yeah/Desktop/数模/第七题/B题/附件（数据）.xlsx"
depth_data, x_coords, y_coords = read_data(file_path)
# 最小二乘法拟合
a, b, c, ssr, r_squared, mse = least_squares_fit(x_coords, y_coords, depth_data)
# 自适应测线布设
all_lines, coverage_map, slope_angle = adaptive_survey(x_coords, y_coords,
depth_data, beam_angle=120)
# 计算评估指标
metrics = calculate_metrics(all_lines, coverage_map, x_coords, y_coords, depth_data)
print(f"\n测线布设评估指标:")
print(f"总测线长度: {metrics['total_length']/1000:.2f} 公里")
print(f"漏测区域占比: {metrics['no_coverage_percent']:.2f}%")
print(f"重叠率超过20%的区域面积: {metrics['overlap_excess_area']/10000:.2f} 平方公里")
print(f"平均重叠率: {metrics['avg_overlap_rate']:.2f}")
print(f"总测线数量: {metrics['num_lines']}")
print(f"总等深线方向: {slope_angle:.2f} 度")
# 可视化覆盖效果
plot_coverage(x_coords, y_coords, depth_data, all_lines, coverage_map)
最新发布
经过几个月的努力，最终将代码行覆盖率从较低水平提升到了50%。此过程不仅涉及单元测试，还融合了集成测试以达到目标。后续将分享所使用的测试框架和技术细节。
扫一扫
628
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
