Numpy库应用实例——GPS定位

背景介绍
 定位系统
 GPS全球定位系统(Global Positioning System)
以GPS系统为例介绍卫星定位的计算方法

GPS定位的基本原理
 GPS定位的基本原理是根据高速运动卫星的瞬间位置作为已知的起算数据，采用空间距离-后方交会的方法，确定待测点的位置。
 假设t时刻在地面待测点上安置GPS接收机，可以测定GPS信号到达接收机的时间△t，再加上接收机所接收到的卫星星历等其它数据，就可以确定一个方程组来对位置信息进行求解。
 假设地球上一个点R，同时收到6颗卫星（S1,S2,…,S6）发射的信号，假设接受信息如下表所示。其中x，y表示卫星的经纬度，z表示卫星的高度。

  由于上述6个卫星和地球在高速运动，从卫星发出的位置信息以光速传输到GPS接收端需要一定的时间。
  假设(x,y,z,t)表示R当前的位置, t是R的相对时间，卫星S1（发出信号时刻）到（当前接收时刻)满足以下关系(其中c是光速)。
 (x-3)^2 + (y-2)^2 + (z-3)^2 = [(10010.00692286 – t)*c]^2，
 该公式表示以(x, y, z,t)为参数的(欧式空间距离)与信号传输距离相等。

对于卫星S1,S2,…,S6，满足方程组:
 ...(1)

其中，光速为常数c=0.299792458km/us，上述方程组是非线性的，但很容易将所有二次项都消去（每个公式减去第一个公式），从而得到：

此时，上述等式变成了A*X=B形式，根据线性代数方法，X=A-1*B，即只需对系数矩阵求逆，再乘以常数矩阵便可以得到方程组的解。

GPS定位的问题建模
 上面给出了GPS的定位原理，如何利用计算机辅助GPS的定位计算呢？
 以6颗卫星为例，GPS定位计算问题的IPO模式----描述如下：
 输入：6颗卫星的欧式坐标和信号时间戳
 处理：GPS定位算法
 输出：GPS接收设备的地理坐标和当前时间
假设第i颗卫星的坐标和时间戳表示为(x_i, y_i ,z_i ,t_i )，结合上述例子，GPS定位算法可以描述为如下公式：

我们下面将使用Numpy函数库实现上述矩阵操作。首预习一下程序中用到的函数：

numpy.dot(a,b)：计算矩阵a与矩阵b的点积
numpy.linalg.inv(a)：求矩阵a 的逆矩阵

GPS定位的程序实现
 Python代码如下：
 其中zeros是NumPy提供的函数，用来建立指定维度的数组,
 zeros用来生成数组x用来存储接受来自外部输入的六颗卫星坐标，
 数组a，b用来存放前面算法中的系数矩阵，
 例程中我们还展示了两种数组的索引方法，最后调用求矩阵逆的函数及点乘操作完成坐标计算。

from numpy import *
def main_GPSLocation():
    i = 1
    c = 0.299792458  # 光速 0.299792458km/us
    x = zeros((6, 4)) #存储6个卫星的（x,y,z,t）参数
    while i<=6:
        print(" %s %d" % ("please input (x,y,z,t) of group",i) )
        temp=input()
        x[i-1]=temp.split()
        j=0
        while j<4:
            x[i-1][j]=float(x[i-1][j])
            j=j+1
        i=i+1
    a=zeros((4,4)) #系数矩阵
    b=zeros((4,1)) #常数项
    j=0
    while j<4:
        a[j][0]=2*(x[5][0]-x[j][0])
        a[j][1]=2*(x[5][1]-x[j][1])
        a[j][2]=2*(x[5][2]-x[j][2])
        a[j][3]=2*c*c*(x[j][3]-x[5][3])
        b[j][0]=x[5][0] * x[5][0] - x[j][0] * x[j][0] + 
                x[5][1] * x[5][1] - x[j][1] * x[j][1] + 
                x[5][2] * x[5][2] - x[j][2] * x[j][2] + 
            c*c*(x[j][3] * x[j][3] - x[5][3] * x[5][3])
        j=j+1
    a_ni=linalg.inv(a) #系数矩阵求逆
    print(dot(a_ni,b))

main_GPSLocation()

运行程序后，依次输入6颗卫星的坐标，运算结果如下: