基线网平差（条件平差，间接平差，附有参数的条件平差，附有限制条件的间接平差）python3.7 自动生成矩阵并迭代平差

总述

先挖个坑 持续更新，目前已经实现了前三个
GayHub 有正常运行的案例

案例详情

一个由24条基线组成的GNSS网，共12个点，4个已知，data.csv如下

1,G01,G02,02912608_G01.zsd-02942607_G02.zsd,-366.6087,-252.4815,72.1291,8.1,11,13.7
2,G01,g03,02912608_G01.zsd-03662608_g03.zsd,-509.2608,22.3781,-426.4979,7.8,10.6,13.3
3,G01,g11,02942630_G01.zsd-03662630_g11.zsd,714.8437,446.3863,-70.7624,4.1,7.7,6.7
4,G02,g03,02942607_G02.zsd-03662608_g03.zsd,-142.6498,274.8577,-498.627,5.4,7.7,9.9
5,G04,g03,02942608_G04.zsd-03662609_g03.zsd,97.926,-192.0434,346.4421,5.3,6.9,5.8
6,G04,g05,02942625_G04.zsd-03662626_g05.zsd,312.0537,204.7549,-46.2603,13.4,38.8,18.6
7,g06,G04,02912631_g06.zsd-02942625_G04.zsd,-701.8827,-464.6907,110.2367,5.7,10.7,6.5
8,g06,g05,02912631_g06.zsd-03662626_g05.zsd,-389.8306,-259.9329,63.9769,5.6,5.9,4.1
9,g06,G07,02912612_g06.zsd-02942612_G07.zsd,348.1171,196.0223,-5.7602,3.3,4.8,3.9
10,g06,g09,02912612_g06.zsd-03662611_g09.zsd,167.4629,-202.2958,413.9846,2.9,4.3,3.3
11,G07,g09,02942612_G07.zsd-03662611_g09.zsd,-180.6533,-398.3194,419.7436,3.6,5.2,4.2
12,G07,g11,02942620_G07.zsd-03662620_g11.zsd,272.0238,-428.7511,818.176,3.7,9.2,4.6
13,G08,g03,02912609_G08.zsd-03662609_g03.zsd,-533.794,-312.1199,27.573,3.4,4.9,3.4
14,G08,G04,02912609_G08.zsd-02942608_G04.zsd,-631.7205,-120.0778,-318.8694,5.1,6.6,5.6
15,G08,g09,02942622_G08.zsd-03662622_g09.zsd,237.6263,142.3337,-15.1263,4.8,7,4.1
16,G10,G07,02912620_G10.zsd-02942620_G07.zsd,-242.2288,144.2423,-388.9749,3.3,6.6,4.1
17,G10,g11,02912620_G10.zsd-03662620_g11.zsd,29.795,-284.5064,429.2001,3.6,9,4.6
18,G12,G01,02912632_G12.zsd-02942630_G01.zsd,-390.41,-259.3411,63.9442,4.1,8.1,6.8
19,G12,G08,02912623_G12.zsd-02942622_G08.zsd,-365.8394,75.1172,-390.1562,6.1,8.8,4.9
20,G12,g09,02912623_G12.zsd-03662622_g09.zsd,-128.2126,217.4485,-405.2823,5.5,8.7,4.3
21,G12,g11,02912632_G12.zsd-03662630_g11.zsd,324.4349,187.0474,-6.82,3.7,8.6,5.7

间接平差

近似坐标手算的，也可以机器算，用的到在github里面的那个求路径的函数，我估计不会去花时间完善了

GNSS控制网类，主要作用是初始化P、X0，以及方便整个网内基线与点的查询

class station(object):
    '''undirected unweighted graph'''

    def __init__(self):
        self.num = int(0)
        self.cat = "未知点"
        self.name = ""
        self.X = 0.0
        self.Y = 0.0
        self.Z = 0.0
        self.N = 0.0
        self.E = 0.0
        self.H = 0.0

        #用来标记是第几个未知点，用来列B矩阵
        self.BeiZhu = int(0)

    def init_station(self,line):
        line = line.rstrip('\n')
        line = line.split(',')

        self.num = int(line[0])
        self.cat = line[1]
        self.name = line[2]

        if self.cat == "已知点":
            self.X = float(line[3])
            self.Y = float(line[4])
            self.Z = float(line[5])
            self.N = float(line[6])
            self.E = float(line[7])
            self.H = float(line[8])


