現場の任意座標を公共座標に直すには? — 既知点で合わせる2次元ヘルマート変換の手順と、残差の読み方

,

現場の任意座標を公共座標に直すには? — 既知点で合わせる2次元ヘルマート変換の手順と、残差の読み方

現場の任意座標を公共座標に直すには? — 既知点で合わせる2次元ヘルマート変換の手順と、残差の読み方

両方の座標を持つ既知点が2点以上あれば直せます。 平面上の「回して・伸縮して・平行移動する」4つのパラメータを既知点から求め、それを全点に適用する——これが2次元ヘルマート変換(平面の相似変換)です。ただし2点で解くと残差がゼロになり、合っているかどうかが分からない。3点以上で解き、残差を読むところまでが手順です。

任意座標とは何だったか

現場で据えたトータルステーションを原点にして、適当な方向を軸にとって観測すると、その現場のなかだけで完結する座標ができます。距離と角度は正しいのに、実世界のどこにあるかという情報だけが抜けている状態です。登記所備付地図データの「任意座標系」も、公図に基づく図面が実世界に対応していないという同じ構図です(地図XMLの公共座標・任意座標の見分け方)。

足りないのは次の4つだけです。

不明なもの パラメータ 意味
どこにあるか ΔX, ΔY 平行移動(2つ)
どちらを向いているか θ 回転(1つ)
長さの尺度は合っているか s 縮尺(1つ)

形は変わらず、向き・大きさ・位置だけが違う——だから相似変換です。角度も、点どうしの距離の比も保たれます。

変換の式

任意座標 (x, y) から公共座標 (X, Y) へは、こう書けます。平面直角座標系はX軸が北、Y軸が東である点に注意してください(XYの軸の向きの落とし穴)。

X = a·x − b·y + ΔX
Y = b·x + a·y + ΔY

  ただし  a = s·cos θ ,  b = s·sin θ
          s = √(a² + b²)  (縮尺)
          θ = atan2(b, a)  (回転角)

未知数は a, b, ΔX, ΔY の4つ。既知点1点につき式が2本(XとY)立つので、2点あれば4本=ちょうど解けるという関係になります。

変換先として選ぶ系番号は、現場の位置で決まります。

平面直角座標系1〜19系の全国区域図

3次元のヘルマート変換とは別物

「ヘルマート変換」という名前は、測地系を変換する話(3パラメータと7パラメータの違い)でも出てきます。あちらは緯度経度をXYZ直交座標に直して3次元で回す話で、未知数は3個か7個。この記事の話は平面上の4パラメータで、目的も使う場面もまったく違います。同じ人名が付いているだけで、混同すると手順を間違えます。

実例で解く

既知点4点を測ったとします。左が現場の任意座標、右がその点の公共座標(平面直角座標系)の成果値です。

点 任意 x 任意 y 公共 X 公共 Y
T1 100.000 100.000 −8513.677 41351.058
T2 350.250 120.400 −8327.593 41519.583
T3 120.800 480.600 −8728.840 41665.647
T4 400.150 455.300 −8491.863 41815.601

未知数4に対して式は8本。余った4本ぶんを最小二乗法で処理すると、次が得られます。

回転   θ = 37.50081°
縮尺   s = 0.9998629
平行移動 ΔX = −8532.134 m
       ΔY =  41210.878 m

そして各点の残差(このパラメータで計算した値と、成果値との差)はこうなります。

点 vX vY |v|
T1 −2.4 mm +12.3 mm 12.6 mm
T2 +3.6 mm −6.2 mm 7.2 mm
T3 −7.2 mm −10.3 mm 12.6 mm
T4 +6.1 mm +4.1 mm 7.4 mm

自由度は 8 − 4 = 4、単位重量の標準偏差は約 10.3 mm。この現場では「1cm前後の食い違いを含んで合っている」と読めます。

残差が教えてくれること・教えてくれないこと

2点では検証できない

未知数4に式4本なので、2点で解くと残差は必ず全部ゼロになります。ぴったり合ったように見えますが、それは合わせたのではなく合わせる余地がなかっただけです。上の例で、T1とT2だけを使って解いたパラメータで残りの点を変換すると、こうなります。

点 ΔX ΔY
T3 −32.9 mm −28.1 mm
T4 −22.4 mm +8.6 mm

残差ゼロだったのに、他の点では3cm食い違う。 2点しか使えないときは、少なくとも「この結果は検証されていない」という認識を持って先へ進む必要があります。

残差は犯人を指さない

もし T3 の公共座標のY値を 12cm 打ち間違えていたら、残差はこうなります。

点 |v|
T1 45.1 mm
T2 5.0 mm
T3 69.8 mm
T4 55.3 mm

