双像解析摄影测量

双像解析摄影测量

立体像对

立体像对是由两个不同摄站对同一区域摄影获得的,具有一定重叠度的像片对。

  • 双像各自的摄站分别记为 S 1 ,   S 2 S_1,\ S_2 S1, S2,用 S 1 O 1 ,   S 2 O 2 S_1O_1, \ S_2O_2 S1O1, S2O2分别代表两个主光轴。摄站的连线记为 B B B,称为摄影基线。
  • 对同一地面点A进行讨论时,A向不同摄站投射的光线 A S 1 AS_1 AS1 A S 2 AS_2 AS2称为同名光线; 同名光线与两像片分别的交点, 即在两像片上的成像点为 a 1 ,   a 2 a_1, \ a_2 a1, a2,称为同名像点。
  • 由摄影基线与地面任意一点组成的平面称为该地面点的核面;过像主点的核面称为主核面,过左片主点的核面称左主核面。
  • 核面交像平面的线称为核线,通过同名像点的核线称为同名核线
  • 摄影基线的延长线与像面的交点称为核点,像平面上的所有核线都交于该像片的核点上。
    立体像对

立体像对空间前方交会

由立体像对的左右两影像的内、外方位元素和同名像点的影像坐标量测值来确定该点的物方空间坐标,称为立体像对的空间前方交会。

基本思路

  1. 使用立体像对上的同名像点,就能得到两条同名光线在空间中的方向,这两条光线在空间中一定相交,其交点就是这两个像点的共同地面点。
  2. 若使用共线方程,对于地面点 A A A,成像点 a 1 a_1 a1,可以列出一组共线方程,其包含A的物方坐标 ( X , Y , Z ) (X,Y,Z) (X,Y,Z)和像方坐标 ( x 1 , y 1 ) (x_1, y_1) (x1,y1);若已知内外方位元素,仅凭两条共线方程是无法解出 ( X , Y , Z ) (X,Y,Z) (X,Y,Z)的。但是使用另一张像片,列出另外两条共线方程,一共就有4条方程,从而可以求解三个未知数。

