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顆衛星的座標，運算結果以下: