Proceedings of the 12th Annual Meeting of Particle Accelerator Society of Japan
August 5-7, 2015, Tsuruga, Japan

PASJ2015 THOM08

相対的な誤差楕円も表示する平面の測量網計算プログラムの開発
ADJUSTMENT PROGRAM OF THE HORIZONTAL SURVEY NETWORK WITHRELATIVE ERROR ELLIPSE
松井佐久夫#, A), 木村洋昭 A), 木内淳 B)
Sakuo Matsui #, A), Hiroaki Kimura A), Jun Kiuchi B)
A)
RIKEN
B)
SPring-8 Service Co.,Ltd.

Abstract
The two-dimensional program for the network survey using Schreiber method is written in FBASIC language. It isuseful for the survey and simulation because also the graphics results is get at the same time. The memory and speed of
PC are improved very much now. Furthermore this program is improved from BASIC language to Excel VBA. Afterthe parameters and survey data are filled in the column of the Excel sheet the calculated results continued in the lower
lines by applying Excel VBA (Visual Basic Application). This program can calculate also the relative error ellipse,
which is much important to the accelerator.

1.

はじめに

2.0

測定個数 角度 1394
距離 1394dr(mm)

1.0
加速器のアライメント関連技術として測量は欠かせないが、土木や建設に比べ加速器特有の部分は少0.0
なく、測量会社に依頼できる。一方、建設前に、最-1.0
2014 SA-USMN
終のアライメントまで必要な精度を得る最適な基準
2014 F-BASIC-2.0
点の位置や数を求めるシミュレーションは欠かせな
0
6
12
18
24
30
36
42
48
いがこの計算は測量会社への依頼は難しい。
Cell number
SPring-8 では 23 年前の蓄積リングの建設時から水Figure 1: Comparison between SA-USNM and FBASIC.
平面の位置計算には測量の教科書にある観測方程式からシュライバーの消去法を用いる最小二乗法の誤差楕円とシミュレーション
BASIC プログラム[1]を使用してきた。実測値と計算 3.
結果からの推定値との差で、入力ミスや測量の不良 3.1 シミュレーション
場所を特定しやすいような部分も加えた。また分散、測定値に似せるために BASIC の一様乱数から下共分散行列を計算するので、誤差楕円[2] 、特に加速記のプログラムで正規乱数に変換している。
器に重要な相対的な誤差楕円、また、既知点 1 点と*TRANS
もう一点は x 又は y の一軸のみ固定する計算や同時' Trans Random Number from flat distribution to Normal oneに図示までできるプログラムを追加した。
RN=RND-.5 : PN=SGN(RN) : IF RN<0 THEN RN=-RN
幸 い 、 富 士 通 の FBASIC は Windows XP や
SUM=0 : XRAN=0 : DXRAN=.01
COMP=RN*SQPI/DXRAN
Windows7 で も 走 る の で そ の 後 の 長 さ 700m のFOR IRAN=1 TO 30000
SACLA や 400m の XSBT の計画段階でのシミュ
SUM=SUM+EXP(-XRAN*XRAN)
レーションや実際の測量計算にも容量やスピードを
IF SUM>COMP THEN RNN=XRAN : GOTO 14210
気にすることなく使用できている。
XRAN=XRAN+DXRAN
レベルは 1 次元でプログラム集[3]のを用いている。 NEXT IRAN
14210 BUNP=PN*RNN*SIGMA/.707107
最近は加速器でレーザートラッカーが多用される
returnようになり、そのコントロールのソフトに測量後の3 次元の計算も含まれるようになってきている。こ 3.2 絶対的な誤差楕円こではこれらの市販のソフトとの比較も試みた。
水平面内での網平均計算には角度と距離測定の測定エラーから推定される重みを取り込んだ最小二乗
2. 市販のソフトとの比較法が用いられる。この重みに計算結果から求まる標
SPring-8 建設時、測量会社の計算にはソフトに保 準偏差を掛けて 2 次元平面で表示したのが誤差楕円証書がついているが同じ測定値を BASIC プログラ である。SACLA 光源棟で角度精度 0.5 秒、距離精度
ムに入力し結果に差がないことを確認した。最近は 10μm+0.4ppm の条件で測量値を決めた値から振らレーザートラッカーのソフト New River Kinematics せ、それを解析し求められた値を 400 セット繰り返
社 の Spatial Analyzer(SA) の USMN(Unified Spatial しプロットしたのが Figure2 である。乱数で振らせMetrology Network)を使用しているので比較し、よく た模擬の測定値の標準偏差と解析時の測定エラーの一致していることがわかった。
値を同じにしたため 2σの赤い楕円内に大部分が___________________________________________
入っているのがわかる。（理論上は 86.5％）# matsui@spring8.or.jp- 140 -

