石家庄铁道大学测量平差课程设计-测角中点三边形-CSDN博客

本文链接：https://blog.csdn.net/m0_60292965/article/details/122092160

刚刚做完平差课设，供学弟学妹们参考，欢迎指正，提出建议。

基于程序目的的，我选择了“测角中点三边形”来作为我的课程设计面向的模型，通过两个已知控制点和观测得到的角度计算未知点的坐标及某些边的长度和坐标方位角。

为了更好的了解所设计的软件，接下来是对软件内容的描述：

数据输入：

考虑到每次作业数据量的不同，数据输入提供了三种方式：手动输入，Excel文件导入，Txt文件导入。在程序主界面有输入框，用户参考我所定义的数据格式即可输入。主界面的菜单栏的数据输入项中可选择“.csv”格式的Excel文件和“.txt”的文本文件进行数据导入。

计算：

软件界面中下方为“计算”按钮，数据输入完成后，点击即可进行计算。计算方法采用的是条件平差模型，根据输入得到的观测值和已知控制点坐标，程序自动计算出改正数和平差值，并以文字数字的形式呈现在主界面的右侧窗口，详细的计算过程请看第三部分。

数据输出：

计算完毕后，用户可根据自己的需求来选择是否进行数据导出。软件提供了两种数据输出的格式“.csv”Excel文件和“.txt”文本文件来实现数据的快速保存，以便保持数据的完备性。

使用说明：

为用户提供软简要的件使用说明书。

测量模型可视化：

为了更直观地展示平差结果的，软件提供了绘图功能，点击“计算按钮”后自动触发该功能。

平差原理

此程序针对的模型为”测角中点三边形”，应用的平差模型为条件平差，具体的计算流程如下：

根据测量模型，列出条件方程： AQAT+W=0

Li=Li+V

根据测角网的几何条件，我们可以根据图形条件，圆周条件，极条件列出五个条件方程。

a1+b1+c1-1800=0

a2+b2+c2-1800=0

a3+b3+c3-1800=0

c1+c2+c3-3600=0

sina1sinb1sina2sinb2sina3sinb3=1

由于根据极条件所列出的方程是非线性的，我们要对其进行线性化处理，最终得到的方程为cota1νa1+cota2νa2+cota3νa3-cotb1νb1-cotb2νb2-cotb3νb3+(1 - sinb1sinb2sinb3sina1sina2sina3)ρ''=0

列出法方程并求解K Naa=AQAT

NaaK+W=0

K=-Naa-1W

将K代入改正数方程，求出V值，并求出平差值

V=QATK

Li=Li+v

用平差值重新列出条件方程进行检验

本程序涉及到的矩阵矩阵均采用Python的科学计算库Numpy进行实现，为避免因为计算而损失精度，计算过程中所有数值保留8位有效数字，输出时根据“四舍六入”进行取舍。

精度评定

σ=νTPvr

φ=f1L1+f2L2+…+fnLn

Qφφ=fTcot-AQfTNaa-1AQf

σφ=σ Qφφ（平差值函数中误差）

'''本软件版作者张斌，于2021年11月17日完成第一次测试。
禁止抄袭我的代码！！！
规定：
（1）数据输入格式：
允许的测角中误差
A,x,y
B,x,y
设站点，瞄准方向，瞄准方向，度.分秒（秒保留三位小数，例如120°50′48.6″“120.50486”）
（2）数据输入和输出的文件均为*.txt和*.csv格式

'''

import math
from tkinter import *
import numpy as np
from tkinter import filedialog
import tkinter.messagebox
import time
import turtle as t

def tranangle(a):  #定义角度转弧度
    x = a.split('.')
    x1 = int(x[1][0:2])
    x2 = "{}.{}".format(x[1][2:4], x[1][4])
    x2 = float(x2)
    x=(x2 / 60 / 60 + x1 / 60 + int(x[0])) / 180 * math.pi
    if x>=2*math.pi:
        x-=2*math.pi
    if x<0.0:
        x+=2*math.pi
    return x
def radtodms(a):
    a = a / math.pi * 180
    b = []
    b.append(int(a))
    b.append(int((a - b[0]) * 60))
    b.append(float("%.1f" % (((a - b[0]) * 60 - int(b[1])) * 60)))
    if (b[2]==60.0):
        b[2]-=60
        b[1]+=1
        if(b[1]==60.0):
            b[1]-=60
            b[0]+=60
   # print(b)
    return b

