激光测量v3.1，使用固定速度拟合圆，每个叶片与对应的圆心分别求桨距角，显示旋转修正后叶片轮廓，排查若干BUG

data_analyse_origin.py

@@ -124,7 +124,7 @@ def data_analyse(path: List[str]):
     rotate_angle = float(path[5])  # 偏航角
     skew_angle_meas = 0  # 偏航角
     noise_reduction = 0.000001  # 如果一个距离值的所有样本量小于总样本量的noise_reduction，则被去掉
-    min_difference = 2  # 如果相邻2个点的距离差大于min_difference，则被注意是否是周期节点
+    min_difference = 1.5  # 如果相邻2个点的距离差大于min_difference，则被注意是否是周期节点
     group_length = [10000, 10000, 5000, 10000]  # 计算叶片轮廓时每个小切片的长度，4个数分别为叶中、叶根、叶尖、轴线计算切片长度
     lift_up_limit = 0  #  如果法兰中的数据大于这个值，则提升5m
     return_list = []
@@ -235,22 +235,24 @@ def data_analyse(path: List[str]):
 
     flange_root_dist = np.sqrt(flange_ava ** 2 + root_ava ** 2 - 2 * flange_ava * root_ava * np.cos(np.deg2rad(angle_flange_new - angle_root)))
 
-    blade_axis = blade_axis_cal(filtered_data_flange, start_flange, end_flange,
-                                   angle_flange + angle_cone + axial_inclination, group_length[3])
-    blade_axis_new, angle_flange_new = flange_coordinate_normalize(blade_axis, angle_flange)
-
     # 对叶片的边界点表做半径计算
     flange_r = radius_cal(border_rows_flange_new, angle_flange_new, dist_cen, angle_cen, axial_inclination, angle_cone)
     root_r = radius_cal(border_rows_root, angle_root, dist_cen, angle_cen, axial_inclination, angle_cone)
     nan_r = radius_cal(border_rows_nan_new, angle_nan_new, dist_cen, angle_cen, axial_inclination, angle_cone)
 
+    blade_axis, blade_axis_tuple, result_line_flange = blade_axis_cal(filtered_data_flange, start_flange, end_flange,
+                                   angle_flange + angle_cone + axial_inclination, group_length[3], flange_r)
+
+    blade_axis_new, angle_flange_new = flange_coordinate_normalize(blade_axis, angle_flange)
+    blade_axis_tuple_new, angle_flange_new = flange_coordinate_normalize(blade_axis_tuple, angle_flange)
+
     if np.abs((root_r - flange_r) - d_radius) > 0.5:
         print(str(root_r - flange_r) + "对比" + str(d_radius))
-        raise ValueError("Radius err1.")
+        root_r = flange_r + d_radius
 
     if np.abs(flange_root_dist - d_radius) > 0.5:
         print(str(flange_root_dist) + "对比" + str(d_radius))
-        raise ValueError("Radius err2.")
+        root_r = flange_r + flange_root_dist
 
     # root_r = d_radius + blade_axis_new["旋转半径"].iloc[0]
     # root_flange_dist = d_radius / np.cos(np.deg2rad(angle_cone))
@@ -259,38 +261,30 @@ def data_analyse(path: List[str]):
     #     * np.cos(np.deg2rad(90 - angle_cone - axial_inclination - angle_flange_new)))
 
     blade_axis_new["中心y"] = blade_axis_new["中心y"] - (flange_ava - root_ava)
+    blade_axis_tuple_new["中心y"] = blade_axis_tuple_new["中心y"] - (flange_ava - root_ava)
 
     # 计算叶片测量位置处的绝对桨距角、相对桨距角、线速度、叶片内部中心点距离
     aero_dist_flange, v_speed_flange, cen_blade_flange = (
         blade_angle_aero_dist(border_rows_flange, flange_r, cycle_len_flange, tower_dist_flange, angle_flange_new))
-    aero_dist_root, v_speed_root, cen_blade_root = (
-        blade_angle_aero_dist(border_rows_root, root_r, cycle_len_root, tower_dist_root, angle_root))
     aero_dist_nan, v_speed_nan, cen_blade_nan = (
         blade_angle_aero_dist(border_rows_nan_new, nan_r, cycle_len_nan, tower_dist_nan, angle_nan_new))
-    pitch_angle_root, v_speed_root = (
-        blade_angle(border_rows_root, blade_axis_new, root_r, cycle_len_root, angle_root + axial_inclination))
+    pitch_angle_root, v_speed_root, blade_axis_tuple_new = (
+        blade_angle(border_rows_root, blade_axis_tuple_new, root_r, cycle_len_root, angle_root + angle_cone + axial_inclination))
 
-    blade_axis_new["中心y"] = blade_axis_new["中心y"]*np.cos(np.deg2rad(angle_root + angle_cone + axial_inclination))
+    blade_axis_tuple_new["中心y"] = blade_axis_tuple_new["中心y"]*np.cos(np.deg2rad(angle_root + angle_cone + axial_inclination))
 
     # 将列表转换为 numpy 数组
-    cen_blade_flange_array = np.array(cen_blade_flange)
     cen_blade_nan_array = np.array(cen_blade_nan)
-    min_flange_array = np.array(min_flange)
     min_nan_array = np.array(min_nan)
 
     # 计算叶片内部中心点距离与叶片最小值之间的差值
-    abs_diff = np.abs(cen_blade_flange_array - min_flange_array)
     abs_diff_nan = np.abs(cen_blade_nan_array - min_nan_array)
-    blade_dist_flange = abs_diff * np.cos(np.deg2rad(angle_flange_new))
     blade_dist_nan = abs_diff_nan * np.cos(np.deg2rad(angle_nan_new))
-    blade_dist_flange.tolist()
     blade_dist_nan.tolist()
 
     # 计算叶片转速-净空散点表
-    dist_distribute = blade_dist_distribute_cal(filtered_data_flange, start_flange, end_flange,
-                                                tower_dist_flange, angle_flange_new, blade_dist_flange)
     dist_distribute_nan = blade_dist_distribute_cal(filtered_data_nan, start_nan, end_nan,
-                                                tower_dist_nan, angle_nan_new, blade_dist_nan)
+                                                tower_dist_nan, angle_nan_new + angle_cone + axial_inclination, blade_dist_nan)
     # dist_distribute = [df.round(5) for df in dist_distribute]
     dist_distribute = [df.round(5) for df in dist_distribute_nan]
 
@@ -314,20 +308,25 @@ def data_analyse(path: List[str]):
         mean_values.append(round((max_values[i] + min_values[i]) / 2, 2))
 
     # 将叶片线表数据乘以线速度，和俯仰角，得到叶片横截面的真实轮廓
-    for df in result_line_flange:
-        first_column = df.iloc[:, 0]
-        sec_column = df.iloc[:, 1]
-        df.iloc[:, 0] = first_column * v_speed_flange
-        min_time = df.iloc[:, 0].min()
-        df.iloc[:, 0] -= min_time
-        df.iloc[:, 1] = sec_column * np.cos(np.deg2rad(angle_flange_new + angle_cone + axial_inclination))
+    # for i in range(3):
+    #     df = result_line_flange[i]
+    #     first_column = df.iloc[:, 0]
+    #     sec_column = df.iloc[:, 1]
+    #     df.iloc[:, 0] = first_column * v_speed_flange
+    #     blade_axis_tuple.iloc[i, 1] = blade_axis_tuple.iloc[i, 1] * v_speed_flange * 5000000
+    #     min_time = df.iloc[:, 0].min()
+    #     df.iloc[:, 0] -= min_time
+    #     blade_axis_tuple.iloc[i, 1] -= min_time
+
 