Proceedings of the 12th Annual Meeting of Particle Accelerator Society of Japan
August 5-7, 2015, Tsuruga, Japan

PASJ2015 THOM08

Figure 2: Simulation results of the SACLA undulator
hall.
3.3

相対的な誤差楕円レーションで達成できる精度を検討する際に測量の実力値を正確に評価しておくことは欠かせない。
Figure 4 は SACLA 光源棟でレーザートラッカーで測量したデータを角度 1 秒、距離は 10μm+距離×
2.5μm/m の条件で計算した時の、測定距離対計算値と測定値の差を青でプロットしたものである。赤はシミュレーションで角度 0.5 秒、距離は 10μm+距離×0.4μm/m の条件で測定値を振らせ、同条件で解析したものである。距離依存性は 10μm+距離×2.5μm/m の直線は立ち上がりも早く、実際とかけ離れており、例えば 0.4μm/m 程度が妥当なことを示している。解析の条件によりこの分布は変化するのでシ
ミュレーションのパラメーターをふり実際のエラーに近い値を見つけることが可能と考えられる。

観測方程式と測定器のエラーから求められる重みの行列から解を求めるために逆行列を計算する。この時求まる分散共分散行列の対角成分は各点の x ま
たはｙ座標の分散のみなので 2 次元で傾きを持つ絶対的な誤差楕円は非対角成分も用いて計算される。
加速器では軌道の滑らかさが重要なので相対的な誤差楕円の方がより重要になる。教科書[4]にある下記の式で求めることができる。プログラムでは簡単のために隣り合う 2 点からのみ計算しているが離れていても共分散の項から求めることができる。
Figure 4: Comparison between measured and simulation.
4.2

Figure 3: Covariance matrix and relative error ellipse.

4.

シミュレーションの応用

4.1

測定エラーの評価とシミュレーション

網平均計算には角度と距離の測定エラーの見積が必要になる。測定器は例えば角度 1 秒、距離は 15μm+距離×6μm/m などメーカーから与えられる。
しかし実際には空気の揺らぎ、ターゲットそのものや台の不安定さ、など他の要因も加わる。シミュ

蓄積リングへの適用

Figure 5 の最大距離 38ｍの測量網で、蓄積リングの測量環境を整えた場合のライカ製レーザートラッ
カーAT402 の実力に近いと考えられる測定エラー(角度 0.5 秒、距離 10μm+0.4ppm)、で反対側は一軸のみ固定でシミュレーションした。固定点と反対側の誤差楕円の半径は Figure6 のように 1σ0.5mm であった。この場合周長のエラーは 0.5mm×2πにはならず、20 回のシミュレーションでは 1436ｍに対して平均値は-0.03mm で分布の rms 値は 0.07mm であった。Figure5 を棒で組みあげられた網のイメージで考えると Figure6 で固定点の反対側で例えば直線の中央から左にずれた場合図の上側の斜めの楕円は中心から左下の方向にずれ結果として周長はほとん
ど変化せずシステマティックなエラーだけが問題になると推測される。この時の相対的な誤差楕円は
Figure7 のように 1σで 20μm 程度で収まっている。
狭い範囲でも密な網と測量機器の精度の結果である。

Figure 5: Control network. There are two tracker positions in one cell.

- 141 -

