原理：空間後方交會是以單幅影像為基礎，從該影像所覆蓋地面範圍內若干控制點的已知地面座標和相應點的像座標量測值出發，根據共線條件方程，解求該影像在航空攝影時刻的外方位元素Xs，Ys，Zs，φ，ω，κ。

演算法：由於每一對像方和物方對應點可列出2個方程，因此若有3個已知地面座標的控制點，則可列出6個方程，解求6個外方位元素的改正數△Xs，△Ys，△Zs，△φ，△ω，△κ。實際應用中為了提高解算精度，常有多餘觀測方程，通常是在影像的四個角上選取4個或均勻地選擇更多的地面控制點，因而要用最小二乘平差方法進行計算。

在不考慮控制點誤差和內方位元素誤差、像點觀測值視為等權的情況下，視、、、φ，ω，κ為待定引數（未知數），列出線性化誤差方程式如下：

 步驟：

1）獲取資料。

2）確定Xs、Ys、Zs、φ，ω，κ的初值。

3）計算旋轉矩陣R。

4）逐點計算像點座標的近似值、。

5）逐點計算誤差方程式的係數和常數項。

6）計演算法方程的係數和常數項。

7）解求未知數Xs、Ys、Zs、φ，ω，κ的改正數。

8）檢查計算是否收斂。

將所求得的外方位元素的改正數與規定的限差比較，通常對φ，ω，κ的改正數△φ，△ω，△κ給予限差，通常為0.000001弧度，當3個改正數均小於0.000001弧度時，迭代結束。否則用近似值加上改正數作為新的近似值，重複3）至7）步的計算。如果迭代次數超過某一規定的次數，則停止迭代，認為不收斂，輸出中間結果和出錯資訊。

9）評定精度，計算單位權中誤差，其中v為觀測值的殘差。

  示例程式碼：

# -*- coding: utf-8 -*-

import numpy as np
import math as m

##示例資料
X=[36589.41,37631.08,39100.97,40426.54]
Y=[25273.32,31324.51,24934.98,30319.81]
Z=[2195.17, 728.69,  2386.50, 757.31]

x=[-86.15,-53.40,-14.78,10.46]
y=[-68.99,82.21,-76.63,64.43]

fo = wo = ko = 0
xo = yo = 0
f = 153.24
Zsum = 0
Xsum=0
Ysum=0
for i in [0,1,2,3]:
    ###Zsum =Zsum + Z[i]
    Xsum =Xsum + X[i]
    Ysum =Ysum + Y[i]

Zos = 7200
Xos = Xsum/4.
Yos = Ysum/4.




def r_mat(f,w,k):
   Rf = np.mat([[m.cos(f), 0, -m.sin(f)],
              [0,         1,          0],
              [m.sin(f), 0,  m.cos(f)]])

   Rw = np.mat([[1,         0,          0],
              [0, m.cos(w), -m.sin(w)],
              [0, m.sin(w),  m.cos(w)]])
    
   Rk = np.mat([[m.cos(k), -m.sin(k), 0],
              [m.sin(k),  m.cos(k), 0],
              [0,         0,          1]])
    
   R = Rf*Rw*Rk 
    
   return R

def xy_approximate(X,Y,Z,x,y,Xs,Ys,Zs,R):
   x_apxm = [0,0,0,0]
   y_apxm = [0,0,0,0]
    
   for i in [0,1,2,3]:
        x_apxm[i] =x[i] - (xo - f*((R[0,0]*(X[i]-Xs)+R[1,0]*(Y[i]-Ys)+R[2,0]*(Z[i]-Zs))
                            /(R[0,2]*(X[i]-Xs)+R[1,2]*(Y[i]-Ys)+R[2,2]*(Z[i]-Zs))))
        y_apxm[i] =y[i] - (yo - f*((R[0,1]*(X[i]-Xs)+R[1,1]*(Y[i]-Ys)+R[2,1]*(Z[i]-Zs))
                            /(R[0,2]*(X[i]-Xs)+R[1,2]*(Y[i]-Ys)+R[2,2]*(Z[i]-Zs))))

   return x_apxm,y_apxm
        