def datainputtxt():
    try:
        root = Tk()
        root.withdraw()

        filepath = filedialog.askopenfilename()

        f = open(filepath, 'r')
        content = f.readlines()
        for i in content:
            userinput.insert('end', i)
        userinput.insert(END,"\n")
        tkinter.messagebox._show('Attention', "数据导入成功!")
    except:
        tkinter.messagebox.showerror('错误', "数据导入失败!")
def datainputexcel():
    try:
        root = Tk()
        root.withdraw()

        filepath = filedialog.askopenfilename()

        f = open(filepath, 'r')
        content = f.readlines()
        userinput.insert('end',content[0][:5])
        userinput.insert('end',"\n")
        userinput.insert('end', content[1][:9])
        userinput.insert('end', "\n")
        userinput.insert('end', content[2][:9])
        userinput.insert('end', "\n")
        for i in range(3,len(content)):
            userinput.insert('end', content[i])
        userinput.insert(END, "\n")
        tkinter.messagebox._show('Attention', "数据导入成功!")
    except:
        tkinter.messagebox.showerror('错误', "数据导入失败!")
def calculatesum(a):
    sum=0
    for i in range(len(a)):
        sum+=tranangle(a[i].split(',')[3])
    return sum
def gettime():
    try:
        timestr = time.strftime("%H:%M:%S")  # 获取当前的时间并转化为字符串
        lb.configure(text=timestr)  # 重新设置标签文本
        root.after(1000, gettime)  # 每隔1s调用函数 gettime 自身获取时间
    except:
        pass
def judgeA(x):
    if x<0.0:
        x+=2*math.pi
    if x>=math.pi*2:
        x-=2*math.pi
    return x
def draw(Loc,p):
    t.screensize(10000, 10000)
    t.title("坐标图")
    t.speed(1000)
    t.write('(0,0)')
    t.goto(-5000, 0)
    t.fd(9000)
    t.write("y")
    t.penup()
    t.goto(0, 4000)
    t.pendown()
    t.write('x')
    t.right(90)
    t.fd(10000)

    for i in range(len(Loc) - 1):
        t.penup()
        t.goto(Loc[3][1], Loc[3][0])
        t.write('D')
        t.pendown()
        t.goto(Loc[i][1], Loc[i][0])
        t.write(p[i])
    for i in range(len(Loc) - 1):
        t.goto(Loc[i][1], Loc[i][0])
    t.penup()
    t.goto(9999, 9999)