Proceedings of the 12th Annual Meeting of Particle Accelerator Society of Japan
August 5-7, 2015, Tsuruga, Japan

PASJ2015 THOM08

入力データの区切り：Figure 8(a)のように数字 9999と上の行と同様にセルを 0 で埋める。各データ間に空白行は入れない。入力が済めば VBA を走らせる。
解析時：時刻と実行しているコメントを出力する。
結果出力部：計算した日時、計算に用いたパラメー
ター、既知点、（半固定点）未知点の座標が出力される。また、角度や距離の測定数、各繰り返しでの標準偏差とその減衰の比、誤差楕円の長軸、短軸の長さと回転角、さらに計算結果から推測される

Figure 6: Absolute error ellipse at SPring-8.

Fig. 8(a)Figure 7: Absolute and relative error ellipse.

5.

Excel VBA プログラム

5.1

Excel シート

Figure8(a),(b),(c)に示すように Excel の 1 枚のシートに必要な数値を入力し、BASIC を Excel VBA に書き直した VBA のプログラムを走らせることで同じシートに計算結果と誤差楕円が出力される。
解析パラメーター：角度と距離の精度、繰返し回数。
既知点：２次元 XY で座標を与える。測量点は番号でなくても名前をそのまま使える。BASIC と異なりVBA では大文字小文字の区別がされる。固定点（既知点）は多くても構わないがこのプログラムでは 1
点だけでは計算できない。
半固定点： 1 点は固定で、もう一点までの正確な距離がわからない場合、その点の X 又は Y 座標だけを観測方程式の行列の対応する”列”だけ０にする
ことで固定できるように改修した。
未知点： XY の近似値を与える。この値が悪いと繰り返しが多くなるが、その収束具合は標準偏差の変化（出力ファイルの Reduction ratio）でわかる。
測定データ部：角度も距離も同時に得られる場合が多いがこのプログラムではまず角度だけ、次に距離だけと別々に入力する。

- 142 -

Fig. 8(b)Proceedings of the 12th Annual Meeting of Particle Accelerator Society of Japan
August 5-7, 2015, Tsuruga, Japan

PASJ2015 THOM08

Fig. 8 (C)Figure 8: Excel worksheet for the survey calculation.
(a) Analysis parameters and measured data.
(b) Calculated Results.
(c) The values and figures of error ellipses.
角度や距離と実際の測定値との差を出力する。この値が大きいところでは、入力のミスや測量が良くな
い、など疑われるのでチェックする。角度だけ、距離だけの平均値も出力される。この２つの平均値は解析する時の角度と距離の精度の設定バランスにも依存する。角度精度を小さくすれば距離の fitting の値は大きくなり逆に距離の精度を小さくすれば角度の fitting の値は大きくなる。
5.2

VBA プログラム

宣言部はわずかだが、データの読み込みと絶対的と相対的誤差楕円の計算まで省かずに以下載せる。
誤差楕円のグラフィックは結果のみ Figure 8(c)に示した。
ガウスの掃き出し法で連立一次方程式を解く時逆行列を計算するが、行数が多くなるため桁落ち防止のため軸の数値がそれより下の行の中で最大値に交換する Pivotting の改造をした。
Sub Macro1()
Dim i As Integer, i0 As Integer, i1 As Integer, ii As Integer
Dim at As Double, b2 As Double, b4 As Double
mp = 500: mn = 2900: mj = mp * 2:
Dim AL(2900) As Double, a(2900, 2900 / 2) As Double,
Dim QK(2900, 3) As String, QP(500) As String,pi = 3.1415926535: sqpi = Sqr(pi): ra = pi / 180#: radsec = 3600# / ra
' ------- Weight Parameters , Iteration times, fixed point--Range("D5").Selectrmt = ActiveCell.Value / radsec: Selection.Offset(1, 0).Selectu1 = ActiveCell.Value * 0.001: Selection.Offset(1, 0).Selectr1 = ActiveCell.Value * 0.000001
jfix = Range("D9"): maxiter = Range("D12")
Range("B16").Select' -------Read X , Y Data------------j=1
Do Until ActiveCell.Value = "9999"
QP(j) = ActiveCell.Value: Selection.Offset(0, 1).SelectXJ(j, 1) = ActiveCell.Value: Selection.Offset(0, 1).SelectYJ(j, 1) = ActiveCell.Value: Selection.Offset(1, -2).Select- 143 -