-    for df in result_line_root:
+    for i in range(3):
+        df = result_line_root[i]
         first_column = df.iloc[:, 0]
         sec_column = df.iloc[:, 1]
         df.iloc[:, 0] = first_column * v_speed_root
         min_time = df.iloc[:, 0].min()
         df.iloc[:, 0] -= min_time
+        blade_axis_tuple_new.iloc[i, 1] -= min_time
         df.iloc[:, 1] = sec_column * np.cos(np.deg2rad(angle_root + angle_cone + axial_inclination))
 
     for df in result_scatter_flange:
@@ -363,6 +362,45 @@ def data_analyse(path: List[str]):
     data_root.iloc[:, 0] = data_root.iloc[:, 0] / 5000000
     lowpass_data.iloc[:, 0] = lowpass_data.iloc[:, 0] / 5000000
 
+    # 将叶片轮廓旋转重合
+    rotated_root = [pd.DataFrame() for _ in range(3)]
+    for i in range(3):
+        angle_rad = np.deg2rad(-pitch_angle_root[i])
+        rotation_matrix = np.array([
+            [np.cos(angle_rad), -np.sin(angle_rad)],
+            [np.sin(angle_rad), np.cos(angle_rad)]
+        ])
+
+        # 获取中心点坐标
+        center_x = blade_axis_tuple_new.iloc[i, 1]
+        center_y = blade_axis_tuple_new.iloc[i, 2]
+
+        # 创建结果DataFrame
+        rotated_points = result_line_root[i].copy()
+
+        # 对每个点进行旋转
+        for idx, row in result_line_root[i].iterrows():
+            # 获取当前点坐标
+            x = row.iloc[0]  # 获取第一列作为x坐标
+            y = row.iloc[1]  # 获取第二列作为y坐标
+
+            # 将点平移到原点
+            translated_x = x - center_x
+            translated_y = y - center_y
+
+            # 应用旋转矩阵
+            rotated = rotation_matrix @ np.array([translated_x, translated_y])
+
+            # 将点平移回原位置
+            final_x = rotated[0] + center_x
+            final_y = rotated[1] + center_y
+
+            # 更新结果DataFrame
+            rotated_points.iloc[idx, 0] = final_x
+            rotated_points.iloc[idx, 1] = final_y
+
+        rotated_root[i % 3] = pd.concat([rotated_root[i % 3], rotated_points])
+
     # 将需要保存到CSV的数据添加到return_list中
     return_list.append(str(time_code))
     return_list.append(str(wind_name))
