平差作業七

平差理論的第七次作業， 計算線性回歸模型，使用全台各處高程取樣， 內差出任意座標的高程值。

線性迴歸內差

給定多組座標與對應高程值與不同座標間高程的共變異數矩陣， 可以用線性迴歸模型內差出任意點的高程。 公式如下：

\hat \gamma = \Sigma _ { \gamma y } ( \Sigma_y + \Sigma_n )^{-1} y

假設距離與高程的共變異數計算如下：

C(d) = 400 exp( \frac {-d^2} {5000^2} )

迴歸模型中 Σy 即為已知點位數的 NxN 矩陣， 其中 Σy(i,j) 為第 i 點與第 j 點間距離 dij 帶入 共變異數函數 C(d) 的結果。 而 Σγy 也類似，只是是 1xN 矩陣， 內容 Σγy(i) 為待求點與第 i 點的距離 di 帶入 C(d) 的結果。

經緯度與距離

給定的 DEM 中座標是使用經緯度， 為了方便計算距離，統一轉換成 twd97 平面座標， 有了平面座標，任二點的距離就能以畢氏定理計算。 使用 proj.4 做轉換過程：

$ proj +init=epsg:3826 <<TWD97
> 121.643707 24.294325
> 121.687935 24.294930
> 121.747860 24.307419
> 121.747704 24.308607
> 121.768066 24.334843
> 121.638969 24.350378
> 121.770851 24.342606
> 121.718277 24.368301
> 121.779121 24.369247
> 121.784531 24.387468
> 121.785973 24.400999
> 121.728401 24.442059
> 121.782791 24.444590
> 121.783516 24.460348
> 121.804192 24.466047
> 121.816399 24.474068
> 121.750992 24.484215
> 121.535110 24.492928
> 121.836288 24.505154
> 121.563942 24.509218
> 121.830872 24.517241
> 121.847576 24.522772
> 121.760681 24.537601
> 121.849899 24.552010
> 121.642738 24.558874
> 121.492306 24.577706
> 121.751938 24.588795
> 121.495644 24.604677
> 121.511536 24.597583
> 121.835426 24.596592
> 121.854218 24.596045
> TWD97

315342.88       2687770.43
319832.29       2687858.90
325908.07       2689273.52
325891.53       2689405.01
327942.29       2692322.12
314833.40       2693976.44
328220.16       2693183.51
322870.44       2696000.86
329042.82       2696138.92
329580.29       2698160.16
329718.07       2699659.68
323854.65       2704175.57
329368.05       2704486.00
329431.68       2706231.78
331524.20       2706875.03
332756.48       2707770.70
326119.96       2708857.01
304234.09       2709720.25
334751.85       2711225.90
307148.95       2711536.12
334194.89       2712561.37
335883.88       2713184.29
327069.50       2714775.47
336099.33       2716424.22
315108.59       2717070.88
299862.27       2719093.94
326152.72       2720440.94
300189.61       2722082.42
301801.80       2721302.58
334603.12       2721353.31
336506.64       2721304.40

共變異數矩陣計算

共變異數矩陣 Σy 是所有已知點的共變異數矩陣， 其中第 i 列 j 行 Σy(i,j) 的值， 為第 i 點與第 j 點間距離 dij 帶入 共變異數函數 C(d) 的結果， 也就是第 i 點與第 j 點的共變異數。

noise 矩陣

如果有高程的標準差，變異數是標準差的平方， 可以計算出 noise 矩陣。 但 DEM 只有給重力異常的標準差， 和高程單位不同，不能相加； 因此沒有使用。

完整計算流程

程式使用 mathjs 計算，語法類似 matlab。