利用点投影系数的空间前方交会法

  1. 坐标定义

    1. 地面测量坐标系定义为 t − X Y Z t-XYZ tXYZ

    2. 立体像对的投影中心(两个摄站)分别为 S 1 ,   S 2 S_1,\ S_2 S1, S2,在 t − X Y Z t-XYZ tXYZ中的坐标分别为 ( X S 1 , Y S 1 , Z S 1 ) (X_{S1},Y_{S1},Z_{S1}) (XS1,YS1,ZS1) ( X S 2 , Y S 2 , Z S 2 ) (X_{S2},Y_{S2},Z_{S2}) (XS2,YS2,ZS2)
      并分别作像空间辅助坐标系 S 1 − X 1 Y 1 Z 1 S_1-X_1Y_1Z_1 S1X1Y1Z1   S 2 − X 2 Y 2 Z 2 \ S_2-X_2Y_2Z_2  S2X2Y2Z2,使坐标轴的方向与地面测量坐标系 t − X Y Z t-XYZ tXYZ一致。

    3. 地面点 A A A在左片的构像为 a 1 a_1 a1,在右片的构像为 a 2 a_2 a2。在地面摄影测量坐标系中的坐标为 ( X A , Y A , Z A ) (X_A,Y_A,Z_A) (XA,YA,ZA),在各自像空间坐标系 S − x y z S-xyz Sxyz 中的坐标分别为 ( x 1 , y 1 , − f ) (x_1, y_1, -f) (x1,y1,f) ( x 2 , y 2 , − f ) (x_2, y_2,-f) (x2,y2,f)
      并由像空系到像空辅系的旋转关系 [ X Y Z ] = R [ x y − f ] \begin{bmatrix}X\\Y\\Z\end{bmatrix}= R\begin{bmatrix}x\\y\\-f\end{bmatrix} XYZ =R xyf ,可以求知构像 a 1 ,   a 2 a_1, \ a_2 a1, a2在像空辅系的坐标 ( X 1 , Y 1 , Z 1 ) (X_1, Y_1, Z_1) (X1,Y1,Z1) ( X 2 , Y 2 , Z 2 ) (X_2, Y_2, Z_2) (X2,Y2,Z2)

      地面摄影测量坐标系:书本上写地面测量坐标系,但它是一个左手系,此处应为右手系,而且在建模工作中,不应该考虑“将模型中的坐标向投影平面转换”这样的问题,这会给建模的工作带来不便。

  2. 在地面测量坐标系中表示摄影基线 B B B
    B B B的各分量 { B X = X S 2 − X S 1 , B Y = Y S 2 − Y S 1 , B Z = Z S 2 − Z S 1 \begin{cases}B_X=X_{S2}-X_{S1}, \\ B_Y=Y_{S2}-Y_{S1}, \\B_Z=Z_{S2}-Z_{S1}\end{cases} BX=XS2XS1,BY=YS2YS1,BZ=ZS2ZS1

  3. 接下来关注由光线 A S 1 AS_1 AS1分别向 Z 1 − S 1 − Y 1 Z_1-S_1-Y_1 Z1S1Y1坐标面、 Z 1 − S 1 − X 1 Z_1-S_1-X_1 Z1S1X1坐标面投影得到的直线与坐标轴围成的两个三角形。
    点投影系数法

    地面点A在 S 1 − X 1 Y 1 Z 1 S_1-X_1Y_1Z_1 S1X1Y1Z1下的坐标为 ( X A − X S 1 , Y A − Y S 1 , Z A − Z S 1 ) (X_A-X_{S1}, Y_A-Y_{S1}, Z_A-Z_{S1}) (XAXS1,YAYS1,ZAZS1),由相似关系可知

    N 1 = d e f S 1 A S 1 a 1 = X A − X S 1 X 1 = Y A − Y S 1 Y 1 = Z A − Z S 1 Z 1 N_1\overset{\mathrm{def}}{=} \frac{S_1 A}{S_1 a_1}=\frac{X_A-X_{S1}}{X_1}=\frac{Y_A-Y_{S1}}{Y_1}=\frac{Z_A-Z_{S1}}{Z_1} N1=defS1a1S1A=X1XAXS1=Y1YAYS1=Z1ZAZS1

    同理对投影光线 A S 2 AS_2 AS2进行讨论可知

    N 2 = d e f S 2 A S 2 a 2 = X A − X S 2 X 2 = Y A − Y S 2 Y 2 = Z A − Z S 2 Z 2 N_2\overset{\mathrm{def}}{=} \frac{S_2 A}{S_2 a_2}=\frac{X_A-X_{S2}}{X_2}=\frac{Y_A-Y_{S2}}{Y_2}=\frac{Z_A-Z_{S2}}{Z_2} N2=defS2a2S2A=X2XAXS2=Y2YAYS2=Z2ZAZS2

    这两个方程实际上是共线条件,在推导共线条件方程时用到过。

    其中的相似比定义为投影系数 N 1 ,   N 2 N_1, \ N_2 N1, N2。一般情况下,不同的点有不同的投影系数,同一点在左右两像片上的投影系数也是不一样的。

  4. 进而可以得到前方交会公式

    [ X A Y A Z A ] = [ N 2 X 2 N 2 Y 2 N 2 Z 2 ] + [ X S 1 Y S 1 Z S 1 ] + [ B X B Y B Z ] = 又 [ N 1 X 1 N 1 Y 1 N 1 Z 1 ] + [ X S 1 Y S 1 Z S 1 ] \begin{bmatrix}X_A\\Y_A\\Z_A\end{bmatrix}=\begin{bmatrix}N_2X_2\\N_2Y_2\\N_2Z_2\end{bmatrix}+\begin{bmatrix}X_{S1}\\Y_{S1}\\Z_{S_1}\end{bmatrix}+\begin{bmatrix}B_X\\B_Y\\B_Z\end{bmatrix}\overset{又}{=}\begin{bmatrix}N_1X_1\\N_1Y_1\\N_1Z_1\end{bmatrix}+\begin{bmatrix}X_{S1}\\Y_{S1}\\Z_{S_1}\end{bmatrix} XAYAZA = N2X2N2Y2N2Z2 + XS1YS1ZS1 + BXBYBZ = N1X1N1Y1N1Z1 + XS1YS1ZS1

    其中又有 { B X = X S 2 − X S 1 , B Y = Y S 2 − Y S 1 , B Z = Z S 2 − Z S 1 \begin{cases}B_X=X_{S2}-X_{S1}, \\ B_Y=Y_{S2}-Y_{S1}, \\B_Z=Z_{S2}-Z_{S1}\end{cases} BX=XS2XS1,BY=YS2YS1,BZ=ZS2ZS1,所以 [ X A Y A Z A ] = [ N 2 X 2 N 2 Y 2 N 2 Z 2 ] + [ X S 1 Y S 1 Z S 1 ] + [ B X B Y B Z ] = 又 [ N 1 X 1 N 1 Y 1 N 1 Z 1 ] + [ X S 1 Y S 1 Z S 1 ] \begin{bmatrix}X_A\\Y_A\\Z_A\end{bmatrix}=\begin{bmatrix}N_2X_2\\N_2Y_2\\N_2Z_2\end{bmatrix}+\begin{bmatrix}X_{S1}\\Y_{S1}\\Z_{S_1}\end{bmatrix}+\begin{bmatrix}B_X\\B_Y\\B_Z\end{bmatrix}\overset{又}{=}\begin{bmatrix}N_1X_1\\N_1Y_1\\N_1Z_1\end{bmatrix}+\begin{bmatrix}X_{S1}\\Y_{S1}\\Z_{S_1}\end{bmatrix} XAYAZA = N2X2N2Y2N2Z2 + XS1YS1ZS1 + BXBYBZ = N1X1N1Y1N1Z1 + XS1YS1ZS1 ,消去 [ X S 1 Y S 1 Z S 1 ] \begin{bmatrix}X_{S1}\\Y_{S1}\\Z_{S_1}\end{bmatrix} XS1YS1ZS1 ,进而有 [ B X B Y B Z ] = [ N 1 X 1 N 1 Y 1 N 1 Z 1 ] − [ N 2 X 2 N 2 Y 2 N 2 Z 2 ] \begin{bmatrix}B_X\\B_Y\\B_Z\end{bmatrix}=\begin{bmatrix}N_1X_1\\N_1Y_1\\N_1Z_1\end{bmatrix}-\begin{bmatrix}N_2X_2\\N_2Y_2\\N_2Z_2\end{bmatrix} BXBYBZ = N1X1N1Y1N1Z1 N2X2N2Y2N2Z2 .
    联立1、3行可以求解投影系数 N 1 = B X Z 2 − B Z X 2 X 1 Z 2 − Z 1 X 2 ,   N 2 = B X Z 1 − B Z X 1 X 1 Z 2 − Z 1 X 2 N_1={{B_X Z_2-B_Z X_2}\over{X_1 Z_2}-Z_1 X_2}, \ N_2=\frac{B_X Z_1-B_Z X_1}{X_1 Z_2-Z_1 X_2} N1=X1Z2Z1X2BXZ2BZX2, N2=X1Z2Z1X2BXZ1BZX1.

  5. 由已知的摄站坐标、投影系数以及像点的像空辅坐标,可以根据前方交会公式求得地面点的地面摄影测量坐标 ( X A , Y A , Z A ) (X_A,Y_A,Z_A) (XA,YA,ZA).

利用共线方程的严格解法

由共线条件方程

x − x 0 = − f a 1 ( X A − X S ) + b 1 ( Y A − Y S ) + c 1 ( Z A − Z S ) a 3 ( X A − X S ) + b 3 ( Y A − Y S ) + c 3 ( Z A − Z S )   y − y 0 = − f a 2 ( X A − X S ) + b 2 ( Y A − Y S ) + c 2 ( Z A − Z S ) a 3 ( X A − X S ) + b 3 ( Y A − Y S ) + c 3 ( Z A − Z S ) x-x_0=-f\frac{a_1(X_A-X_S)+b_1(Y_A-Y_S)+c_1(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)}\\ \ \\ y-y_0=-f\frac{a_2(X_A-X_S)+b_2(Y_A-Y_S)+c_2(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)} xx0=fa3(XAXS)+b3(YAYS)+c3(ZAZS)a1(XAXS)+b1(YAYS)+c1(ZAZS) yy0=fa3(XAXS)+b3(YAYS)+c3(ZAZS)a2(XAXS)+b2(YAYS)+c2(ZAZS)

X A ,   Y A ,   Z A X_A,\ Y_A,\ Z_A XA, YA, ZA分别求偏导得

x = x ( A 0 ) + ∂ x ∂ X A Δ X A + ∂ x ∂ Y A Δ Y A + ∂ x ∂ Z A Δ Z A   y = y ( A 0 ) + ∂ y ∂ X A Δ X A + ∂ y ∂ Y A Δ Y A + ∂ y ∂ Z A Δ Z A x=x(A_0)+ \frac{\partial x}{\partial X_A}\Delta X_A+ \frac{\partial x}{\partial Y_A}\Delta Y_A+ \frac{\partial x}{\partial Z_A}\Delta Z_A \\ \ \\ y=y(A_0)+ \frac{\partial y}{\partial X_A}\Delta X_A+ \frac{\partial y}{\partial Y_A}\Delta Y_A+ \frac{\partial y}{\partial Z_A}\Delta Z_A x=x(A0)+XAxΔXA+YAxΔYA+ZAxΔZA y=y(A0)+XAyΔXA+YAyΔYA+ZAyΔZA

其中, A 0 A_0 A0是地面点A的近似点。对于左右影像上的一对同名像点,可以列出4个上述线性方程,而未知数个数为3,可以用最小二乘法求解。若以n个影像中含有同一个空间点,则可由共2n个方程求解 Δ X , Δ Y , Δ Z \Delta X, \Delta Y, \Delta Z ΔX,ΔY,ΔZ。进而计算地面点A的精确位置,这是一种严格的空间前方交会物方坐标解算方法,需要物方空间坐标的初值。


相对定向

相对定向元素

  • 相对定向:
    利用立体像对中摄影时存在的同名光线对应相交的关系,通过量测像点坐标,以解析计算的方法(不需野外控制点),求解两像片的像对方位元素的过程称为解析像对定向

    相对定向的目的是建立一个与被摄物体相似的几何模型,以确定模型点的三维坐标。

  • 相对定向元素:
    是确定相邻像片相对位置和姿态的要素

    • 连续法相对定向元素:
      将像空间辅助坐标系的原点取在左摄站点上,其坐标轴系保持与立体像对中左片的像空间坐标系的轴系分别重合,则左片相对于像空辅系的外方位元素中 角元素均为零;
      即:以左像为基准,求出右片相对于左片的外方位元素
      两像片的外方位元素相对差为:
      B X = X S 2 − X S 1 ,   B Y = Y S 2 − Y S 1 ,   B Z = Z S 2 − Z S 1 ,   Δ φ = φ 2 − φ 1 = φ 2 ,   Δ ω = ω 2 ,   Δ κ = κ 2 B_X=X_{S2}-X_{S1}, \ B_Y=Y_{S2}-Y_{S1}, \ B_Z=Z_{S2}-Z_{S1}, \ \Delta\varphi=\varphi_2-\varphi_1=\varphi_2,\ \Delta\omega=\omega_2, \ \Delta\kappa =\kappa_2 BX=XS2XS1, BY=YS2YS1, BZ=ZS2ZS1, Δφ=φ2φ1=φ2, Δω=ω2, Δκ=κ2
      以上 B Y , B Z , φ 2 , ω 2 , κ 2 B_Y,B_Z,\varphi_2,\omega_2,\kappa_2 BY,BZ,φ2,ω2,κ2称为相对定向元素, B X B_X BX只决定立体模型的大小,不影响相对方位,采用连续法相对定向元素的相对定向方法,称为连续法相对定向
      连续法相对定向
      单独法相对定向
    • 单独法相对定向元素:
      将像空辅系的原点取在左摄站点上,坐标系的X轴保持与摄影基线B的方向重合,并使坐标系的Z轴落在像对左片的主核面内,使得两像片的外方位元素保持 Y S 2 = Y S 1 ,   Z S 2 = Z S 1 , ω 1 = 0 Y_{S2}=Y_{S1},\ Z_{S2}=Z_{S1},\omega_1=0 YS2=YS1, ZS2=ZS1,ω1=0的相对关系,这时候相对变化的元素 φ 1 , κ 1 , φ 2 , ω 2 , κ 2 \varphi_1,\kappa_1,\varphi_2,\omega_2,\kappa_2 φ1,κ1,φ2,ω2,κ2称为相对定向元素,B只决定模型比例尺,不改变相对方位

共面条件

如图为相对定向示意图, m 1 ,   m 2 m_1, \ m_2 m1, m2为模型点 M M M在左、右两幅影像上的构像, S 1 m 1 ,   S 2 m 2 S_1m_1,\ S_2m_2 S1m1, S2m2表示一对同名光线,它们与空间基线 S 1 S 2 S_1S_2 S1S2共面,用三个向量 R 1 ,   R 2 ,   B \bold R_1, \ \bold R_2, \ \bold B R1, R2, B的混合积表示这个平面,
B ⋅ ( R 1 × R 2 ) = 0 \bold B \cdot (\bold R_1\times \bold R_2)=0 B(R1×R2)=0

共面条件

向量混合积:

向量 a ⃗ ,   b ⃗ ,   c ⃗ \vec a, \ \vec b, \ \vec c a , b , c 的混合积为 a ⃗ ⋅ ( b ⃗ × c ⃗ ) \vec a\cdot(\vec b \times \vec c) a (b ×c ),表示三个向量形成的平行六面体的体积

共面条件方程

使用 a ⃗ ⋅ ( b ⃗ × c ⃗ ) = 0 \vec a\cdot(\vec b \times \vec c)=0 a (b ×c )=0 约束,使得平行六面体的体积为0,即可得到三个向量共面的条件方程

a ⃗ × b ⃗ = ∣ i ⃗      j ⃗     k ⃗ a x   a y   a z b x   b y   b z ∣ = A 11 i ⃗ + A 12 j ⃗ + A 13 k ⃗ \vec a \times \vec b=\begin{vmatrix}\vec i\ \ \ \ \vec j\ \ \ \vec k \\a_x \ a_y \ a_z \\b_x \ b_y \ b_z \\\end{vmatrix}=A_{11}\vec i+A_{12}\vec j+A_{13} \vec k a ×b = i     j    k ax ay azbx by bz =A11i +A12j +A13k ,得 a ⃗ ⋅ ( b ⃗ × c ⃗ ) = a ⃗ ⋅ [ A 11   A 12   A 13 ] T = A 11 a x + A 12 a y + A 13 a z = ∣   a x   a y   a z b x   b y   b z c x   c y   c z ∣ \vec a \cdot (\vec b \times \vec c)=\vec a \cdot [A_{11}\ A_{12} \ A_{13}]^\mathrm T=A_{11}a_x+A_{12}a_y+A_{13}a_z=\begin{vmatrix}\ a_x\ a_y\ a_z\\b_x \ b_y \ b_z \\c_x \ c_y \ c_z \\\end{vmatrix} a (b ×c )=a [A11 A12 A13]T=A11ax+A12ay+A13az=  ax ay azbx by bzcx cy cz

进一步 B ⋅ ( R 1 × R 2 ) = ∣   B X   B Y   B Z X 1    Y 1     Z 1 X 2    Y 2     Z 2 ∣ = 0 \bold B\cdot(\bold R_1\times \bold R_2)=\begin{vmatrix}\ B_X\ B_Y\ B_Z\\X_1 \ \ Y_1 \ \ \ Z_1 \\X_2 \ \ Y_2 \ \ \ Z_2 \\\end{vmatrix}=0 B(R1×R2)=  BX BY BZX1  Y1   Z1X2  Y2   Z2 =0,其中 [ X 1 Y 1 Z 1 ] = R [ x 1 y 1 − f 1 ] ,   [ X 2 Y 2 Z 2 ] = R ′ [ x 2 y 2 − f 2 ] \begin{bmatrix}X_1\\Y_1\\Z_1\end{bmatrix}= R\begin{bmatrix}x_1\\y_1\\-f_1\end{bmatrix}, \ \begin{bmatrix}X_2\\Y_2\\Z_2\end{bmatrix}= R'\begin{bmatrix}x_2\\y_2\\-f_2\end{bmatrix} X1Y1Z1 =R x1y1f1 ,  X2Y2Z2 =R x2y2f2 R R R中涉及的角元素应都为0,而 R ′ R' R中的角元素为 Δ φ ,   Δ ω ,   Δ κ \Delta \varphi, \ \Delta\omega, \ \Delta\kappa Δφ, Δω, Δκ,而 B X ,   B Y B_X, \ B_Y BX, BY又可由 μ ,   ν \mu, \ \nu μ, ν求解

B Y = tan ⁡ μ ⋅ B X ,   B Z = B X cos ⁡ μ ⋅ tan ⁡ ν B_Y=\tan \mu \cdot B_X, \ B_Z={B_X\over\cos\mu}\cdot\tan\nu BY=tanμBX, BZ=cosμBXtanν

当角元素为小角度时, B Y ≈ B X ⋅ μ ,   B Z ≈ B X ⋅ ν B_Y\approx B_X\cdot\mu, \ B_Z\approx B_X \cdot \nu BYBXμ, BZBXν

对于单独法相对定向系统,其 Y S 2 = Y S 1 ,   Z S 2 = Z S 1 , ω 1 = 0 Y_{S2}=Y_{S1},\ Z_{S2}=Z_{S1},\omega_1=0 YS2=YS1, ZS2=ZS1,ω1=0,所以 B Y = 0 ,   B Z = 0 B_Y=0,\ B_Z=0 BY=0, BZ=0,则共面条件方程为 ∣ B X    0       0 X 1    Y 1     Z 1 X 2    Y 2     Z 2 ∣ = 0 \begin{vmatrix}B_X\ \ 0\ \ \ \ \ 0\\X_1 \ \ Y_1 \ \ \ Z_1 \\X_2 \ \ Y_2 \ \ \ Z_2 \\\end{vmatrix}=0 BX  0     0X1  Y1   Z1X2  Y2   Z2 =0,进一步简化得 ∣ Y 1   Z 1 Y 2   Z 2 ∣ = 0 \begin{vmatrix}Y_1\ Z_1\\Y_2\ Z_2\end{vmatrix}=0 Y1 Z1Y2 Z2 =0

共面条件方程的线性化求解

共面条件方程是关于 μ ,   ν ,   φ 2 ,   ω 2 ,   κ 2 \mu,\ \nu, \ \varphi_2, \ \omega_2, \ \kappa_2 μ, ν, φ2, ω2, κ2 的方程, F ( μ ,   ν ,   φ 2 ,   ω 2 ,   κ 2 ) = 0 F(\mu,\ \nu, \ \varphi_2, \ \omega_2, \ \kappa_2) = 0 F(μ, ν, φ2, ω2, κ2)=0,将其泰勒一阶展开后得

F = F 0 + ∂ F ∂ φ d φ + ∂ F ∂ ω d ω + ∂ F ∂ κ d κ + ∂ F ∂ μ d μ + ∂ F ∂ ν d ν F=F_0+{\partial F \over \partial\varphi}\mathrm d \varphi+{\partial F \over \partial\omega}\mathrm d \omega+{\partial F \over \partial\kappa}\mathrm d \kappa+{\partial F \over \partial\mu}\mathrm d \mu+{\partial F \over \partial\nu}\mathrm d \nu F=F0+φFdφ+ωFdω+κFdκ+μFdμ+νFdν

其中, F 0 F_0 F0是根据相对定向元素近似值得到的 F F F的初值, d φ = φ 2 − φ 2 0 ,   ⋯ \mathrm d\varphi=\varphi_2-\varphi_2^0,\ \cdots dφ=φ2φ20, 

进一步计算各项偏导数,其中 Q Q Q是模型点的上下视差,当一个立体像对完成相对定向时, Q = 0 Q=0 Q=0;当一个立体像对未完成相对定向,即同名光线不相交, Q ≠ 0 Q\neq0 Q=0

Q = F 0 N 2 B X Z 2 = N 1 Y 1 − N 2 Y 2 − B Y Q=\frac{F_0N_2}{B_XZ_2}=N_1Y_1-N_2Y_2-B_Y Q=BXZ2F0N2=N1Y1N2Y2BY

在解算过程中,将 Q Q Q视为观测值,求解 d μ ,   d ν ,   d φ 2 ,   d ω 2 ,   d κ 2 \mathrm d \mu,\ \mathrm d \nu, \ \mathrm d \varphi_2, \ \mathrm d \omega_2, \ \mathrm d \kappa_2 dμ, dν, dφ2, dω2, dκ2,需观测5对以上同名像点的像点坐标,列立误差方程,按间接平差求解。实际计算中用迭代法。

解算步骤

模型点坐标的计算

在相对定向元素正确求得之后,立体模型就得以建立起来。由于求解相对定向元素的过程中,基线分量 B X B_X BX是任意选取的,因此模型的比例尺也是任意的。计算单个模型中各模型点的坐标可按照立体像对前方交会公式进行


核线影像

重叠影像上的同名像点必定位于同名核线上,确定了同名核线,就可以在同名核线上搜索同名像点

确定同名核线的方法

  • 基于影像几何纠正的核线解析关系
    将两像片本互不平行的核线投影到一个平行于摄影基线的水平上,投影核线相互平行
    在倾斜影像上采样核线影像的灰度,赋值给水平核线,这样就得到了水平影像对上的同名核线
  • 基于共面条件的同名核线几何关系
    已知左影像上具有构像p,确定左影像上通过该点的核线 l 与右片上的同名核线 l’
    该问题转化为等价问题:确定左片上另一点q,右片同名核线上的两个点p’, q’ (p, p’; q, q’)不一定为同名像点
    利用共面条件来解决问题 B ⋅ ( S p × S q ) = 0 B\cdot (Sp\times Sq)=0 B(Sp×Sq)=0

解析绝对定向

在相对定向后,得到了立体模型,但其在地面坐标系中的方位和大小是不确定的。确定立体模型在地面坐标系中的方位和大小,则需把模型坐标转换为地面坐标,称绝对定向。

其目的是将建立的模型纳入到地面坐标系,并归化为规定的比例尺。

坐标变换

一般情况下,模型坐标是属于右手空间直角的摄影测量坐标系,而地面坐标为左手空间直角坐标系。因此在进行模型的绝对定向之前,需要做空间直角坐标系的转换,即将地面测量坐标系转换为地面摄影测量坐标系

  1. 地面坐标点一般给到的是大地坐标 ( B , L , H ) (B,L,H) (B,L,H),首先转换到投影平面坐标(高斯-克吕格或UTM投影)
  2. 再从高斯-克吕格坐标转换到摄影测量坐标

在绝对定向完成,模型点地面坐标需要归化到大地坐标时,则从摄影测量坐标转换到高斯坐标再转回大地坐标。

在归化过程中,需要对摄影测量所得高程引入地球曲率改正。

可以看出,以下的讨论的地面坐标系是地面摄影测量坐标系,以某个控制点为原点。

解析绝对定向元素

在基本的立体像对中,每个影像具有6个外方位元素,则立体像对具有12个外方位元素。通过像对定向,可以求得5个相对定向元素,确定模型的绝对位置和方位,另需7个绝对定向元素,包括平移、旋转和缩放。

其中平移元素为 Δ X ,   Δ Y ,   Δ Z \Delta X, \ \Delta Y, \ \Delta Z ΔX, ΔY, ΔZ,旋转元素为 Φ ,   Ω ,   K \Phi, \ \Omega, \ \Kappa Φ, Ω, K,缩放元素为 λ \lambda λ。其中旋转元素可以进一步计算其方向余弦,得到旋转矩阵 R R R。则模型点与地面点之间的空间相似变换(布尔莎模型)为:

[ X t p Y t p Z t p ] = λ R [ X P Y P Z P ] + [ Δ X Δ Y Δ Z ] = λ [ a 1 a 2 a 3 b 1 b 2 b 3 c 1 c 2 c 3 ] [ X P Y P Z P ] + [ Δ X Δ Y Δ Z ] \begin{bmatrix} X_{tp}\\Y_{tp}\\Z_{tp} \end{bmatrix}=\lambda R \begin{bmatrix} X_P\\Y_P\\Z_P\\ \end{bmatrix}+\begin{bmatrix} \Delta X\\ \Delta Y\\ \Delta Z\\ \end{bmatrix} =\lambda\begin{bmatrix} a_1&a_2&a_3\\ b_1&b_2&b_3\\ c_1&c_2&c_3 \end{bmatrix} \begin{bmatrix} X_P\\Y_P\\Z_P\\ \end{bmatrix}+\begin{bmatrix} \Delta X\\ \Delta Y\\ \Delta Z\\ \end{bmatrix} XtpYtpZtp =λR XPYPZP + ΔXΔYΔZ =λ a1b1c1a2b2c2a3b3c3 XPYPZP + ΔXΔYΔZ

其中 ( X t p , Y t p , Z t p ) (X_{tp},Y_{tp}, Z_{tp}) (Xtp,Ytp,Ztp)为地面控制点的地面摄影测量坐标, ( X P , Y P Z P ) (X_P,Y_PZ_P) (XP,YPZP)为对应模型点的摄影测量坐标。

以上方程揭示了模型点与地面点的关系,其中包含7个未知参数。可以通过建立已知的地面控制点与同名的模型点之间的关系,来反解这7个参数。

将上述方程表示成关于7个未知参数的方程得 F ( Δ X , Δ Y , Δ Z , Φ , Ω , K , λ ) = 0 F(\Delta X, \Delta Y, \Delta Z, \Phi, \Omega, \Kappa, \lambda)=0 F(ΔX,ΔY,ΔZ,Φ,Ω,K,λ)=0,并对其进行一阶泰勒展开得线性形式:

F = F 0 + ∂ F ∂ Δ X d Δ X + ∂ F ∂ Δ Y d Δ Y + ∂ F ∂ Δ Z d Δ Z + ∂ F ∂ Φ d Φ + ∂ F ∂ Ω d Ω + ∂ F ∂ K d K + ∂ F ∂ λ d λ F=F_0+{\partial F \over \partial\Delta X}\mathrm d \Delta X+{\partial F \over \partial\Delta Y}\mathrm d \Delta Y+{\partial F \over \partial\Delta Z}\mathrm d \Delta Z+{\partial F \over \partial\Phi}\mathrm d \Phi+{\partial F \over \partial\Omega}\mathrm d \Omega+{\partial F \over \partial\Kappa}\mathrm d \Kappa+{\partial F \over \partial\lambda}\mathrm d \lambda F=F0+ΔXFdΔX+ΔYFdΔY+ΔZFdΔZ+ΦFdΦ+ΩFdΩ+KFdK+λFdλ

需注意,上式中的 F F F事实上是一个函数向量, F ⃗ = [ X t p Y t p Z t p ] \vec F=\begin{bmatrix}X_{tp}\\Y_{tp}\\Z_{tp}\end{bmatrix} F = XtpYtpZtp ,针对其每一个分量,都可以展开得到一条单独的方程,进一步可写出误差项 v x ,   v y ,   v z v_x, \ v_y, \ v_z vx, vy, vz

v x = ∂ X ∂ Δ X d Δ X + ∂ X ∂ Φ d Φ + ∂ X ∂ Ω d Ω + ∂ X ∂ K d K + ∂ X ∂ λ d λ − l X v y = ∂ Y ∂ Δ Y d Δ Y + ∂ Y ∂ Φ d Φ + ∂ Y ∂ Ω d Ω + ∂ Y ∂ K d K + ∂ Y ∂ λ d λ − l Y v z = ∂ Z ∂ Δ Z d Δ Z + ∂ Z ∂ Φ d Φ + ∂ Z ∂ Ω d Ω + ∂ Z ∂ K d K + ∂ Z ∂ λ d λ − l Z l X = X t p − Δ X − λ X ′ l Y = Y t p − Δ Y − λ Y ′ l Z = Z φ − Δ Z − λ Z ′ ( X ′ Y ′ Z ′ ) = ( a 1 a 2 a 3 b 1 b 2 b 3 c 1 c 2 c 3 ) ( X P Y P Z P ) \begin{gathered}v_{x}=\frac{\partial X}{\partial\Delta X}\mathrm{d}\Delta X+\frac{\partial X}{\partial{\Phi}}\mathrm{d}{\Phi}+\frac{\partial X}{\partial{\Omega}}\mathrm{d}{\Omega}+\frac{\partial X}{\partial K}\mathrm{d}K+\frac{\partial X}{\partial\lambda}\mathrm{d}\lambda-l_{X} \\v_{y}=\frac{\partial Y}{\partial\Delta Y}\mathrm{d}\Delta Y+\frac{\partial Y}{\partial{\Phi}}\mathrm{d}{\Phi}+\frac{\partial Y}{\partial\Omega}\mathrm{d}{\Omega}+\frac{\partial Y}{\partial K}\mathrm{d}K+\frac{\partial Y}{\partial\lambda}\mathrm{d}\lambda-l_{Y} \\v_{z}=\frac{\partial Z}{\partial\Delta Z}\mathrm{d}\Delta Z+\frac{\partial Z}{\partial\Phi}\mathrm{d}\Phi+\frac{\partial Z}{\partial\Omega}\mathrm{d}\Omega+\frac{\partial Z}{\partial K}\mathrm{d}K+\frac{\partial Z}{\partial\lambda}\mathrm{d}\lambda-l_{Z} \\\begin{aligned}l_{X}&=X_{tp}-\Delta X-\lambda X^{\prime}\\l_{Y}&=Y_{tp}-\Delta Y-\lambda Y^{\prime}\\l_{Z}&=Z_{\varphi}-\Delta Z-\lambda Z^{\prime}\end{aligned}\quad\begin{pmatrix}X^{\prime}\\Y^{\prime}\\Z^{\prime}\end{pmatrix}=\begin{pmatrix}a_{1}&a_{2}&a_{3}\\b_{1}&b_{2}&b_{3}\\c_{1}&c_{2}&c_{3}\end{pmatrix}\begin{pmatrix}X_{P}\\Y_{P}\\Z_{P}\end{pmatrix} \end{gathered} vx=ΔXXdΔX+ΦXdΦ+ΩXdΩ+KXdK+λXdλlXvy=ΔYYdΔY+ΦYdΦ+ΩYdΩ+KYdK+λYdλlYvz=ΔZZdΔZ+ΦZdΦ+ΩZdΩ+KZdK+λZdλlZlXlYlZ=XtpΔXλX=YtpΔYλY=ZφΔZλZ XYZ = a1b1c1a2b2c2a3b3c3 XPYPZP

不加推导地给出误差方程的矩阵形式为

[ v x v y v z ] = [ 1 0 0 X ′ − Z ′ 0 − Y 0 1 0 Y ′ 0 − Z ′ X 0 0 1 Z ′ X ′ Y ′ 0 ] [ d Δ X d Δ Y d Δ Z d λ d Φ d Ω d K ] − [ l x l y l z ] \begin{bmatrix} v_x \\ v_y \\ v_z \end{bmatrix}=\begin{bmatrix} 1 & 0 & 0 & X' & -Z' & 0 & -Y \\ 0 & 1 & 0 & Y' & 0 & -Z'& X \\ 0 & 0 & 1 & Z' & X' & Y' & 0\\ \end{bmatrix}\begin{bmatrix} \mathrm d \Delta X\\ \mathrm d \Delta Y\\ \mathrm d \Delta Z\\ \mathrm d \lambda\\ \mathrm d \Phi\\ \mathrm d \Omega\\ \mathrm d \Kappa\\ \end{bmatrix} -\begin{bmatrix} l_x\\l_y\\l_z \end{bmatrix} vxvyvz = 100010001XYZZ0X0ZYYX0 dΔXdΔYdΔZdλdΦdΩdK lxlylz

式中 X ′ ,   Y ′ ,   Z ′ X', \ Y', \ Z' X, Y, Z是模型点 ( X P , Y P , Z P ) (X_P, Y_P, Z_P) (XP,YP,ZP)经旋转变换 R R R后得到的坐标,上式是单元模型空间相似变换的误差方程。将 ( X P , Y P , Z P ) (X_P, Y_P, Z_P) (XP,YP,ZP) 作为观测值, l x , l y , l z l_x, ly,l_z lx,ly,lz为方程的常数项, d Δ X , d Δ Y , d Δ Z , d Φ , d Ω , d K , d λ \mathrm d \Delta X,\mathrm d \Delta Y,\mathrm d \Delta Z,\mathrm d \Phi,\mathrm d \Omega,\mathrm d \Kappa,\mathrm d \lambda dΔX,dΔY,dΔZ,dΦ,dΩ,dK,dλ 作为待定参数近似值的改正数。

坐标重心化

坐标重心化是区域网平差中经常采用的一种数据预处理方法

  1. 能够减少模型点坐标在计算过程中的有效位数,保证计算的精度
  2. 采用了重心化坐标,可以使得法方程系数简化,个别项的数值变成0,部分未知数可以分开求解,从而提高了计算速度。

取单元模型中全部控制点的摄影测量坐标系和地面摄影测量坐标系,计算它们的重心坐标(即计算平均坐标)

X t p g = ∑ X t p n , Y t p g = ∑ Y t p n , Z t p g = ∑ Z t p n , X g = ∑ X P n , Y g = ∑ Y P n , Z g = ∑ Z P n \begin{aligned}X_{tpg}&=\frac{\sum X_{tp}}{n},Y_{tpg}&=\frac{\sum Y_{tp}}{n},Z_{tpg}&=\frac{\sum Z_{tp}}{n},\\X_{g}&=\frac{\sum X_{P}}{n},Y_{g}&=\frac{\sum Y_{P}}{n},Z_{g}&=\frac{\sum Z_{P}}{n}\end{aligned} XtpgXg=nXtp,Ytpg=nXP,Yg=nYtp,Ztpg=nYP,Zg=nZtp,=nZP

其中 X t p g , Y t p g , Z t p g X_{tpg},Y_{tpg},Z_{tpg} Xtpg,Ytpg,Ztpg为各控制点的地面摄影测量坐标系重心坐标, X g , Y g , Z g X_g,Y_g,Z_g Xg,Yg,Zg为各控制点的摄影测量坐标系重心坐标, n n n为控制点个数。

在重心化后,各模型点坐标应表示为相对于 由各控制点所得的重心 的偏移 { X ‾ = X P − X g Y ‾ = Y p − Y g Z ‾ = Z p − Z g \begin{cases}\begin{aligned}\overline{X}&=X_P-X_g\\\overline{Y}&=Y_p-Y_g\\\overline{Z}&=Z_p-Z_g\end{aligned}\end{cases} XYZ=XPXg=YpYg=ZpZg,各控制点的地面摄影测量坐标也应表示为相对于其重心的偏移 { X ‾ t p = X t p − X t p g Y ‾ t p = Y t p − Y t p g Z ‾ t p = Z t p − Z t p g \begin{cases}\begin{aligned}\overline{X}_{tp}&=X_{tp}-X_{tpg}\\\overline{Y}_{tp}&=Y_{tp}-Y_{tpg}\\\overline{Z}_{tp}&=Z_{tp}-Z_{tpg}\end{aligned}\end{cases} XtpYtpZtp=XtpXtpg=YtpYtpg=ZtpZtpg

对于误差方程式,若直接使用重心化坐标表示,则化为以下形式

( v X v Y v Z ) = ( 1 0 0 X ‾ − Z ‾ 0 − Y ‾ 0 1 0 Y ‾ 0 − Z ‾ X ‾ 0 0 1 Z ‾ X ‾ Y ‾ 0 ) ( d Δ X d Δ Y d Δ Z d λ d Φ d Φ d K ) − ( l X l Y l z ) 其中, ( l X l Y l Z ) = ( X ‾ p Y ‾ p Z ‾ p ) − λ R ( X ‾ Y ‾ Z ‾ ) − ( Δ X Δ Y Δ Z ) \begin{aligned}\begin{pmatrix}v_X\\v_Y\\v_Z\end{pmatrix}=\begin{pmatrix}1&0&0&\overline{X}&-\overline{Z}&0&-\overline{Y}\\0&1&0&\overline{Y}&0&-\overline{Z}&\overline{X}\\0&0&1&\overline{Z}&\overline{X}&\overline{Y}&0\end{pmatrix}\begin{pmatrix}d\Delta X\\d\Delta Y\\d\Delta Z\\d\lambda\\d\Phi\\d\Phi\\dK\end{pmatrix}-\begin{pmatrix}l_X\\l_Y\\l_z\end{pmatrix}\\其中,\begin{pmatrix}l_X\\l_Y\\l_Z\end{pmatrix}=\begin{pmatrix}\overline{X}_p\\\overline{Y}_p\\\overline{Z}_p\end{pmatrix}-\lambda R\begin{pmatrix}\overline{X}\\\overline{Y}\\\overline{Z}\end{pmatrix}-\begin{pmatrix}\Delta X\\\Delta Y\\\Delta Z\end{pmatrix}\end{aligned} vXvYvZ = 100010001XYZZ0X0ZYYX0 dΔXdΔYdΔZdλdΦdΦdK lXlYlz 其中, lXlYlZ = XpYpZp λR XYZ ΔXΔYΔZ

坐标重心化并没有省去旋转,只是把旋转部分转移到了常数项

绝对定向元素的解算

绝对定向元素的解算过程实际上就是列立至少7个误差方程,求解7个待定参数。在控制点最少的情况下,需要两个平面高程控制点和一个高程控制点(一共可以列出7个方程)。在具有多余控制点的情况下,按照间接平差原理进行计算, V = A X − L V=AX-L V=AXL X = [ d Δ X , d Δ Y , d Δ Z , d λ , d Φ , d Ω , d K ] T X=[\mathrm d \Delta X, \mathrm d \Delta Y, \mathrm d \Delta Z, \mathrm d \lambda,\mathrm d \Phi, \mathrm d \Omega, \mathrm d \Kappa]^\mathrm T X=[dΔX,dΔY,dΔZ,dλ,dΦ,dΩ,dK]T,给出的解为 X = ( A T A ) − 1 A T L X=(A^\mathrm T A)^{-1}A^\mathrm TL X=(ATA)1ATL

采用重心化坐标,模型点和地面点的坐标系原点都移到各自的重心处,使得 A T A A^\mathrm TA ATA以及 A T L A^TL ATL中关于平移量的元素为0,此时 d Δ X , d Δ Y , d Δ Z \mathrm d \Delta X, \mathrm d \Delta Y, \mathrm d \Delta Z dΔX,dΔY,dΔZ均为0,因而免去了求解

因为 d Δ X , d Δ Y , d Δ Z \mathrm d \Delta X, \mathrm d \Delta Y, \mathrm d \Delta Z dΔX,dΔY,dΔZ是在各控制点的附近进行展开的改正数,本来各控制点在地面摄影测量坐标系下的坐标是相对某一控制点而言,而在摄影测量坐标系下的坐标是相对模型的原点而言,进行重心化后,双方都是以各自的重心为原点,因此各模型点的平移量就为0

d λ \mathrm d \lambda dλ可以直接求出,这样实际上法方程及法方程的常数项极大简化,待求参数只剩下角元素改正数 d Φ , d Ω , d K d \Phi, \mathrm d \Omega, \mathrm d \Kappa dΦ,dΩ,dK

模型比例尺归化,是坐标重心化后,选择两个相距最远的地面控制点和对应的模型点计算模型比例因子 λ ′ \lambda' λ,并用 λ ′ \lambda' λ逐个乘以模型点坐标,从而使模型与地面的比例基本一致。

在角元素解算完成后,线元素只是两个重心的差值 { Δ X = X g − X t p g Δ Y = Y g − Y t p g Δ Z = Z g − Z t p g \begin{cases}\begin{aligned}\Delta{X}&=X_g-X_{tpg}\\\Delta{Y}&=Y_g-Y_{tpg}\\\Delta{Z}&=Z_g-Z_{tpg}\end{aligned}\end{cases} ΔXΔYΔZ=XgXtpg=YgYtpg=ZgZtpg

角元素的具体解算流程就按照迭代法进行

  • 确定待定参数的初值: Φ 0 = Ω 0 = K 0 = 0 ,   λ 0 = 1 , Δ X = Δ Y = Δ Z = 0 \Phi^0=\Omega^0=\Kappa^0=0,\ \lambda^0=1 , \Delta X=\Delta Y=\Delta Z = 0 Φ0=Ω0=K0=0, λ0=1,ΔX=ΔY=ΔZ=0
  • 计算地面摄影测量坐标系的重心坐标和控制点重心化坐标
  • 计算摄影测量坐标系的重心坐标和控制点重心化坐标
  • 计算旋转矩阵 R R R以及法方程系数阵、常数项
  • 法方程求逆,给出参数的解
  • 计算参数新值(初值+解)
  • 判断改正数 d Φ , d Ω , d K d \Phi, \mathrm d \Omega, \mathrm d \Kappa dΦ,dΩ,dK是否收敛,若收敛则完成计算,否则以新值作为新的初值进行下一轮迭代

参考文献

  • 《摄影测量学》王双亭 测绘出版社 2017
  • 《摄影测量学(第三版)》潘励 武汉大学出版社 2023
  • 21
    点赞
  • 32
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值