AUTUMN的自动化技术分享
0%

动态解算打分

Posted on Edited on

使用背景

在动态解算中，GNSS基站通过ntrip传输差分数据给GNSS测站，测站通过接收到的差分数据在板卡上进行前端解算，并通过串口输出解算后的GGA数据。本脚本通过输入GNSS测站输出的POS文件，可以计算出基线的内附和误差以及真误差。

调用库

1
2
3
4
from tkinter import filedialog
import xlwt
import math
import numpy

角度格式转化函数

将“度度分分.分分分分”的数据格式转化为“度.度度度度”。

1
2
3
4
5
6
7
def NEtoangle(blh):
    """数据格式为ddmm.mmmm"""
    data = float(blh)
    du = data // 100
    fen = data - du * 100
    result = du + fen / 60
    return result

大地坐标转化空间直角坐标函数

使用WGS84椭球参数计算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def BLH2XYZ(B, L, H):
    """
     该函数把某点的大地坐标(B, L, H)转换为空间直角坐标（X, Y, Z).
    :param B:  大地纬度, 角度制, 规定南纬为负，范围为[-90, 90]
    :param L:  大地经度, 角度制, 规定西经为负, 范围为[-180, 180]
    :param H:  海拔，大地高, 单位 m
    :param a:  地球长半轴，即赤道半径，单位 m
    :param b:  地球短半轴，即大地坐标系原点到两级的距离, 单位 m
    :return:   X, Y, Z, 空间直角坐标, 单位 m
    """
    a = 6378137.0  # 参考椭球的长半轴, 单位 m
    b = 6356752.31414  # 参考椭球的短半轴, 单位 m

    sqrt = math.sqrt
    sin = math.sin
    cos = math.cos

    B = angle2rad(B)  # 角度转为弧度
    L = angle2rad(L)  # 角度转为弧度

    e = sqrt((a ** 2 - b ** 2) / (a ** 2))  # 参考椭球的第一偏心率
    N = a / sqrt(1 - e * e * sin(B) * sin(B))  # 卯酉圈半径, 单位 m

    X = (N + H) * cos(B) * cos(L)
    Y = (N + H) * cos(B) * sin(L)
    Z = (N * (1 - e * e) + H) * sin(B)
    return X, Y, Z  # 返回空间直角坐标(X, Y, Z)

角度转弧度

1
2
3
4
5
6
7
8
def angle2rad(a):
    """
    该函数可以实现角度到弧度的转换.
    :param a:  角度
    :return:  r, 对应的弧度
    """
    r = a * math.pi / 180.0
    return r

真误差计算函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
def real_solution(b, l, h, b1, l1, h1):
    level_1_accuracy = 0
    level_2_accuracy = 0
    level_3_accuracy = 0
    level_4_accuracy = 0
    pass_group = 0
    alldatafirst, alldatasecond, alldatathird = [], [], []
    # 基准站实际坐标
    X, Y, Z = BLH2XYZ(b, l, h)
    # 测站实际坐标
    X1, Y1, Z1 = BLH2XYZ(b1, l1, h1)

    # 基线真实长度(mm)
    real_distance = math.sqrt((X - X1) ** 2 + (Y - Y1) ** 2 + (Z - Z1) ** 2) * 1000
    # 基线在水平方向上的长度(mm)
    real_distance_horizon = math.sqrt((X - X1) ** 2 + (Y - Y1) ** 2) * 1000

    path = filedialog.askopenfilename()
    with open(path, 'r') as f:
        lines = f.readlines()
        line = [line for line in lines if line.strip() and "RECV ASCII" not in line]
        total = len(line) - 19
        # 对文件每一行进行分析
        for index in range(len(line) - 19):
            success = 0
            pass_status = False
            for j in range(20):
                numbers = line[index + j].strip().split(',')
                if numbers[6] == '4':  # 根据固定位的值填写
                    success += 1
            if success >= 19:  # 19
                # 通过的掩膜数量
                pass_group += 1
                # 判断当前掩膜是否通过
                pass_status = True
            if pass_status and index % 20 == 0:
                for k in range(20):
                    numbers = line[index + k].strip().split(',')
                    alldatafirst.append(NEtoangle(numbers[2]))
                    alldatasecond.append(NEtoangle(numbers[4]))
                    alldatathird.append(float(numbers[9]) + float(numbers[11]))
    # 计算坐标
    xyz = [BLH2XYZ(b, l, h) for b, l, h in zip(alldatafirst, alldatasecond, alldatathird)]
    # 计算测量出的基线长度
    measure_distance = [math.sqrt((x - X) ** 2 + (y - Y) ** 2 + (z - Z) ** 2) * 1000 for x, y, z in xyz]
    # 计算测量出的基线在水平方向的长度
    measure_distance_horizon = [math.sqrt((x - X) ** 2 + (y - Y) ** 2) * 1000 for x, y, z in xyz]

    # 测量出的基线真误差
    delta_distance = [x - real_distance for x in measure_distance]
    # 测量出的基线在水平方向上的真误差
    delta_distance_horizon = [x - real_distance_horizon for x in measure_distance_horizon]

    # 测量出的点和基站的高程真误差
    deltaH = [(z - Z1) * 1000 for x, y, z in xyz]

    pass_rate = pass_group / total
    # 创建一个新的Excel工作簿
    workbook = xlwt.Workbook(encoding='utf-8')
    # 新建sheet
    sheet1 = workbook.add_sheet("真误差")
    sheet1.write(0, 1, "掩膜通过率")
    sheet1.write(0, 2, f"{pass_rate * 100}%")
    sheet1.write(1, 1, "基线长真误差（mm）")
    sheet1.write(1, 2, "基线长水平方向真误差（mm）")
    sheet1.write(1, 3, "高程真误差（mm）")
    for i in range(1, len(delta_distance) + 1):
        sheet1.write(i + 1, 1, delta_distance[i - 1])
        sheet1.write(i + 1, 2, delta_distance_horizon[i - 1])
        sheet1.write(i + 1, 3, deltaH[i - 1])
    save_path = r"D:\9999868真误差报告.xls"
    workbook.save(save_path)
    print(f"解算成功!\n文件已经保存为：{save_path}")

