激光测量v3，允许小角度偏斜测量，计算法兰处拟合圆，通过旋转中心而非前缘确定角度偏差；添加偏斜角度计算但未使用；拟合圆半径过大的问题需解决

data_analyse_origin.py

@@ -11,13 +11,16 @@ import time
 import sys
 import frequency_filter as ff
 from datetime import datetime
+from scipy.optimize import least_squares, differential_evolution
+from scipy.signal import savgol_filter
+
 
 
 warnings.filterwarnings("ignore", category=FutureWarning) # 忽略特定警告
 plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
 plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
 
-
+# TODO: 将偏斜角度还原叶片数据的代码转移到最后的步骤
 
 def result_main():
 
@@ -107,33 +110,90 @@ def data_analyse(path: List[str]):
 
     """
     创建data目录，把分析数据保存到历史记录中，同时返回全量分析数据
+    locate_file_0：非中心线数据，1通道下方，2通道上方
+    locate_file：轮毂中心数据，1通道叶尖，2通道轮毂中心
+    measure_file：测量数据，1通道法兰，2通道叶片最宽处
     """
 
     # 基础配置参数
-    locate_file = path[0]
-    measure_file = path[1]
-    angle_cone = float(path[2])  # 锥角
-    axial_inclination = float(path[3])  # 轴向倾角
-    skew_angle = 5  # 偏航角
+    locate_file_0 = path[0]
+    locate_file = path[1]
+    measure_file = path[2]
+    angle_cone = float(path[3])  # 锥角
+    axial_inclination = float(path[4])  # 轴向倾角
+    rotate_angle = float(path[5])  # 偏航角
+    skew_angle_meas = 0  # 偏航角
     noise_reduction = 0.000001  # 如果一个距离值的所有样本量小于总样本量的noise_reduction，则被去掉
-    min_difference = 1  # 如果相邻2个点的距离差大于min_difference，则被注意是否是周期节点
-    group_length = [10000, 10000, 5000]  # 计算叶片轮廓时每个小切片的长度，三个数分别为叶中、叶根、叶尖切片长度
+    min_difference = 2  # 如果相邻2个点的距离差大于min_difference，则被注意是否是周期节点
+    group_length = [10000, 10000, 5000, 10000]  # 计算叶片轮廓时每个小切片的长度，4个数分别为叶中、叶根、叶尖、轴线计算切片长度
+    lift_up_limit = 0  #  如果法兰中的数据大于这个值，则提升5m
     return_list = []
 
     # 读取文件信息，包括风场名、风机编号、时间、采样频率、2个通道俯仰角
+    wind_name_0, turbine_code_0, time_code_0, sampling_fq_0, angle_locate_1, angle_locate_2 = find_param(locate_file_0)
     wind_name, turbine_code, time_code, sampling_fq, angle_nan, angle_cen = find_param(locate_file)
-    wind_name_1, turbine_code_1, time_code_1, sampling_fq_1, angle_tip, angle_root = find_param(measure_file)
+    wind_name_1, turbine_code_1, time_code_1, sampling_fq_1, angle_flange, angle_root = find_param(measure_file)
+    sampling_fq_0 = sampling_fq_0 * 1000
     sampling_fq_1 = sampling_fq_1 * 1000
     sampling_fq = sampling_fq * 1000
+    print(wind_name_0, turbine_code_0, time_code_0, sampling_fq_0, angle_locate_1, angle_locate_2)
     print(wind_name, turbine_code, time_code, sampling_fq, angle_nan, angle_cen)
-    print(wind_name_1, turbine_code_1, time_code, sampling_fq_1, angle_tip, angle_root)
+    print(wind_name_1, turbine_code_1, time_code, sampling_fq_1, angle_flange, angle_root)
 
     # 读取数据，并检查是否有时间序列异常，分离2通道数据
-    data_nan, data_cen = process_data(locate_file)
-    data_tip, data_root = process_data(measure_file)
+    data_locate1, data_locate2 = process_data(locate_file_0, 0)  # 定位数据，1是下方，2是上方
+    data_locate0, data_locate_nan = process_data(locate_file, 0)
+    data_nan, data_cen = process_data(locate_file, skew_angle_meas)
+    data_flange, data_root = process_data(measure_file, skew_angle_meas)
+
+    # 偏斜计算部分
+    locate1, data_loc1 = locate_filter(data_locate1)
+    locate2, data_loc2 = locate_filter(data_locate2)
+    locate0_filter = tower_filter(data_locate0, noise_reduction)
+    locate0 = locate0_filter['distance'].mean()
+    start_loc1, end_loc1, filtered_data_loc1 = cycle_calculate(data_loc1, noise_reduction, min_difference)
+    start_loc2, end_loc2, filtered_data_loc2 = cycle_calculate(data_loc2, noise_reduction, min_difference)
+    if np.abs(start_loc2.iloc[1, 1] - start_loc1.iloc[0, 1]) < np.abs(start_loc2.iloc[0, 1] - start_loc1.iloc[0, 1]):
+        start_loc2 = start_loc2.iloc[1:]
+        end_loc2 = end_loc2.iloc[1:]
+
+    angle_list = []
+    for i in range(2, min(len(start_loc2), len(start_loc1)) - 2):
+        time_period = start_loc2.iloc[i, 0] - start_loc1.iloc[i, 0]  # 2个测量点的时间差
+        time_angle1 = start_loc1.iloc[i + 1, 0] - start_loc1.iloc[i, 0]  # 以下方测量点计一个周期
+        time_angle2 = start_loc2.iloc[i, 0] - start_loc2.iloc[i - 1, 0]  # 以上方测量点计一个周期
+        if np.abs(time_angle1 - time_angle2) < 50000:
+            angle_diff =  360 * time_period / (time_angle1 * 3)
+            angle_list.append(angle_diff)
+        else:
+            pass
+    angle_avg = np.mean(angle_list)  # 计算2个测量点的平均相位角度
+    locate0_coordinate = [locate0 * np.cos(np.deg2rad(angle_cen)), 0, locate0 * np.sin(np.deg2rad(angle_cen))]
+    locate1_coordinate = [locate1 * np.cos(np.deg2rad(angle_locate_1)) * np.cos(np.deg2rad(rotate_angle)),
+                          locate1 * np.cos(np.deg2rad(angle_locate_1)) * np.sin(np.deg2rad(rotate_angle)),
+                          locate1 * np.sin(np.deg2rad(angle_locate_1))]
+    locate2_coordinate = [locate2 * np.cos(np.deg2rad(angle_locate_2)) * np.cos(np.deg2rad(rotate_angle)),
+                          locate2 * np.cos(np.deg2rad(angle_locate_2)) * np.sin(np.deg2rad(rotate_angle)),
+                          locate2 * np.sin(np.deg2rad(angle_locate_2))]
+    equation_list = [locate0_coordinate, locate1_coordinate, locate2_coordinate,
+                     axial_inclination, angle_avg, angle_cone]
+
+    print("\033[34m偏斜方程计算输入：\033[0m", equation_list)
+    # skew_angle = ff.process_equations(equation_list)
+
+    if lift_up_limit >= 0.1:
+        discrete_values = np.arange(0, 0.101, 0.001)
+        condition = data_flange['distance'] > lift_up_limit
+        n = condition.sum()
+        random_discrete = np.random.choice(discrete_values, size=n)
+        data_flange.loc[condition, 'distance'] = lift_up_limit + 3 + random_discrete
+    elif np.abs(lift_up_limit) < 0.1:
+        pass
+    else:
+        raise ValueError("lift_up_limit error.")
 
     # 全部数据进行降噪、去除异常点处理，叶根叶尖数据计算叶片扫掠起始、结束点，轮毂中心数据计算距离均值
-    start_tip, end_tip, filtered_data_tip = cycle_calculate(data_tip, noise_reduction, min_difference)
+    start_flange, end_flange, filtered_data_flange = cycle_calculate(data_flange, noise_reduction, min_difference)
     start_root, end_root, filtered_data_root = cycle_calculate(data_root, noise_reduction, min_difference)
     start_nan, end_nan, filtered_data_nan = cycle_calculate(data_nan, noise_reduction, min_difference)
 
@@ -142,67 +202,93 @@ def data_analyse(path: List[str]):
     dist_cen = np.mean(filtered_data_cen.iloc[:, 1].tolist())
     filtered_data_cen.iloc[:, 1] = filtered_data_cen.iloc[:, 1] * np.cos(np.deg2rad(angle_cen + axial_inclination))
 