σ₀ は 10.3mm から 50.0mm へ跳ね上がり、異常は一目で分かります。しかし誤りは T3 だけなのに、残差は無関係な T1・T4 にも回ってしまう。最小二乗法は誤差を全体へ均すので、残差の最大値が犯人とは限らないのです。σ₀ が急に大きくなったら、1点ずつ外して解き直し、どれを外したときに σ₀ が急落するかを見るのが確実な探し方になります。

縮尺 s は「1」ではない

ここが実務で一番つまずくところです。トータルステーションで測った任意座標は地上の実長ですが、公共座標は投影面上の長さです。両者は一致しません。

要素 内容
投影の縮尺係数 平面直角座標系は原点で 0.9999。原点子午線から離れるほど大きくなる
標高補正 標高 H の地上距離を楕円体面へ落とすと (1 − H/R) 倍になる

この2つを掛けた合成縮尺は、現場の場所と標高で決まります。

原点子午線からの距離 標高 合成縮尺 1kmあたり
0 km 0 m 0.99990000 −100 mm
30 km 100 m 0.99989539 −105 mm
60 km 300 m 0.99989726 −103 mm
90 km 500 m 0.99992131 −79 mm
120 km 1000 m 0.99992043 −80 mm

1kmあたり8〜10cm。これを無視して s = 1 と決め打つと、基点からの距離に比例した系統誤差が出ます。

基点からの距離 出てしまう誤差
100 m 10 mm
200 m 20 mm
400 m 40 mm
1000 m 100 mm

実際、先ほどの4点で縮尺を1に固定して(回転と平行移動だけで)解き直すと、残差は 20.3 / 33.2 / 36.5 / 36.6 mm へ悪化します。4パラメータ版の最大 12.6mm と比べて3倍近い。残差が距離に比例して外側ほど大きくなっていたら、まず縮尺を疑うのが定石です。縮尺係数そのものの話は平面直角座標系の縮尺係数と距離にまとめています。

ただし、小さい現場では縮尺は決めにくい

逆の注意もあります。上の例で推定された s = 0.9998629 は、理論上の合成縮尺(0.99990前後)と 3〜4×10⁻⁵ ずれています。これは既知点の座標が1cm前後ばらついているためで、現場の広がり400m × 1cmの誤差 = 2.5×10⁻⁵ という不確かさが縮尺にそのまま乗る計算です。

つまり数百m規模の現場では、推定した縮尺と理論値を区別できません。このときは s を自由に推定させるより、その現場の合成縮尺を計算して固定してしまうほうが素直です。s を推定させると、既知点の誤差を縮尺が吸収してしまい、現場の外側へ外挿したときに効いてきます。

手順のまとめ

  1. 両方の座標を持つ既知点を3点以上そろえる(4点以上が望ましい)。2点では検証できない
  2. 既知点は現場の外周に散らす。片側に寄せると、反対側が外挿になって精度が落ちる
  3. 4パラメータの最小二乗解を求め、σ₀ と各点の残差を見る
  4. σ₀ が想定より大きければ、1点ずつ外して解き直し、σ₀ が急落する点を探す
  5. 残差が外側ほど大きい傾向なら縮尺を疑う。現場の合成縮尺を計算して固定する判断もある
  6. 変換後、変換に使わなかった既知点でチェックする(これができるように1点残しておくとよい)
  7. どの既知点で・どのパラメータで合わせたかを記録に残す。後から成果を説明できるかどうかは、この記録で決まる

逆に、できないこと

相似変換は形を変えません。もし現場の任意座標そのものに歪み(測距の系統誤差、閉合差の無理な配分など)があると、その歪みは残差として残るだけで、吸収されません。残差に規則的な傾向(片側だけ外へ膨らむ、ねじれているなど)が見えるときは、パラメータを増やす(アフィン変換にする)前に、元の観測を疑うのが順序です。パラメータを増やせば残差は必ず小さくなりますが、それは合ったのではなく、誤差を数式に吸わせただけになりがちです。

出典

  • 作業規程の準則 — 国土地理院: https://www.gsi.go.jp/be/jyunsoku.html
  • 平面直角座標系(平成十四年国土交通省告示第九号) — 国土地理院: https://www.gsi.go.jp/LAW/heimencho.html
  • 測量計算(座標変換・距離と方位角の計算) — 国土地理院: https://vldb.gsi.go.jp/sokuchi/surveycalc/main.html
  • 登記所備付地図データ — 法務省/G空間情報センター: https://www.geospatial.jp/ckan/dataset/legalaffairsbureau-fudemap

次に確認する


開発者より: アプリ・Kindle本・公開プロジェクトの一覧は GitHub: amru195704 にまとめています。


お願い
本記事の情報は参考目的で掲載しており、正確性・完全性を保証するものではありません。誤記・不正確な情報がございましたら、コメント欄よりご指摘いただければ、確認のうえ修正いたします。


コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です

トップへ戻る