平差作業八

gnss 基線平差練習，分別使用極大權約制， 與使用條件平差二種作法，比較成果差異。

基線觀測方程式

基線是三維空間中連接二個點的向量， 而三維空間中的一點由三個維度的座標決定。 所以一條基線方程式應該在三個維度上展開， 成為三條觀測方程式； 對應的未知點也應該在三維上展開，成為三個未知數。

程式中函數 expand_3d_design_matrix 可以將針對向量的設計矩陣， 在三維上展開成 x y z 座標的設計矩陣。

間接平差模型

程式

% constrain adjustment

design_vector_matrix = [
  -1   1   0   0   0
   0  -1   1   0   0
   0   0  -1   1   0
   0   0   0  -1   1
  -1   0   1   0   0
   0  -1   0   1   0
   0   0  -1   0   1
];

function m = expand_3d_design_matrix(mv)
  %% expand [1 0 -1 0 0]
  %% to [1   0   0   0   0   0  -1   0   0   0   0   0   0   0   0
  %%     0   1   0   0   0   0   0  -1   0   0   0   0   0   0   0
  %%     0   0   1   0   0   0   0   0  -1   0   0   0   0   0   0]
  m = zeros(size(mv)*3);

  for i = 1:rows(m)
    for j = 1:columns(m)
      if (mod(j,3) == mod(i,3))
        m(i,j) = mv(ceil(i/3), ceil(j/3));
      else
        m(i,j) = 0;
      end
    end
  end
end

design_magrix = expand_3d_design_matrix(design_vector_matrix);

baseline_vector_vector = [
 11488.973 -13224.866 19.446
 -5605.975 -6047.976 -7.814
 13505.271 28936.647 5.247
 -14911.055 -35571.005 -6.257
 5882.998 -19272.842 11.632
 7899.296 22888.671 -2.567
 -1405.784 -6634.358 -1.010
];

baseline_vector = reshape(
                      baseline_vector_vector',
                      1, prod(size(baseline_vector_vector)))';

constrain_matrix = expand_3d_design_matrix([1 0 0 0 0]);
constrain_value = [177597.421 2639738.678 25.375]';

priority = 100;

normal_matrix = (
  design_magrix'*design_magrix +
  priority * constrain_matrix' * constrain_matrix
);

point = inv(normal_matrix) * (
          design_magrix' * baseline_vector +
          priority * constrain_matrix' * constrain_value
        );

point_vector = reshape(point, 3, rows(point)/3)';

observation_residule = baseline_vector - design_magrix * point;
variance_post = observation_residule' * observation_residule / (
    rows(baseline_vector) - rows(point) + rows(constrain_value)
);
point_covariance_matrix = variance_post * inv(normal_matrix);

結果

用 matlab 計算結果誤差極小。 若照題目看來，實際上是沒有誤差的， 各基線在題目給定的有效位數上閉合差為 0。 但用 matlab 計算會產生誤差， 實際上是因為二進位浮點數無法精確表示十進位小數。 例如許多程式語言中，0.1+0.2-0.3≠0。 因此這些誤差是沒有意義的，實際上都能視作 0。

\begin{bmatrix}
A \\ B \\ C \\ D \\ E
\end{bmatrix}
=
\begin{bmatrix}
1.7760 \times 10^{+05} & 2.6397 \times 10^{+06} & 2.5375 \times 10^{+01} \\
1.8909 \times 10^{+05} & 2.6265 \times 10^{+06} & 4.4821 \times 10^{+01} \\
1.8348 \times 10^{+05} & 2.6205 \times 10^{+06} & 3.7007 \times 10^{+01} \\
1.9699 \times 10^{+05} & 2.6494 \times 10^{+06} & 4.2254 \times 10^{+01} \\
1.8207 \times 10^{+05} & 2.6138 \times 10^{+06} & 3.5997 \times 10^{+01}
\end{bmatrix}

\\

\hat {\sigma_0} = 4.4228 \times 10 ^ {-10}

\\

\begin{bmatrix}
\sigma_{A_x}^2 & \sigma_{A_y}^2 & \sigma_{A_z}^2 \\
\vdots \\
\sigma_{E_x}^2 & \sigma_{E_y}^2 & \sigma_{E_z}^2
\end{bmatrix}
=
\begin{bmatrix}
1.9561 \times 10^{-23} & 1.9561 \times 10^{-23} & 1.9561 \times 10^{-23} \\
1.2111 \times 10^{-19} & 1.2111 \times 10^{-19} & 1.2111 \times 10^{-19} \\
1.2111 \times 10^{-19} & 1.2111 \times 10^{-19} & 1.2111 \times 10^{-19} \\
1.7700 \times 10^{-19} & 1.7700 \times 10^{-19} & 1.7700 \times 10^{-19} \\
2.2357 \times 10^{-19} & 2.2357 \times 10^{-19} & 2.2357 \times 10^{-19}
\end{bmatrix}

條件觀測平差

條件觀測平差是計算量較少，但較難以程式計算的平差方法。 需要針對觀測量，列出所有不相關的條件式。

題目中有 5 未知數，7 觀測方程式，1 條件， 應列出 7-5+1=3 條條件式。 各條件式寫成向量如下，純量形式一樣是在三維上展開。