#已经算过近似坐标了
    def init_station2(self,line):
        line = line.rstrip('\n')
        line = line.split(',')

        self.num = int(line[0])
        self.cat = line[1]
        self.name = line[2]

        #
        self.X = float(line[3])
        self.Y = float(line[4])
        self.Z = float(line[5])
        self.N = float(line[6])
        self.E = float(line[7])
        self.H = float(line[8])



    def station_inf(self):
        return [ self.num ,self.name, self.X, self.Y , self.Z , self.cat , self.BeiZhu]

    def station_match(self, point_name):
        if self.name == point_name:
            return True
        else :
            return False

    def cat_match(self):
        if self.cat == "已知点":
            return True
        elif self.cat == "未知点":
            return False
        else:
            return None

    def beizhu(self,i):
        self.BeiZhu = i

    def unknow_num(self):
        return self.BeiZhu




class Baseline(object):
    '''undirected unweighted graph'''

    def __init__(self):
        self.num = int(0)
        self.origin = ""
        self.target = ""
        self.name = ""
        self.DX = 0.0
        self.DY = 0.0
        self.DZ = 0.0
        self.sigema_DX = 0.0
        self.sigema_DY = 0.0
        self.sigema_DZ = 0.0

    def init_baseline(self,line):
        self.num = int(line[0])
        self.origin = line[1]
        self.target = line[2]
        self.name = line[3]
        self.DX = float(line[4])
        self.DY = float(line[5])
        self.DZ = float(line[6])
        self.sigema_DX = float(line[7])
        self.sigema_DY = float(line[8])
        self.sigema_DZ = float(line[9])

    def baseline_inf(self):
        return [self.num,self.origin,self.target,self.name ,self.DX,self.DY,self.DZ,self.sigema_DX,self.sigema_DY,self.sigema_DZ,]
        
    def baseline_match(self,baseline):
        if baseline[0] == self.origin and baseline[1] == self.target:
            return True,1,self.num
        elif baseline[0] == self.target and baseline[1] == self.origin:
            return True,-1,self.num
        else:
            return False,0,-1

class GNSSNet(object):
    '''undirected unweighted graph'''

    def __init__(self):

        self.BaselineSet = []

        self.StationSet = []

        self.ControlPts = []
        self.UnknowPts = []

    def insert_Baseline(self,baseline):
        self.BaselineSet.append(baseline)

#初始化站点信息
    def insert_Station(self, station):
        self.StationSet.append(station)

#已知点
    def insert_ControlPoint(self, station):
        self.StationSet.append(station)
#未知点
    def insert_UnknowPoint(self, station):
        self.StationSet.append(station)

    def Net_baseline_match(self,baseline):
        for Baseline in self.BaselineSet:
            a,towards,number = Baseline.baseline_match(baseline)
            if a:
                return towards,self.BaselineSet[number-1].baseline_inf()

    def station_match(self,station_name):
        for station in self.StationSet:
            if station.station_match(station_name):
                return True,station.station_inf()
        return False,None

    def Station_Cat_Match(self,station_name):
        for station in self.StationSet:
            if station.station_match(station_name):
                return True,station.cat_match(),station.station_inf()
        return False,None,None


    def init_P(self,P):
        for Baseline in self.BaselineSet:
            P[Baseline.num*3-3][Baseline.num*3-3] = 1.0 / (Baseline.sigema_DX * Baseline.sigema_DX)
            P[Baseline.num*3-2][Baseline.num*3-2] = 1.0 / (Baseline.sigema_DY * Baseline.sigema_DY)
            P[Baseline.num*3-1][Baseline.num*3-1] = 1.0 / (Baseline.sigema_DZ * Baseline.sigema_DZ)
        return P


    def init_L(self,L):
        for Baseline in self.BaselineSet:
            L[Baseline.num*3-3][0] = Baseline.DX 
            L[Baseline.num*3-2][0] = Baseline.DY
            L[Baseline.num*3-1][0] = Baseline.DZ
        return L

    def init_X0(self,X0):
        for station in self.StationSet:
            if not station.cat_match():
                num = int(station.unknow_num())
                X0[num*3-3][0] = station.X 
                X0[num*3-2][0] = station.Y
                X0[num*3-1][0] = station.Z
                print(station.station_inf()[5] , station.station_inf()[6])
        return X0  

非常直白的平差过程

'''
间接平差
main.py
'''

import pandas as pd
import numpy as np
import sys
np.set_printoptions(threshold=sys.maxsize) #print显示完整array
from gnssnet import *
import math

def Save2Excel(mats,name):
        data = pd.DataFrame(mats)
        writer = pd.ExcelWriter("C:\\Users\\sheld\\Desktop\\result3\\"+ name + ".xlsx")
        data.to_excel(writer, "page_1", float_format = '%.6f')#浮点数，精确到6位小数
        writer.save()
        writer.close()

def init_GNSSNet(G):

#读入所有的基线坐标信息
    with open("data.csv",'r') as f:
        for line in f.readlines():
            b= Baseline()
            b.init_baseline(line.split(','))
            G.insert_Baseline(b)

        f.close()

#已经手算过近似坐标了
    with open("zhandian.csv",'r',encoding='UTF-8-sig') as f:
        i = 1
        for line in f.readlines():
                
                p = station()
                p.init_station2(line)

                if not p.cat_match():
                        p.beizhu(i)
                        i = i + 1
                G.insert_Station(p)
        f.close()