-    # 检查起始结束点顺序，确保叶根叶尖测点同步开始、结束
-    if end_tip.iloc[0, 0] < start_root.iloc[0, 0]:
-        start_tip = start_tip.drop(start_tip.index[0])
-        end_tip = end_tip.drop(end_tip.index[0])
-    if start_root.iloc[0, 0] < start_tip.iloc[0, 0] < end_tip.iloc[0, 0] < end_root.iloc[0, 0]:
-        pass
-    else:
-        raise ValueError("The elements are not in the expected order.")
-
     # 计算叶根、叶中、叶尖处的塔筒距离，对轮毂中心做FFT分析
-    tower_dist_tip = ff.tower_cal(filtered_data_tip, start_tip, end_tip, sampling_fq_1)
+    tower_dist_flange = ff.tower_cal(filtered_data_flange, start_flange, end_flange, sampling_fq_1)
     tower_dist_root = ff.tower_cal(filtered_data_root, start_root, end_root, sampling_fq_1)
     tower_dist_nan = ff.tower_cal(filtered_data_nan, start_nan, end_nan, sampling_fq)
     lowpass_data, fft_x, fft_y, tower_freq, tower_max= ff.process_fft(filtered_data_cen, sampling_fq)
 
     # 根据起始结束点，对叶根、对叶片数据进行归一化处理，计算每个叶片的散点表、线表、边界点表、标准循环周期长度、每个叶片平均最小值
-    result_line_tip, result_scatter_tip, border_rows_tip, cycle_len_tip, min_tip \
-        = data_normalize(filtered_data_tip, start_tip, end_tip, group_length[0])
+    result_line_flange, result_scatter_flange, border_rows_flange, cycle_len_flange, min_flange \
+        = data_normalize(filtered_data_flange, start_flange, end_flange, group_length[0])
     result_line_root, result_scatter_root, border_rows_root, cycle_len_root, min_root \
         = data_normalize(filtered_data_root, start_root, end_root, group_length[1])
     result_line_nan, result_scatter_nan, border_rows_nan, cycle_len_nan, min_nan \
         = data_normalize(filtered_data_nan, start_nan, end_nan, group_length[2])
 
     # 计算3个叶片的平均轮廓，3个叶片的形状差
-    result_avg_tip, result_diff_tip = blade_shape(result_line_tip)
+    result_avg_flange, result_diff_flange = blade_shape(result_line_flange)
     result_avg_root, result_diff_root = blade_shape(result_line_root)
 
     # 对叶尖的边界点表和俯仰角做坐标归一化处理
-    border_rows_tip_new, angle_tip_new = coordinate_normalize(border_rows_tip, angle_tip)
+    border_rows_flange_new, angle_flange_new = coordinate_normalize(border_rows_flange, angle_flange)
     border_rows_nan_new, angle_nan_new = coordinate_normalize(border_rows_nan, angle_nan)
 
+    # 将所有 DataFrame 合并成一个并求平均
+    flange_ava = pd.concat([df['distance'] for df in border_rows_flange_new]).mean(numeric_only=True).mean()
+    root_ava = pd.concat([df['distance'] for df in border_rows_root]).mean(numeric_only=True).mean()
+
+    d_radius = np.abs((flange_ava * np.cos(np.deg2rad(angle_flange_new))
+               - root_ava * np.cos(np.deg2rad(angle_root))) * np.sin(np.deg2rad(axial_inclination))
+               + (flange_ava * np.sin(np.deg2rad(angle_flange_new))
+               - root_ava * np.sin(np.deg2rad(angle_root))) * np.cos(np.deg2rad(axial_inclination)))
+
+    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)
+
     # 对叶片的边界点表做半径计算
-    tip_r = radius_cal(border_rows_tip_new, angle_tip_new, dist_cen, angle_cen, axial_inclination, angle_cone)
+    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)
 
+    if np.abs((root_r - flange_r) - d_radius) > 0.5:
+        print(str(root_r - flange_r) + "对比" + str(d_radius))
+        raise ValueError("Radius err1.")
+
+    if np.abs(flange_root_dist - d_radius) > 0.5:
+        print(str(flange_root_dist) + "对比" + str(d_radius))
+        raise ValueError("Radius err2.")
+
+    # root_r = d_radius + blade_axis_new["旋转半径"].iloc[0]
+    # root_flange_dist = d_radius / np.cos(np.deg2rad(angle_cone))
+    # blade_axis_new["中心y"] = np.sqrt(blade_axis_new["中心y"] ** 2 + root_flange_dist ** 2 - 2 *
+    #     root_flange_dist * blade_axis_new["中心y"]
+    #     * np.cos(np.deg2rad(90 - angle_cone - axial_inclination - angle_flange_new)))
+
+    blade_axis_new["中心y"] = blade_axis_new["中心y"] - (flange_ava - root_ava)
+
     # 计算叶片测量位置处的绝对桨距角、相对桨距角、线速度、叶片内部中心点距离
-    pitch_angle_tip, aero_dist_tip, v_speed_tip, cen_blade_tip = (
-        blade_angle_aero_dist(border_rows_tip, tip_r, cycle_len_tip, tower_dist_tip, angle_tip_new))
-    pitch_angle_root, aero_dist_root, v_speed_root, cen_blade_root = (
+    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))
-    pitch_angle_nan, aero_dist_nan, v_speed_nan, cen_blade_nan = (
+    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))
+
+    blade_axis_new["中心y"] = blade_axis_new["中心y"]*np.cos(np.deg2rad(angle_root + angle_cone + axial_inclination))
 
     # 将列表转换为 numpy 数组
-    cen_blade_tip_array = np.array(cen_blade_tip)
+    cen_blade_flange_array = np.array(cen_blade_flange)
     cen_blade_nan_array = np.array(cen_blade_nan)
-    min_tip_array = np.array(min_tip)
+    min_flange_array = np.array(min_flange)
     min_nan_array = np.array(min_nan)
 
     # 计算叶片内部中心点距离与叶片最小值之间的差值
-    abs_diff = np.abs(cen_blade_tip_array - min_tip_array)
+    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_tip = abs_diff * np.cos(np.deg2rad(angle_tip_new))
+    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_tip.tolist()
+    blade_dist_flange.tolist()
     blade_dist_nan.tolist()
 
     # 计算叶片转速-净空散点表
-    dist_distribute = blade_dist_distribute_cal(filtered_data_tip, start_tip, end_tip,
-                                                tower_dist_tip, angle_tip_new, blade_dist_tip)
+    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)
     # dist_distribute = [df.round(5) for df in dist_distribute]
@@ -228,23 +314,27 @@ def data_analyse(path: List[str]):
         mean_values.append(round((max_values[i] + min_values[i]) / 2, 2))
 
     # 将叶片线表数据乘以线速度，和俯仰角，得到叶片横截面的真实轮廓
-    for df in result_line_tip:
+    for df in result_line_flange:
         first_column = df.iloc[:, 0]
         sec_column = df.iloc[:, 1]
-        df.iloc[:, 0] = first_column * v_speed_tip
-        df.iloc[:, 1] = sec_column * np.cos(np.deg2rad(angle_tip_new + angle_cone + axial_inclination))
+        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:
         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
         df.iloc[:, 1] = sec_column * np.cos(np.deg2rad(angle_root + angle_cone + axial_inclination))
 
-    for df in result_scatter_tip:
+    for df in result_scatter_flange:
         first_column = df.iloc[:, 0]
         sec_column = df.iloc[:, 1]
-        df.iloc[:, 0] = first_column * v_speed_tip
-        df.iloc[:, 1] = sec_column * np.cos(np.deg2rad(angle_tip_new + angle_cone + axial_inclination))
+        df.iloc[:, 0] = first_column * v_speed_flange
+        df.iloc[:, 1] = sec_column * np.cos(np.deg2rad(angle_flange_new + angle_cone + axial_inclination))
 
     for df in result_scatter_root:
         first_column = df.iloc[:, 0]
@@ -253,22 +343,23 @@ def data_analyse(path: List[str]):
         df.iloc[:, 1] = sec_column * np.cos(np.deg2rad(angle_root + angle_cone + axial_inclination))
 
     # 将叶片平均轮廓数据乘以线速度，得到实际叶片长度
-    avg_tip = result_avg_tip.iloc[:, 0]
-    result_avg_tip.iloc[:, 0] = avg_tip * v_speed_tip
+    avg_flange = result_avg_flange.iloc[:, 0]
+    result_avg_flange.iloc[:, 0] = avg_flange * v_speed_flange
     avg_root = result_avg_root.iloc[:, 0]
     result_avg_root.iloc[:, 0] = avg_root * v_speed_root
 
     # 计算叶片扭转角度