j=j+1
Loop
nr = j - 1: ' nr=number of known points
For ii = 1 To nr
For iter = 2 To maxiter: XJ(ii, iter) = XJ(ii, 1): YJ(ii,iter) = YJ(ii, 1): Next iter
Next ii
Selection.Offset(1, 0).Select' ---------------- Assumed X Y data ------------------Do Until ActiveCell.Value = "9999"
QP(j) = ActiveCell.Value: Selection.Offset(0, 1).SelectXJ(j, 1) = ActiveCell.Value: Selection.Offset(0, 1).SelectYJ(j, 1) = ActiveCell.Value: Selection.Offset(1, -2).Select
j=j+1
Loop
Selection.Offset(1, 0).Selectnp = j - 1: nq = np - nr: 'np=number of total points nq= unknown points
For iter = 1 To maxiter
For i = 1 To np
For j = 1 To np: IAN(i, j) = 0: Next j
Next i'----- Initialize Matrix A(j,k) , AK(i)-------n = nq * 2: ' N=twice of unknown points
For i = 1 To mn
For j = 1 To mn / 2: a(i, j) = 0#: Next j
P(i) = 0#: AL(i) = 0#: AK(i) = 0#
Next i' -------Read angle data----------i = 0: l = 0: sigma = rmt
Do
l=l+1
If iter = 1 Then
QK(l, 1) = ActiveCell.Value: Selection.Offset(0, 1).SelectQK(l, 2) = ActiveCell.Value: Selection.Offset(0, 1).SelectQK(l, 3) = ActiveCell.Value: Selection.Offset(0, 1).SelectANG(l) = ActiveCell.Value: Selection.Offset(1, -3).SelectIf ANG(l) < 0 Then ANG(l) = ANG(l) + 360#
End IfIf QK(l, 1) = "9999" Thennang = l - 1
Else
For j = 1 To np
If QK(l, 1) = QP(j) Then j1 = j
If QK(l, 2) = QP(j) Then j0 = j
If QK(l, 3) = QP(j) Then j2 = j
Next j'-- Calculation of Coefficients of observation eq.
c1 = XJ(j0, iter) - XJ(j1, iter): d1 = YJ(j0, iter) - YJ(j1, iter)
If (c1 = 0) * (d1 = 0) ThenElse: e1 = d1 / (c1 ^ 2 + d1 ^ 2): f1 = -c1 / (c1 ^ 2 + d1 ^ 2)
End Ifc2 = XJ(j0, iter) - XJ(j2, iter): d2 = YJ(j0, iter) - YJ(j2, iter)
If (c2 = 0) * (d2 = 0) ThenElse: e2 = d2 / (c2 ^ 2 + d2 ^ 2): f2 = -c2 / (c2 ^ 2 + d2 ^ 2)
End Ift1 = arctan(-c1, -d1)
t2 = arctan(-c2, -d2)
at = ANG(l) * raIf t2 < t1 Then t2 = t2 + 2 * pi
ANGL(l, iter) = (t2 - t1) / ra: BUN(l) = ANG(l) - ANGL(l, iter)b2 = (t2 - t1) - at
If IAN(j1, j0) = 1 Then GoTo 10820: ' J1=reference angle

Proceedings of the 12th Annual Meeting of Particle Accelerator Society of Japan
August 5-7, 2015, Tsuruga, Japan

PASJ2015 THOM08