def work():
    try:
        a = userinput.get('1.0', END)
        content = a.split('\n')
        for i in range(len(content)):  # 消除字符串末尾的换行符
            content[i] = content[i].replace("\n", '')
        A1 = []
        A2 = []
        Loc=[]
        f = content[0]
        A1.append(content[1])
        A2.append(content[2])


        ABD = []
        BCD = []
        ACD = []
        w = []  # 单位是″

        for i in range(len(content)):  # 按规定图形分成三个三角形
            if sorted(content[i].split(','))[1:] == ['A', 'B', 'D']:
                ABD.append(content[i])
            if sorted(content[i].split(','))[1:] == ['B', 'C', 'D']:
                BCD.append(content[i])
            if sorted(content[i].split(','))[1:] == ['A', 'C', 'D']:
                ACD.append(content[i])

        ABD.sort()
        ACD.sort()  # 将每个三角形的内角按规定顺序排列，便于后期计算
        BCD.sort()
        w.append(radtodms(calculatesum(ABD) - math.pi)[2])
        w.append(radtodms(calculatesum(BCD) - math.pi)[2])
        w.append(radtodms(calculatesum(ACD) - math.pi)[2])

        w.append(radtodms(tranangle(ABD[2].split(',')[3]) + tranangle(ACD[2].split(',')[3]) + tranangle(
            BCD[2].split(',')[3]) - 2 * math.pi)[2])  # ρ=206264.80264
        # 极条件的常数项
        w.append((1 - math.sin(tranangle(ABD[1].split(',')[3])) * math.sin(tranangle(ACD[0].split(',')[3])) * math.sin(
            tranangle(BCD[1].split(',')[3])) / math.sin(tranangle(ABD[0].split(',')[3])) / math.sin(
            tranangle(ACD[1].split(',')[3])) / math.sin(tranangle(BCD[0].split(',')[3]))) * 206264.80624)
        # 给W常数项项矩阵赋值
        W = np.ones((5, 1))
        for i in range(5):
            W[i][0] = w[i]
        # 协因数矩阵
        Q = np.eye(9)
        # 系数矩阵
        A = np.eye(5, 9)
        for i in range(5):
            A[i][i] = 0
        for i in range(3):  # 图形条件和圆周方程的系数
            A[0][i] = 1  # ABD
            A[1][i + 3] = 1  # BCD
            A[2][i + 6] = 1  # ACD
            A[3][3 * i + 2] = 1
        # 极条件的系数
        A[4][0] = 1.0 / math.tan(tranangle(ABD[0].split(',')[3]))
        A[4][3] = 1.0 / math.tan(tranangle(BCD[0].split(',')[3]))
        A[4][6] = 1.0 / math.tan(tranangle(ACD[1].split(',')[3]))
        A[4][1] = -1.0 / math.tan(tranangle(ABD[1].split(',')[3]))
        A[4][4] = -1.0 / math.tan(tranangle(BCD[1].split(',')[3]))
        A[4][7] = -1.0 / math.tan(tranangle(ACD[0].split(',')[3]))

        '''条件平差计算公式       精度评定
           Naa=A*Q*At           单位权中误差fz=(Vt*P*V)/r
           Naa*K+W=0
           V=Q*At*K'''
        At = A.T
        Naa = A @ Q @ At
        K = -np.linalg.inv(Naa) @ W
        V = Q @ At @ K
        #print(ABD, BCD, ACD)
        # 改正角度观测值
        newABD = []
        newBCD = []
        newACD = []
        # 精度评定
        fz = math.sqrt((V.T @ np.linalg.inv(Q) @ V) / 5)
        #print(fz)
        for i in range(3):
            newABD.append("{}{}".format(ABD[i][:6], radtodms(tranangle(ABD[i][6:]) + V[i][0] / 3600 / 180 * math.pi)))
            newBCD.append(
                "{}{}".format(BCD[i][:6], radtodms(tranangle(BCD[i][6:]) + V[i + 3][0] / 3600 / 180 * math.pi)))
            newACD.append(
                "{}{}".format(ACD[i][:6], radtodms(tranangle(ACD[i][6:]) + V[i + 6][0] / 3600 / 180 * math.pi)))

        # 计算坐标
        if A1[0].split(',')[0] in "ABD" and A2[0].split(',')[0] in "ABD":


            Loc.append((float(A1[0].split(',')[1]),float(A1[0].split(',')[2])))
            Loc.append((float(A2[0].split(',')[1]),float(A2[0].split(',')[2])))


            fx = float(A2[0].split(',')[1])-float(A1[0].split(',')[1])
            fy = float(A2[0].split(',')[2])-float(A1[0].split(',')[2])


            s1 = math.sqrt(fx ** 2 + fy ** 2)

            sbd=math.sin(tranangle(ABD[0].split(',')[3]))*s1/math.sin(tranangle(ABD[2].split(',')[3]))
            #print(s1)
            #print(sbd)
            if fy>0.0:
                A1A2 = math.atan(fx / fy)
            elif fy<0.0:
                A1A2 = math.atan(fx / fy)+math.pi
            else:

                if fx>0.0:
                    A1A2=math.pi/2
                if fx<0.0:
                    A1A2=math.pi


            A2A3 = judgeA(A1A2 - float(tranangle(ABD[1].split(',')[3])) + math.pi)
            yd = float(A2[0].split(',')[2]) + sbd * math.cos(A2A3)
            xd = float((A2[0].split(',')[1])) + sbd * math.sin(A2A3)
            Scd=s1*math.sin(tranangle(ABD[0].split(',')[3]))*math.sin(tranangle(BCD[0].split(',')[3]
                ))/math.sin(tranangle(ABD[2].split(',')[3]))/math.sin(tranangle(BCD[1].split(',')[3]))
            Acd=radtodms(judgeA(A2A3+float(tranangle(BCD[2].split(',')[3]))-math.pi))

            xc=xd+Scd*math.sin(judgeA(A2A3+float(tranangle(BCD[2].split(',')[3]))-math.pi))
            yc=yd+Scd*math.cos(judgeA(A2A3+float(tranangle(BCD[2].split(',')[3]))-math.pi))

            Loc.append((xc,yc))
            Loc.append((xd,yd))

            #精度评定
            fy=np.zeros((9,1))
            fy[0][0]=math.cos(tranangle(ABD[0].split(',')[3]))*math.cos(A2A3)/math.sin(tranangle(ABD[2].split(',')[3]))
            fy[1][0]=math.sin(A2A3)*math.sin(tranangle(ABD[0].split(',')[3]))/math.sin(tranangle(ABD[2].split(',')[3]))
            fy[2][0]=-math.sin(tranangle(ABD[0].split(',')[3]))*math.cos(A2A3)*math.cos(tranangle(ABD[2].split(',')[3]))/math.sin(tranangle(ABD[2].split(',')[3]))
            AQfy=A@Q@fy
            AQfyT=AQfy.T
            Qyd=(fy.T)@Q@fy-AQfyT@(np.linalg.inv(Naa))@A@Q@fy
            fyd=fz*math.sqrt(Qyd[0])*s1/206265  #单位是m

            fx=np.zeros((9,1))
            fx[0][0]=math.cos(tranangle(ABD[0].split(',')[3]))*math.sin(A2A3)/math.sin(tranangle(ABD[2].split(',')[3]))
            fx[1][0]=-math.cos(A2A3)*math.sin(tranangle(ABD[0].split(',')[3]))/math.sin(tranangle(ABD[2].split(',')[3]))
            fx[2][0]=-math.sin(A2A3)*math.sin(tranangle(ABD[0].split(',')[3]))*math.cos(tranangle(ABD[2].split(',')[3]))/math.sin(tranangle(ABD[2].split(',')[3]))
            AQfx = A @ Q @ fx
            AQfxT=AQfx.T
            Qxd = (fx.T) @ Q @ fx - AQfxT @ (np.linalg.inv(Naa)) @ A @ Q @ fx
            fxd = fz * math.sqrt(Qxd[0]) * s1 / 206265  # 单位是m
            #点位中误差
            fd=math.sqrt(fxd**2+fyd**2)

            #绘图模块
            p=['A','B','C']
            try:
                draw(Loc,p)
            except:
                draw(Loc,p)
            #print(Loc,type(Loc[0]))
            #print(radtodms(A1A2+tranangle(ABD[0].split(',')[3])+tranangle(ABD[2].split(',')[3])+tranangle(BCD[2].split(',')[3])-math.pi))
            #print(Acd)
            #print(xc,yc)
            #print(xd, yd)

        # 输出改正值
        endingoutput.insert('end', '-——————改正后的角度观测值及精度评定——————')
        endingoutput.insert('end', '\n\n')
        for i in range(len(newABD)):
            str = "{0}={1}°{2}′{3}″]".format(newABD[i][:5], newABD[i].split(',')[3], newABD[i].split(',')[4],
                                             newABD[i].split(',')[5][1:-1])
            endingoutput.insert('end', str)
            endingoutput.insert('end', '\n')

        for i in range(len(newBCD)):
            str = "{0}={1}°{2}′{3}″]".format(newBCD[i][:5], newBCD[i].split(',')[3], newBCD[i].split(',')[4],
                                             newBCD[i].split(',')[5][1:-1])
            endingoutput.insert('end', str)
            endingoutput.insert('end', '\n')
        for i in range(len(newACD)):
            str = "{0}={1}°{2}′{3}″]".format(newACD[i][:5], newACD[i].split(',')[3], newACD[i].split(',')[4],
                                             newACD[i].split(',')[5][1:-1])
            endingoutput.insert('end', str)
            endingoutput.insert('end', '\n')
        endingoutput.insert('end', '\n' )
        endingoutput.insert('end', "平差后计算的D点坐标为 ({:.4f},{:.4f})".format(xd, yd))
        endingoutput.insert('end', '\n' * 2)
        endingoutput.insert('end', "平差后计算的c点坐标为 ({:.4f},{:.4f})".format(xc, yc))
        endingoutput.insert('end', '\n' * 2)
        endingoutput.insert('end', "BD的距离为{:.3f}".format(sbd))
        endingoutput.insert('end', '\n' * 2)
        endingoutput.insert('end',"CD的距离为{:.3f}".format(Scd))
        endingoutput.insert('end', '\n' * 2)
        endingoutput.insert('end', "CD的坐标方位角为{}°{}′{}″".format(Acd[0],Acd[1],Acd[2]))
        endingoutput.insert('end', '\n' * 2)
        endingoutput.insert('end', "经平差后求得测角中误差fz={:.1f}″".format(fz))
        endingoutput.insert('end', '\n' * 2)
        endingoutput.insert('end', "经平差后求得D点x坐标的中误差fxd={:.1f}mm".format(fxd*1000))  #单位为mm
        endingoutput.insert('end', '\n' * 2)
        endingoutput.insert('end', "经平差后求得D点y坐标的中误差fyd={:.1f}mm".format(fyd * 1000))  # 单位为mm
        endingoutput.insert('end', '\n' * 2)
        endingoutput.insert('end', "经平差后求得D的点位中误差fd={:.1f}mm".format(fd*1000))  # 单位为mm
        endingoutput.insert('end', '\n'*2)
        try:
            if fz <= float(f.split(',')[0]):
                endingoutput.insert('end', '测角精度符合要求   ')
            else:
                endingoutput.insert('end', '测角精度超限   ')
            if fd <= float(f.split(',')[1]):
                endingoutput.insert('end', '点位中误差精度符合要求')
            else:
                endingoutput.insert('end', '点位中误差精度超限')
            endingoutput.insert('end', '\n' * 3)
        except:
            pass

    except:
        tkinter.messagebox.showerror('错误', "请输入正确的数据格式！")