-    twist_1 = round(np.abs(pitch_angle_root[0] - pitch_angle_tip[0]), 2)
-    twist_2 = round(np.abs(pitch_angle_root[1] - pitch_angle_tip[1]), 2)
-    twist_3 = round(np.abs(pitch_angle_root[2] - pitch_angle_tip[2]), 2)
+    pitch_angle_flange = [0, 0, 0]
+    twist_1 = round(np.abs(pitch_angle_root[0] - pitch_angle_flange[0]), 2)
+    twist_2 = round(np.abs(pitch_angle_root[1] - pitch_angle_flange[1]), 2)
+    twist_3 = round(np.abs(pitch_angle_root[2] - pitch_angle_flange[2]), 2)
     twist_avg = round((twist_1 + twist_2 + twist_3) / 3, 2)
 
     # 降低给数据采样频率，降低接口负担
     sampling_num = int(0.015 * sampling_fq_1)
 
     # 将原始数据的时间列由计时时钟转换为实际时间
-    data_tip.iloc[:, 0] = data_tip.iloc[:, 0] / 5000000
+    data_flange.iloc[:, 0] = data_flange.iloc[:, 0] / 5000000
     data_root.iloc[:, 0] = data_root.iloc[:, 0] / 5000000
     lowpass_data.iloc[:, 0] = lowpass_data.iloc[:, 0] / 5000000
 