point
If i = 0 Then GoTo 10710: ' J0=survey instrument position
End If
For j = 1 To np * 2: a(i, j) = CS(j): Next j
AL(i) = cl: P(i) = -1# / i1
If QK$(l, 1) = "9999" Then nt = i: Exit Do: ' NT=Angles number
10710: For j = 1 To np * 2: CS(j) = 0#: Next j
cl = 0#: i1 = 0: IAN(j1, j0) = 1
i=i+1a(i, j1 * 2 - 1) = e1: a(i, j1 * 2) = f1
a(i, j0 * 2 - 1) = -e1: a(i, j0 * 2) = -f1
AL(i) = 0#: P(i) = 1#
For j = 1 To 2 * np: CS(j) = CS(j) + a(i, j): Next j
cl = cl + AL(i): i1 = i1 + 1i=i+110820: a(i, j2 * 2 - 1) = e2: a(i, j2 * 2) = f2
a(i, j0 * 2 - 1) = -e2: a(i, j0 * 2) = -f2
AL(i) = b2: P(i) = 1#
For j = 1 To np * 2: CS(j) = CS(j) + a(i, j): Next j
cl = cl + AL(i): i1 = i1 + 1
i=i+1
Loop
'------- Read Distance data -----i = 1: l = 0
Do
l=l+1
If iter = 1 Then
QS(l, 1) = ActiveCell.Value: Selection.Offset(0, 1).SelectQS(l, 2) = ActiveCell.Value: Selection.Offset(0, 1).SelectSS(l) = ActiveCell.Value: Selection.Offset(1, -2).SelectEnd IfIf QS(l, 1) = "9999" Then ns = i - 1: ndis = l - 1: Exit Do: 'NS=Distance number
i0 = i + nt
For j = 1 To np
If QS(l, 1) = QP(j) Then j1 = j
If QS(l, 2) = QP(j) Then j2 = j
Next jc1 = XJ(j2, iter) - XJ(j1, iter): d1 = YJ(j2, iter) - YJ(j1, iter)s2 = c1 ^ 2 + d1 ^ 2: s1 = Sqr(s2)e1 = c1 / s2: f1 = d1 / s2a(i0, 2 * j1 - 1) = -e1: a(i0, 2 * j1) = -f1
a(i0, 2 * j2 - 1) = e1: a(i0, 2 * j2) = f1
b4 = (s1 - SS(l)) / s1: SSR(l, iter) = s1: SBUN(l) = SS(l) SSR(l, iter)' ----- Weight ( Distance ) -----------pd = rmt ^ 2 * s2 / (u1 ^ 2 + r1 ^ 2 * s2)
AL(i0) = b4: P(i0) = pd
i=i+1
Loop
'--- components arrange of A
If iter = maxiter Then Exit For
For i = 1 To nt + ns
For j = 1 To nq * 2: sss = a(i, j + 2 * nr): a(i, j) = sss: Next j
Next iIf (jfix = 1) Or (jfix = 2) Then
For i = 1 To mn: a(i, jfix) = 0#: Next i: ' fix X or Y coordinate
End If'--------- Normal Equations A ----------n = nq * 2: m = nt + ns
For i = 1 To n
For j = 1 To n: S(i, j) = 0#: Next j: AK(i) = 0#
Next i