v_1 + v_2 - v_5 = \vec{A B} + \vec{B C} - \vec{A C} \\
v_2 + v_3 - v_6 = \vec{B C} + \vec{C D} - \vec{B D} \\
v_3 + v_4 - v_7 = \vec{C D} + \vec{D E} - \vec{C E}

解法則是要先解出聯系數矩陣 K， 才能解殘差 V。 因為是條件式，解出來只有觀測量的殘差， 要再將殘差修正後的觀測量組合才能得出絕對點座標。 而本次作業等權，所以 P 矩陣為單位矩陣，可以忽略。

A V = W \\

A = \begin{bmatrix}
1 & 1 & 0 & 0 & -1 & 0 & 0 \\
0 & 1 & 1 & 0 & 0 & -1 & 0 \\
0 & 0 & 1 & 1 & 0 & 0 & -1
\end{bmatrix}

\\

W = \begin{bmatrix}
\vec{A B} + \vec{B C} - \vec{A C} \\
\vec{B C} + \vec{C D} - \vec{B D} \\
\vec{C D} + \vec{D E} - \vec{C E}
\end{bmatrix}

\\

K = N^{-1} W = (A P^{-1} A^T)^{-1} W \\
V = P^{-1} A^T K

程式

condition_vector_equation = [
1 1 0 0 -1 0 0 
0 1 1 0 0 -1 0 
0 0 1 1 0 0 -1
];
condition_equation = expand_3d_design_matrix(condition_vector_equation);

condition_vector_value = [
  baseline_vector_vector(1,:) + baseline_vector_vector(2,:) - baseline_vector_vector(5,:);
  baseline_vector_vector(2,:) + baseline_vector_vector(3,:) - baseline_vector_vector(6,:);
  baseline_vector_vector(3,:) + baseline_vector_vector(4,:) - baseline_vector_vector(7,:);
];

condition_value = reshape(
                      condition_vector_value',
                      prod(size(condition_vector_value)), 1);

condition_normal = condition_equation * condition_equation';
condition_K = inv(condition_normal) * condition_value;
condition_residue = condition_equation' * condition_K;

condition_variance_post = (
  condition_residue' * condition_residue / rows(condition_equation)
);

結果

條件平差解出來的結果是各觀測量的殘差， 所以需要把殘差加回觀測量， 再把觀測量組合，求出各點的絕對座標。

condition_fix_observation = baseline_vector - condition_residue;
condition_fix_vector = reshape(condition_fix_observation, 3, rows(condition_fix_observation) / 3)';

condition_point = zeros(columns(design_vector_matrix), 3);
condition_point(1,:) = constrain_value';
condition_point(2,:) = condition_point(1,:) + condition_fix_vector(1,:);
condition_point(3,:) = condition_point(1,:) + condition_fix_vector(5,:);
condition_point(4,:) = condition_point(3,:) + condition_fix_vector(3,:);
condition_point(5,:) = condition_point(4,:) + condition_fix_vector(4,:);

結果應該要與間接平差的結果相同。 但已知題目的基線是無誤差的， 在間接平差時只有浮點數誤差，條件平差也有。 但這二者的誤差，又因為計算不同，浮點數誤差也不同。

\begin{bmatrix}
A \\ B \\ C \\ D \\ E
\end{bmatrix}
=
\begin{bmatrix}
1.7760 \times 10^{+05} & 2.6397 \times 10^{+06} & 2.5375 \times 10^{+01} \\
1.8909 \times 10^{+05} & 2.6265 \times 10^{+06} & 4.4821 \times 10^{+01} \\
1.8348 \times 10^{+05} & 2.6205 \times 10^{+06} & 3.7007 \times 10^{+01} \\
1.9699 \times 10^{+05} & 2.6494 \times 10^{+06} & 4.2254 \times 10^{+01} \\
1.8207 \times 10^{+05} & 2.6138 \times 10^{+06} & 3.5997 \times 10^{+01}
\end{bmatrix}

比較

浮點數誤差

條件平差結果應該會與間接觀測平差相同， 但程式計算數值上，二者的結果不同。 單純是因為計算過程不同，造成浮點數誤差不同的結果。

例如以 AB BC AC 三條基線舉例， 這三條基線應該是無閉合差的， 但在 matlab 中計算後， 三條基線卻在 z 方向出現了 1.7764e-15 的誤差。

>> ( [11488.973 -13224.866 19.446] + 
     [-5605.975 -6047.976 -7.814] -
     [5882.998 -19272.842 11.632] )

ans =

   0.0000e+00   0.0000e+00   1.7764e-15

如果在以無閉合差計算，條件平差的 W 應為 0 向量， 連帶連系數矩陣 K 與殘差 V 也都會是 0。 但因為浮點數誤差，如同上面的計算 W 不為 0， 解出來的殘差 V 也就不為 0。

條件式的列式

相比間接平差直接列出所有觀測方程即可， 條件平差需要列出所有獨立的條件式。 我一開始還列不出來，把 ABCDE 連成一圈， 之後就找不到剩下的二條了。 後來用最小閉合圈去想才找到三條獨立的條件式。

而且所謂計算量較小，也是相對間接平差而已。 在計算基線平差時，由於是三維向量的平差， 一條向量式就是三條純量式，要手算也是很辛苦。 只有傳統的三角網，只有角度未知數會少一點， 用條件式可能較好算而已。