# 上面計算的平面座標，欄位如下
#   [經度,      緯度,      E,         N,          高程,     重力異常標準差]
data = [
    121.643707, 24.294325, 315342.88, 2687770.43, 178.699, 0.04;
    121.687935, 24.294930, 319832.29, 2687858.90, 119.084, 0.05;
    121.747860, 24.307419, 325908.07, 2689273.52, 25.441, 0.03;
    121.747704, 24.308607, 325891.53, 2689405.01, 23.657, 0.04;
    121.768066, 24.334843, 327942.29, 2692322.12, 25.600, 0.03;
    121.638969, 24.350378, 314833.40, 2693976.44, 61.407, 0.02;
    121.770851, 24.342606, 328220.16, 2693183.51, 17.983, 0.01;
    121.718277, 24.368301, 322870.44, 2696000.86, 135.044, 0.04;
    121.779121, 24.369247, 329042.82, 2696138.92, 10.943, 0.05;
    121.784531, 24.387468, 329580.29, 2698160.16, 5.072, 0.06;
    121.785973, 24.400999, 329718.07, 2699659.68, 5.371, 0.04;
    121.728401, 24.442059, 323854.65, 2704175.57, 81.512, 0.04;
    121.782791, 24.444590, 329368.05, 2704486.00, 13.584, 0.06;
    121.783516, 24.460348, 329431.68, 2706231.78, 12.202, 0.07;
    121.804192, 24.466047, 331524.20, 2706875.03, 10.725, 0.05;
    121.816399, 24.474068, 332756.48, 2707770.70, 3.144, 0.06;
    121.750992, 24.484215, 326119.96, 2708857.01, 99.407, 0.09;
    121.535110, 24.492928, 304234.09, 2709720.25, 168.485, 0.05;
    121.836288, 24.505154, 334751.85, 2711225.90, 17.825, 0.04;
    121.563942, 24.509218, 307148.95, 2711536.12, 190.186, 0.05;
    121.830872, 24.517241, 334194.89, 2712561.37, 1.243, 0.06;
    121.847576, 24.522772, 335883.88, 2713184.29, 13.435, 0.05;
    121.760681, 24.537601, 327069.50, 2714775.47, 122.804, 0.06;
    121.849899, 24.552010, 336099.33, 2716424.22, 15.764, 0.07;
    121.642738, 24.558874, 315108.59, 2717070.88, 113.688, 0.08;
    121.492306, 24.577706, 299862.27, 2719093.94, 6.237, 0.07;
    121.751938, 24.588795, 326152.72, 2720440.94, 70.245, 0.06;
    121.495644, 24.604677, 300189.61, 2722082.42, 40.148, 0.06;
    121.511536, 24.597583, 301801.80, 2721302.58, 8.103, 0.08;
    121.835426, 24.596592, 334603.12, 2721353.31, 10.748, 0.07;
    121.854218, 24.596045, 336506.64, 2721304.40, 8.320, 0.09
]

# 平均高程
data_mean_height = mean(data[:,5])

# 第 i 個點的座標
data_point(i) = [data[i,3], data[i,4]]
# 第 i 個點的高程
data_height(i) = data[i,5]
# 第 i 個點的重力異常標準差
data_std(i) = data[i,6]


# 平面二點距離公式
distance(p1, p2) = norm(p1 - p2)

# 根據距離計算共變異數
covariance_distance(d) = 400*exp(-d^2 / 5000^2)
covariance_point(p1,p2) = covariance_distance(distance(p1,p2))

rows(A) = size(A)[1]
columns(A) = size(A)[2]

# 計算一點與所有已知點的 Σγy 矩陣公式
Sigma_ph(p) = map(1:rows(data), cf(i) = covariance_point(p, data_point(i)))

# 所有高程點的共變異數矩陣
Sigma_hh = map(zeros(rows(data), rows(data)),
               cf(x, index) = covariance_point(data_point(index[1]),
                                               data_point(index[2])))

# 如果有給定高程的標準差，就能算出 noise 矩陣
# Sigma_n = diag(map(1:rows(data), data_std))

# 高程殘差
height_residual = data_height(1:rows(data)) - data_mean_height

# 由整體計算特定一點的內差
interpolate_point(p) = (
    Sigma_ph(p) * inv(Sigma_hh) * height_residual + data_mean_height
)

p1 = [326039.96, 2705066.8] # [124.45, 121.75]
p2 = [346165.75, 2727344.37] # [24.65, 121.95]

interpolate_point(p1) # 40.32176363428874
interpolate_point(p2) # 51.80235472239709

結果