四组数据算七参数

1.布尔莎模型转换原理

（1）模型公式：

 

（2）七参数求解过程：

首先，展开模型，得到：

 

其次，构建误差方程，选用间接平差方式，其形式为：

 

写成矩阵形式为：

 

对其进行求解并默认为独立同精度观测,有：

 

对其进行精度评定，有

 

求其方差，有

 

对其开算术平方根，即可得到绝对精度中误差。

 

2.数据源及数据信息

表1 四个公共点两套坐标值

点号/坐标	1	2	3	4
/m	-2722720.343	-2723130.106	-2726120.504	-2725958.007
/m	4429306.301	4431260.780	4427716.358	4430477.961
/m	3682136.625	3679550.261	3681514.511	3678382.473
/m	-2722708.864	-2723118.632	-2726109.020	-2725946.525
/m	4429315.210	4431269.685	4427725.263	4430486.870
/m	3682132.936	3679546.589	3681510.823	3678378.793

来源于《应用大地测量学》.张海华著

3.程序设计

# -*- coding: utf-8 -*-
import numpy as np

from numpy import *

list1 = [-2722720.343, 4429306.301, 3682136.625]
list2 = [-2722708.864, 4429315.210, 3682132.936]
list3 = [-2723130.106, 4431260.780, 3679550.261]
list4 = [-2723118.632, 4431269.685, 3679546.589]
list5 = [-2726120.504, 4427716.358, 3681514.511]
list6 = [-2726109.020, 4427725.263, 3681510.823]
list7 = [-2725958.007, 4430477.961, 3678382.473]
list8 = [-2725946.528, 4430486.870, 3678378.793]
# 七参数转换法
B = np.mat([[list1[0], 1, 0, 0, 0, -list1[2], list1[1]],
            [list1[1], 0, 1, 0, list1[2], 0, -list1[0]],
            [list1[2], 0, 0, 1, -list1[1], list1[0], 0],
            [list3[0], 1, 0, 0, 0, -list3[2], list3[1]],
            [list3[1], 0, 1, 0, list3[2], 0, -list3[0]],
            [list3[2], 0, 0, 1, -list3[1], list3[0], 0],
            [list5[0], 1, 0, 0, 0, -list5[2], list5[1]],
            [list5[1], 0, 1, 0, list5[2], 0, -list5[0]],
            [list5[2], 0, 0, 1, -list5[1], list5[0], 0],
            [list7[0], 1, 0, 0, 0, -list7[2], list7[1]],
            [list7[1], 0, 1, 0, list7[2], 0, -list7[0]],
            [list7[2], 0, 0, 1, -list7[1], list7[0], 0]])  # B为系数矩阵

L = np.mat([[list2[0]],
            [list2[1]],
            [list2[2]],
            [list4[0]],
            [list4[1]],
            [list4[2]],
            [list6[0]],
            [list6[1]],
            [list6[2]],
            [list8[0]],
            [list8[1]],
            [list8[2]]])  # L为常数列

d = np.mat([[list1[0]],
            [list1[1]],
            [list1[2]],
            [list3[0]],
            [list3[1]],
            [list3[2]],
            [list5[0]],
            [list5[1]],
            [list5[2]],
            [list7[0]],
            [list7[1]],
            [list7[2]]])  # d为常数列

X = (B.T * B).I * B.T * (d - L)

n = 12
t = 7
r = n - t
v = B * X - (d - L)

mid0 = (v.T * v) / r

Qxx = (B.T * B).I

Dm = sqrt(Qxx[0, 0] * mid0)
DX0 = sqrt(Qxx[1, 1] * mid0)
DY0 = sqrt(Qxx[2, 2] * mid0)
DZ0 = sqrt(Qxx[3, 3] * mid0)
Dx = sqrt(Qxx[4, 4] * mid0)
Dy = sqrt(Qxx[5, 5] * mid0)
Dz = sqrt(Qxx[6, 6] * mid0)