For j = 1 To n
For k = 1 To n: s3 = 0#
For i = 1 To m: s3 = s3 + P(i) * a(i, j) * a(i, k): Next i
S(j, k) = s3
Next kNext j
For j = 1 To n
For i = 1 To m: AK(j) = AK(j) + P(i) * a(i, j) * AL(i): Next i
Next j'---------------------Inverse Matrix ------------For k = 1 To n
If k <> jfix Then
smax1 = S(k, k)
For i = k To n: If Abs(S(k, i)) >= smax1 Then smax1 = Abs(S(k, i)): ich = i
Next i
For j = 1 To ntemp(j) = S(k, j): S(k, j) = S(ich, j): S(ich, j) = temp(j)Next jtempak = AK(k): AK(k) = AK(ich): AK(ich) = tempak
t = 1# / S(k, k): S(k, k) = 1#
For j = 1 To n: S(k, j) = S(k, j) * t: Next j
For i = 1 To n
If i <> k Thenq = S(i, k): S(i, k) = 0#
For j = 1 To n: S(i, j) = S(i, j) - q * S(k, j): Next j
End IfNext iEnd IfNext k'-----------------Calculate X------------For i = 1 To n: XX(i, iter) = 0#
For j = 1 To n: XX(i, iter) = XX(i, iter) + S(i, j) * AK(j): Next j
Next ij=1
For i = nr + 1 To np: XJ(i, iter + 1) = XJ(i, iter) - XX(j,iter): YJ(i, iter + 1) = YJ(i, iter) - XX(j + 1, iter)
j=j+2
Next i' -------Standard Deviation ------------------For i = 1 To m: V(i) = 0
For j = 1 To n: V(i) = V(i) + a(i, j) * XX(j, iter): Next j
V(i) = V(i) + AL(i): 'AL(I) is lack in the textbook program
Next ivv = 0#For i = 1 To m: vv = vv + P(i) * V(i) ^ 2: Next iSD(iter) = Sqr(vv / (m - n))sdifang = 0#: sdifdis = 0#
For i = 1 To nang: sdifang = sdifang + (ANGL(i, iter) ANG(i)) ^ 2: Next i: Angfit(iter) = Sqr(sdifang / nang) * 3600
For i = 1 To ndis: sdifdis = sdifdis + (SSR(i, iter) - SS(i)) ^2: Next i: Disfit(iter) = Sqr(sdifdis / ndis) * 1000#
Next iter'Print "ITER= "; iter
For i = 2 To maxiter: Reduction(i) = Sqr(SD(i) / SD(i - 1)): Next i' ----------- Error ellipse (Absolute)----------------For i = 1 To nq
SN_abs(i) = SD(maxiter - 1) * Sqr(Abs(S(i * 2 - 1, i * 2 1))) * 1000: EW_abs(i) = SD(maxiter - 1) * Sqr(Abs(S(i * 2, i *2))) * 1000
SNEW_abs(i) = SD(maxiter - 1) * Sqr(Abs(S(i * 2 - 1, i *2))) * 1000 * Sgn(S(i * 2 - 1, i * 2))
Angle_abs(i) = 0.5 / ra * Atn(2 * SNEW_abs(i) ^ 2 /- 144 -

Proceedings of the 12th Annual Meeting of Particle Accelerator Society of Japan
August 5-7, 2015, Tsuruga, Japan

PASJ2015 THOM08(SN_abs(i) ^ 2 - EW_abs(i) ^ 2)) * Sgn(S(i * 2 - 1, i * 2))
SMAX_abs(i) = Sqr(SN_abs(i) ^ 2 + EW_abs(i) ^ 2 +Sqr((SN_abs(i) ^ 2 - EW_abs(i) ^ 2) ^ 2 + 4 * SNEW_abs(i) ^4)) / Sqr(2)SMIN_abs(i) = Sqr(Abs(SN_abs(i) ^ 2 + EW_abs(i) ^ 2 Sqr((SN_abs(i) ^ 2 - EW_abs(i) ^ 2) ^ 2 + 4 * SNEW_abs(i) ^4))) / Sqr(2)Next i' -------------Error ellipse (Relative) -----------------For i = 1 To nq - 1
SN_rela(i) = SD(maxiter - 1) * Sqr(Abs(S(i * 2 - 1, i * 2 1) + S(i * 2 + 1, i * 2 + 1) - 2 * S(i * 2 - 1, i * 2 + 1))) * 1000
EW_rela(i) = SD(maxiter - 1) * Sqr(Abs(S(i * 2, i * 2) +S(i * 2 + 2, i * 2 + 2) - 2 * S(i * 2, i * 2 + 2))) * 1000
SNEW_rela(i) = SD(maxiter - 1) * Sqr(Abs(S(i * 2 - 1, i *2) + S(i * 2 + 1, i * 2 + 2) - S(i * 2, i * 2 + 1) - S(i * 2 - 1, i * 2 +2))) * 1000 * Sgn(S(i * 2 - 1, i * 2) + S(i * 2 + 1, i * 2 + 2) - S(i* 2, i * 2 + 1) - S(i * 2 - 1, i * 2 + 2))
Angle_rela(i) = 0.5 / ra * Atn(2 * SNEW_rela(i) ^ 2 /(SN_rela(i) ^ 2 - EW_rela(i) ^ 2)) * Sgn(S(i * 2 - 1, i * 2) + S(i *2 + 1, i * 2 + 2) - S(i * 2, i * 2 + 1) - S(i * 2 - 1, i * 2 + 2)) +(Sgn(EW_rela(i) ^ 2 - SN_rela(i) ^ 2) + 1) * 45 * ra
SMAX_rela(i) = Sqr(SN_rela(i) ^ 2 + EW_rela(i) ^ 2 +Sqr((SN_rela(i) ^ 2 - EW_rela(i) ^ 2) ^ 2 + 4 * SNEW_rela(i) ^4)) / Sqr(2)SMIN_rela(i) = Sqr(Abs(SN_rela(i) ^ 2 + EW_rela(i) ^ 2 Sqr((SN_rela(i) ^ 2 - EW_rela(i) ^ 2) ^ 2 + 4 * SNEW_rela(i) ^4))) / Sqr(2)Next i
Function arctan(x As Double, y As Double) As Double' arc tangent
Dim pi As Double
pi = 3.141592653
arctan = Atn(y / x)
If x < 0 Then
arctan = arctan + pi
ElseIf y < 0 Then
arctan = arctan + pi * 2
End If
End Function