内附和精度计算函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
def inertia_solution():
    level_1_accuracy = 0
    level_2_accuracy = 0
    level_3_accuracy = 0
    level_4_accuracy = 0
    alldatafirst, alldatasecond, alldatathird = [], [], []
    path = filedialog.askopenfilename()
    with open(path, 'r') as f:
        lines = f.readlines()
        line = [line for line in lines if line.strip() and "RECV ASCII" not in line]
        # 对文件每一行进行分析
        for index in range(len(line)):
            numbers = line[index].strip().split(',')
            alldatafirst.append(NEtoangle(numbers[2]))
            alldatasecond.append(NEtoangle(numbers[4]))
            alldatathird.append(float(numbers[9]) + float(numbers[11]))

    averagefirst = numpy.mean(alldatafirst)
    averagesecond = numpy.mean(alldatasecond)
    averagethird = numpy.mean(alldatathird)
    # 基准站实际坐标
    X, Y, Z = BLH2XYZ(31.159529544, 121.178237368, 41.326345761)
    # 测站实际坐标(内附和)
    X1, Y1, Z1 = BLH2XYZ(averagefirst, averagesecond, averagethird)

    # 基线真实长度(mm)
    real_distance = math.sqrt((X - X1) ** 2 + (Y - Y1) ** 2 + (Z - Z1) ** 2) * 1000
    # 基线在水平方向上的长度(mm)
    real_distance_horizon = math.sqrt((X - X1) ** 2 + (Y - Y1) ** 2) * 1000

    # 计算坐标
    xyz = [BLH2XYZ(b, l, h) for b, l, h in zip(alldatafirst, alldatasecond, alldatathird)]
    # 计算测量出的基线长度
    measure_distance = [math.sqrt((x - X) ** 2 + (y - Y) ** 2 + (z - Z) ** 2) * 1000 for x, y, z in xyz]
    # 计算测量出的基线在水平方向的长度
    measure_distance_horizon = [math.sqrt((x - X) ** 2 + (y - Y) ** 2) * 1000 for x, y, z in xyz]
    # 测量出的基线真误差
    delta_distance = [x - real_distance for x in measure_distance]
    # 测量出的基线在水平方向上的真误差
    delta_distance_horizon = [x - real_distance_horizon for x in measure_distance_horizon]

    # 测量出的点和基站的高程真误差
    deltaH = [(z - Z1) * 1000 for x, y, z in xyz]

    # 创建一个新的Excel工作簿
    workbook = xlwt.Workbook(encoding='utf-8')
    # 新建sheet
    sheet1 = workbook.add_sheet("真误差")
    sheet1.write(1, 1, "基线长真误差（mm）")
    sheet1.write(1, 2, "基线长水平方向真误差（mm）")
    sheet1.write(1, 3, "高程真误差（mm）")
    for i in range(1, len(delta_distance) + 1):
        sheet1.write(i + 1, 1, delta_distance[i - 1])
        sheet1.write(i + 1, 2, delta_distance_horizon[i - 1])
        sheet1.write(i + 1, 3, deltaH[i - 1])
    save_path = r"D:\966内附和报告.xls"
    workbook.save(save_path)
    print(f"解算成功!\n文件已经保存为：{save_path}")

主函数

将data中的1~4号点的坐标修改为实际点的空间直角坐标。其中高程为椭球高。

1
2
3
4
5
6
7
8
9
10
11
12
13
data = {
    "1": (31.159529544, 121.178237368, 41.326345761),
    "2": (31.159482390, 121.178156434, 41.381853455),
    "3": (31.159473582, 121.178158438, 41.375060047),
    "4": (31.159427102, 121.178077343, 41.336077423)
}
base_id = input("请输入基站编号：")
rover_id = input("请输入测站编号：")
base_coords = data[base_id]
rover_coords = data[rover_id]

real_solution(base_coords[0], base_coords[1], base_coords[2], rover_coords[0], rover_coords[1], rover_coords[2])
# inertia_solution()