@@ -303,8 +394,8 @@ def data_analyse(path: List[str]):
     json_output = {
         'original_plot': {
             'blade_tip': {
-                'xdata': data_tip.iloc[:, 0].tolist()[::sampling_num],
-                'ydata': data_tip.iloc[:, 1].tolist()[::sampling_num]
+                'xdata': data_flange.iloc[:, 0].tolist()[::sampling_num],
+                'ydata': data_flange.iloc[:, 1].tolist()[::sampling_num]
             },
             'blade_root': {
                 'xdata': data_root.iloc[:, 0].tolist()[::sampling_num],
@@ -331,20 +422,20 @@ def data_analyse(path: List[str]):
         },
         'blade_tip': {
             'first_blade': {
-                'xdata': result_line_tip[0].iloc[:, 0].tolist(),
-                'ydata': result_line_tip[0].iloc[:, 1].tolist()
+                'xdata': result_line_flange[0].iloc[:, 0].tolist(),
+                'ydata': result_line_flange[0].iloc[:, 1].tolist()
             },
             'second_blade': {
-                'xdata': result_line_tip[1].iloc[:, 0].tolist(),
-                'ydata': result_line_tip[1].iloc[:, 1].tolist()
+                'xdata': result_line_flange[1].iloc[:, 0].tolist(),
+                'ydata': result_line_flange[1].iloc[:, 1].tolist()
             },
             'third_blade': {
-                'xdata': result_line_tip[2].iloc[:, 0].tolist(),
-                'ydata': result_line_tip[2].iloc[:, 1].tolist()
+                'xdata': result_line_flange[2].iloc[:, 0].tolist(),
+                'ydata': result_line_flange[2].iloc[:, 1].tolist()
             },
             'avg_blade': {
-                'xdata': result_avg_tip.iloc[:, 0].tolist(),
-                'ydata': result_avg_tip.iloc[:, 1].tolist()
+                'xdata': result_avg_flange.iloc[:, 0].tolist(),
+                'ydata': result_avg_flange.iloc[:, 1].tolist()
             }
         },
         'blade_root': {
@@ -452,42 +543,47 @@ def data_analyse(path: List[str]):
     with open(json_file_path, 'w') as json_file:
         json.dump(json_output, json_file, indent=4)
 
+    print('-' * 50)
     print('csv文件路径' + str(csv_file_path))
-    # print(result_line_tip[0].iloc[:, 0])
-    # print(result_line_root[0].iloc[:, 0])
+    # print(result_line_flange[0].iloc[:, 0])
     print('振动主频' + str(tower_freq))
     print('振动幅值' + str(tower_max))
     print('净空最小值', min_values)
     print('最小值对应的键', min_keys)
     print('净空最大值', max_values)
     print('最大值对应的键', max_keys)
-    print('叶尖速度' + str(v_speed_tip), '叶根速度' + str(v_speed_root))
-    print('新俯仰角' + str(angle_tip_new))
+    print('叶尖速度' + str(v_speed_flange), '叶根速度' + str(v_speed_root))
+    print('新俯仰角' + str(angle_flange_new))
     print('轮毂中心距离' + str(dist_cen))
     print('叶根原始数据采样时间长度' + str(data_root.iloc[-1, 0]))
     print('-' * 50)
-    print(json.dumps(json_output, indent=4, ensure_ascii=False))
 
-    # plot_data(result_line_tip, 'line', 'data1')
-    # plot_data(result_diff_tip, 'line', 'data_diff_1')
-    # plot_data(result_scatter_tip, 'scatter', 'data1')
-    plot_data(result_line_root, 'line', 'data2')
+    # plot_data(result_line_flange, '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_diff_root, 'line', 'data_diff_2')
-    plot_data(result_scatter_root, 'scatter', 'data2')
+    # plot_data(result_scatter_root, blade_axis_new,  'scatter', 'data2')
     # plot_data(dist_distribute, 'scatter', 'dist_distribute')
 
     return json_output
 
 
-def process_data(file_path):
+def process_data(file_path, skew_angle):
 
     """
     打开、解决时间重置、按时间清洗异常值、分列数据
     """
 
+    skew_angle = 90 - skew_angle
     # 读取第2、4、9列的数据
     data = pd.read_csv(file_path, usecols=[1, 3, 4, 8, 9], header=None, engine='c')
     data = data.head(int(len(data) * 0.95))
+    # n = len(data)
+    # start = int(n * 0.35)  # 开始位置在35%处，这样可以取中间30%
+    # end = int(n * 0.65)  # 结束位置在65%处，这样总共是30%的数据
+    # data = data.iloc[start:end]
+
     print('原始数据长度' + str(len(data)))
 
     '''
@@ -540,10 +636,48 @@ def process_data(file_path):
     data_1.columns = ['time', 'distance', 'grey']
     data_2.columns = ['time', 'distance', 'grey']
 
+    data_1['time'] = data_1['time'] + data_1['distance'] * np.cos(np.deg2rad(skew_angle))
+    data_1['distance'] = data_1['distance'] * np.sin(np.deg2rad(skew_angle))
+    data_2['time'] = data_2['time'] + data_2['distance'] * np.cos(np.deg2rad(skew_angle))
+    data_2['distance'] = data_2['distance'] * np.sin(np.deg2rad(skew_angle))
+    data_1['time'] = data_1['time'].round().astype(int)
+    data_2['time'] = data_2['time'].round().astype(int)
 
     return data_1, data_2
 
 
+def locate_filter(data_group: pd.DataFrame):
+
+    """
+    求定位数据的叶片距离均值，并将数据标准化，以便放入cycle_calculate，用作求2组定位点的旋转角度
+    :param data_group: process_data计算完成后的数据。
+    :return: distance_mean：测量距离值
+    :return: filtered_data：w后的数据
+    """
+
+    # 创建新的DataFrame以避免修改原始数据
+    filtered_data = data_group.copy()
+
+    # 去掉distance为0的行
+    data_group = data_group[data_group['distance'] != 0]
+
+    # 计算grey列20%-80%分位数的均值
+    grey_mean = data_group['grey'].quantile(0.2)
+    grey_upper = data_group['grey'].quantile(0.8)
+    grey_threshold = data_group[(data_group['grey'] >= grey_mean) & (data_group['grey'] <= grey_upper)]['grey'].mean()
+
+    # 去掉grey小于均值*0.8的行
+    data_group = data_group[data_group['grey'] >= grey_threshold * 0.8]
+
+    # 计算剩余行distance列的均值
+    distance_mean = data_group['distance'].mean()
+
+    # 将distance为0的值替换为distance_mean+10
+    filtered_data.loc[filtered_data['distance'] == 0, 'distance'] = distance_mean + 10
+
+    return distance_mean, filtered_data
+
+
 def tower_filter(data_group: pd.DataFrame, noise_threshold: float):
 
     """
@@ -600,15 +734,19 @@ 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
 
-    # 选择频率最大的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-30) | (
-                data_group['distance'] > mean_values*1.1), 'distance'] = np.nan
+    if 0 in distance_counts.nlargest(3).index:
+        pass
+    else:
+        # 选择频率最大的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-30) | (
+                    data_group['distance'] > mean_values*1.1), 'distance'] = np.nan
+        print(f"中值是：{mean_values}")
 
     nan_count = data_group['distance'].isna().sum()
     all_count = data_group.shape[0]
-    print(f"中值是：{mean_values}，替换为NaN的distance异常值的数量是: {nan_count}, 总数量是: {all_count},"
+    print(f"替换为NaN的distance异常值的数量是: {nan_count}, 总数量是: {all_count},"
           f"占比: {nan_count / all_count * 100:.2f}%")
 
     # 前向填充
@@ -629,7 +767,7 @@ def cycle_calculate(data_group: pd.DataFrame, noise_threshold: float, min_distan
         # 获取当前行的 distance 值
         current_distance = filtered_data.loc[idx, 'distance']
 
-        next_rows_large = filtered_data.loc[idx - 200: idx - 1]
+        next_rows_large = filtered_data.loc[idx - 500: idx - 1]
 
         # 检查是否任意 distance 的值小于 current_distance - 2
         if next_rows_large['distance'].le(current_distance - min_distance).all():
@@ -640,13 +778,15 @@ def cycle_calculate(data_group: pd.DataFrame, noise_threshold: float, min_distan
         # 获取当前行的 distance 值
         current_distance = filtered_data.loc[idx - 1, 'distance']
 
-        next_rows_small = filtered_data.iloc[idx: idx + 200]
+        next_rows_small = filtered_data.iloc[idx: idx + 500]
 
         # 检查是否任意 distance 的值小于 current_distance - 2
         if next_rows_small['distance'].le(current_distance - min_distance).all():
             # 如果都小于，则将当前行和下一行添加到 special_points 中
             start_points = pd.concat([start_points, filtered_data.loc[[idx]]])
 
+    if 0 in distance_counts.nlargest(3).index:
+        end_points, start_points = start_points, end_points  # 互换
 
     if end_points.iloc[0, 0] < start_points.iloc[0, 0]:
         end_points = end_points.drop(end_points.index[0])
@@ -686,7 +826,7 @@ def data_normalize(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
 
     # 将 start_points 中的时间点转换为列表
     start_times = combined_df_sorted['time'].tolist()
-    print('本次测量风机完整旋转圈数：'+ str(len(start_times) / 2))
+    print('本次测量风机完整旋转圈数：'+ str(int(len(start_times) / 6)))
     time.sleep(1)
 
 
@@ -741,8 +881,6 @@ def data_normalize(data_group: pd.DataFrame, start_points: pd.DataFrame, end_poi
 
         # 分组，时间列每1000为一组（每40个时间点一组）
         bins = list(range(int(first_time), int(turbine_sorted['time'].max()), group_len))
-        # 原始代码
-        # bins = list(range(int(first_time), int(turbine_sorted['time'].max()) + len(start_times), int(fs / 50)))
         grouped = turbine_sorted.groupby(pd.cut(turbine_sorted['time'], bins=bins, right=False))
 
         # 初始化一个空的 DataFrame 用于存储处理后的数据
@@ -909,10 +1047,10 @@ def blade_shape(turbines_processed: List[pd.DataFrame]):
 def coordinate_normalize(tip_border_rows: List[pd.DataFrame], tip_angle):
 
     """
-    将叶尖测量数据和叶根、轮毂中心的测量原点归一化。
+    将叶尖、法兰测量数据和叶根、轮毂中心的测量原点归一化。
     :param tip_border_rows: 3个叶尖边缘数据
     :param tip_angle: 叶尖测量俯仰角
-    :return: 归一化后叶尖数据，叶尖俯仰角
+    :return: 归一化后叶尖数据，新俯仰角
     """
 
     tip_angle1 = np.deg2rad(tip_angle)
@@ -929,11 +1067,248 @@ def coordinate_normalize(tip_border_rows: List[pd.DataFrame], tip_angle):
 
     tip_angle_new = float(np.mean(tip_angle_list))
     tip_angle_new1 = np.rad2deg(tip_angle_new)
-    print('坐标转换后的新叶尖俯仰角: ' + str(tip_angle_new1))
+    print('叶尖俯仰角: ' + str(tip_angle) + '坐标转换后的新叶尖俯仰角: ' + str(tip_angle_new1))
     
     return tip_border_rows, tip_angle_new1
 
 
+def flange_coordinate_normalize(flange_cen_row: pd.DataFrame, flange_angle):
+    """
+    将法兰中心计算数据测量原点归一化。
+    :param flange_cen_row: 法兰数据
+    :param flange_angle: 法兰测量俯仰角
+    :return: 归一化后法兰数据，新俯仰角
+    """
+    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']))
+    flange_angle_cal = np.arctan(flange_angle_cal0)
+
+    # 更新中心y列
+    flange_cen_row['中心y'] = (flange_cen_row['中心y'] ** 2 + 0.0057881664 -
+                                     0.15216 * flange_cen_row['中心y'] * np.sin(flange_angle1)) ** 0.5
+
+    # 计算新的俯仰角（由于现在只有一个值，直接使用计算出的值）
+    flange_angle_new = float(flange_angle_cal)
+    flange_angle_new1 = np.rad2deg(flange_angle_new)
+    print('坐标转换后的新法兰俯仰角: ' + str(flange_angle_new1))
+
+    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,
+                   group_len: int) \
+        -> pd.DataFrame:
+
+    """
+    每个叶片周期计算旋转中心及半径
+    :param data_group: cycle_calculate计算完成后的叶根法兰数据。
+    :param start_points: 所有每个周期开始点，叶片前缘突变点。
+    :param end_points: 叶片后缘突变点。
+    :param horizon_angle: 叶根法兰水平角度。
+    :param group_len: 每数据长度。
+    :return: 叶根法兰圆轮廓的圆心、测量点的旋转半径
+    """
+    print('正在进行叶片轴线计算......')
+
+    def fit_circle(df, v_bounds=(2, 10.0), top_k=5, prefilter=True):
+        """
+        Savitzky–Golay 平滑 + 差分进化 + least_squares 精修
+        Savitzky–Golay 平滑在每个局部窗口里用多项式拟合来代替原始点，从而减少噪声影响，同时保留曲线的整体趋势和细节
+        """
+
+        # ========== Savitzky–Golay 平滑 ==========
+        def smooth_savgol(y, window_length=101, polyorder=3):
+            wl = min(window_length, len(y) if len(y) % 2 == 1 else len(y) - 1)
+            if wl < 3:
+                return y
+            if wl % 2 == 0:
+                wl -= 1
+            return savgol_filter(y, wl, polyorder)
+
+        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)]
+
+        def residuals_sq(params):
+            v, xc, yc, R = params
+            if v <= 0 or 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,
+            strategy='rand2bin',
+            mutation=(0.8, 1.2),
+            recombination=0.8,
+            popsize=30,
+            maxiter=1000,
+            polish=False,
+            seed=42,
+            workers=1
+        )
+
+        # 多候选点精修
+        pop = result.population
+        energies = result.population_energies
+        idx = np.argsort(energies)[:top_k]
+        candidates = pop[idx]
+
+        best_rmse = np.inf
+        best_result = None
+
+        for cand in candidates:
+            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]),
+                method='trf',
+                loss='linear',
+                max_nfev=50000,
+                xtol=1e-12,
+                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)
+            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]
+
+        result_df = pd.DataFrame([best_result],
+                                 columns=['旋转半径', '中心x', '中心y', '圆半径', '几何RMSE'])  # 旋转半径本身为测量点线速度
+        return result_df
+
+    def plot_circle_fit(df, fit_result):
+        """
+        df: 原始数据，包含 'time' 和 'distance'
+        fit_result: fit_circle 输出的 DataFrame
+        """
+        t = np.asarray(df['time'])
+        d_raw = np.asarray(df['distance'])
+
+        # 最优参数
+        v, xc, yc, R = fit_result[['旋转半径', '中心x', '中心y', '圆半径']].values[0]
+
+        # x 轴乘速度后的原始散点
+        x_all = v * t
+
+        plt.figure(figsize=(8, 8))
+
+        # 原始散点
+        plt.scatter(x_all, d_raw, s=10, color='blue', alpha=0.5, label='原始散点 (v*t)')
+
+        # 拟合圆
+        theta = np.linspace(0, 2 * np.pi, 500)
+        x_circle = xc + R * np.cos(theta)
+        y_circle = yc + R * np.sin(theta)
+        plt.plot(x_circle, y_circle, color='red', lw=2, label='拟合圆')
+
+        # 圆心
+        plt.scatter([xc], [yc], color='green', s=50, label='圆心')
+
+        plt.xlabel('x = v * t')
+        plt.ylabel('distance')
+        plt.axis('equal')
+        plt.legend()
+        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]):
+        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)
+    angle_speed = (np.pi / full_cycle) * 5000000
+    turbines = [pd.DataFrame() for _ in range(3)]
+
+    for i in range(0, len(start_times), 2):
+
+
+        start_time = start_times[i]
+        end_time = start_times[i + 1]
+        segment = data_group[(data_group['time'] > start_time) & (data_group['time'] <= end_time)]
+        if segment is None or segment.empty:
+            raise ValueError("Segment is empty")
+
+        segment = segment.copy()
+        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_processed = []
+    turbines_scattered = []
+    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))
+
+        process_df = pd.DataFrame()
+        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)
+        processed_df = process_df[(process_df['time'] >= lower_bound) & (process_df['time'] <= upper_bound)]
+        blade_cen_est = fit_circle(processed_df)
+
+        plot_circle_fit(processed_df, blade_cen_est)
+        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["中心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
+
 
 def radius_cal(border_rows, meas_angle, cen_dist, cen_angle, angle_main, angle_rotate):
 
@@ -951,7 +1326,6 @@ def radius_cal(border_rows, meas_angle, cen_dist, cen_angle, angle_main, angle_r
     aero_dist = (pd.concat([df['distance'] for df in border_rows]).mean())
     radius = np.abs(aero_dist * np.sin(np.deg2rad(meas_angle - angle_main))
                     - cen_dist * np.sin(np.deg2rad(cen_angle - angle_main)))
-    print('测量点旋转半径：' + str(radius))
 
     return radius
 
@@ -967,43 +1341,89 @@ def blade_angle_aero_dist(border_rows: List[pd.DataFrame], radius: float, full_c
     :param full_cycle: 全周期
     :param tower_dist: 塔筒距离
     :param v_angle: 俯仰角度
-    :return: 绝对桨距角，净空距离，叶片线速度
+    :return: 净空距离，叶片线速度，叶片中心距离
     """
 
-    print('正在进行相对桨距角和叶片净空距离计算......')
+    print('正在进行叶片净空距离计算......')
     v_speed = 2 * np.pi * radius / full_cycle  # 叶片线速度m/(1计时器单位）
-    pitch_angle_list = []
     aero_dist_list = []
     cen_blade = []
     for turbine in border_rows:
 
-        diff_time = turbine.iloc[1, 0] - turbine.iloc[0, 0]
-
-        diff_len = (turbine.iloc[1, 1] - turbine.iloc[0, 1]) * np.cos(np.deg2rad(v_angle))
         mean_col2 = (turbine.iloc[1, 1] + turbine.iloc[0, 1]) / 2
         aero_dist = abs(mean_col2 - tower_dist) * np.cos(np.deg2rad(v_angle))
 
-        pitch_angle = np.degrees(np.arctan(diff_len / (diff_time * v_speed)))
-        print('单个叶片绝对桨距角' + str(pitch_angle))
-        pitch_angle_list.append(pitch_angle)
         aero_dist_list.append(aero_dist)
         cen_blade.append(mean_col2)
+
+    aero_dist_list.append(np.mean(aero_dist_list))
+    aero_dist_list = [round(num, 2) for num in aero_dist_list]
+
+    return aero_dist_list, v_speed, cen_blade
+
+
+def blade_angle(border_rows: List[pd.DataFrame], cen_data: pd.DataFrame, radius: float, full_cycle: int,
+                                   v_angle: float):
+
+    """
+    计算叶片相对桨距角。
+    :param border_rows: 三个叶片的边界
+    :param cen_data: 旋转中心数据
+    :param radius: 旋转半径
+    :param full_cycle: 全周期
+    :param v_angle: 俯仰角度
+    :return: 绝对桨距角，叶片线速度
+    """
+
+    print('正在进行相对桨距角计算......')
+    print('叶根半径：' + str(radius))
+    v_speed = 2 * np.pi * radius / full_cycle  # 叶片线速度m/(1计时器单位）
+
+    values = []
+    for df in border_rows:
+        if df.shape[0] >= 2 and df.shape[1] >= 2:
+            values.append(df.iloc[0, 1])
+            values.append(df.iloc[1, 1])
+
+    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])
+
+    pitch_angle_list = []
+    for idx, turbine in enumerate(border_rows, start=1):
+
+        # 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))
+        pitch_angle_list.append(pitch_angle)
+
     pitch_mean = np.mean(pitch_angle_list)
     pitch_angle_list = [angle - pitch_mean for angle in pitch_angle_list]
     pitch_angle_list.append(max(pitch_angle_list) - min(pitch_angle_list))
-    aero_dist_list.append(np.mean(aero_dist_list))
     pitch_angle_list = [round(num, 2) for num in pitch_angle_list]
-    aero_dist_list = [round(num, 2) for num in aero_dist_list]
-
-    return pitch_angle_list, aero_dist_list, v_speed, cen_blade
+    print("\033[34m叶片相对角度偏差：\033[0m", pitch_angle_list)
 
+    return pitch_angle_list, v_speed
 
 
-def plot_data(data, plot_type: str, data_name: str):
+def plot_data(data, point, plot_type: str, data_name: str):
 
     """
     绘制数据图表并保存为文件。
     :param data: 数据列表，每个元素是一个 DataFrame。
+    :param point: 测量点。
     :param plot_type: 图表类型，'line' 或 'scatter'。
     :param data_name: 数据名称，用于生成文件名。
     """
@@ -1023,6 +1443,7 @@ def plot_data(data, 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='测量点')
     axy = plt.gca()  # 获取当前坐标轴对象
     plt.grid(which='both', linewidth=2)  # 设置网格线宽度为2
     axy.xaxis.set_major_locator(MaxNLocator(nbins=200))  # 设置x轴主刻度的最大数量为10
@@ -1129,31 +1550,98 @@ def blade_dist_distribute_cal(data_group: pd.DataFrame, start_points: pd.DataFra
 
     return tower_clearance
 
-
-# locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tf-20_20250630223600_20_13.03_23.32.csv"
-# measure_path= "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tf-20_20250630223849_20_17.89_21.07.csv"
-
-# locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tf-50_20250630223358_50_13.03_23.32.csv"
-# measure_path= "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tf-50_20250630224408_50_17.89_21.07.csv"
-
-# locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tf-100_20250630222752_100_13.03_23.32.csv"
-# measure_path= "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tf-100_20250630225119_100_17.89_21.07.csv"
-
-# locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tff-20_20250630231223_20_12.51_20.06.csv"
-# measure_path= "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tff-20_20250630232052_20_15.36_18.17.csv"
-
-# locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tff-50_20250630231417_50_12.51_20.06.csv"
-# measure_path= "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tff-50_20250630233420_50_15.35_18.16.csv"
-
-# locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tff-100_20250630231610_100_12.51_20.06.csv"
-# measure_path= "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250728/gy_10-tff-100_20250630234012_100_15.35_18.16.csv"
-
-locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/gytest/测试数据/gy_18-RF-1_20250701154647_50_23.70_40.01.csv"
-measure_path= "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/gytest/测试数据/gy_18-RF-2_20250701155057_50_29.30_36.78.csv"
-
-start_t = time.time()  # 记录开始时间
-data_path = [locate_path, measure_path, 5, 6]
-list_1 = data_analyse(data_path)
-# print(list_1)
-print(f"耗时: {time.time() - start_t:.2f} 秒")
+if __name__ == "__main__":
+    # 左偏
+    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250919怀来/1/hl_F14-L_20250819134200_50_11.25_13.42.csv"
+    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250919怀来/1/hl_F14-L_20250819133535_50_7.55_12.10.csv"
+    # measure_path= "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250919怀来/1/hl_F14-L_20250819132731_50_11.71_10.34.csv"
+    # 右偏
+    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250919怀来/1/hl_F14-last_20250819144135_50_9.36_11.02.csv"
+    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250919怀来/1/hl_F14-last_20250819143634_50_7.44_10.04.csv"
+    # measure_path= "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250919怀来/1/hl_F14-last_20250819143011_50_9.81_8.57.csv"
+    # 左微偏
+    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250919怀来/2/hl_F14-z_20250819124400_50_11.32_13.07.csv"
+    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250919怀来/2/hl_F14-z_20250819124108_50_7.51_12.05.csv"
+    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20250919怀来/2/hl_F14-z_20250819123507_50_11.23_10.51.csv"
+    # 第一组正对
+    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-z_20251011100836_50_26.95_30.34.csv"
+    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-z_20251011095956_50_20.23_28.22.csv"
+    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-z_20251011095339_50_27.64_26.44.csv"
+    # 第二组正对
+    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-z_20251011100836_50_26.95_30.34.csv"
+    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-z_20251011102933_50_20.48_28.28.csv"
+    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-z_20251011102339_50_27.72_26.33.csv"
+    # 第三组左侧
+    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-l_20251011105053_50_25.67_27.83.csv"
+    # locate_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-l_20251011104654_50_18.84_26.34.csv"
+    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-l_20251011104100_50_25.84_24.66.csv"
+    # 第四组右侧
+    # 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_20-rr_20251011114224_50_22.36_30.22.csv"
+    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-r_20251011111012_50_29.57_28.26.csv"
+    # 第五组左侧
+    # 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_20-ll_20251011121305_50_25.22_34.48.csv"
+    # measure_path = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-ll_20251011120724_50_33.79_32.64.csv"
+    # 第六组正对
+    # locate_path0 = "C:/Users/laiwe/Desktop/风电/激光测量/测试数据/20251011沽源/验证/gy_20-zz_20251011123416_50_31.48_33.93.csv"
+    # 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_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"
+
+    start_t = time.time()  # 记录开始时间
+    data_path = [locate_path0, locate_path, measure_path, 5, 5, 1.54]  # 偏斜测量数据、轮毂数据、叶根数据、锥角、轴向倾角、偏航角
+    list_1 = data_analyse(data_path)
+    print(f"耗时: {time.time() - start_t:.2f} 秒")
 

data_clean.py

@@ -8,6 +8,9 @@ import warnings
 import sys
 import frequency_filter as ff
 from datetime import datetime
+from scipy.optimize import least_squares, differential_evolution
+from scipy.signal import savgol_filter
+
 
 warnings.filterwarnings("ignore", category=FutureWarning)  
 plt.rcParams['font.sans-serif'] = ['SimHei']  
@@ -89,76 +92,106 @@ def data_analyse(path: List[str]):
     min_difference = 1.5  
     angle_cone = float(path[2])
     axial_inclination = float(path[3])
+    lift_up_limit = float(path[4])
     group_length = [10000, 20000, 5000]
     return_list = []
 
     wind_name, turbine_code, time_code, sampling_fq, angle_nan, angle_cen = find_param(locate_file)
-    wind_name_1, turbine_code_1, time_code_1, sampling_fq_1, angle_tip, angle_root = find_param(measure_file)
+    wind_name_1, turbine_code_1, time_code_1, sampling_fq_1, angle_flange, angle_root = find_param(measure_file)
 
     sampling_fq_1 = sampling_fq_1 * 1000
     sampling_fq = sampling_fq * 1000
 
     data_nan, data_cen = process_data(locate_file)
-    data_tip, data_root = process_data(measure_file)
+    data_flange, data_root = process_data(measure_file)
+
+    if lift_up_limit >= 0.1:
+        discrete_values = np.arange(0, 0.101, 0.001)
+        condition = data_flange['distance'] > lift_up_limit
+        n = condition.sum()
+        random_discrete = np.random.choice(discrete_values, size=n)
+        data_flange.loc[condition, 'distance'] = lift_up_limit + 3 + random_discrete
+    elif np.abs(lift_up_limit) < 0.1:
+        pass
+    else:
+        raise ValueError("lift_up_limit error.")
 
-    start_tip, end_tip, filtered_data_tip = cycle_calculate(data_tip, noise_reduction, min_difference)
+    start_flange, end_flange, filtered_data_flange = cycle_calculate(data_flange, noise_reduction, min_difference)
     start_root, end_root, filtered_data_root = cycle_calculate(data_root, noise_reduction, min_difference)
     start_nan, end_nan, filtered_data_nan = cycle_calculate(data_nan, noise_reduction, min_difference)
     filtered_data_cen = tower_filter(data_cen, noise_reduction)
     dist_cen = np.mean(filtered_data_cen.iloc[:, 1].tolist())
     filtered_data_cen.iloc[:, 1] = filtered_data_cen.iloc[:, 1] * np.cos(np.deg2rad(angle_cen + axial_inclination))
 
-    if end_tip.iloc[0, 0] < start_root.iloc[0, 0]:
-        start_tip = start_tip.drop(start_tip.index[0])
-        end_tip = end_tip.drop(end_tip.index[0])
-    if start_root.iloc[0, 0] < start_tip.iloc[0, 0] < end_tip.iloc[0, 0] < end_root.iloc[0, 0]:
-        pass
-    else:
-        raise ValueError("The elements are not in the expected order.")
-
-    tower_dist_tip = ff.tower_cal(filtered_data_tip, start_tip, end_tip, sampling_fq_1)
+    tower_dist_flange = ff.tower_cal(filtered_data_flange, start_flange, end_flange, sampling_fq_1)
     tower_dist_root = ff.tower_cal(filtered_data_root, start_root, end_root, sampling_fq_1)
     tower_dist_nan = ff.tower_cal(filtered_data_nan, start_nan, end_nan, sampling_fq)
     lowpass_data, fft_x, fft_y, tower_freq, tower_max = ff.process_fft(filtered_data_cen, sampling_fq)
 
-    result_line_tip, result_scatter_tip, border_rows_tip, cycle_len_tip, min_tip \
-        = data_normalize(filtered_data_tip, start_tip, end_tip, group_length[0])
+    result_line_flange, result_scatter_flange, border_rows_flange, cycle_len_flange, min_flange \
+        = data_normalize(filtered_data_flange, start_flange, end_flange, group_length[0])
     result_line_root, result_scatter_root, border_rows_root, cycle_len_root, min_root \
         = data_normalize(filtered_data_root, start_root, end_root, group_length[1])
     result_line_nan, result_scatter_nan, border_rows_nan, cycle_len_nan, min_nan \
         = data_normalize(filtered_data_nan, start_nan, end_nan, group_length[2])
 
-    result_avg_tip, result_diff_tip = blade_shape(result_line_tip)
+    result_avg_flange, result_diff_flange = blade_shape(result_line_flange)
     result_avg_root, result_diff_root = blade_shape(result_line_root)
 
-    border_rows_tip_new, angle_tip_new = coordinate_normalize(border_rows_tip, angle_tip)
+    border_rows_flange_new, angle_flange_new = coordinate_normalize(border_rows_flange, angle_flange)
     border_rows_nan_new, angle_nan_new = coordinate_normalize(border_rows_nan, angle_nan)
 
-    tip_r = radius_cal(border_rows_tip_new, angle_tip_new, dist_cen, angle_cen, axial_inclination, angle_cone)
+    flange_ava = pd.concat([df['distance'] for df in border_rows_flange_new]).mean(numeric_only=True).mean()
+    root_ava = pd.concat([df['distance'] for df in border_rows_root]).mean(numeric_only=True).mean()
+
+    d_radius = np.abs((flange_ava * np.cos(np.deg2rad(angle_flange_new))
+               - root_ava * np.cos(np.deg2rad(angle_root))) * np.sin(np.deg2rad(axial_inclination))
+               + (flange_ava * np.sin(np.deg2rad(angle_flange_new))
+               - root_ava * np.sin(np.deg2rad(angle_root))) * np.cos(np.deg2rad(axial_inclination)))
+
+    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)
 
-    pitch_angle_tip, aero_dist_tip, v_speed_tip, cen_blade_tip = (
-        blade_angle_aero_dist(border_rows_tip, tip_r, cycle_len_tip, tower_dist_tip, angle_tip_new))
-    pitch_angle_root, aero_dist_root, v_speed_root, cen_blade_root = (
+    if np.abs((root_r - flange_r) - d_radius) > 0.5:
+        raise ValueError("Radius err1.")
+
+    if np.abs(flange_root_dist - d_radius) > 0.5:
+        raise ValueError("Radius err2.")
+
+    blade_axis_new["中心y"] = blade_axis_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))
-    pitch_angle_nan, aero_dist_nan, v_speed_nan, cen_blade_nan = (
+    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))
 
+    blade_axis_new["中心y"] = blade_axis_new["中心y"]*np.cos(np.deg2rad(angle_root + angle_cone + axial_inclination))
     
-    cen_blade_tip_array = np.array(cen_blade_tip)
+    cen_blade_flange_array = np.array(cen_blade_flange)
     cen_blade_nan_array = np.array(cen_blade_nan)
-    min_tip_array = np.array(min_tip)
+    min_flange_array = np.array(min_flange)
     min_nan_array = np.array(min_nan)
-    abs_diff = np.abs(cen_blade_tip_array - min_tip_array)
+    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_tip = abs_diff * np.cos(np.deg2rad(angle_tip_new))
+    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_tip.tolist()
+    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, blade_dist_nan)
     dist_distribute = [df.round(5) for df in dist_distribute_nan]
 
     
@@ -180,30 +213,35 @@ 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_tip:
+    for df in result_line_flange:
         first_column = df.iloc[:, 0]
         sec_column = df.iloc[:, 1]
-        df.iloc[:, 0] = first_column * v_speed_tip
-        df.iloc[:, 1] = sec_column * np.cos(np.deg2rad(angle_tip_new + angle_cone + axial_inclination))
+        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:
         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
         df.iloc[:, 1] = sec_column * np.cos(np.deg2rad(angle_root + angle_cone + axial_inclination))
 
-    avg_tip = result_avg_tip.iloc[:, 0]
-    result_avg_tip.iloc[:, 0] = avg_tip * v_speed_tip
+    avg_flange = result_avg_flange.iloc[:, 0]
+    result_avg_flange.iloc[:, 0] = avg_flange * v_speed_flange
     avg_root = result_avg_root.iloc[:, 0]
     result_avg_root.iloc[:, 0] = avg_root * v_speed_root
 
-    twist_1 = round(np.abs(pitch_angle_root[0] - pitch_angle_tip[0]), 2)
-    twist_2 = round(np.abs(pitch_angle_root[1] - pitch_angle_tip[1]), 2)
-    twist_3 = round(np.abs(pitch_angle_root[2] - pitch_angle_tip[2]), 2)
+    pitch_angle_flange = [0, 0, 0]
+    twist_1 = round(np.abs(pitch_angle_root[0] - pitch_angle_flange[0]), 2)
+    twist_2 = round(np.abs(pitch_angle_root[1] - pitch_angle_flange[1]), 2)
+    twist_3 = round(np.abs(pitch_angle_root[2] - pitch_angle_flange[2]), 2)
     twist_avg = round((twist_1 + twist_2 + twist_3) / 3, 2)
 
     sampling_num = int(0.015 * sampling_fq_1)
-    data_tip.iloc[:, 0] = data_tip.iloc[:, 0] / 5000000
+    data_flange.iloc[:, 0] = data_flange.iloc[:, 0] / 5000000
     data_root.iloc[:, 0] = data_root.iloc[:, 0] / 5000000
     lowpass_data.iloc[:, 0] = lowpass_data.iloc[:, 0] / 5000000
 
@@ -238,8 +276,8 @@ def data_analyse(path: List[str]):
     json_output = {
         'original_plot': {
             'blade_tip': {
-                'xdata': data_tip.iloc[:, 0].tolist()[::sampling_num],
-                'ydata': data_tip.iloc[:, 1].tolist()[::sampling_num]
+                'xdata': data_flange.iloc[:, 0].tolist()[::sampling_num],
+                'ydata': data_flange.iloc[:, 1].tolist()[::sampling_num]
             },
             'blade_root': {
                 'xdata': data_root.iloc[:, 0].tolist()[::sampling_num],
@@ -266,20 +304,20 @@ def data_analyse(path: List[str]):
         },
         'blade_tip': {
             'first_blade': {
-                'xdata': result_line_tip[0].iloc[:, 0].tolist(),
-                'ydata': result_line_tip[0].iloc[:, 1].tolist()
+                'xdata': result_line_flange[0].iloc[:, 0].tolist(),
+                'ydata': result_line_flange[0].iloc[:, 1].tolist()
             },
             'second_blade': {
-                'xdata': result_line_tip[1].iloc[:, 0].tolist(),
-                'ydata': result_line_tip[1].iloc[:, 1].tolist()
+                'xdata': result_line_flange[1].iloc[:, 0].tolist(),
+                'ydata': result_line_flange[1].iloc[:, 1].tolist()
             },
             'third_blade': {
-                'xdata': result_line_tip[2].iloc[:, 0].tolist(),
-                'ydata': result_line_tip[2].iloc[:, 1].tolist()
+                'xdata': result_line_flange[2].iloc[:, 0].tolist(),
+                'ydata': result_line_flange[2].iloc[:, 1].tolist()
             },
             'avg_blade': {
-                'xdata': result_avg_tip.iloc[:, 0].tolist(),
-                'ydata': result_avg_tip.iloc[:, 1].tolist()
+                'xdata': result_avg_flange.iloc[:, 0].tolist(),
+                'ydata': result_avg_flange.iloc[:, 1].tolist()
             }
         },
         'blade_root': {
@@ -455,17 +493,17 @@ 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
 
-    
-    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 - 31) | (
-            data_group['distance'] > mean_values * 1.1), 'distance'] = np.nan
-
-
+    if 0 in distance_counts.nlargest(3).index:
+        pass
+    else:
+        # 选择频率最大的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 - 30) | (
+                data_group['distance'] > mean_values * 1.1), 'distance'] = np.nan
     
     data_group['distance'] = data_group['distance'].fillna(method='ffill')
     filtered_data = data_group
-
     
     filtered_data['distance_diff'] = filtered_data['distance'].diff()
     large_diff_indices = filtered_data[filtered_data['distance_diff'] > min_distance].index
@@ -480,7 +518,7 @@ def cycle_calculate(data_group: pd.DataFrame, noise_threshold: float, min_distan
         
         current_distance = filtered_data.loc[idx, 'distance']
 
-        next_rows_large = filtered_data.loc[idx - 201: idx - 1]
+        next_rows_large = filtered_data.loc[idx - 500: idx - 1]
 
         
         if next_rows_large['distance'].le(current_distance - min_distance).all():
@@ -491,12 +529,15 @@ def cycle_calculate(data_group: pd.DataFrame, noise_threshold: float, min_distan
         
         current_distance = filtered_data.loc[idx - 1, 'distance']
 
-        next_rows_small = filtered_data.iloc[idx: idx + 201]
+        next_rows_small = filtered_data.iloc[idx: idx + 500]
 
         
         if next_rows_small['distance'].le(current_distance - min_distance).all():
             
             start_points = pd.concat([start_points, filtered_data.loc[[idx]]])
+            
+    if 0 in distance_counts.nlargest(3).index:
+        end_points, start_points = start_points, end_points
 
     if end_points.iloc[0, 0] < start_points.iloc[0, 0]:
         end_points = end_points.drop(end_points.index[0])
@@ -693,7 +734,7 @@ def coordinate_normalize(tip_border_rows: List[pd.DataFrame], tip_angle):
     tip_angle_list = []
     for turbine in tip_border_rows:
         tip_angle_cal0 = ((np.sin(tip_angle1) * turbine['distance'] - 0.07608) /
-                          (np.cos(tip_angle1) * turbine['distance']))
+                             (np.cos(tip_angle1) * turbine['distance']))
         tip_angle_cal = np.arctan(tip_angle_cal0)
         turbine['distance'] = (turbine['distance'] ** 2 + 0.0057881664 -
                                0.15216 * turbine['distance'] * np.sin(tip_angle1)) ** 0.5
@@ -706,6 +747,181 @@ def coordinate_normalize(tip_border_rows: List[pd.DataFrame], tip_angle):
     return tip_border_rows, tip_angle_new1
 
 
+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']))
+    flange_angle_cal = np.arctan(flange_angle_cal0)
+
+    flange_cen_row['中心y'] = (flange_cen_row['中心y'] ** 2 + 0.0057881664 -
+                                     0.15216 * flange_cen_row['中心y'] * np.sin(flange_angle1)) ** 0.5
+
+    flange_angle_new = float(flange_angle_cal)
+    flange_angle_new1 = np.rad2deg(flange_angle_new)
+
+    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 fit_circle(df, v_bounds=(2, 10), 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)
+            if wl < 3:
+                return y
+            if wl % 2 == 0:
+                wl -= 1
+            return savgol_filter(y, wl, polyorder)
+
+        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)]
+
+        def residuals_sq(params):
+            v, xc, yc, R = params
+            if v <= 0 or 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,
+            strategy='rand2bin',
+            mutation=(0.8, 1.2),
+            recombination=0.8,
+            popsize=30,
+            maxiter=1000,
+            polish=False,
+            seed=42,
+            workers=1
+        )
+
+        # 多候选点精修
+        pop = result.population
+        energies = result.population_energies
+        idx = np.argsort(energies)[:top_k]
+        candidates = pop[idx]
+
+        best_rmse = np.inf
+        best_result = None
+
+        for cand in candidates:
+            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]),
+                method='trf',
+                loss='linear',
+                max_nfev=50000,
+                xtol=1e-12,
+                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)
+            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]
+
+        result_df = pd.DataFrame([best_result],
+                                 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)
+    angle_speed = (np.pi / full_cycle) * 5000000
+    turbines = [pd.DataFrame() for _ in range(3)]
+
+    for i in range(0, len(start_times), 2):
+
+
+        start_time = start_times[i]
+        end_time = start_times[i + 1]
+        segment = data_group[(data_group['time'] > start_time) & (data_group['time'] <= end_time)]
+        if segment is None or segment.empty:
+            raise ValueError("Segment is empty")
+
+        segment = segment.copy()
+        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_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))
+
+        process_df = pd.DataFrame()
+        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)
+        processed_df = process_df[(process_df['time'] >= lower_bound) & (process_df['time'] <= upper_bound)]
+        blade_cen_est = fit_circle(processed_df)
+
+        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["中心y"] = result_df["中心y"] / np.cos(np.deg2rad(horizon_angle))
+
+    return result_df
+
+
+
 def radius_cal(border_rows, meas_angle, cen_dist, cen_angle, angle_main, angle_rotate):
 
     aero_dist = (pd.concat([df['distance'] for df in border_rows]).mean())
@@ -715,32 +931,58 @@ 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):
+                                   tower_dist: float, v_angle: float):
 
-
-    v_speed = 2 * np.pi * radius / full_cycle  
-    pitch_angle_list = []
+    v_speed = 2 * np.pi * radius / full_cycle  # 叶片线速度m/(1计时器单位）
     aero_dist_list = []
     cen_blade = []
     for turbine in border_rows:
-        diff_time = turbine.iloc[1, 0] - turbine.iloc[0, 0]
 
-        diff_len = (turbine.iloc[1, 1] - turbine.iloc[0, 1]) * np.cos(np.deg2rad(v_angle))
         mean_col2 = (turbine.iloc[1, 1] + turbine.iloc[0, 1]) / 2
         aero_dist = abs(mean_col2 - tower_dist) * np.cos(np.deg2rad(v_angle))
 
-        pitch_angle = np.degrees(np.arctan(diff_len / (diff_time * v_speed)))
-        pitch_angle_list.append(pitch_angle)
         aero_dist_list.append(aero_dist)
         cen_blade.append(mean_col2)
+
+    aero_dist_list.append(np.mean(aero_dist_list))
+    aero_dist_list = [round(num, 2) for num in aero_dist_list]
+
+    return aero_dist_list, v_speed, cen_blade
+
+def blade_angle(border_rows: List[pd.DataFrame], cen_data: pd.DataFrame, radius: float, full_cycle: int,
+                                   v_angle: float):
+
+    v_speed = 2 * np.pi * radius / full_cycle
+
+    values = []
+    for df in border_rows:
+        if df.shape[0] >= 2 and df.shape[1] >= 2:
+            values.append(df.iloc[0, 1])
+            values.append(df.iloc[1, 1])
+
+    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
+
+    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)))
+        pitch_angle = np.degrees(np.arctan(diff_len / diff_time))
+        pitch_angle_list.append(pitch_angle)
+
     pitch_mean = np.mean(pitch_angle_list)
     pitch_angle_list = [angle - pitch_mean for angle in pitch_angle_list]
     pitch_angle_list.append(max(pitch_angle_list) - min(pitch_angle_list))
-    aero_dist_list.append(np.mean(aero_dist_list))
     pitch_angle_list = [round(num, 2) for num in pitch_angle_list]
-    aero_dist_list = [round(num, 2) for num in aero_dist_list]
 
-    return pitch_angle_list, aero_dist_list, v_speed, cen_blade
+    return pitch_angle_list, v_speed
 
 
 def find_param(path: str):

frequency_filter.py

@@ -6,6 +6,10 @@ from scipy.signal import butter, filtfilt
 from scipy.signal import welch
 from scipy.fft import fft, fftfreq
 import numpy as np
+from sympy import symbols, solve, Eq, simplify, nsolve
+import itertools
+import multiprocessing
+from concurrent.futures import ThreadPoolExecutor
 
 
 warnings.filterwarnings("ignore", category=SettingWithCopyWarning)
@@ -147,3 +151,151 @@ def process_fft(data, fs):
     '''
 
     return segment, fft_x, fft_y_scaled, round(fft_x[max_index], 2), round(max_value, 2)
+
+
+def process_equations(list_param):
+    """
+    处理方程组并求解
+
+    参数:
+    a, b, p: 坐标点列表，P是轮毂中心
+    d, e, f: 轴向倾角、相位角、锥角
+
+    返回:
+    solutions: 包含所有有效解的列表
+    positive_angles: 包含所有正角度值的列表
+    """
+
+    a, b, p, d, e, f = [list_param[i] for i in range(6)]
+
+    # 定义符号
+    ya, yb, za, zb, x, y = symbols('ya yb za zb x y', real=True)
+
+    # 计算全局表达式
+    za_expr = solve(Eq(a[0] - p[0] + ya * (a[1] - p[1]) + za * (a[2] - p[2]), 0), za)[0]
+    zb_expr = solve(Eq(b[0] - p[0] + yb * (b[1] - p[1]) + zb * (b[2] - p[2]), 0), zb)[0]
+
+    # 解线性方程组得到 x_expr, y_expr
+    sol_xy = solve([Eq(x + ya * y + za_expr * d, 0), Eq(x + yb * y + zb_expr * d, 0)], [x, y], dict=True)
+    if not sol_xy:
+        raise RuntimeError("解线性方程组求 x,y 失败（意外），请检查输入符号。")
+    x_expr = simplify(sol_xy[0][x])
+    y_expr = simplify(sol_xy[0][y])
+
+    # 最终方程
+    f1 = simplify(Eq(1 + ya * yb + za_expr * zb_expr + e * (1 + ya ** 2 + za_expr ** 2) ** 0.5 * (
+                1 + yb ** 2 + zb_expr ** 2) ** 0.5, 0).lhs)
+    f2 = simplify(Eq(x_expr ** 2 + y_expr ** 2 + d ** 2 - 1, 0).lhs)
+
+    # 存储正角度值的列表
+    positive_angles = []
+
+    # 处理单个解的函数
+    def process_solution(args):
+        guess, index, total, solutions, positive_angles = args
+        try:
+            # 尝试求解
+            solution = nsolve([f1, f2], [ya, yb], guess, dict=True, prec=6)
+
+            # 获取解的值
+            ya_val = solution[0][ya]
+            yb_val = solution[0][yb]
+
+            # 计算方程残差（误差）
+            f1_error = abs(f1.subs({ya: ya_val, yb: yb_val}))
+            f2_error = abs(f2.subs({ya: ya_val, yb: yb_val}))
+            err_tol = 1e-1
+
+            is_new = True
+            if f1_error > err_tol or f2_error > err_tol:
+                is_new = False
+
+            # 检查解是否已经存在（允许小的数值误差）
+            for existing_sol in solutions:
+                if (abs(solution[0][ya] - existing_sol[ya]) < 1e-3 and
+                        abs(solution[0][yb] - existing_sol[yb]) < 1e-3):
+                    is_new = False
+                    break
+
+            if is_new:
+                solution_dict = solution[0]
+                solutions.append(solution_dict)
+
+                # 计算其他值
+                ya_val = simplify(solution_dict[ya])
+                yb_val = simplify(solution_dict[yb])
+                x_val = simplify(x_expr.subs({ya: ya_val, yb: yb_val}))
+                y_val = simplify(y_expr.subs({ya: ya_val, yb: yb_val}))
+                za_val = simplify(za_expr.subs({ya: ya_val}))
+                zb_val = simplify(zb_expr.subs({yb: yb_val}))
+
+                angle = np.degrees(np.arctan(float(y_val) / float(x_val)))
+
+                # 检查角度是否为正
+                if angle > 0:
+                    positive_angles.append(angle)
+
+                print(f"初始值 {guess} 的解:")
+                print(f"  x  = {x_val}")
+                print(f"  y  = {y_val}")
+                print(f"  方程 f1 残差: {f1_error}")
+                print(f"  方程 f2 残差: {f2_error}")
+                print(f"  angle = {angle}°")
+                print("-" * 50)
+                return 1
+            else:
+                return 0
+        except Exception as e:
+            if index % 1000 == 0:
+                print(f"已完成 {index} 个")
+            return 0
+
+
+    # 设置Windows下的多进程启动方法
+    multiprocessing.freeze_support()
+
+    # 存储结果的列表
+    manager = multiprocessing.Manager()
+    solutions = manager.list()
+    positive_angles = manager.list()
+
+    # 初始猜测值
+    initial_values = np.arange(-10, 10, 0.1)
+    initial_guesses = list(itertools.product(initial_values, repeat=2))
+
+    # 使用进程池并行处理
+    print(f"总共尝试 {len(initial_guesses)} 个初始组合...")
+    num_processes = int(multiprocessing.cpu_count() - 1)
+
+    with ThreadPoolExecutor(max_workers=num_processes) as executor:
+        args_list = [(guess, i, len(initial_guesses), solutions, positive_angles)
+                    for i, guess in enumerate(initial_guesses)]
+        results = list(executor.map(process_solution, args_list))
+
+    # 统计结果
+    success_count = sum(results)
+    print(f"\n计算完成:")
+    print(f"- 成功找到解: {success_count} 个")
+    print(f"- 总共找到 {len(solutions)} 个不同的解")
+    print(f"- 找到 {len(positive_angles)} 个正角度值")
+
+    # 打印符号求解得到的 (ya,yb) 解集
+    if solutions:
+        print("\n符号求解得到的 (ya,yb) 解集：")
+        for idx, sol in enumerate(solutions, 1):
+            ya_val = simplify(sol[ya])
+            yb_val = simplify(sol[yb])
+            x_val = simplify(x_expr.subs({ya: ya_val, yb: yb_val}))
+            y_val = simplify(y_expr.subs({ya: ya_val, yb: yb_val}))
+            za_val = simplify(za_expr.subs({ya: ya_val}))
+            zb_val = simplify(zb_expr.subs({yb: yb_val}))
+            angle = np.degrees(np.arctan(float(y_val) / float(x_val)))
+
+            print(f"解 {idx}:")
+            print(f"  x  = {x_val}")
+            print(f"  y  = {y_val}")
+            print(f"  angle = {angle}°")
+            print("-" * 50)
+
+    return list(positive_angles)
+