def dataoutput():
    try:
        root = Tk()
        root.withdraw()
        filepath = filedialog.asksaveasfilename()
        f = open(filepath, 'w')
        content = endingoutput.get('1.0', END)
        content = content.split('\n')
        for i in content:
            f.write(i)
            f.write('\n')
        f.close()
        tkinter.messagebox._show('Attention', "数据导出成功!")
    except:
        tkinter.messagebox.showerror('错误','数据导出失败！')
def show():
    show=Tk()
    show.title("使用说明")
    show.geometry('520x600')
    s1=Label(show,text="规定：(Txt数据输入用英文逗号','隔开)\n"
                       "（1）数据输入格式：\n"
                       "    允许的测角中误差,允许的点位中误差\n"
                       "    A,x,y\n"
                       "    B,x,y\n"
                       "    设站点，瞄准方向，瞄准方向，度.分秒\n（秒保留三位小数，例如120°50′48.6″“120.50486”）\n"
                       "（2）数据输入和输出的文件均为*.txt和*.csv格式)\n"
                       "（3）严格遵守输入规定程序才能正常运行\n"
                       " (4) 使用本软件前需要搭载Pyhon 3.x的编译环境,\n"
                       "  并且提前安装tkinter,numpy,filedialog库\n"
                       "  from tkinter import *\n"
                       "  import numpy as np\n"
                       "  from tkinter import filedialog\n"
                       "  import tkinter.messagebox\n",fg='black',font="宋体",justify="left")
    s1.pack()