@@ -436,6 +474,10 @@ def data_analyse(path: List[str]):
             'avg_blade': {
                 'xdata': result_avg_flange.iloc[:, 0].tolist(),
                 'ydata': result_avg_flange.iloc[:, 1].tolist()
+            },
+            'blade_center': {
+                'xdata': blade_axis_tuple.iloc[:, 1].tolist(),
+                'ydata': blade_axis_tuple.iloc[:, 2].tolist()
             }
         },
         'blade_root': {
@@ -454,6 +496,22 @@ def data_analyse(path: List[str]):
             'avg_blade': {
                 'xdata': result_avg_root.iloc[:, 0].tolist(),
                 'ydata': result_avg_root.iloc[:, 1].tolist()
+            },
+            'first_rotate_blade': {
+                'xdata': rotated_root[0].iloc[:, 0].tolist(),
+                'ydata': rotated_root[0].iloc[:, 1].tolist()
+            },
+            'second_rotate_blade': {
+                'xdata': rotated_root[1].iloc[:, 0].tolist(),
+                'ydata': rotated_root[1].iloc[:, 1].tolist()
+            },
+            'third_rotate_blade': {
+                'xdata': rotated_root[2].iloc[:, 0].tolist(),
+                'ydata': rotated_root[2].iloc[:, 1].tolist()
+            },
+            'blade_center': {
+                'xdata': blade_axis_tuple_new.iloc[:, 1].tolist(),
+                'ydata': blade_axis_tuple_new.iloc[:, 2].tolist()
             }
         },
         'dist_distribution': {
@@ -558,12 +616,13 @@ def data_analyse(path: List[str]):
     print('叶根原始数据采样时间长度' + str(data_root.iloc[-1, 0]))
     print('-' * 50)
 
-    # plot_data(result_line_flange, 'line', 'data1')
+    plot_data(result_line_flange, blade_axis_tuple, 'line', 'data1')
     # plot_data(result_diff_flange, 'line', 'data_diff_1')
     # plot_data(result_scatter_flange, 'scatter', 'data1')
-    plot_data(result_line_root, blade_axis_new, 'line', 'data2')
+    plot_data(result_line_root, blade_axis_tuple_new, 'line', 'data2')
+    plot_data(rotated_root, blade_axis_tuple_new, 'line', 'data3')
     # plot_data(result_diff_root, 'line', 'data_diff_2')
-    # plot_data(result_scatter_root, blade_axis_new,  'scatter', 'data2')
+    # plot_data(result_scatter_root, blade_axis_tuple_new,  'scatter', 'data2')
     # plot_data(dist_distribute, 'scatter', 'dist_distribute')
 
     return json_output
@@ -699,7 +758,7 @@ def tower_filter(data_group: pd.DataFrame, noise_threshold: float):
     # 选择频率最大的5个值
     top_5_distances = distance_counts.head(5).index
     mean_values = data_group[data_group['distance'].isin(top_5_distances)]['distance'].mean()
-    data_group.loc[(data_group['distance'] < mean_values-20) | (
+    data_group.loc[(data_group['distance'] < mean_values*0.9) | (
                 data_group['distance'] > mean_values*1.1), 'distance'] = np.nan
 
     nan_count = data_group['distance'].isna().sum()
@@ -733,17 +792,21 @@ def cycle_calculate(data_group: pd.DataFrame, noise_threshold: float, min_distan
     noise_distance_threshold = distance_counts[distance_counts < noise_threshold].index
     noise_indices = data_group[data_group['distance'].isin(noise_distance_threshold)].index
     data_group.loc[noise_indices, 'distance'] = np.nan
+    print(distance_counts.get(0, 0))
 
-    if 0 in distance_counts.nlargest(3).index:
-        pass
+    if distance_counts.get(0, 0) >= 0.1:
+        top_5_distances = distance_counts[distance_counts.index != 0].head(5).index
+        mean_values = data_group[data_group['distance'].isin(top_5_distances)]['distance'].mean()
+        data_group.loc[((data_group['distance'] > 0) & (data_group['distance'] < mean_values - 30)) |
+                       (data_group['distance'] > mean_values * 1.1), 'distance'] = np.nan
     else:
         # 选择频率最大的5个值
-        top_5_distances = distance_counts.head(5).index
+        top_5_distances = distance_counts[distance_counts.index != 0].head(5).index
         mean_values = data_group[data_group['distance'].isin(top_5_distances)]['distance'].mean()
         data_group.loc[(data_group['distance'] < mean_values-30) | (
                     data_group['distance'] > mean_values*1.1), 'distance'] = np.nan
-        print(f"中值是：{mean_values}")
 
+    print(f"中值是：{mean_values}")
     nan_count = data_group['distance'].isna().sum()
     all_count = data_group.shape[0]
     print(f"替换为NaN的distance异常值的数量是: {nan_count}, 总数量是: {all_count},"
@@ -830,7 +893,7 @@ def data_normalize(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
     time.sleep(1)
 
 
-    normalize_cycle = start_times[1] - start_times[0]
+    normalize_cycle = int(start_times[1] - start_times[0])
     full_cycle = int((start_times[2] - start_times[0]) * 3)
     turbines = [pd.DataFrame() for _ in range(3)]
 
@@ -844,16 +907,14 @@ def data_normalize(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
         # 根据当前起始时间点和结束时间点对数据进行分段
         segment = data_group[(data_group['time'] > start_time) & (data_group['time'] <= end_time)]
 
-        if segment is None:
-            pass
-        else:
+        if not segment.empty:
         # 周期归一化
             ratio = (end_time - start_time) / normalize_cycle
             segment.loc[:, 'time'] = (segment['time'] - start_time) / ratio
             # segment.loc[:, 'distance'] = ff.butter_lowpass_filter(segment['distance'], cutoff_low, fs)
 
             # 将结果添加到相应的 turbine 数据框中
-            turbines[i % 3] = pd.concat([turbines[i % 3], segment])
+            turbines[int(i / 2) % 3] = pd.concat([turbines[int(i / 2) % 3], segment])
 
 
     # 数据分组清洗、求平均
@@ -873,8 +934,18 @@ def data_normalize(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
         grey_start_index = int(len(turbine_sorted) * 0.1)
         grey_end_index = int(len(turbine_sorted) * 0.9)
         subset_grey = turbine_sorted[grey_start_index:grey_end_index]
+        turbine_sorted = ff.median_filter_correction(turbine_sorted, 2, 10)
         mean_grey = subset_grey['grey'].mean() * 0.8  # 0.8
-        turbine_sorted = turbine_sorted[turbine_sorted['grey'] > mean_grey]
+        # 计算分界点
+        n = len(turbine_sorted)
+        n_10 = int(0.1 * n)
+
+        # 创建筛选掩码
+        is_extreme = (turbine_sorted.index < n_10) | (turbine_sorted.index >= len(turbine_sorted) - n_10)
+        meets_condition = turbine_sorted['grey'] > mean_grey
+
+        # 筛选：如果是前10%或后10%，需要满足条件；如果是中间80%，直接保留
+        turbine_sorted = turbine_sorted[~is_extreme | (is_extreme & meets_condition)]
 
         # 找到time列的第一个值
         first_time = turbine_sorted['time'].iloc[0]
@@ -1080,10 +1151,11 @@ def flange_coordinate_normalize(flange_cen_row: pd.DataFrame, flange_angle):
     :return: 归一化后法兰数据，新俯仰角
     """
     flange_angle1 = np.deg2rad(flange_angle)
+    center_y_mean = flange_cen_row['中心y'].mean()
 
     # 计算新的俯仰角
-    flange_angle_cal0 = ((np.sin(flange_angle1) * flange_cen_row['中心y'] - 0.07608) /
-                         (np.cos(flange_angle1) * flange_cen_row['中心y']))
+    flange_angle_cal0 = ((np.sin(flange_angle1) * center_y_mean - 0.07608) /
+                         (np.cos(flange_angle1) * center_y_mean))
     flange_angle_cal = np.arctan(flange_angle_cal0)
 
     # 更新中心y列
@@ -1099,8 +1171,7 @@ def flange_coordinate_normalize(flange_cen_row: pd.DataFrame, flange_angle):
 
 
 def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_points: pd.DataFrame, horizon_angle: float,
-                   group_len: int) \
-        -> pd.DataFrame:
+                   group_len: int, radius_blade: float):
 
     """
     每个叶片周期计算旋转中心及半径
@@ -1109,14 +1180,15 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
     :param end_points: 叶片后缘突变点。
     :param horizon_angle: 叶根法兰水平角度。
     :param group_len: 每数据长度。
+    :param radius_blade: 半径。
     :return: 叶根法兰圆轮廓的圆心、测量点的旋转半径
     """
     print('正在进行叶片轴线计算......')
 
-    def fit_circle(df, v_bounds=(2, 10.0), top_k=5, prefilter=True):
+    def fit_circle(df, v_fixed, top_k=5, prefilter=True):
         """
+        修改版：速度v作为已知量，只优化圆心和半径
         Savitzky–Golay 平滑 + 差分进化 + least_squares 精修
-        Savitzky–Golay 平滑在每个局部窗口里用多项式拟合来代替原始点，从而减少噪声影响，同时保留曲线的整体趋势和细节
         """
 
         # ========== Savitzky–Golay 平滑 ==========
@@ -1134,13 +1206,19 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
         # ---------- 平滑 ----------
         d_smooth = smooth_savgol(d_raw, window_length=101, polyorder=3) if prefilter else d_raw
 
-        bounds = [v_bounds, (0, 5), (200, 201), (1, 3)]
+        # 计算x坐标（已知速度v）
+        x = v_fixed * t
+
+        # 新的参数边界：只优化圆心坐标(xc, yc)和半径R
+        bounds = [(min(x) - 5, max(x) + 5),  # xc范围：x值范围附近
+                  (min(d_smooth), max(d_smooth) + 10),  # yc范围：y值范围附近
+                  (0.5, 10)]  # R范围
 
         def residuals_sq(params):
-            v, xc, yc, R = params
-            if v <= 0 or R <= 0:
+            # 参数现在只有三个：xc, yc, R
+            xc, yc, R = params
+            if R <= 0:
                 return 1e6 * np.ones_like(t)
-            x = v * t
             return (x - xc) ** 2 + (d_smooth - yc) ** 2 - R ** 2
 
         def objective_mean_sq(params):
@@ -1174,8 +1252,8 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
             res = least_squares(
                 residuals_sq,
                 x0=cand,
-                bounds=([v_bounds[0], -np.inf, -np.inf, 1e-6],
-                        [v_bounds[1], np.inf, np.inf, np.inf]),
+                bounds=([-np.inf, -np.inf, 1e-6],  # 下界
+                        [np.inf, np.inf, np.inf]),  # 上界
                 method='trf',
                 loss='linear',
                 max_nfev=50000,
@@ -1183,17 +1261,16 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
                 ftol=1e-12,
                 gtol=1e-12
             )
-            v_opt, xc_opt, yc_opt, R_opt = res.x
-            x_all = v_opt * t
-            Ri_all = np.sqrt((x_all - xc_opt) ** 2 + (d_smooth - yc_opt) ** 2)
+            xc_opt, yc_opt, R_opt = res.x
+            Ri_all = np.sqrt((x - xc_opt) ** 2 + (d_smooth - yc_opt) ** 2)
             geo_rmse = np.sqrt(np.mean((Ri_all - R_opt) ** 2))
 
             if geo_rmse < best_rmse:
                 best_rmse = geo_rmse
-                best_result = [v_opt, xc_opt, yc_opt, R_opt, geo_rmse]
+                best_result = [v_fixed, xc_opt, yc_opt, R_opt, geo_rmse]
 
         result_df = pd.DataFrame([best_result],
-                                 columns=['旋转半径', '中心x', '中心y', '圆半径', '几何RMSE'])  # 旋转半径本身为测量点线速度
+                                 columns=['旋转半径', '中心x', '中心y', '圆半径', '几何RMSE'])
         return result_df
 
     def plot_circle_fit(df, fit_result):
@@ -1231,7 +1308,6 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
         plt.title('速度乘时间后的散点与拟合圆')
         plt.show()
 
-    group_len = group_len
     combined_df_sorted = pd.concat([start_points, end_points]).sort_values(by='time')
     # 检查排序后的数据从start开始，end结束
     if combined_df_sorted.iloc[0].equals(end_points.iloc[0]):
@@ -1246,6 +1322,7 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
     data_group['distance'] = data_group['distance'] * np.cos(np.deg2rad(horizon_angle))
     normalize_cycle = start_times[1] - start_times[0]
     full_cycle = int((start_times[2] - start_times[0]) * 3)
+    v_blade = 10000000 * np.pi * radius_blade / full_cycle
     angle_speed = (np.pi / full_cycle) * 5000000
     turbines = [pd.DataFrame() for _ in range(3)]
 
@@ -1262,7 +1339,7 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
         ratio = (end_time - start_time) / normalize_cycle
         segment.loc[:, 'time'] = (segment['time'] - start_time) / ratio
 
-        turbines[i % 3] = pd.concat([turbines[i % 3], segment])
+        turbines[int(i / 2) % 3] = pd.concat([turbines[int(i / 2) % 3], segment])
 
     turbines_processed = []
     turbines_scattered = []
@@ -1290,24 +1367,26 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
 
         process_df['time'] = process_df['time'] / 5000000
         lower_bound = process_df['time'].quantile(0.2)
-        upper_bound = process_df['time'].quantile(0.7)
+        upper_bound = process_df['time'].quantile(0.8)
         processed_df = process_df[(process_df['time'] >= lower_bound) & (process_df['time'] <= upper_bound)]
-        blade_cen_est = fit_circle(processed_df)
+        blade_cen_est = fit_circle(processed_df, v_blade)
 
         plot_circle_fit(processed_df, blade_cen_est)
+        processed_df['time'] = processed_df['time'] * v_blade
         turbines_processed.append(processed_df)
         turbines_scattered.append(turbine)
         result_df = pd.concat([result_df, blade_cen_est], ignore_index=True)
         if blade_cen_est['几何RMSE'].iloc[0] >= 0.1:
             raise ValueError("叶片几何误差过大")
 
-    result_df = result_df.mean(numeric_only=True).to_frame().T
+    result_df_mean = result_df.mean(numeric_only=True).to_frame().T
+    result_df_mean["中心y"] = result_df_mean["中心y"] / np.cos(np.deg2rad(horizon_angle))
     result_df["中心y"] = result_df["中心y"] / np.cos(np.deg2rad(horizon_angle))
     # turbines_mean["中心x"] = turbines_mean["中心x"] * 5000000
     # turbines_mean["圆半径"] = turbines_mean["圆半径"] / np.cos(np.deg2rad(horizon_angle))
     # turbines_mean['旋转半径'] = turbines_mean['旋转半径'] / angle_speed
     print(result_df)
-    return result_df
+    return result_df_mean, result_df, turbines_processed
 
 
 def radius_cal(border_rows, meas_angle, cen_dist, cen_angle, angle_main, angle_rotate):
@@ -1387,26 +1466,32 @@ def blade_angle(border_rows: List[pd.DataFrame], cen_data: pd.DataFrame, radius:
 
     mean_value = sum(values) / len(values) if values else float('nan')
 
-    if np.abs(cen_data['中心y'].iloc[0] - mean_value) > 0.3:
-        cen_data['中心y'].iloc[0] = mean_value
-        print('y_change')
-    if cen_data['中心x'].iloc[0] > 1.5:
-        cen_data['中心x'].iloc[0] = 1.5
-        print('x_change')
-    if cen_data['中心x'].iloc[0] < 0.75:
-        cen_data['中心x'].iloc[0] = 0.75
-        print('x_change')
-    print(cen_data['中心x'].iloc[0])
+    for i in [0, 1, 2]:
+        if np.abs(cen_data['中心y'].iloc[i] - mean_value) > 0.5:
+            print('原本：' + str(cen_data['中心y'].iloc[i]) + '标准：' + str(mean_value))
+            cen_data['中心y'].iloc[i] = mean_value
+            print('y_change')
+        if cen_data['中心x'].iloc[i] > 1.5:
+            cen_data['中心x'].iloc[i] = 1.5
+            print('x_change')
+        if cen_data['中心x'].iloc[i] < 0.75:
+            cen_data['中心x'].iloc[i] = 0.75
+            print('x_change')
+        print(cen_data['中心x'].iloc[i])
 
     pitch_angle_list = []
     for idx, turbine in enumerate(border_rows, start=1):
+        # 使用圆拟合的情况（注释掉的是备选方案）
+        diff_time = np.abs(cen_data['中心x'].iloc[idx - 1] - turbine.iloc[1, 0] * v_speed)
+        # diff_time = np.abs((turbine.iloc[0, 0] - turbine.iloc[1, 0]) * 0.66 * v_speed)
+
+        # 计算长度差，使用对应的cen_data行
+        diff_len = np.abs((cen_data['中心y'].iloc[idx - 1] - turbine.iloc[1, 1]) * np.cos(np.deg2rad(v_angle)))
 
-        # diff_time = np.abs(cen_data['中心x'].iloc[0] - turbine.iloc[1, 0] * v_speed)
-        # TODO 如果圆拟合成功，则使用上面代码，删除下面的
-        diff_time = np.abs((turbine.iloc[0, 0] - turbine.iloc[1, 0]) * 0.66 * v_speed)
-        diff_len = np.abs((cen_data['中心y'].iloc[0] - turbine.iloc[1, 1]) * np.cos(np.deg2rad(v_angle)))
+        # 计算桨距角
         pitch_angle = np.degrees(np.arctan(diff_len / diff_time))
-        print('单个叶片绝对桨距角' + str(pitch_angle))
+
+        print(f'第{idx}个叶片绝对桨距角: {pitch_angle}')
         pitch_angle_list.append(pitch_angle)
 
     pitch_mean = np.mean(pitch_angle_list)
@@ -1415,7 +1500,7 @@ def blade_angle(border_rows: List[pd.DataFrame], cen_data: pd.DataFrame, radius:
     pitch_angle_list = [round(num, 2) for num in pitch_angle_list]
     print("\033[34m叶片相对角度偏差：\033[0m", pitch_angle_list)
 
-    return pitch_angle_list, v_speed
+    return pitch_angle_list, v_speed, cen_data
 
 
 def plot_data(data, point, plot_type: str, data_name: str):
@@ -1443,7 +1528,7 @@ def plot_data(data, point, plot_type: str, data_name: str):
     else:
         raise ValueError("plot_type must be either 'line' or 'scatter'")
 
-    plt.scatter(point['中心x'].iloc[0], point['中心y'].iloc[0], color='black', s=1000, marker='o', label='测量点')
+    plt.scatter(point['中心x'], point['中心y'], color='black', s=1000, marker='o', label='测量点')
     axy = plt.gca()  # 获取当前坐标轴对象
     plt.grid(which='both', linewidth=2)  # 设置网格线宽度为2
     axy.xaxis.set_major_locator(MaxNLocator(nbins=200))  # 设置x轴主刻度的最大数量为10
@@ -1588,60 +1673,20 @@ if __name__ == "__main__":
     # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-zz_20251011123102_50_25.23_32.98.csv"
     # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-zz_20251011122457_50_32.24_31.02.csv"
     # 第七组右侧
-    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-r-z_20251011124902_50_31.19_34.81.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-r-z_20251011124620_50_24.85_32.90.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-r-z_20251011124036_50_32.24_31.14.csv"
-    # 21-l
     locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-r-z_20251011124902_50_31.19_34.81.csv"
-    locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_21-l_20251013154240_50_26.12_35.64.csv"
-    measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_21-l_20251013153649_50_34.96_32.84.csv"
-    # 10-170
+    locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-r-z_20251011124620_50_24.85_32.90.csv"
+    measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-r-z_20251011124036_50_32.24_31.14.csv"
+    # dq-8
     # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-r-z_20251011124902_50_31.19_34.81.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_10_20251012110528_50_23.73_38.14.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_10_20251012105925_50_36.92_34.99.csv"
-    # 17
-    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-rr_20251011114639_50_28.64_31.72.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试第3天/gy_17_20251014101746_50_17.86_33.37.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试第3天/gy_17_20251014101154_50_31.96_29.80.csv"
-    # 17-1
-    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-rr_20251011114639_50_28.64_31.72.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试第3天/gy_17-1_20251014103313_50_16.24_30.50.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试第3天/gy_17-1_20251014102718_50_29.27_27.02.csv"
-    # 21-292
-    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-rr_20251011114639_50_28.64_31.72.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_21_20251013145823_50_23.04_31.70.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_21_20251013145148_50_31.24_29.51.csv"
-    # 5-232
-    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-rr_20251011114639_50_28.64_31.72.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_5_20251012115014_50_18.06_29.47.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_5_20251012114417_50_28.88_26.92.csv"
-    # 5右侧-231.8
-    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-rr_20251011114639_50_28.64_31.72.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_5_20251012115014_50_18.06_29.47.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_5-1_20251012120929_50_28.9_26.98.csv"
-    # 5左侧
-    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-rr_20251011114639_50_28.64_31.72.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_5_20251012115014_50_18.06_29.47.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_5-l_20251012102212_50_26.60_24.66.csv"
-    # 1--198
-    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-rr_20251011114639_50_28.64_31.72.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_1-6_20251013114932_50_20.26_34.14.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_1-6_20251013114336_50_33.39_31.41.csv"
-    # 2--154
-    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-rr_20251011114639_50_28.64_31.72.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_2-2_20251013112737_50_26.43_42.10.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_2-2_20251013112018_50_40.99_38.86.csv"
-    # 7--198.2
-    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-rr_20251011114639_50_28.64_31.72.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_7_20251012123914_50_20.36_32.06.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_7_20251012123253_50_31.15_28.97.csv"
-    # 8-219
-    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-rr_20251011114639_50_28.64_31.72.csv"
-    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_8_20251013122350_50_14.79_26.51.csv"
-    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/测试/gy_8_20251013121806_50_26.02_24.14.csv"
+    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251023大庆/第三天/dh_8_20251109102949_50_17.22_27.45.csv"
+    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251023大庆/第三天/dh_8_20251109102255_50_26.18_24.06.csv"
+    # zb-9
+    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-r-z_20251011124902_50_31.19_34.81.csv"
+    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251121淄博/数据/zb_9_20251122110132_50_14.19_19.83.csv"
+    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251121淄博/数据/zb_9_20251122105534_50_19.37_17.73.csv"
 
     start_t = time.time()  # 记录开始时间
-    data_path = [locate_path0, locate_path, measure_path, 5, 5, 1.54]  # 偏斜测量数据、轮毂数据、叶根数据、锥角、轴向倾角、偏航角
+    data_path = [locate_path0, locate_path, measure_path, 3.5, 5, 1.54]  # 偏斜测量数据、轮毂数据、叶根数据、锥角、轴向倾角、偏航角
     list_1 = data_analyse(data_path)
     print(f"耗时: {time.time() - start_t:.2f} 秒")
 

data_clean.py

@@ -93,7 +93,7 @@ def data_analyse(path: List[str]):
     angle_cone = float(path[2])
     axial_inclination = float(path[3])
     lift_up_limit = float(path[4])
-    group_length = [10000, 20000, 5000]
+    group_length = [10000, 10000, 5000, 10000]
     return_list = []
 
     wind_name, turbine_code, time_code, sampling_fq, angle_nan, angle_cen = find_param(locate_file)
@@ -151,47 +151,43 @@ def data_analyse(path: List[str]):
 
     flange_root_dist = np.sqrt(flange_ava ** 2 + root_ava ** 2 - 2 * flange_ava * root_ava * np.cos(np.deg2rad(angle_flange_new - angle_root)))
 
-    blade_axis = blade_axis_cal(filtered_data_flange, start_flange, end_flange,
-                                   angle_flange + angle_cone + axial_inclination)
-    blade_axis_new, angle_flange_new = flange_coordinate_normalize(blade_axis, angle_flange)
-
     flange_r = radius_cal(border_rows_flange_new, angle_flange_new, dist_cen, angle_cen, axial_inclination, angle_cone)
     root_r = radius_cal(border_rows_root, angle_root, dist_cen, angle_cen, axial_inclination, angle_cone)
     nan_r = radius_cal(border_rows_nan_new, angle_nan_new, dist_cen, angle_cen, axial_inclination, angle_cone)
 
+    blade_axis, blade_axis_tuple, result_line_flange = blade_axis_cal(filtered_data_flange, start_flange, end_flange,
+                                   angle_flange + angle_cone + axial_inclination, group_length[3], flange_r)
+
+    blade_axis_new, angle_flange_new = flange_coordinate_normalize(blade_axis, angle_flange)
+    blade_axis_tuple_new, angle_flange_new = flange_coordinate_normalize(blade_axis_tuple, angle_flange)
+
     if np.abs((root_r - flange_r) - d_radius) > 0.5:
-        raise ValueError("Radius err1.")
+        root_r = flange_r + d_radius
 
     if np.abs(flange_root_dist - d_radius) > 0.5:
-        raise ValueError("Radius err2.")
+        root_r = flange_r + flange_root_dist
 
     blade_axis_new["中心y"] = blade_axis_new["中心y"] - (flange_ava - root_ava)
+    blade_axis_tuple_new["中心y"] = blade_axis_tuple_new["中心y"] - (flange_ava - root_ava)
 
     aero_dist_flange, v_speed_flange, cen_blade_flange = (
         blade_angle_aero_dist(border_rows_flange, flange_r, cycle_len_flange, tower_dist_flange, angle_flange_new))
-    aero_dist_root, v_speed_root, cen_blade_root = (
-        blade_angle_aero_dist(border_rows_root, root_r, cycle_len_root, tower_dist_root, angle_root))
     aero_dist_nan, v_speed_nan, cen_blade_nan = (
         blade_angle_aero_dist(border_rows_nan_new, nan_r, cycle_len_nan, tower_dist_nan, angle_nan_new))
-    pitch_angle_root, v_speed_root = (
-        blade_angle(border_rows_root, blade_axis_new, root_r, cycle_len_root, angle_root + axial_inclination))
+    pitch_angle_root, v_speed_root, blade_axis_tuple_new = (
+        blade_angle(border_rows_root, blade_axis_tuple_new, root_r, cycle_len_root, angle_root + angle_cone + axial_inclination))
+
+    blade_axis_tuple_new["中心y"] = blade_axis_tuple_new["中心y"]*np.cos(np.deg2rad(angle_root + angle_cone + axial_inclination))
 
-    blade_axis_new["中心y"] = blade_axis_new["中心y"]*np.cos(np.deg2rad(angle_root + angle_cone + axial_inclination))
-    
-    cen_blade_flange_array = np.array(cen_blade_flange)
     cen_blade_nan_array = np.array(cen_blade_nan)
-    min_flange_array = np.array(min_flange)
     min_nan_array = np.array(min_nan)
-    abs_diff = np.abs(cen_blade_flange_array - min_flange_array)
     abs_diff_nan = np.abs(cen_blade_nan_array - min_nan_array)
-    blade_dist_flange = abs_diff * np.cos(np.deg2rad(angle_flange_new))
     blade_dist_nan = abs_diff_nan * np.cos(np.deg2rad(angle_nan_new))
-    blade_dist_flange.tolist()
     blade_dist_nan.tolist()
 
 
     dist_distribute_nan = blade_dist_distribute_cal(filtered_data_nan, start_nan, end_nan,
-                                                tower_dist_nan, angle_nan_new, blade_dist_nan)
+                                                tower_dist_nan, angle_nan_new + angle_cone + axial_inclination, blade_dist_nan)
     dist_distribute = [df.round(5) for df in dist_distribute_nan]
 
     
@@ -213,20 +209,14 @@ def data_analyse(path: List[str]):
     for i in range(3):
         mean_values.append(round((max_values[i] + min_values[i]) / 2, 2))
 
-    for df in result_line_flange:
-        first_column = df.iloc[:, 0]
-        sec_column = df.iloc[:, 1]
-        df.iloc[:, 0] = first_column * v_speed_flange
-        min_time = df.iloc[:, 0].min()
-        df.iloc[:, 0] -= min_time
-        df.iloc[:, 1] = sec_column * np.cos(np.deg2rad(angle_flange_new + angle_cone + axial_inclination))
-
-    for df in result_line_root:
+    for i in range(3):
+        df = result_line_root[i]
         first_column = df.iloc[:, 0]
         sec_column = df.iloc[:, 1]
         df.iloc[:, 0] = first_column * v_speed_root
         min_time = df.iloc[:, 0].min()
         df.iloc[:, 0] -= min_time
+        blade_axis_tuple_new.iloc[i, 1] -= min_time
         df.iloc[:, 1] = sec_column * np.cos(np.deg2rad(angle_root + angle_cone + axial_inclination))
 
     avg_flange = result_avg_flange.iloc[:, 0]
@@ -246,6 +236,37 @@ def data_analyse(path: List[str]):
     lowpass_data.iloc[:, 0] = lowpass_data.iloc[:, 0] / 5000000
 
 
+    rotated_root = [pd.DataFrame() for _ in range(3)]
+    for i in range(3):
+        angle_rad = np.deg2rad(-pitch_angle_root[i])
+        rotation_matrix = np.array([
+            [np.cos(angle_rad), -np.sin(angle_rad)],
+            [np.sin(angle_rad), np.cos(angle_rad)]
+        ])
+
+        center_x = blade_axis_tuple_new.iloc[i, 1]
+        center_y = blade_axis_tuple_new.iloc[i, 2]
+
+        rotated_points = result_line_root[i].copy()
+
+        for idx, row in result_line_root[i].iterrows():
+            x = row.iloc[0]
+            y = row.iloc[1]
+
+            translated_x = x - center_x
+            translated_y = y - center_y
+
+            rotated = rotation_matrix @ np.array([translated_x, translated_y])
+
+            final_x = rotated[0] + center_x
+            final_y = rotated[1] + center_y
+
+            rotated_points.iloc[idx, 0] = final_x
+            rotated_points.iloc[idx, 1] = final_y
+
+        rotated_root[i % 3] = pd.concat([rotated_root[i % 3], rotated_points])
+
+
     return_list.append(str(time_code))
     return_list.append(str(wind_name))
     return_list.append(str(turbine_code))
@@ -318,6 +339,10 @@ def data_analyse(path: List[str]):
             'avg_blade': {
                 'xdata': result_avg_flange.iloc[:, 0].tolist(),
                 'ydata': result_avg_flange.iloc[:, 1].tolist()
+            },
+            'blade_center': {
+                'xdata': blade_axis_tuple.iloc[:, 1].tolist(),
+                'ydata': blade_axis_tuple.iloc[:, 2].tolist()
             }
         },
         'blade_root': {
@@ -336,6 +361,22 @@ def data_analyse(path: List[str]):
             'avg_blade': {
                 'xdata': result_avg_root.iloc[:, 0].tolist(),
                 'ydata': result_avg_root.iloc[:, 1].tolist()
+            },
+            'first_rotate_blade': {
+                'xdata': rotated_root[0].iloc[:, 0].tolist(),
+                'ydata': rotated_root[0].iloc[:, 1].tolist()
+            },
+            'second_rotate_blade': {
+                'xdata': rotated_root[1].iloc[:, 0].tolist(),
+                'ydata': rotated_root[1].iloc[:, 1].tolist()
+            },
+            'third_rotate_blade': {
+                'xdata': rotated_root[2].iloc[:, 0].tolist(),
+                'ydata': rotated_root[2].iloc[:, 1].tolist()
+            },
+            'blade_center': {
+                'xdata': blade_axis_tuple_new.iloc[:, 1].tolist(),
+                'ydata': blade_axis_tuple_new.iloc[:, 2].tolist()
             }
         },
         'dist_distribution': {
@@ -475,7 +516,7 @@ def tower_filter(data_group: pd.DataFrame, noise_threshold: float):
     
     top_5_distances = distance_counts.head(5).index
     mean_values = data_group[data_group['distance'].isin(top_5_distances)]['distance'].mean()
-    data_group.loc[(data_group['distance'] < mean_values - 20) | (
+    data_group.loc[(data_group['distance'] < mean_values * 0.9) | (
             data_group['distance'] > mean_values * 1.1), 'distance'] = np.nan
 
     
@@ -493,14 +534,16 @@ def cycle_calculate(data_group: pd.DataFrame, noise_threshold: float, min_distan
     noise_indices = data_group[data_group['distance'].isin(noise_distance_threshold)].index
     data_group.loc[noise_indices, 'distance'] = np.nan
 
-    if 0 in distance_counts.nlargest(3).index:
-        pass
+    if distance_counts.get(0, 0) >= 0.1:
+        top_5_distances = distance_counts[distance_counts.index != 0].head(5).index
+        mean_values = data_group[data_group['distance'].isin(top_5_distances)]['distance'].mean()
+        data_group.loc[((data_group['distance'] > 0) & (data_group['distance'] < mean_values - 30)) |
+                       (data_group['distance'] > mean_values * 1.1), 'distance'] = np.nan
     else:
-        # 选择频率最大的5个值
-        top_5_distances = distance_counts.head(5).index
+        top_5_distances = distance_counts[distance_counts.index != 0].head(5).index
         mean_values = data_group[data_group['distance'].isin(top_5_distances)]['distance'].mean()
-        data_group.loc[(data_group['distance'] < mean_values - 30) | (
-                data_group['distance'] > mean_values * 1.1), 'distance'] = np.nan
+        data_group.loc[(data_group['distance'] < mean_values-30) | (
+                    data_group['distance'] > mean_values*1.1), 'distance'] = np.nan
     
     data_group['distance'] = data_group['distance'].fillna(method='ffill')
     filtered_data = data_group
@@ -564,7 +607,7 @@ def data_normalize(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
     
     start_times = combined_df_sorted['time'].tolist()
 
-    normalize_cycle = start_times[1] - start_times[0]
+    normalize_cycle = int(start_times[1] - start_times[0])
     full_cycle = int((start_times[2] - start_times[0]) * 3)
     turbines = [pd.DataFrame() for _ in range(3)]
 
@@ -578,15 +621,13 @@ def data_normalize(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
         
         segment = data_group[(data_group['time'] > start_time) & (data_group['time'] <= end_time)]
 
-        if segment is None:
-            pass
-        else:
+        if not segment.empty:
             
             ratio = (end_time - start_time) / normalize_cycle
             segment.loc[:, 'time'] = (segment['time'] - start_time) / ratio
 
             
-            turbines[i % 3] = pd.concat([turbines[i % 3], segment])
+            turbines[int(i / 2) % 3] = pd.concat([turbines[int(i / 2) % 3], segment])
 
     
     turbines_processed = []
@@ -602,8 +643,13 @@ def data_normalize(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
         grey_start_index = int(len(turbine_sorted) * 0.1)
         grey_end_index = int(len(turbine_sorted) * 0.9)
         subset_grey = turbine_sorted[grey_start_index:grey_end_index]
+        turbine_sorted = ff.median_filter_correction(turbine_sorted, 2, 10)
         mean_grey = subset_grey['grey'].mean() * 0.8
-        turbine_sorted = turbine_sorted[turbine_sorted['grey'] > mean_grey]
+        n = len(turbine_sorted)
+        n_10 = int(0.1 * n)
+        is_extreme = (turbine_sorted.index < n_10) | (turbine_sorted.index >= len(turbine_sorted) - n_10)
+        meets_condition = turbine_sorted['grey'] > mean_grey
+        turbine_sorted = turbine_sorted[~is_extreme | (is_extreme & meets_condition)]
         
         first_time = turbine_sorted['time'].iloc[0]
 
@@ -751,8 +797,10 @@ def flange_coordinate_normalize(flange_cen_row: pd.DataFrame, flange_angle):
 
     flange_angle1 = np.deg2rad(flange_angle)
 
-    flange_angle_cal0 = ((np.sin(flange_angle1) * flange_cen_row['中心y'] - 0.07608) /
-                         (np.cos(flange_angle1) * flange_cen_row['中心y']))
+    center_y_mean = flange_cen_row['中心y'].mean()
+
+    flange_angle_cal0 = ((np.sin(flange_angle1) * center_y_mean - 0.07608) /
+                         (np.cos(flange_angle1) * center_y_mean))
     flange_angle_cal = np.arctan(flange_angle_cal0)
 
     flange_cen_row['中心y'] = (flange_cen_row['中心y'] ** 2 + 0.0057881664 -
@@ -764,10 +812,10 @@ def flange_coordinate_normalize(flange_cen_row: pd.DataFrame, flange_angle):
     return flange_cen_row, flange_angle_new1
 
 
-def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_points: pd.DataFrame, horizon_angle: float) \
-        -> pd.DataFrame:
+def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_points: pd.DataFrame, horizon_angle: float,
+                   group_len: int, radius_blade: float):
 
-    def fit_circle(df, v_bounds=(2, 10), top_k=5, prefilter=True):
+    def fit_circle(df, v_fixed, top_k=5, prefilter=True):
 
         def smooth_savgol(y, window_length=101, polyorder=3):
             wl = min(window_length, len(y) if len(y) % 2 == 1 else len(y) - 1)
@@ -780,23 +828,23 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
         t = np.asarray(df['time'])
         d_raw = np.asarray(df['distance'])
 
-        # ---------- 平滑 ----------
         d_smooth = smooth_savgol(d_raw, window_length=101, polyorder=3) if prefilter else d_raw
 
-        bounds = [v_bounds, (0, 5), (200, 201), (1, 3)]
+        x = v_fixed * t
+        bounds = [(min(x) - 5, max(x) + 5),
+                  (min(d_smooth), max(d_smooth) + 10),
+                  (0.5, 10)]
 
         def residuals_sq(params):
-            v, xc, yc, R = params
-            if v <= 0 or R <= 0:
+            xc, yc, R = params
+            if R <= 0:
                 return 1e6 * np.ones_like(t)
-            x = v * t
             return (x - xc) ** 2 + (d_smooth - yc) ** 2 - R ** 2
 
         def objective_mean_sq(params):
             res = residuals_sq(params)
             return np.mean(res ** 2)
 
-        # ---------- 差分进化 ----------
         result = differential_evolution(
             objective_mean_sq,
             bounds,
@@ -810,7 +858,6 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
             workers=1
         )
 
-        # 多候选点精修
         pop = result.population
         energies = result.population_energies
         idx = np.argsort(energies)[:top_k]
@@ -823,8 +870,8 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
             res = least_squares(
                 residuals_sq,
                 x0=cand,
-                bounds=([v_bounds[0], -np.inf, -np.inf, 1e-6],
-                        [v_bounds[1], np.inf, np.inf, np.inf]),
+                bounds=([-np.inf, -np.inf, 1e-6],
+                        [np.inf, np.inf, np.inf]),
                 method='trf',
                 loss='linear',
                 max_nfev=50000,
@@ -832,34 +879,31 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
                 ftol=1e-12,
                 gtol=1e-12
             )
-            v_opt, xc_opt, yc_opt, R_opt = res.x
-            x_all = v_opt * t
-            Ri_all = np.sqrt((x_all - xc_opt) ** 2 + (d_smooth - yc_opt) ** 2)
+            xc_opt, yc_opt, R_opt = res.x
+            Ri_all = np.sqrt((x - xc_opt) ** 2 + (d_smooth - yc_opt) ** 2)
             geo_rmse = np.sqrt(np.mean((Ri_all - R_opt) ** 2))
 
             if geo_rmse < best_rmse:
                 best_rmse = geo_rmse
-                best_result = [v_opt, xc_opt, yc_opt, R_opt, geo_rmse]
+                best_result = [v_fixed, xc_opt, yc_opt, R_opt, geo_rmse]
 
         result_df = pd.DataFrame([best_result],
-                                 columns=['旋转半径', '中心x', '中心y', '圆半径', '几何RMSE'])  # 旋转半径本身为测量点线速度
+                                 columns=['旋转半径', '中心x', '中心y', '圆半径', '几何RMSE'])
         return result_df
 
-    group_len = 10000
     combined_df_sorted = pd.concat([start_points, end_points]).sort_values(by='time')
-    # 检查排序后的数据从start开始，end结束
     if combined_df_sorted.iloc[0].equals(end_points.iloc[0]):
         combined_df_sorted = combined_df_sorted.iloc[1:]
     if combined_df_sorted.iloc[-1].equals(start_points.iloc[-1]):
         combined_df_sorted = combined_df_sorted.iloc[:-1]
     combined_df_sorted.reset_index(drop=True, inplace=True)
 
-    # 将 start_points 中的时间点转换为列表
     start_times = combined_df_sorted['time'].tolist()
 
     data_group['distance'] = data_group['distance'] * np.cos(np.deg2rad(horizon_angle))
     normalize_cycle = start_times[1] - start_times[0]
     full_cycle = int((start_times[2] - start_times[0]) * 3)
+    v_blade = 10000000 * np.pi * radius_blade / full_cycle
     angle_speed = (np.pi / full_cycle) * 5000000
     turbines = [pd.DataFrame() for _ in range(3)]
 
@@ -876,21 +920,17 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
         ratio = (end_time - start_time) / normalize_cycle
         segment.loc[:, 'time'] = (segment['time'] - start_time) / ratio
 
-        turbines[i % 3] = pd.concat([turbines[i % 3], segment])
+        turbines[int(i / 2) % 3] = pd.concat([turbines[int(i / 2) % 3], segment])
 
     turbines_processed = []
     turbines_scattered = []
-    turbines_result = []
     result_df = pd.DataFrame()
     time_list = list(range(0, normalize_cycle, group_len))
     for turbine in turbines:
-        # 按时间排序
         turbine_sorted = turbine.sort_values(by='time').reset_index(drop=True)
 
-        # 找到time列的第一个值
         first_time = turbine_sorted['time'].iloc[0]
 
-        # 分组，时间列每1000为一组（每40个时间点一组）
         bins = list(range(int(first_time), int(turbine_sorted['time'].max()), group_len))
         grouped = turbine_sorted.groupby(pd.cut(turbine_sorted['time'], bins=bins, right=False))
 
@@ -898,27 +938,28 @@ def blade_axis_cal(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
         for index, (bin, group) in enumerate(grouped):
 
             mid_point = group.mean()
-            # 将中点转换为 DataFrame 并添加到处理后的 DataFrame 中
+
             mid_point_df = pd.DataFrame([mid_point])
             mid_point_df.iloc[0, 0] = time_list[index]
             process_df = pd.concat([process_df, mid_point_df], ignore_index=True)
 
         process_df['time'] = process_df['time'] / 5000000
         lower_bound = process_df['time'].quantile(0.2)
-        upper_bound = process_df['time'].quantile(0.7)
+        upper_bound = process_df['time'].quantile(0.8)
         processed_df = process_df[(process_df['time'] >= lower_bound) & (process_df['time'] <= upper_bound)]
-        blade_cen_est = fit_circle(processed_df)
-
+        blade_cen_est = fit_circle(processed_df, v_blade)
+        processed_df['time'] = processed_df['time'] * v_blade
         turbines_processed.append(processed_df)
         turbines_scattered.append(turbine)
         result_df = pd.concat([result_df, blade_cen_est], ignore_index=True)
         if blade_cen_est['几何RMSE'].iloc[0] >= 0.1:
             raise ValueError("叶片几何误差过大")
 
-    result_df = result_df.mean(numeric_only=True).to_frame().T
+    result_df_mean = result_df.mean(numeric_only=True).to_frame().T
+    result_df_mean["中心y"] = result_df_mean["中心y"] / np.cos(np.deg2rad(horizon_angle))
     result_df["中心y"] = result_df["中心y"] / np.cos(np.deg2rad(horizon_angle))
 
-    return result_df
+    return result_df_mean, result_df, turbines_processed
 
 
 
@@ -933,7 +974,7 @@ def radius_cal(border_rows, meas_angle, cen_dist, cen_angle, angle_main, angle_r
 def blade_angle_aero_dist(border_rows: List[pd.DataFrame], radius: float, full_cycle: int,
                                    tower_dist: float, v_angle: float):
 
-    v_speed = 2 * np.pi * radius / full_cycle  # 叶片线速度m/(1计时器单位）
+    v_speed = 2 * np.pi * radius / full_cycle
     aero_dist_list = []
     cen_blade = []
     for turbine in border_rows:
@@ -962,18 +1003,24 @@ def blade_angle(border_rows: List[pd.DataFrame], cen_data: pd.DataFrame, radius:
 
     mean_value = sum(values) / len(values) if values else float('nan')
 
-    if np.abs(cen_data['中心y'].iloc[0] - mean_value) > 0.3:
-        cen_data['中心y'].iloc[0] = mean_value
-    if cen_data['中心x'].iloc[0] > 1.5:
-        cen_data['中心x'].iloc[0] = 1.5
-    if cen_data['中心x'].iloc[0] < 0.75:
-        cen_data['中心x'].iloc[0] = 0.75
+    for i in [0, 1, 2]:
+        if np.abs(cen_data['中心y'].iloc[i] - mean_value) > 0.5:
+            print('原本：' + str(cen_data['中心y'].iloc[i]) + '标准：' + str(mean_value))
+            cen_data['中心y'].iloc[i] = mean_value
+            print('y_change')
+        if cen_data['中心x'].iloc[i] > 1.5:
+            cen_data['中心x'].iloc[i] = 1.5
+            print('x_change')
+        if cen_data['中心x'].iloc[i] < 0.75:
+            cen_data['中心x'].iloc[i] = 0.75
+            print('x_change')
+        print(cen_data['中心x'].iloc[i])
 
     pitch_angle_list = []
     for idx, turbine in enumerate(border_rows, start=1):
 
-        diff_time = np.abs((turbine.iloc[0, 0] - turbine.iloc[1, 0]) * 0.66 * v_speed)
-        diff_len = np.abs((cen_data['中心y'].iloc[0] - turbine.iloc[1, 1]) * np.cos(np.deg2rad(v_angle)))
+        diff_time = np.abs(cen_data['中心x'].iloc[idx - 1] - turbine.iloc[1, 0] * v_speed)
+        diff_len = np.abs((cen_data['中心y'].iloc[idx - 1] - turbine.iloc[1, 1]) * np.cos(np.deg2rad(v_angle)))
         pitch_angle = np.degrees(np.arctan(diff_len / diff_time))
         pitch_angle_list.append(pitch_angle)
 
@@ -982,7 +1029,7 @@ def blade_angle(border_rows: List[pd.DataFrame], cen_data: pd.DataFrame, radius:
     pitch_angle_list.append(max(pitch_angle_list) - min(pitch_angle_list))
     pitch_angle_list = [round(num, 2) for num in pitch_angle_list]
 
-    return pitch_angle_list, v_speed
+    return pitch_angle_list, v_speed, cen_data
 
 
 def find_param(path: str):

frequency_filter.py

@@ -103,7 +103,7 @@ def process_fft(data, fs):
     """
 
     cutoff_low = 0.001 * fs  # 设置低通滤波截止频率
-    segment = data.head(int(len(data) * 0.95))
+    segment = data.head(int(len(data) * 1))
 
     # 确保segment的长度是2的n次幂
     desired_length = 2**((len(segment) - 1).bit_length() - 1)
@@ -122,7 +122,7 @@ def process_fft(data, fs):
     fft_y_scaled = [0] + fft_y_scaled
     fft_x = xf_filter_truncated.tolist()
     fft_x = [0] + fft_x
-    fft_y_scaled = [0 if a_val < 0.1 else b_val for a_val, b_val in zip(fft_x, fft_y_scaled)]
+    fft_y_scaled = [0 if a_val < 0.04 else b_val for a_val, b_val in zip(fft_x, fft_y_scaled)]
     max_value = max(fft_y_scaled)
     max_index = fft_y_scaled.index(max_value)
 
@@ -299,3 +299,29 @@ def process_equations(list_param):
 
     return list(positive_angles)
 
+def median_filter_correction(df, column_index, window_size):
+    """
+    对指定列进行中值滤波修正
+
+    参数:
+    df: pandas DataFrame，包含数据的DataFrame
+    column_index: int，要处理的列索引（从0开始）
+    window_size: int，窗口大小
+
+    返回:
+    修正后的DataFrame
+    """
+    df_corrected = df.copy()
+    gray_values = df.iloc[:, column_index].values
+    corrected_values = gray_values.copy()
+
+    # 对每window_size个点计算中值，并赋值给第一个点
+    # 最后window_size-1个点不做处理
+    for i in range(len(gray_values) - window_size + 1):
+        window_data = gray_values[i:i + window_size]
+        median_value = np.median(window_data)
+        corrected_values[i] = median_value
+
+    df_corrected.iloc[:, column_index] = corrected_values
+
+    return df_corrected