print("""七参数转化法求得转换参数为
m:{}
平移量X0:{}m
平移量Y0:{}m
平移量Z0:{}m
旋转量x：{}"
旋转量y：{}"
旋转量z：{}"
---------------------------------------------
精度评定中误差为：
单位权中误差为{}
缩放值m中误差为{}
平移量X0的中误差为{}m
平移量Y0的中误差为{}m
平移量Z0的中误差为{}m
旋转量x的中误差为{}
旋转量y的中误差为{}
旋转量z的中误差为{}
---------------------------------------------""".format(round(X[0, 0], 11), round(X[1, 0], 4), round(X[2, 0], 4), round(X[3, 0], 4),
                      round(np.rad2deg(X[4, 0]) * 3600, 4), round(np.rad2deg(X[5, 0]) * 3600, 4),
                      round(np.rad2deg(X[6, 0]) * 3600, 4), mid0, Dm[0, 0], DX0[0, 0], DY0[0, 0], DZ0[0, 0], Dx[0, 0],
                      Dy[0, 0], Dz[0, 0]))

C = np.mat([[list8[0], 1, 0, 0, 0, -list8[2], list8[1]],
            [list8[1], 0, 1, 0, list8[2], 0, -list8[0]],
            [list8[2], 0, 0, 1, -list8[1], list8[0], 0]])

D = np.mat([[list8[0]],
            [list8[1]],
            [list8[2]]])

print(round(X[4, 0], 30))

X1 = np.mat([
    [round(X[0, 0], 11)],
    [round(X[1, 0], 4)],
    [round(X[2, 0], 4)],
    [round(X[3, 0], 4)],
    [round(X[4, 0],9)],
    [round(X[5, 0], 9)],
    [round(X[6, 0], 9)]
])

# [round(np.rad2deg(X[4, 0]) * 3600, 4)],
# [round(np.rad2deg(X[5, 0]) * 3600, 4)],
# [round(np.rad2deg(X[6, 0]) * 3600, 4)]

E = C * X1 + D

print("""将点号4的X2,Y2,Z2带回重新进行验证，所得结果如下：
点号4第一套坐标({}，{}，{})
""".format(round(E[0, 0], 3), round(E[1, 0], 3), round(E[2, 0], 3)))

4.输出结果

1. 七参数转化法求得转换参数为
2. m:1.72153e-06
3. 平移量X0:-13.2685m
4. 平移量Y0:-29.2514m
5. 平移量Z0:7.8575m
6. 旋转量x：0.4814"
7. 旋转量y：0.013"
8. 旋转量z：0.3126"
9. ---------------------------------------------
10. 精度评定中误差为：
11. 单位权中误差为[[1.4454205e-05]]
12. 缩放值m中误差为7.459597287807743e-07
13. 平移量X0的中误差为6.6354744837286885m
14. 平移量Y0的中误差为5.78210708382355m
15. 平移量Z0的中误差为5.998230759476042m
16. 旋转量x的中误差为9.730780488227572e-07
17. 旋转量y的中误差为9.367198813454598e-07
18. 旋转量z的中误差为9.85159398180071e-07
19. ---------------------------------------------
20. 将点号4的X2,Y2,Z2带回重新进行验证，所得结果如下：
21. 点号4第一套坐标(-2725958.008，4430477.962，3678382.47)

5.结果与分析

数据采用课本中的案例，使用了四个联测点8套坐标全部参与7参数求解，所得结果与课本结果一致，即此程序对误差方程中X的平差值求解无误，对七参数求解无误。当把点号4中的第二套坐标带回到此七参数下进行反算时，出现了x坐标差1mm,y坐标差1mm，z坐标差3mm的情况。我认为是在对欧勒角的小数点后所取位数上和课本存在不一致的现象，我们无法得知课本中对弧度下的欧勒角所取位数，只知道在角度下其取到了小数点后四位，我尝试按照弧度依然取四位有效数字进行演算，经验证，猜想基本正确，测试结果如下：

旋转角x、y、z取到有效数字第四位，则输出结果为：(-2725958.007，4430477.961，3678382.47)，即此时只存在z坐标差3mm的情况，x、y坐标不存在差异。

即我们如果知道了课本中对弧度制下欧勒角的具体取舍方式，便可得到与课本完全一致的结果。