root = Tk()
root.title('张斌的平差课设')
root.geometry('800x610')
#创建菜单
menubar = Menu(root)  #顶级菜单

filemenu = Menu(menubar, tearoff=False)
filemenu.add_command(label='导入Excel文件', command=datainputexcel)
filemenu.add_command(label='导入Txt文件', command=datainputtxt)
menubar.add_cascade(label='数据输入', menu=filemenu)  # 创建级联菜单，menu选项指定下一级的菜单是什么

outmenu = Menu(menubar, tearoff=False)
outmenu.add_command(label='导出Excel文件', command=dataoutput)
outmenu.add_command(label='导出Txt文件', command=dataoutput)
menubar.add_cascade(label='数据输出', menu=outmenu)
menubar.add_cascade(label="使用说明",command=show)
root.config(menu=menubar)


#定义标签
label=Label(root,text='欢迎使用！！！！！\n\n土1903-4 张斌 20190788',fg='red',font=('宋体',15))
label.pack()

label1=Label(root,text="输入窗口")
label1.place(relx=0.23,rely=0.20)
label2=Label(root,text="输出窗口")
label2.place(relx=0.73,rely=0.20)
#定义输入窗口
userinput=Text(root)
userinput.place(relx=0.01,rely=0.25,width=388,height=400)
#定义计算按钮
B=Button(root,text='计算',width=6,height=2,command=work)
B.place(relx=0.47,rely=0.91)
#定义结果输出窗口
endingoutput=Text(root)
endingoutput.place(relx=0.505,rely=0.25,width=388,height=400)


try:
    lb = tkinter.Label(root, text='', fg='blue', font=("黑体", 20))
    lb.place(relx=0.8, rely=0.02)
    gettime()
except:
    pass


root.mainloop()