n = 63
t = 24
r = n - t

global G
G = GNSSNet()
init_GNSSNet(G)

P = np.zeros([n,n], dtype = float)
L = np.zeros([n,1], dtype = float)
l = np.zeros([n,1], dtype = float)

'''
L^ = BX^ + d
l = L - (BX0 + d)
V = Bx^ - l
'''
B =  np.zeros([n,t], dtype = float)
d = np.zeros([n,1], dtype = float)

X0 = np.zeros([t,1], dtype = float)
x = np.zeros([t,1], dtype = float)

#写入B,d矩阵
for baseline in G.BaselineSet:
        bl_inf = baseline.baseline_inf()
        num = bl_inf[0]
        origin = bl_inf[1]
        target = bl_inf[2]

        #起点
        found,ifknow,stat_inf  = G.Station_Cat_Match(origin)
        if ifknow:
                d[3 * num - 3][0] = d[3 * num - 3][0] - stat_inf[2]
                d[3 * num - 2][0] = d[3 * num - 2][0] - stat_inf[3]
                d[3 * num - 1][0] = d[3 * num - 1][0] - stat_inf[4]
##                print(num,bl_inf,"起点是已知点",stat_inf)
        elif not ifknow:
                B[3 * num - 3][3 * stat_inf[6] - 3] = -1
                B[3 * num - 2][3 * stat_inf[6] - 2] = -1
                B[3 * num - 1][3 * stat_inf[6] - 1] = -1
##                print(num,bl_inf,"起点是未知点",stat_inf)
                
        else:
                print("有错误的站点名")
        
        #末尾
        found,ifknow,stat_inf  = G.Station_Cat_Match(target)
        if ifknow:
                d[3 * num - 3][0] = d[3 * num - 3][0] + stat_inf[2]
                d[3 * num - 2][0] = d[3 * num - 2][0] + stat_inf[3]
                d[3 * num - 1][0] = d[3 * num - 1][0] + stat_inf[4]
##                print(num,bl_inf,"终点是已知点",stat_inf)
                
        elif not ifknow:
                B[3 * num - 3][3 * stat_inf[6] - 3] = 1
                B[3 * num - 2][3 * stat_inf[6] - 2] = 1
                B[3 * num - 1][3 * stat_inf[6] - 1] = 1
##                print(num,bl_inf,"终点是未知点",stat_inf)
        else:
                print("有错误的站点名")

##        print(B[3 * num - 3])
##        print(B[3 * num - 2])
##        print(B[3 * num - 1])
##        print(d[3 * num - 3])
##        print(d[3 * num - 2])
##        print(d[3 * num - 1])
            
#定权
P = G.init_P(P)
L = G.init_L(L)
x0 = G.init_X0(X0)

P = np.matrix(P)
L = np.matrix(L)
X0 = np.matrix(X0)
d = np.matrix(d)
B = np.matrix(B)

Q = P.I

print("X0:\n",X0)
print("L:\n",L)

print("d:\n",d)
print("B:\n",B)

l = L - np.dot(B,X0) - d
print("l:\n",l)

#化成毫米单位
l = np.matrix(np.dot(l, 1000))

#Nbbx^ - W = 0
Nbb = np.dot(np.dot(B.T, P), B)
W = np.dot(np.dot(B.T, P), l)

#矩阵的秩
print(np.linalg.matrix_rank(B, tol=None, hermitian=False))
print(np.linalg.matrix_rank(l, tol=None, hermitian=False))

#x^
x =  np.dot(Nbb.I,W)
V = np.dot(B,x) - l

print('x',x,'\n')
print('V',V,'\n')

time =  1
V_total = V
x_total = x
##
##while  abs(V.max()) > 0.00001 and time < 10000:
##        L = L + V/1000
##        X0 = X0 + x/1000
##
##        l = L - np.dot(B,X0) - d
##
##        l = np.matrix(np.dot(l, 1000))
##        x =  np.dot(Nbb.I,W)
##        V = np.dot(B,x) - l
##        
##        V_total = V_total + V
##        x_total = x_total + x
##        time = time + 1
##
##

L = L + V/1000
X0 = X0 + x/1000

print('L^',L)
print("time:",time,'\n')
print('V_total',V_total,'\n')
print('x_total',x_total,'\n')

sigema02 = np.dot(np.dot(V_total.T,P),V_total)/r
sigema0 = math.sqrt(sigema02)
print(sigema0)

Save2Excel(B,"B")
Save2Excel(d,"d")
Save2Excel(L,"L^")
Save2Excel(l,"l")
Save2Excel(X0,"X^")
Save2Excel(V_total,"V_total")
Save2Excel(x_total,"x_total")
Save2Excel(P,"P")