def a_parameter(X,Y,Z,Xs,Ys,Zs,x,y,w,k,R):
    parameter = np.zeros((2,6))
    mean = np.zeros((3,1))
    minus = np.zeros((3,1))
    
    minus = np.array([[X-Xs],
                     [Y-Ys],
                     [Z-Zs]])
    mean = R.T * np.mat(minus)
    
    parameter[0][0] = (R[0,0]*f+R[0,2]*(x-xo))/mean[2]
    parameter[0][1] = (R[1,0]*f+R[1,2]*(x-xo))/mean[2]
    parameter[0][2] = (R[2,0]*f+R[2,2]*(x-xo))/mean[2]
    parameter[1][0] = (R[0,1]*f+R[0,2]*(y-yo))/mean[2]
    parameter[1][1] = (R[1,1]*f+R[1,2]*(y-yo))/mean[2]
    parameter[1][2] = (R[2,1]*f+R[2,2]*(y-yo))/mean[2]
    
    parameter[0][3] = (y-yo)*m.sin(w)-(((x-xo)/f)*((x-xo)*m.cos(k)-(y-yo)*m.sin(k))+f*m.cos(k))*m.cos(w)
    parameter[0][4] = -f*m.sin(k)-((x-xo)/f)*((x-xo)*m.sin(k)+(y-yo)*m.cos(k))
    parameter[0][5] = y-yo
    parameter[1][3] = -(x-xo)*m.sin(w)-(((y-yo)/f)*((x-xo)*m.cos(k)-(y-yo)*m.sin(k))-f*m.cos(k))*m.cos(w)
    parameter[1][4] = -f*m.cos(k)-((y-yo)/f)*((x-xo)*m.sin(k)+(y-yo)*m.cos(k))
    parameter[1][5] = -(x-xo)
    
    return parameter
  

x_apxm = [0,0,0,0]
y_apxm = [0,0,0,0]

R = np.mat(np.zeros((3,3)))
L = np.mat(np.zeros((8,1)))
A = np.mat(np.zeros((8,6)))

f_cor=w_cor=k_cor = 1
flag = 0
             
while (abs(f_cor)>0.000001) | (abs(w_cor)>0.000001) | (abs(k_cor)>0.000001):
    R = r_mat(fo,wo,ko)
    x_apxm,y_apxm = xy_approximate(X,Y,Z,x,y,Xos,Yos,Zos,R)
    
    for i in [0,1,2,3]:
        L[2*i] = x_apxm[i]
        L[2*i+1] = y_apxm[i]
    
    for i in [0,1,2,3]:
        A[2*i:2*i+2,:] = a_parameter(X[i],Y[i],Z[i],Xos,Yos,Zos,x[i],y[i],wo,ko,R)
        
    X_mat = np.mat(np.zeros((6,1)))
    X_mat = (A.T*A).I*A.T*L
    
    f_cor = X_mat[3,0]
    w_cor = X_mat[4,0]
    k_cor = X_mat[5,0]
    
    Xos = Xos + X_mat[0,0]
    Yos = Yos + X_mat[1,0]
    Zos = Zos + X_mat[2,0]
    fo = fo + X_mat[3,0]
    wo = wo + X_mat[4,0]
    ko = ko + X_mat[5,0]
    
    flag += 1
    if flag <= 100:
        print("第 %d 次迭代：f_cor = %f,w_cor = %f,k_cor = %f" %(flag,f_cor,w_cor,k_cor))
    else :
        print("級數不收斂，中間結果為：")
        break

if flag <= 100:
    print("級數收斂，最終結果為：")
print(" Xs=%f,\n Ys=%f,\n Zs=%f,\n f=%f,\n w=%f,\n k=%f" %(Xos,Yos,Zos,fo,wo,ko))
print(" R = ", end='')
print(R)

V = np.mat(np.zeros((8,1)))
V = A*X_mat - L
errorValue = m.sqrt((V.T*V)/(2*4-6))
print("單位權中誤差：%f" %errorValue)

      

&n