GPS原理与应用实验
题 目： GPS单点定位
专 业： 测绘工程
班 级： 12-01
学 号：
姓 名：
指导教师：
时间：20##.11
目录
一、实验目的......................................................3
二、实验原理………..........................................3
三、实验内容……..............................................3
四、实验效果图..................................................9
五、实验总结......................................................9
一． 实验目的
1.深入了解单点定位的计算过程；
2.加强单点定位基本公式和误差方程式，法线方程式的记忆；
3.通过上机调试程序加强动手能力的培养。
二．实验原理
一个接收机接受三个火三个以上卫星信号，得出卫星坐标和伪距，利用间接平差计算接收机的坐标。
三．实验内容
1.程序流程图
2、实验数据
3、实验程序代码
Private Sub Command1_Click()
CommonDialog1.Filter = "TXT files|*.txt|"
CommonDialog1.FilterIndex = 1
CommonDialog1.ShowOpen
Open Me.CommonDialog1.FileName For Input As #1
Do While Not EOF(1)
Line Input #1, Text
textbuff = textbuff + Text + vbCrLf
Loop
Close #1
kk = MSFlexGrid1.Rows - 1
Dim a
ReDim a(kk - 1)
a = Split(textbuff, vbCrLf)
For j = 1 To kk
For i = 1 To 5
MSFlexGrid1.TextMatrix(j, i) = a(j - 1 + 5 * (i - 1))
Next i
Next j
For k = 1 To kk
MSFlexGrid1.TextMatrix(k, 0) = "第" & k & "个点"
Next k
MSFlexGrid1.TextMatrix(0, 1) = "X"
MSFlexGrid1.TextMatrix(0, 2) = "Y"
MSFlexGrid1.TextMatrix(0, 3) = "Z"
MSFlexGrid1.TextMatrix(0, 4) = "伪距"
MSFlexGrid1.TextMatrix(0, 5) = "钟差"
End Sub
Private Sub Command2_Click()
kk = MSFlexGrid1.Rows - 1
X0 = 0: Y0 = 0: Z0 = 0
c = 299792458
Dim a()
ReDim a(kk - 1, 3)
Dim ll()
ReDim ll(kk - 1, 0)
For ii = 1 To 100
For i = 1 To kk
l = (MSFlexGrid1.TextMatrix(i, 1) - X0) / Sqr((MSFlexGrid1.TextMatrix(i, 1) - X0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 2) - Y0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 3) - Z0) ^ 2)
m = (MSFlexGrid1.TextMatrix(i, 2) - Y0) / Sqr((MSFlexGrid1.TextMatrix(i, 1) - X0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 2) - Y0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 3) - Z0) ^ 2)
n = (MSFlexGrid1.TextMatrix(i, 3) - Z0) / Sqr((MSFlexGrid1.TextMatrix(i, 1) - X0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 2) - Y0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 3) - Z0) ^ 2)
a(i - 1, 0) = l
a(i - 1, 1) = m
a(i - 1, 2) = n
a(i - 1, 3) = -1
lk = MSFlexGrid1.TextMatrix(i, 4) - Sqr((MSFlexGrid1.TextMatrix(i, 1) - X0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 2) - Y0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 3) - Z0) ^ 2) + c * MSFlexGrid1.TextMatrix(i, 5)
ll(i - 1, 0) = lk
Next i
gzs = xc(qiuni(xc(zz(a), a)), xc(zz(a), ll))
X0 = X0 - gzs(0, 0)
Y0 = Y0 - gzs(1, 0)
Z0 = Z0 - gzs(2, 0)
j = j + 1
Next ii
Text2.Text = "X=" & X0 & vbCrLf & vbCrLf & "Y=" & Y0 & vbCrLf & vbCrLf & "Z=" & Z0
V = jian(ll, xc(a, gzs))
zjl = xc(zz(V), V)
σ0 = Sqr(zjl(0, 0)) / (kk - 3)
Qx = qiuni(xc(zz(a), a))
Text3.Text = "σX=" & σ0 * Sqr(Qx(0, 0)) & vbCrLf & vbCrLf & "σY=" & σ0 * Sqr(Qx(1, 1)) & vbCrLf & vbCrLf & "σZ=" & σ0 * Sqr(Qx(2, 2))
End Sub
Private Sub Form_Load()
MSFlexGrid1.ColWidth(1) = 1300
MSFlexGrid1.ColWidth(2) = 1300
MSFlexGrid1.ColWidth(3) = 1300
MSFlexGrid1.ColWidth(4) = 1300
Text2.Text = ""
Text3.Text = ""
End Sub
'矩阵相减
Public Function jian(m, n)
Dim i, j As Integer
If UBound(m, 1) ＜＞ UBound(n, 1) Or UBound(m, 2) ＜＞ UBound(n, 2) Then
MsgBox ("请确认输入数组是否可以相减！")
Else
Dim c()
ReDim c(UBound(m, 1), UBound(n, 2))
For i = 0 To UBound(c, 1)
For j = 0 To UBound(c, 2)
c(i, j) = m(i, j) - n(i, j)
Next j
Next i
jian = c
End If
End Function
'矩阵的转置
Public Function zz(a)
Dim i As Integer, j As Integer, t As Integer, b()
If UBound(a, 1) = UBound(a, 2) Then
For i = 0 To UBound(a, 1)
For j = 0 To UBound(a, 2)
If i ＜ j Then
t = a(i, j)
a(i, j) = a(j, i)
a(j, i) = t
End If
Next j
Next i
zz = a
Else
ReDim b(UBound(a, 2), UBound(a, 1))
For i = 0 To UBound(a, 2)
For j = 0 To UBound(a, 1)
b(i, j) = a(j, i)
Next j
Next i
zz = b
End If
End Function
'两矩阵相乘
Public Function xc(a, b)
Dim i As Integer, j As Integer, k As Integer
If UBound(a, 2) ＜＞ UBound(b, 1) Then
MsgBox ("这两个矩阵不能够相乘")
Exit Function
End If
ReDim sd(UBound(a, 1), UBound(b, 2))
For i = 0 To UBound(a, 1)
For j = 0 To UBound(b, 2)
For k = 0 To UBound(b, 1)
sd(i, j) = sd(i, j) + a(i, k) * b(k, j)
Next k
Next j
Next i
xc = sd
End Function
Public Function qiuni(a)
Dim c, m%, n%, p#, l%, i%, j%, ab#
m = UBound(a, 1)
n = UBound(a, 2)
If m ＜＞ n Then
MsgBox ("该矩阵不可逆！！！")
Exit Function
End If
ReDim c(m, 2 * n + 1)
For i = 0 To m
For j = 0 To n
c(i, j) = a(i, j)
Next j
Next i
For i = 0 To m
For j = m + 1 To 2 * m + 1
c(i, j) = 0
Next j
Next i
i = 0
For j = m + 1 To 2 * m + 1
c(i, j) = 1
i = i + 1
Next j
For k = 0 To n
If c(k, k) = 0 Then
For i = k + 1 To n
If c(i, k) ＜＞ 0 Then
GoTo this
End If
Next i
If i = n + 1 Then
MsgBox ("该矩阵不可逆!!!")
Exit Function
End If
this:
For j = 0 To 2 * m + 1
p = c(k, j)
c(k, j) = c(i, j)
c(i, j) = p
Next j
End If
ab = 1# / c(k, k)
For j = 0 To 2 * m + 1
c(k, j) = c(k, j) * ab
Next j
For i = 0 To n
If i ＜＞ k Then
For j = 0 To 2 * m + 1
If j ＜＞ k Then
c(i, j) = c(i, j) - c(i, k) * c(k, j)
End If
Next j
c(i, k) = 0
End If
Next i
Next k
For i = 0 To m
For j = 0 To m
a(i, j) = c(i, j + n + 1)
a(i, j) = Round(a(i, j), 4)
Next j
Next i
qiuni = a
End Function
四．实验结果图
五．实验总结
此次实验让我深入了解单点定位的计算过程，加强了对单点定位基本公式和误差方程式，法线方程式的记忆。并通过上机调试程序加强了自己动手能力的培养。通过此次实验为以后的GPS的学习打下了基础。
第二篇：声源定位和GPS模拟实验报告
声源定位和GPS模拟实验报告
一．实验原理
1. 声源定位(二维)
如图所示,3个接收传感器S 0、S 1、S 2的坐标分别是（X 0，Y 0）、（X 1,Y 1）、（X 2,Y 2），当平面上某处（X,Y）发出超声波时，该信号将先后被3个传感器所接收，设时间分别是t 0,t 1,t 2,实验只能测出它们到达各个传感器的时间差△t 1=t 1-t 0, △t 2=t 2-t 0。设声波沿媒质表面的传播度为c，对换能器S 0、S 1而言，声源发生的位置应该在到该两点的距离差为c△t 1的曲线上，这是一条双曲线。显然，利用△t 1、△t 2可以得到两条双曲线，它们的交点即是声源
将S0设为原点，声源位置为（X,Y），用极坐标表示，则满足：
化简整理的
令
引入
,
。
则可得
至此，声源位置已通过极坐标给出。
2. GPS模拟
本实验对GPS过程的声学模拟是在一个二维的平面上进行的，位置（Xi，Yi）(i=1,2,……)已知的发送换能器（传感器做发送用，模拟“导航卫星“）发出声波（模拟卫星发出的电磁波），被位置（X,Y）的待求接收传感器（模拟用户）接收，它们之间有关系：
（i=1,2，……）
式中c为波的传播速度。显然，对一个二维的定位问题（确定X和Y），如果传播速度已知，要算出X和Y可以归结为一个求解两个变量的代数方程问题，也就是说原则上只要有两颗不同位置的模拟卫星就行了。实际的GPS定位，则至少要对四颗卫星同时进行测量，才能确定地球坐标系中的三维坐标和因接收机时钟不同步所造成的种差修正。在本实验中，为了减小时差不准对定位精度的影响，应当获取来自多个“卫星“的位置和时差信息，并通过最小二乘法求得”用户“的位置和声波的传播速度。即当n大于3时，
F(X,Y,c)=
=min
由此可获得X,Y,c应满足的一组代数方程：
这是一个三元的非线性代数方程组，可通过计算方法求得数值解。若c未知，而且需要考虑钟差修正，则获得一组四元代数方程。可由牛顿迭代法求数值解。
二．实验仪器
传播媒质、模拟源（铅笔芯折断）、压力传感器、接收放大器、时差测定装置、计算机、隔离放大器、单脉冲发生器。
实验电路：
三．数据处理
1. 声源定位
选择坐标（50,200）单位：时间us 位置mm
实验测得的结果（49.8,201.9）跟实际值相比相差很少，笔芯折断的位置和网格的误差都有可能引起实验结果与实际值的差异。
2. GPS模拟
计算的P3的坐标为X=-0.00170731 mm Y=-0.00144067mm
实验误差较小
主要误差来源：发生器的位置读数有误差，传感器与媒质界面的接触不均匀也可能带来误差。
思考题
1.1从测量原理、仪器组成和信息特点来讨论声源定位和GPS仿真实验的相同点和不同点。
答：测量原理：前者是已知接收器位置求声源位置；后者是已知发生器位置求接收器位置。
仪器组成：前者：传播媒质、电压传感器、接收放大器、时差测定仪、模拟源。
后者：传播媒质、电压传感器，接收放大器，时差测定仪、隔离放大器、单脉冲发生器。
信息特点：都是位置坐标和时间差。
2考虑到钢媒质的表面波传播速度约为3000 m/s,传播距离最长约为54cm,那么时差测定仪最大的读数应该不会大于多少？
答：t=d/c=180us
2.1声源定位的发射源位置如何确定？是否用两个时差就可以了？
答：用时间差所确定的双曲线的交点坐标即发射源坐标。两个时间差就可以了。
2GPS仿真实验用最小二乘法求数值解时，一般会涉及4个待测参数，除用户接收机的X,Y坐标和传播速度C以外，另一个参数是什么？有无具体的物理意义？
答：另一个参数为钟差，它是卫星与接收器时间同步性的一个衡量指标。
1、 声源定位：
#include
#include
using namespace std;
#define c 2982
#define pai 3.1415926535
int main()
{
double x0,x1,x2,y0,y1,y2,t0,t1,t2,t3,A,B,D,angle1,angle2,dt1,dt2,r;
x0=0,x1=0,x2=0.3;
y0=0,y1=0.45,y2=0.45;
int i = 8;
while(i--)
{
cout ＜＜"请输入第"＜＜8-i＜＜"组时间数据"＜＜ endl;
cin ＞＞ t1 ＞＞ t2＞＞ t3;
t1 = t1 / 1000000;
t2 = 0;
t3 = t3 / 1000000;
dt1 = t1 - t2;
dt2 = t3 - t2;
A = x2*(pow(x1,2)+pow(y1,2)-pow(c*dt1,2)) - x1*(pow(x2,2)+pow(y2,2)-pow(c*dt2,2));
B = y2*(pow(x1,2)+pow(y1,2)-pow(c*dt1,2)) - y1*(pow(x2,2)+pow(y2,2)-pow(c*dt2,2));
D = c*dt1*(pow(x2,2)+pow(y2,2)-pow(c*dt2,2)) - c*dt2*(pow(x1,2)+pow(y1,2)-pow(c*dt1,2));
angle2 = atan(B/A);
if(acos(D/sqrt(pow(A,2)+pow(B,2)))+angle2 ＜ pai/2)
angle1 = acos(D/sqrt(pow(A,2)+pow(B,2)))+angle2;
else
angle1 = angle2 - acos(D/sqrt(pow(A,2)+pow(B,2)));
r = (pow(x1,2)+pow(y1,2)-pow(c*dt1,2))/(2*(x1*cos(angle1)+y1*sin(angle1)+c*dt1));
cout ＜＜ "声源坐标为：（"＜
} // end while;
return 0;
} // end main;
2、 GPS模拟
#include
#include
#include
using namespace std;
#define c 2982
#define min 0.0000000001
void GPS(double x[10] , double y[10] , double t[10] , double x1 , double y1 , int num)
{
double a1,a2,a3,a4,b1,b2,mult,D=0,E=0;
int i;
while(1)
{
a1 = a2 = a3 = a4 = b1 = b2 = 0;
for(i = 0; i ＜ 10; i ++)
{
a1 = a1-(3*pow((x[i]-x1),2)+pow((y[i]/1000-y1),2)-pow(c*t[i]/1000000,2));
a2 = a2+2*(x[i]/1000-x1)*(y1-y[i]/1000);
a3 = a2;
a4 = a4-(3*pow((y[i]-y1),2)+pow((x[i]/1000-x1),2)-pow(c*t[i]/1000000,2));
b1 = b1+(pow((x[i]/1000-x1),2)+pow((y[i]/1000-y1),2)-pow(c*t[i]/1000000,2))*(x1-x[i]/1000);
b2 = b2+(pow((x[i]/1000-x1),2)+pow((y[i]/1000-y1),2)-pow(c*t[i]/1000000,2))*(y1-y[i]/1000);
} //end for;
mult = a3/a1;
a4 -= a2*mult;
b2 -= b1*mult;
E = b2/a4;
D = (b1-a2*E)/a1;
x1 += D;
y1 += E;
if((pow(D,2)+pow(E,2)) ＜ min)
break;
} // end while;
x1 = x1*1000;
y1 = y1*1000;
cout ＜＜"点P"＜＜ num ＜＜"定位于: X = "＜＜ fixed ＜＜ setprecision(8)＜＜ x1;
cout ＜＜"(mm) , Y = "＜＜ fixed ＜＜ setprecision(8)＜＜ y1 ＜＜"(mm)"＜
} // end GPS;
int main()
{
int yy;
double t1[10] , t2[10] , t3[10] , x[10] , y[10];
cout ＜＜ "请输入时间数据(us)与坐标数据(mm):" ＜＜ endl;
int i;
for(i = 0; i ＜= 9; i ++)
{
cout ＜＜ "第" ＜＜ i+1 ＜＜ "组:" ＜＜ endl;
cin ＞＞ t1[i] ＞＞ t2[i] ＞＞ t3[i] ＞＞ x[i] ＞＞ y[i];
} // end for;
GPS(x,y,t1,0,0.45,1);
GPS(x,y,t2,0.3,0.45,2);
GPS(x,y,t3,0,0,3);
cin＞＞yy;
return 0;
} // end main;