6.

まとめ

6.1

BASIC 系ソフト （Excel VBA 含む）•

BASIC のプログラムから Excel VBA への変換
は同じ文法が使える等スムーズに実現できた。
Excel （VBA）で解析でき使い勝手が向上した。
シミュレーションや相対的な誤差楕円の計算など目的や設定が決まっている場合は測量器に付属しているソフトに比べ使いやすい。
VBA でも計算と同時に誤差楕円も図示できる。
今回の計算時間は FBASIC より VBA の方が速かった。
外部からダウンロードできるところに VBA のソフトを置き、使用者が使いやすい様に改良していってもらえればと考えている。（それまでは必要ならメールで問い合わせてください。）

•
•
•
•
•

6.2

ライカ製 AT402 の実力を十分発揮できれば、蓄積リングの架台間の半径方向のエラーは 20μm(1σ)程度で収まり、また周長のエラーは極わずかでシステ
マティックエラーのみによることがわかった。
6.3

参考文献
[1] 細野武庸他, “測量叢書 第 1 巻 改訂版 基準点測量”, 日本測量協会 (1992).
[2] 原田健久, “わかりやすい測量厳密計算法”, 鹿島出版会(1992).
[3] 石川甲子男他, “パーソナルコンピュータによる測量計
算プログラム”, 山海堂 (1991).
[4] TJM Kennie and G Petrie, “Engineering Surveying
Technology”, p163.

Table 1: Comparison of Calculation Time between
Application FBASIC and That of Excel VBA
角度
111
1300

距離
133
1400

FBASIC<1 秒
26 分(36 分)

測定値のエラーの評価

空気などの環境やターゲット、測量器の実力から測定値のエラーの見積もりはシミュレーションと実測値の比較からある程度評価が可能だと考えられる。

5.3 計算時間の比較
FBASIC のコンパイラーと VBA で同様のプログラ
ム を Windows7 の ノ ー ト PC （ Core (TM)i3 CPU
1.33GH ｚ ） で の 1 回 計 算 時 間 （ () 内 の 時 間 はWindows XP Pentium(R)D CPU 3.4GHz 使用）は Table
1 の通り VBA の方がかなり速いことがわかった。

場所
測量点
光源棟 50
蓄積リ 486
ング

シミュレーション結果

VBA
<1 秒
10 分
(19 分)

測量点数が 500、測量の数が角度、距離含め 3000程度なら観測方程式の行列の要素数は 3000 行×1000 列=300 万にのぼるが 10 分程度で計算できるのでもう少し多くても実用的な範囲内と考えられる。

- 145 -