1/1,000수치지형도 좌표계변환

표준 작업지침(Ver 1.0)

2005.  10.

			건 설 교 통 부
		국토지리정보원

목  차

I. 지침의 개요	1

1. 제정배경 및 목적	1

가. 제정배경	1

나. 제정목적	1

2. 지침의 제정근거 및 내용 등	2

가. 지침의 제정 근거	2

나. 지침의 주요내용	2

다. 관련 기준 등	3

3. 기대 효과	4

II. 1/1,000 수치지형도 좌표변환 지침	5

1. 총칙	5

가. 지침의 목적	5

나. 용어의 정의	5

다. 적용범위	6

2. 표준 작업절차	7

가. 세계측지계로의 변환 개념	7

나. 작업절차	7

3. 작업 단계별 적용기준	8

가. 준비 단계	8

나. 공통점의 확보	10

다. 세계측지계로 기본 변환	13

라. 왜곡량의 Modeling 및 보정	14

마. 변환된 1/1,000수치지형도의 정확도 평가	17

III. 변환 S/W의 활용 및 비용산정	20

1. 변환 S/W의 활용	20

가. 좌표변환 S/W(1000GDKtrans)	20

나. 왜곡량 계산 S/W개발 및 활용	25

2. 변환비용의 산정	25

가. 적용원칙	25

나. 비용산정 기준	26

<별첨 1> 세계측지계로의 전환 개념도	29

<별첨 2> 세계측지계로의 변환작업 흐름도	30

<별첨 3> 세계측지계로의 기본 변환 기준수식	31

<별첨 4> 왜곡량의 계산 및 모델링 작업 예시	34

<별첨 5> 세계측지계 변환 비용산정 예시	62

<첨부 1> 세계측지계 변환에 필요한 투영식	64

지침의 개요

제정배경 및 목적

제정배경

측량법(□01.12.19 개정)이 개정됨에 따라 2003년 1월 1일부터는 모든 측량의 기준이 세계측지기준계를 이용하여야 함

다만 측량법(측량에 관한 경과조치)에 의한 고시(국토지리정보원 고시 2002-304호, □02.12.31)에 따라 기존의 지역측지계와 세계측지계를 성과를 2006년 12월 31일까지 병행하여 사용 가능

세계측지계는 종전의 지역측지계에 비하여 공간데이터의 획득, 유지, 관리 및 활용이 용이하며, 높은 정확도 확보가 가능

2007년 1월 1일부터는 1/1,000수치지형도를 비롯한 지도 등 모든 측량성과는 세계측지계를 기준으로 사용하여야 함

지역측지계로 제작된 1/1,000수치지형도 등의 성과는 세계측지계를 기준한 성과로 변환하거나 또는 새로이 제작하여 사용하여야 함

따라서 기존의 지역측지계를 기준하여 제작된 1/1,000 수치지형도의 효율적인 변환을 수행할 수 있는 지침이 필요

측지기준의 전환에 따른 혼란과 오차를 최소화하기 위해서 세계측지계로의 변환문제에 대한 통일되고, 표준화된 전략이 요구됨

제정목적

측량법 개정(□01.12.19)으로 측량의 기준이 세계측지기준계로 변경됨에 따라, 1/1,000 수치지형도를 세계좌표계로의 전환 필요

1/1,000수치지형도는 지방자치단체의 각종 UIS업무 등에 기본이 되는 자료로 변환에 따른 혼란을 방지하기 위해서는 지자체 지원 필요

지방자치단체에서 세계측지계로의 전환이 효율적으로 추진될 수 있도록 지원하고, 세계측지계로의 좌표변환 작업이 일관성 있고 정확도가 확보되도록 전환에 필요한 통일된 기준 필요

변환에 필요한 기술적 기본방향 및 기준이 없이 개별적으로 추진될 경우 1/1,000수치지형도의 정확도가 저하될 우려가 있음

변환에 필요한 연구의 중복 실시 및 변환기술의 수준 차이로 인한 일정 수준의 정확도가 확보되지 않아 활용이 저하될 우려가 있음

수치지형도의 효율적이고, 일관성 있는 변환을 위하여 표준적인 단계별 변환 작업지침을 제시함으로써 체계적 변환 도모

1/1,000수치지형도를 세계측지계로 변환하기 위한 표준 작업지침을 제시함으로써 변환에 따른 혼란방지 및 효율성을 증대

지침의 제정근거 및 내용 등

지침의 제정 근거

측량법 제5조(측량의 기준) 및 제6조(측량기술의 연구개발등)

측량법시행령 제2조의5(세계측지계 등)

측량법시행령 제3조제3항

건설교통부장관은 측량에 관한 작업방법 및 기준 등에 관하여 필요한 사항을 정할 수 있다

지침의 주요내용

지침의 제정목적 및 적용기준 등 지침의 개요

단계별 세계측지계 변환 방법

변환 작업검토 및 공통점의 확보

세계측지계로의 기본 변환 및 왜곡량 보정

변환된 1/1,000수치지형도의 정확도 평가

소요비용 등 세계측지계 변환 관련 기준

변환 S/W의 활용 매뉴얼

변환비용의 산정 기준

기타 관련 참고자료

세계측지계로의 전환 개념 및 작업 흐름도

기본 변환에 필요한 기준수식 및 변환에 필요한 투영식

왜곡량 계산 및 변환비용의 산정 예시

MATLAB을 이용한 왜곡량 모델링 프로그램 예시

관련 기준 등

관련기준

고시자료

측량의 사용기준 고시 : 건설교통부 제2002-304호(□02.12.31)

국가좌표변환 계수 : 국토지리정보원 제2002-433(□02.12.30)

국가좌표변환 계수 변경 : 국토지리정보원 제2003-497호(□03.12.24)

관련 지침

세계측지계 기술지침서(국토지리정보원)

지침의 이력

이력관리(Ver)	제□개정일	개정판 내용
제1판	2005년 10월 18일	제정

기대 효과

1/1,000 수치지형도를 새로운 국가측지계에 기준한 수치지형도로 변환하여 사용함으로써 수치지형도의 활용성 증대

GPS 측량성과가 변환 과정 없이 호환되므로 측량 비용의 감소와 함께 위치 정확도 향상에 따른 측량 작업의 효율성 증대

전환이 용이하도록 표준 변환 지침을 제공함으로써 1/1,000수치지형도를 비롯한 관련 DB의 세계측지계로의 변환촉진 도모

GPS 측량의 활성화와 함께 높은 위치 정보의 실시간 연동이 가능하게 되므로 GIS와 GPS 산업발전의 촉매제 역할 기대

국가에서 통일된 변환 지침을 제공함으로써 각 지자체별 대축척 수치지형도 변환 작업의 효율성을 극대화하고 일관성 있는 지리정보의 변환을 통한 국가 지리정보의 신뢰성 향상

측지기준의 전환에 따른 사용자 혼란과 오차를 최소화

국제기구(ICAO, IMO, ITU, IHO 등)에서 권고하고 있는 세계측지계를 사용함으로써의 국제협력 및 국제 표준화 사업에 기여

1/1,000 수치지형도 좌표변환 지침

총칙

지침의 목적

본 지침은 지역측지계를 기준으로 제작된 1/1,000수치지형도를 세계측지계로 변환함에 있어 표준 작업방법 등을 제시함으로써 변환의 효율성을 확보하고, 변환된 성과가 필요한 정확도를 확보할 수 있도록 하는 것을 목적으로 한다

용어의 정의

□국가기준점□이란 건설교통부 국토지리정보원에서 절대적인 위치를 결정한 기준점(삼각점)을 말하는 것으로, 측량성과는 국가기준점이 정비되기 이전에 고시된 성과를 □구성과□, 기준점이 정비된 후 고시된 성과를 □신성과□로 구분한다

□잠정성과□란 정확한 수치지형도 제작을 위해서는 정확도가 확보된 기준점(삼각점)의 사용이 불가피함에 따라 국가기준점(삼각점) 정비과정에서 산출된 것으로 조정하여 고시되지 않은 측량성과를 잠정성과라 말한다

□항측기준점□이란 항공사진측량 방법으로 1/1,000수치지형도를 제작하기 위하여 국가기준점(잠정성과 포함)을 기초(절대적인 값으로 고정)로 대상 지점의 절대좌표를 얻기 위하여 측량된 지상기준점(평면기준점)을 말한다

□공통점□이란 수치지형도 제작 당시에 사용된 기준점(국가기준점 및 항측기준점 등)에 대하여 2개의 측지계(지역측지계, 세계측지계) 성과(평면직각좌표값)를 모두 가지고 있는 지점을 말한다

□왜곡량□이란 공통점의 지역측지계 성과를 좌표변환 과정을 거쳐 변환된 세계측지계 기준의 평면직각좌표와 GPS 관측을 통해 세계측지계 기준으로 계산한 평면직각좌표와의 차이를 말한다

적용범위

적용원칙

본 표준작업 지침은 지역측지계로 제작된 1/1,000수치지형도를 세계측지계로 변환할 경우에 적용하되, 세계측지계로의 변환과 함께 해당 성과를 수정·갱신하는 경우에는 본 표준지침을 준용하여 사용하거나 또는 별도의 기준을 마련하여 적용할 수 있다

공공측량계획기관에서 1/1,000수치지형도를 세계측지계로 변환하고자 할 경우에는 공공측량 작업절차에 따라 시행하여야 한다

관련 DB 변환에 활용

본 표준작업 지침은 1/1,000수치지형도를 기반으로 사용하는 관련 DB(도로 및 지하시설물 등)를 세계측지계로 변환할 경우에도 이를 준용하여 활용할 수 있다

관련 DB를 세계측지계로 변환할 경우에는 1/1,000수치지형도 변환시 사용하였던 변환 기준을 동일하게 사용하여야 한다

예외 적용기준

1/1,000수치지형도를 세계측지계로 변환함에 있어 잔여 왜곡량이 발생하지 않아 표준지침을 적용할 필요가 없거나 기관 또는 사업의 특성상 본 지침을 적용하기가 곤란한 경우에는 공공측량작업 절차에 따라 작업규정을 승인 받아 예외 적용이 가능하다

표준 작업절차

세계측지계로의 변환 개념

좌표변환의 정확성과 신뢰성의 증진을 위해서 국가에서 고시된 변환계수를 적용하여 기본 변환을 수행함과 동시에 각 지역별 왜곡량 모델링을 통하여 잔여 왜곡량 보정이 이루어져야 하며 기본적인 전환 개념은 별첨1과 같다

1/1,000수치지형도의 좌표 변환 후에는 소요의 정확도 확보가 매우 중요하므로 세계측지계로 변환 후 성과검증을 실시하여야 하며, 변환된 성과는 요구 정확도 범위 이내이어야 한다

지침에 따라 세계측지계로 변환했음에도 불구하고 요구되는 정확도가 확보되지 않는 경우에는 1/1,000수치지형도를 신규로 제작하는 등 다른 방안을 강구하여야 한다

작업절차

1/1,000수치지형도를 세계측지계로 변환하기 위해서는, 첫째 기준점성과 등의 분석 및 변환에 필요한 공통점의 확보 등의 준비단계, 둘째 세계측지계로의 기본 변환, 셋째 정확도 확보를 위한 왜곡량의 Modeling 및 보정, 마지막으로 변환된 1/1,000수치지형도의 정확도 평가 과정을 걸쳐 수행되어야 하며 변환작업의 흐름도는 별첨 2와 같다

준비단계

수치지형도 제작 당시의 각종 현황을 파악하여 공통점을 확보하고 지침에 따라 세계측지계로 변환할 경우 정확도 확보가 가능한지 여부 등을 사전에 검토한 후 변환 방법을 선정하고 변환에 필요한 작업계획을 수립하는 과정을 말한다

세계측지계로의 기본 변환

1/5,000이상 축척의 측량성과(수치지형도 등)를 세계측지계로 변환하기 위하여 고시(국토지리정보원 고시 제2003-497호, □03.12.24)된 기준에 의거 세계측지계로 변환하는 과정을 말한다

세계측지계로의 기본 변환은 지역측지계의 평면직각좌표를 지역측지계의 경위도로 변환하고, 지역측지계의 경위도를 3차원 직각지심좌표로 변환하고, 지역측지계의 3차원 직각지심좌표를 세계측지계의 3차원 직각좌표로 변환한 후, 세계측지계의 경위도를 세계측지계의 평면직각좌표로 변환하여야 한다

왜곡량의 모델링(Modeling) 및 보정

1/1,000수치지형도는 제작 당시에 사용한 기준점성과의 특성 등이 있어 기본 변환만으로는 정확도 확보가 불가능하므로 위치정확도를 확보하기 위해 세계측지계로 기본 변환한 후 잔여 왜곡량 모델링을 통하여 결정된 왜곡량을 보정하여야 한다

정확도 평가

세계측지계로 변환된 1/1,000수치지형도가 요구되는 정확도가 확보하는지 검토하기 위해 변환된 수치지형도상에 검사점을 선정하고 측량을 실시 한 후 정확도를 검증하는 과정을 말한다

작업 단계별 적용기준

준비 단계

수치지형도제작 성과의 분석

세계측지계로 변환하고자 하는 1/1,000수치지형도에 대하여 제작시 제출된 용역결과보고서 등을 이용하여 사용된 기준점성과 등을 분석하여야 한다

항공사진측량의 표정을 위하여 사용된 국가기준점의 성과를 구성과, 신성과 및 잠정성과로 구분하고 해당 지역의 수치지형도 제작 당시에 사용된 기준점성과가 현재 보유여부 등을 파악하여야 한다

변환작업 방법의 결정

공통점 등 각종 현황 및 상태를 분석한 후 일부 지역을 시범 변환하여 요구되는 정확도로의 변환이 가능한지 여부 등을 사전에 검토한 후 변환 방법을 결정하여야 한다

표준지침에 따라 변환함에도 불구하고 요구되는 정확도가 확보되지 않을 경우에는 1/1,000수치지형도를 신규로 제작하는 등의 방안을 강구하여야 한다

검토결과 잔여 왜곡량 등이 없어 표준지침을 적용할 필요가 없거나 수정을 병행하는 등 사업의 특성상 본 지침을 적용하기가 곤란한 경우에는 다른 변환 방법을 모색하여야 한다

작업계획의 수립

1/1,000수치지형도를 세계측지계로 전환하기 위한 시범변환 결과 표준 작업지침 적용에 문제가 없을 경우에 한하여 본 지침에 의거 작업계획을 수립하여야 한다

세계측지계로의 전환 계획에는 좌표변환 작업방법, 사용하는 주요한 기기, 프로그램, 작업자, 추진일정 등에 대하여 적절한 작업계획을 수립하여 시행하여야 한다

공통점의 확보

공통점의 확보기준

확보기준

공통점은 수치지형도 제작 당시에 사용한 국가기준점 및 항측기준점(표고기준점은 제외)을 대상으로 하여야 하며, 공통점은 지역측지계 및 세계측지계에 대한 평면직각좌표를 가지고 있어야 한다

공통점은 도시별로 확보하되 변환하여야 할 대상지역에는 반듯이 일정 수량 이상의 공통점이 있어야 한다

공통점에 대해 지역측지계와 세계측지계 성과가 모두 존재할 경우에는 1/1,000수치지형도 제작시 사용된 모든 평면기준점을 활용하는 것을 원칙으로 한다

공통점에 대한 성과가 존재하지 않거나 또는 일부 존재할 경우에는 제작시 측량된 평면기준점의 80% 이상 또는 일정 수량이 균일하게 분포되도록 공통점 확보를 위한 별도의 측량을 실시하여야 한다

· 공통점의 수량 = 10점+(도엽수×0.25□/1□)점

공통점의 확보지역

도시별로 공통점을 확보하여야 함에도 불구하고 인접된 도시가 있고 변환할 1/1,000수치지형도가 존재하는 경우에는 인접 지역은 공통점을 추가로 확보하여 공동으로 사용하여야 한다

도심지역이라 할지라도 1/1,000수치지형도가 제작되지 않은 지역은 공통점을 확보하기 위한 측량을 수행하지 않아도 된다

공통점의 확보방법

공통점의 선정조건

기존의 지역측지계를 기준하여 제작된 수치지형도를 세계측지계에 기준한 수치지형도로 변환하기 위해서는 공통점의 확보가 필요하다

공통점은 수치지형도제작 용역결과보고서 등에서 확보가 가능하며, 3가지의 경우로 나뉘어진다

경우 ① : 용역완성결과보고서 등에서 두 측지계에 기준하여 계산된 성과를 얻을 수 있는 경우

경우 ② : 용역완성결과보고서 등에서 지역측지계에 기준하여 계산된 성과만 얻을 수 있는 경우

경우 ③ : 용역완성결과보고서 등에서 측량성과가 없는 경우

경우 ①과 ②는 다음과 같이 세부적으로 나뉘어진다

경우 ⓐⓓ : 계산이 완료된 성과와 기준점에 대한 GPS 관측 데이터가 존재하는 경우

경우 ⓑⓔ : 계산이 완료된 성과만 존재하는 경우

경우 ⓒⓕ : 기준점에 대한 GPS 관측 데이터만 존재하는 경우

조건별 공통점의 확보방법

㉮ 방법의 추진내용

기준점에 대한 GPS관측 데이터를 인접한 1등기준점(GPS상시관측소)의 세계측지계에 기준한 절대 좌표를 고정하여 기선해석을 실시한다

기선해석된 공통점(국가기준점 및 항측기준점)의 세계측지계 성과와 용역완성결과보고서의 세계측지계 성과와의 차이를 구한다

차이가 평균 3cm를 벗어나는 경우에는 공통점의 세계측지계 성과는 기선해석 하여 계산된 성과를 사용하여야 한다

평균 3cm의 한계를 규정하는 이유는 『정밀2차기준점측량 작업규정의 제3조(좌표의 정확도)』를 준수하기 위함이다

㉯ 방법의 추진내용

용역완성결과보고서의 2개의 측지계에 기준한 성과를 확보한다

2개의 측지계 성과를 공통점으로 이용한다

㉰ 방법의 추진내용

기준점에 대한 GPS 관측 데이터를 인접한 1등기준점(GPS상시관측소)의 두 측지계에 기준한 절대 좌표를 각각 고정하여 기선해석을 실시한다

2개의 측지계 성과를 공통점으로 이용한다

㉱ 방법의 추진내용

기준점에 대한 GPS 관측 데이터를 인접한 1등기준점(GPS상시관측소)의 두 측지계에 기준한 절대 좌표를 각각 고정하여 기선해석을 실시한다

재계산한 공통점(국가기준점 및 항측기준점)의 지역측지계 성과와 용역완성결과보고서의 지역측지계 성과와의 차이를 구한다

차이가 평균 3cm를 벗어나는 경우에는 공통점의 지역측지계 성과는 용역완성결과보고서의 성과를 사용하도록 한다

㉲ 방법의 추진내용

수치지형도 제작 당시에 이용된 기준점에 대한 점의 조서를 확보하고 조서상의 측점 위치에 현장 GPS 관측(정적관측)을 실시한다

현장 GPS 관측은 항측기준점 측량 방법에 의한다

점의 조서상의 기준점에 대한 GPS 관측 데이터를 인접한 1등기준점(GPS상시관측소)의 두 측지계에 기준한 절대 좌표를 각각 고정하여 기선해석을 실시한다

2개의 측지계 성과를 공통점으로 이용한다

세계측지계로 기본 변환

변환의 개념

세계측지계로의 기본 변환은 지역측지계의 평면직각좌표를 세계측지계 기준의 평면직각좌표로 변환하는 것으로써 아래의 단계를 걸쳐 변환하여야 한다

지역측지계의 평면직각좌표를 지역측지계의 경위도로 변환

지역측지계의 경위도를 지역측지계의 3차원 지심직각좌표로 변환

지역측지계의 3차원 지심직각좌표를 세계측지계의 3차원 지심직각좌표로 변환

세계측지계의 3차원 지심직각좌표를 세계측지계의 경위도로 변환

세계측지계의 경위도를 세계측지계의 평면직각좌표로 변환

지역측지계의 경위도를 세계측지계의 경위도로 변환시 고시된 내용에는 1/5,000축척을 초과하는 경우에는 사용할 수 없으나 본 지침에 의거 변환하는 것을 전제로 1/1,000수치지형도에도 적용할 수 있다

단계별 변환 방법

지역측지계로 제작된 수치지형도상의 평면직각좌표를 TM 역투영식에 의하여 Bessel1841 타원체면 상의 경위도 좌표(측지좌표 또는 구면좌표)로 변환하는 기준은 첨부 1의 『평면직각좌표를 타원체면상의 경위도 좌표로 변환』기준을 사용하여 변환하여야 한다

지역측지계상의 경위도 좌표를 세계측지계에 기준한 3차원 지심직각좌표로의 변환은 측량법 제21조제1항의 규정에 의거 기존 1/5,000 이하 축척에 사용하기 위하여 고시된 국가좌표변환계수(국토지리정보원 고시 제2003-497호, 2003.12.23)를 사용하여 변환하여야 한다

세계측지계로 변환된 3차원 지심직각좌표([eqn: X_{ g }],[eqn: Y_{ g }],[eqn: Z_{ g }])를 GRS80 타원체면 상의 경위도 좌표(측지좌표 또는 구면좌표)로 변환하는 기준은 별첨 3의 『변환된 3차원 지심직각좌표를 평면직각좌표로 변환』기준을 사용하여 변환하여야 한다

왜곡량의 모델링(Modeling) 및 보정

작업순서

왜곡량의 모델링(Modeling) 및 보정은 기준점에 포함된 우연오차와 좌표변환에 포함된 수학적인 한계에 의하여 발생하는 오차 등을 소거하여 세계측지계 기준으로 변환한 수치지형도의 정확도를 향상시키기 위하여 수행하는 작업이다

왜곡량 모델링 및 왜곡량의 보정은 다음의 순서에 따라 수행하여야 하며 왜곡량의 계산 및 모델링에 대한 예시는 별첨 4와 같다

왜곡량의 계산

비상사성인 공통점의 제거

편중된 공통점 제거

왜곡량의 경향분석 및 왜곡량 계산

해석적 공분산 함수 결정

최소 제곱 콜로케이션법에 의한 격자화

왜곡량 격자파일 제작

도곽 왜곡량을 이용한 내삽

작업방법

왜곡량의 계산

공통점의 지역측지계 성과를 좌표변환 과정을 거쳐 변환된 세계측지계 기준의 평면직각좌표와 측량을 통해 세계측지계 기준으로 계산한 평면직각좌표와의 차이를 계산하여야 한다

왜곡량의 계산은 세계측지계 기준의 평면직각좌표와 지역측지계의 성과를 기본 변환한 평면직각좌표계상의 X축(North)과 Y축(East)의 좌표 차이 성분인 왜곡량(△X, △Y)을 따로 계산하여야 한다

비상사성인 공통점의 제거

왜곡량 모델링 결과의 신뢰성을 확보하고 일관성을 유지하기 위하여 공통점의 왜곡량이 모든 공통점과 비교하여 비상사성인 공통점은 제거하여야 한다

비상사성인 공통점의 제거는 축별(X축, Y축)로 분석하여야 하며, 한 개의 축이라도 기준 이상을 초과하는 경우에는 이를 제거하여야 한다

비상사성인 공통점의 제거 기준은 통계학적으로 신뢰구간을 설정하여 왜곡량이 신뢰구간을 벗어난 경우이거나 또한 공통점으로 사용하기가 곤란한 공통점은 제거하여야 한다

평면직각좌표계의 각 축별 왜곡량이 표준편차(±□)보다 3배(±3□)이상의 차이를 나타나는 경우

공통점이라 할지라도 1/1,000수치지형도 제작 또는 갱신을 위한 사진기준점측량 작업시 오차가 많이 발생하여 소거된 지점

편중된 공통점 제거(Data Thinning)

고르게 분포되지 않은 공통점중에서 특정점으로 인하여 변환시 중대한 문제가 발생되지 않도록 편중성을 분석하여 편중된 데이터를 제거하여야 한다

공통점과 공통점과의 최소거리는 100m로 하는 것을 원칙으로 하나, 불가피한 경우에는 100m이하로 하거나 제거하지 않을 수 있다

왜곡량의 경향분석 및 왜곡량 계산

비상사성과 편중성이 제거된 공통점에 대하여 X축, Y축의 왜곡량에 대한 경향(trend)을 분석하여 경향방정식을 만들어야 한다

선형 회귀식을 사용하여 최소제곱법에 의한 X축, Y축의 경향방정식을 결정하거나 또는 부등각사상변환(affine변환) 식을 이용하여 미지수를 결정한다

경향방정식 또는 부등각사상변환을 이용하여 공통점의 위치에서 변환값을 계산한 뒤 이 값을 공통점의 왜곡량에서 가감하여 X축, Y축의 잔여 왜곡량을 계산하여야 한다

왜곡량의 모델링

X축, Y축의 잔여 왜곡량에 대한 경험적 공분산 함수를 계산하며, 경험적 공분산 함수에 최소제곱 곡선접합법(Least Square Curve Fitting)을 사용하여 경험적 공분산함수들에 대한 적절한 해석적 모델을 선택하여 2가지 변환계수인 초기공분산(분산)값과 상관거리(Correlation Length)를 계산한다

해석적 공분산 함수식은 다음의 가우시안(Gaussian) 공분산 함수를 사용하도록 한다

[eqn: C(S)=C(o)e^{ -S^{ 2 } A^{ 2 } }]

공분산 함수식에서 [eqn: S]는 거리, [eqn: A]는 상관거리, [eqn: C(o)]:[eqn: S=0]  일 때에 [eqn: C(S)]  로써 취한 값인 초기 분산값이다

왜곡량의 보정기준

인접도시와의 관계

인접지역에 1/1,000수치지형도가 제작된 도심지역이 존재할 경우 인접 1/1,000수치지형도와의 인접을 위하여 상호 인접되는 지점에 대하여 공동의 공통점을 반영하여야 한다

왜곡량의 보정간격 등

왜곡량 보정을 위하여 산정되는 보정량의 간격은 1/1,000수치지형도 도곽 간격(18˙∀0.005¸)으로 설정하고 계산되어야 한다

세계측지계에서 규정하는 수치지형도의 표고 기준은 기존의 국가수준면(평균해수면)을 준용하므로 표고 방향 왜곡량의 계산과 보정은 하지 않는 것을 원칙으로 한다

변환된 1/1,000수치지형도의 정확도 평가

정확도 평가 방법

현장 검증측량 방법

변환된 수치지형도의 변환 정확도를 평가하기 위해서는 변환된 수치지형도상에서 명확하게 확인이 되는 위치를 선정하여 현장 측량을 실시하여야 한다

검사점은 차도와 인도 경계부, 도로 구획선, 중앙 분리대, 맨홀 등과 같이 위치확인이 용이한 지점을 선정하여야 한다

현장 검증측량 지점은 변환하고자 하는 대상 지역이 균일하게 분포되도록 배치하고 현장검증 측량은 4급기준점 측량 방법을 적용하여야 한다

GPS 측량방법으로 검사점을 측량할 경우에는 정적 GPS 측량을 실시하되 소정의 정확도를 확보하기 위하여 GPS 위성신호 수신 시간은 최소 30분 이상을 하도록 한다

현장 검증 측량시 기준점은 세계측지계 성과가 존재하는 삼각점 또는 상시관측소의 세계측지계 기준의 절대 좌표를 기준으로 측량하여야 한다

GPS로 측량시 GPS 데이터 처리는 세계측지계 기준의 절대 좌표를 고정점으로 하여 기선 해석을 실시하여야 한다

현장 검증측량 수량

현장 검증측량은 도엽당(약0.25□) 1점을 기준으로 하되 최소 10점 이상이어야 하며, 지역이 광범위할 경우에는 현장 검증측량 수량을 작업량과 관련하여 조정하여 적용할 수 있다

작업량(도엽)	20이상 ∼50미만	50이상 ∼200미만	200이상 ∼1000미만	1000이상
검사수량	70%	50% (35점 이상)	30% (100점 이상)	10% (300점 이상)

정확도 평가방법

현장에서 측정한 평면직각좌표와 변환한 수치지형도상의 검사점에 해당되는 평면직각좌표와의 차이를 분석하여 정확도를 평가한다

정확도 검사를 위한 통계적 비교 기준은 RMSE(Root Mean Square Error)를 구하여 평가하도록 한다. RMSE를 구하기 위한 계산식은 다음과 같다

[eqn: "RMSE"= SQRT { { SMALLSUM (C_i -C_t )^2} over {n-1 } }]

여기서, [eqn: C_i] : 지도좌표, [eqn: C_t] : 측량 좌표, [eqn: n] : 검사점의 수

정확도 평가기준

변환 수치지형도의 정확도 평가 기준은『공공측량의 작업규정 세부 기준』제98조에서 항공사진에 의한 수치지형도 제작시 수치지형도의 정확도는 제47조의 지도 정확도를 준용하도록 하고 있으므로 이를 평가 기준으로 한다

변환 수치지형도의 정확도 기준은 다음과 같다.

항 목			축 척		비 고
			1/500 이상	1/1,000 이하
표준편차	평면위치		0.5mm 이내	0.7mm 이내	도상거리
	표고	표고점	△h/4 이내	△h/3 이내	△h는 주곡선의 간격
		등고선	△h/2 이내

변환 S/W의 활용 및 비용산정

변환 S/W의 활용

좌표변환 S/W(1000GDKtrans)

개발개요

개발목적

복잡한 계산과정을 프로그램을 이용하여 손쉽게 변환할 수 있도록 지원함으로써 세계측지계로의 변환에 대한 효율성 도모

복잡한 수식에 의하여 변환하는 과정을 S/W를 이용하여 단순화

변환작업 공정을 자동화하여 변환 작업에 대한 효율성 도모

변환에 따른 일관성을 도모하여 측량성과의 정확도 확보

S/W를 이용하여 일괄변환 함으로써 측량성과 균질성 확보

사용조건

1000GDKtrans는 1/1,000 수치지형도를 세계측지계로 변환하기 위해서만 사용이 가능하다

지역측지계(Bessel 1841타원체) 기준으로 작성된 1/1,000수치지형도를 변환하기 위하여 개발된 프로그램으로 남한의 육지(제주도 포함) 및 연안에 해당하는 수치지형도의 변환에만 사용이 가능

1000GDKtrans S/W의 소유권은 국토지리정보원에 있으며, 공공 목적을 위해서는 국토지리정보원 홈페이지(http://www.ngii.go.kr)에서 S/W를 다운받아 무료로 활용할 수 있다

왜곡량 계산이 완료된 지역으로써 왜곡량 모델링을 통하여 왜곡량 값이 입력된 지역의 수치지형도만 세계측지계로 변환 가능

입력된 이외의 지역은 왜곡량이 포함되어 있지 않아 왜곡량 보정이 불가능하며, 변환될지라도 단순히 국가좌표변환계수만을 이용한 변환만 가능

수치지형도 좌표변환 도구의 설치

설치 환경 및 S/W의 구성

1000GDKtrans S/W의 설치 및 구동환경은 다음과 같다.

O/S : Windows-9X/NT 이상

기종 : IBM PC 호환기종

CPU : Pentium Ⅱ 300Mhz 이상

RAM : 기본메모리이상 (64M 이상)

HDD : 200M 이상 (변환하고자 하는 자료량에 따라 달라짐)

1/1,000 수치지형도의 좌표변환을 위한 프로그램은 2개의 폴더(conf 및 distortion)와 하나의 실행파일(1000GDKtrans.exe)로 구성이 되어있으며, conf 폴더는 왜곡량 정보파일(infor.dat)이 distortion 폴더에는 왜곡량 격자파일(dx.txt 및 dy.txt)로 구성되어 있다.

1000GDKtrans.exe : 1/1,000수치지형도를 변환하는 실행프로그램

conf 폴더 : infor.dat 파일이 있으며, 워드패드나 메모장 등에서 작성이 가능하며, distortion 폴더에 있는 자료의 격자점의 갯수 등을 정하고 있는 파일

distortion 폴더 : dx.txt와 dy.txt 파일로 구성되어 있어 워드패드나 메모장 등에서 작성 및 수정이 가능하며, X·Y축의 왜곡량을 표시하는 파일

1000GDKtrans는 지역측지계 성과를 세계측지계로 정변환하거나 또는 세계측지계 성과를 지역측지계로 역변환할 경우 사용이 가능하며, DXF 또는 NGI Format에 한하여 사용할 수 있다

지역별 왜곡량 작성방법

conf 폴더

infor.dat는 distortion 폴더에 있는 자료의 격자 라인수, 왜곡량 격자의 범위 및 격자간격을 정하고 있는 파일이다

X축으로 3개 및 Y축으로 4개의 격자 라인이고 변환하고자하는 구역의 좌상단의 위도가 36.5¸이고 경도가 128.1¸인 경우 작성 방법은 아래의 예시와 같다

conf 폴더 예시		내 용
[Distortion Info]	∀	파일의 헤더
NX=3	∀	X축 방향(N)의 왜곡량 격자 라인수(예 : 3 라인)
NY=4	∀	Y축 방향(E)의 왜곡량 격자 라인수(예 : 4 라인)
NXE=36.5	∀	왜곡량 격자의 위도 범위(상단 위도, 단위 : °)
NYS=128.1	∀	왜곡량 격자의 경도 범위(좌단 경도, 단위 : °)
Interval=0.005	∀	왜곡량 격자 간격(단위 : °)

distortion 폴더

dx.txt는 X축 방향에 대한 왜곡량값을 표시하며, dy.txt는 Y축 방향에 대한 왜곡량값을 격자 라인수별로 명시한 2개의 아스키 파일로 구성된다

*.txt Data의 왜곡량값은 최소한 소수점 5째 자리 이상으로 구성하여야 하며, 격자별 Data의 구분은 공백(Space)으로 구별한다

dx.txt 및 dy.txt 작성 기준은 아래와 같다

· 첫째행 : 왜곡량의 최소값과 최대값을 명시

· 둘째행 이후 : Y축 격자 라인수에 해당하는 격자점별 X 또는 Y방향의 왜곡량을 공백(space)으로 구분하여 입력하되 Y축 격자 라인수가 10개 이상을 초과할 경우에는 한행에 10개씩 줄바꿈(□)으로 구분하여 입력하고, 모든 격자값이 표기되어 X축 라인이 변경될 경우에는 빈줄을 삽입

X축으로 3개 및 Y축으로 4개의 격자를 가지는 경우 dx.txt 및 dy.txt 파일의 작성 형태는 다음과 같다

dx.txt 작성예시		내 용
0.073046939661912 0.16485397133959	∀	Min, Max
0.14084978566105 0.13998660363273 0.13917240596561 0.13843999225511	∀	4개 Data 줄바꿈
0.13954686062261 0.14035773847061 0.1411846768414 0.14199530710432	∀	4개 Data 줄바꿈
0.14773886018632 0.14849753809762 0.14919162408192 0.14983941223018	∀	4개 Data 줄바꿈

dy.txt 작성예시		내용
-0.12414107799183 -0.043316168218771	∀	Min, Max
-0.1169500057748 -0.11441009787301 -0.11187620999873 -0.10934784166526	∀	4개 Data 줄바꿈
-0.11655302385944 -0.11411856269966 -0.11168671171767 -0.1092620970904	∀	4개 Data 줄바꿈
-0.11616446581696 -0.11383042978425 -0.11150067117134 -0.10917829165958	∀	4개 Data 줄바꿈

변환S/W의 활용

사용전 준비사항

왜곡량의 격자 정보를 1000GDKtrans의 conf 폴더의 infor.dat를 열어 수정하고, 작성된 왜곡량의 격자 파일을 distortion 폴더의 dx.txt와 dy.txt와 교체한 뒤 1000GDKtrans.exe를 실행하여야 하며, 왜곡량 정보가 포함된 지역에 한하여 수치지형도의 변환이 가능하다

변환하고자 하는 지역의 왜곡량 정보를 입력하고 지역측지계 성과를 세계측지계로 변환하고자 하는 수치지형도를 확보하여 변환을 실시하여야 한다

변환도구 사용방법

S/W의 구성

1000GDKtrans S/W는  주메뉴·파일관리·Layer관리 창으로 구성되어 있으며, 주요 기능은 주메뉴의 창에서 구현되도록 되어 있다

S/W의 주기능을 보완하는 파일관리·Layer관리 및 메인화면에 마우스의 위치를 표시하는 좌표값보기(↙↔)로 구성되어 있다

S/W의 활용방법

1000GDKtrans의 주된 기능은 주메뉴 창의 메뉴로 구성되어 있으며, 주메뉴의 기능은『파일열기·닫기·저장·출력 기능』과 메인화면을 컨트롤하는『화면확대·축소·이동·전체화면 보기·화면의 바탕색을 조정하는 White』와 Layer 구성 상태를 표시해주는 Layer 관리창으로 구성되어 있다

정변환은 지역측지계 성과를 세계측지계로 변환하기 위하여 사용하는 기능으로 세계측지계 변환 파라미터에 의한 변환이 이루어진 후 conf 및 distortion 폴더에 주어진 왜곡량을 이용하여 세계측지계로 일괄 변환한다

역변환은 세계측지계 성과를 지역측지계로 변환하기 위하여 사용하는 기능으로 세계측지계 변환 역파라미터에 의한 변환이 이루어진 후 conf 및 distortion 폴더에 주어진 역왜곡량을 이용하여 지역측지계로 일괄 변환한다

왜곡량 계산 S/W개발 및 활용

S/W의 개발 및 보급

1/1,000수치지형도를 세계측지계로 전환하기 위한 도시별 왜곡량을 자동 산출하는 프로그램을 개발하여 활용할 경우에는 국토지리정보원과 협의하여야 한다

개발하는 S/W는 동 지침에서 제시하고 있는 표준작업지침을 사용하여 개발되어야 한다

왜곡량 계산 S/W의 인증

왜곡량 계산 S/W를 개발하여 보급 할 경우에는 국토지리정보원에 왜곡량 계산에 대한 검증 및 인증을 받아야 한다

인증되지 않은 S/W는 세계측지계변환에 활용할 수 없으며, 이를 보급하여서는 안된다

변환비용의 산정

적용원칙

1/1,000수치지형도 또는 각종 GIS Data에 대하여 지역측지계 성과를 세계측지계로 변환하고자 할 경우에는 이에 소요되는 비용을 지급하여야 한다

용역대가의 적용은 측량용역대가의 기준을 적용하여야 한다

비용산정 기준

1/1,000 수치지형도 좌표변환

산정기준

1/1,000수치지형도에 대하여 지역측지계 성과를 세계측지계로 정변환하거나 또는 세계측지계 성과를 지역측지계로 역변환하는 경우에 사용하여야 한다

1/1,000수치지형도의 좌표계 변환에 소요되는 비용은 1도엽(약 0.25□) 단위를 기준으로 하며, 이에 소요되는 비용은 아래와 같으며 비용산정 방법에 대한 예시는 별첨 5와 같다

(단위 : 인/도엽당)

구 분	특 급 기술자	고 급 기술자	초 급 기술자	정보처 리기사	중 급 기능사	비 고
현황조사 및 성과분석	0.05	0.05
왜곡량 계산 및 모델링	0.05	0.05
변환프로그램 제작		0.05		0.05
세계측지계 변환		0.05	0.1		0.5
GPS 현지측량						품셈적용
정확도 분석	0.05	0.05
합 계	0.15	0.25	0.1	0.05	0.5

※ 중급기능사는 지도제작기능사를 적용

적용기준

1/1,000수치지형도를 세계측지계로 변환하기 위한 비용의 산정기준은 1/1,000수치지형도 1도엽(약 0.25□)을 1작업단위로 한다

작업량에 따른 증감계수는 다음 계수를 곱하여 계상한다

작업량 (도엽)	10이상 ∼50미만	50이상 ∼200미만	200이상 ∼1000미만	1000이상	비고
작업량에 따른 증감계수	1.0	0.8	0.6	0.5

변환하고자 하는 지역에 대하여 공통점이 확보된 지역을 기준으로 하며, 공통점 확보를 위한 측량이 별도로 필요한 경우에는 별도 산정한다

공통점에 대한 측량 및 변환성과의 정확도 분석을 위한 현지측량 작업은 건설표준품셈(21-3-4, 4급기준점 측량)을 적용한다

수치지형도는 DXF 또는 NGI Format을 기준으로 도곽변경 없이 변환하는 것으로써 다른 Format을 변환하거나 또는 도곽 변경 등이 필요한 경우에는 별도의 품을 적용할 수 있다

각종 GIS 데이터의 좌표변환

산정기준

1/1,000수치지형도를 기반으로 사용되는 각종 GIS Data로써 변환에 필요한 왜곡량이 이미 계산되어 1/1,000수치지형도의 좌표변환이 선행된 후 연계 시행하는 작업에 한하여 적용할 수 있다

(단위 : 인/도엽당)

구 분	고급 기술자	초급 기술자	중급기능사 (지도제작)	비 고
각종 GIS Data 변환	0.1	0.2	0.6

적용기준

각종 GIS Data의 좌표계 변환에 소요되는 비용은 1도엽(약 0.25□) 단위를 기준으로 한다

작업량에 따른 증감은『1/1,000 수치지형도 좌표변환』증감계수를 준용하여 계상할 수 있다

1/1,000수치지형도 변환없이 각종 GIS Data만 변환하는 경우에는 GIS Data 변환비용 이외에 1/1,000 수치지형도 변환에 소요되는 조사 및 성과분석, 왜곡량 모델링 및 정확도 검증 비용 등은 별도로 산정하여야 한다

<별첨 1>

세계측지계로의 전환 개념도

<별첨 2>

세계측지계로의 변환작업 흐름도

<별첨 3>

세계측지계로의 기본 변환 기준수식

수치지형도의 평면직각좌표를 3차원 직각좌표로 변환

TM 투영식은 <첨부 1>를 사용하도록 하며, Bessel1841 타원체의 제원과 투영 원점의 가산값은 다음과 같다

장반경 a = 6 377 397.155 m

편평률 f = 1/299.1528128

축척계수 [eqn: k_0]=1.0000

원점 가산값 □X(N) = 500,000 m (단, 제주 지역은 550,000 m)

□Y(E) = 200,000 m

Bessel1841 타원체면상의 경위도 좌표를 지역측지계상의 3차원 직각좌표로 변환하는 식은 다음과 같다

[eqn: bmatrix { X_b # Y_b # Z_b } = bmatrix { (N+h_b )``cos □_b ``cos □_b # (N+h_b )``cos □_b ```sin □_b # (N{(1-e^{ 2 } )+h_b )} sin □_b `` }]

여기서,

[eqn: &N(묘유선의~ 곡률~반경)= { a } over { sqrt { 1-e^ 2 sin^{ 2 }□ } } ##&e^2~ :~ 타원체의~ 이심률 ##&□_b ~:~rmBessel ~타원체상의 ~위도, ~□_b ~:~rmBessel ~타원체상의~경도 ##&h_b ~:~rmBessel ~타원체고]

지역측지계상의 3차원 직각좌표를 세계측지계로 변환

지역측지계상의 3차원 직각좌표를 세계측지계에 기준한 3차원 지심직각좌표로 변환을 수행한다

3차원 좌표변환에 필요한 국가좌표변환계수와 변환식은『국토지리정보원 고시 제2003-497호』에서 제시한 Molodensky-Badekas 모델에 의하여 결정된 값과 식을 사용한다

고시된 국가좌표변환 계수와 변환식은 다음과 같다

구분	평행이동량(m)			회전량(□)			축척변화량(ppm)
	□x	□y	□z	Rx	Ry	Rz	□
변 환 계 수	-145.907	505.034	685.756	-1.162	2.347	1.592	6.342

[eqn: bmatrix { X_{ g } # Y_{ g } # Z_{ g } } = bmatrix { X_{ o } # Y_{ o } # Z_{ o } } + bmatrix { □X # □Y # □Z } +(1+□) bmatrix { 1 & R_{ z } & -R_{ y } # -R_{ z } & 1 & R_{ x } # R_{ y } & -R_{ x } & 1 } bmatrix { X_{ b } -X_{ o } # Y_{ b } -Y_{ o } # Z_{ b } -Z_{ o } }]

단, [eqn: X_{ g }],[eqn: Y_{ g }],[eqn: Z_{ g }] : GRS80 타원체 상의 3차원 지심직각좌표(세계측지계)

[eqn: X_{ b }],[eqn: Y_{ b }],[eqn: Z_{ b }] : Bessel 타원체 상의 3차원 직각좌표(지역측지계)

[eqn: X_{ o }],[eqn: Y_{ o }],[eqn: Z_{ o }] : 좌표변환 기준좌표

구 분	기 준 좌 표		비 고
	Bessel => GRS80	GRS80 => Bessel
[eqn: X_{ o }] [eqn: Y_{ o }] [eqn: Z_{ o }]	-3,159,521.31 4,068,151.32 3,748,113.85	-3,159,666.86 4,068,655.70 3,748,799.65

변환된 3차원 지심직각좌표를 평면직각좌표로 변환

변환된 3차원 지심직각좌표([eqn: X_{ g }],[eqn: Y_{ g }],[eqn: Z_{ g }])를 GRS80 타원체면 상의 경위도 좌표(측지좌표 또는 구면좌표)로 변환한다

변환된 3차원 지심직각좌표를 측지 좌표로 변환하는 식은 다음을 이용한다

[eqn: &□_g = tan ^-1 {{Z_g +e'^2 ·b· {sin^3}□ }over{p-a·e^2 {cos^3}□}}## &□_g = tan ^-1 {Y_g over X_g}## &h_g =p·cos□+Z_g sin□_g- { a^2} over { N}]

여기서,

[eqn: &a ~:~타원체의~ 장반경,~b~:~타원체의~ 단반경 ## &p=sqrt{X_g^2 +Y_g^2 } `, ~□~=~ tan ^-1 { Z_g ·a} over {p·b } ## &e^2 (제1이심률) = {a^2 - b^2} over a^2 ## &e'^2 (제2이심률) = {a^2 - b^2} over b^2]

GRS80 타원체의 제원은 다음과 같다

장반경 a = 6 378 137 m

편평률 f = 1/298.257222101

GRS80 타원체면 상의 경위도 좌표를 <별첨 1>의 TM 투영식을 사용하여 평면직각좌표로 변환하며, 투영원점의 축척계수와 원점 가산값은 별첨 3의 값과 같다

<별첨 4>

왜곡량의 계산 및 모델링 작업 예시

왜곡량의 계산

공통점에 대한 차이 계산

확보한 공통점의 성과 중에서 지역측지계 기준의 평면직각좌표를 세계측지계로의 기초변환 단계를 거쳐 세계측지계 기준의 평면직각좌표를 도출한다

확보한 공통점의 성과 중에서 세계측지계 기준의 평면직각좌표와 변환된 평면직각좌표간의 차이(X축, Y축)를 계산한다

비상사성인 공통점 제거

비상사성이란 공통점별로 각 축(X축, Y축)의 왜곡량의 크기가 공통점 대부분의 왜곡량보다 현저히 크게 발생된 것으로 비상사성인 공통점을 제거하도록 한다

왜곡량 모델링 결과의 신뢰성을 확보하고, 일관성을 유지하기 위하여 평면직각좌표계의 각 축별 왜곡량의 표준편차(±□)보다 3배(±3□)이상의 차이를 보이는 공통점을 비상사성으로 판단하고 제거하도록 한다

편중된 공통점 제거(Data Thinning)

편중성이란 공통점의 분포가 고르게 이루어지지 않고 일부 집중되어 모여있는 것으로 편중된 공통점을 제거하도록 한다

고르게 분포되지 않은 공통점중에서 특정점의 예측치가 편기되는 문제를 발생시켜 불규칙한 데이터 분포로 인한 불완전한 왜곡 모델링의 수행을 방지하기 위하여 일정한 공간적 분포를 얻기 위한 데이터의 편중성을 분석하여 편중된 데이터를 그림과 같이 제거한다

왜곡량 경향분석 및 왜곡량 계산

비상사성과 편중성이 제거된 공통점의 X축, Y축의 왜곡량에 대한 경향(trend)을 분석하고 선형 회귀식을 사용하여 최소제곱법에 의한 X축, Y축의 경향방정식을 결정하거나 또는 부등각사상변환(affine변환) 식을 이용 미지수를 결정한다

결정된 경향방정식 또는 부등각사상변환을 사용하여 공통점의 위치에서 경향값을 계산한 뒤 이 경향값을 공통점의 왜곡량에서 가감하여 X축, Y축의 잔여 왜곡량을 계산한다

왜곡량의 모델링

왜곡량 모델링

X축, Y축의 잔여 왜곡량에 대한 경험적 공분산 함수를 계산하되 경험적 공분산 함수에 최소제곱 곡선접합법(Least Square Curve Fitting)을 사용한다. 경험적 공분산함수들에 대한 적절한 해석적 모델을 선택하여 2가지 변환계수인 초기공분산(분산)값과 상관거리(Correlation Length)를 계산한다

해석적 공분산 함수식에는 다음의 가우시안 공분산 함수를 사용하도록 한다

[eqn: C(S)=C(o)e^{ -S^{ 2 } A^{ 2 } }]

공분산 함수식에서 [eqn: S]는 거리, [eqn: A]는 상관거리, [eqn: C(o)]:[eqn: S=0]일 때에 [eqn: C(S)]로써 취한 값인 초기 분산값이다

격자 간격의 왜곡량 산정

각 공통점에서의 X축, Y축 왜곡량으로부터 격자 간격의 왜곡량을 생성하고 결정하고자 하는 왜곡량 모델의 격자파일의 각 격자점에서 X축, Y축 왜곡량 성분을 최소제곱콜로케이션법 및 보간법들을 사용하여 계산한다

최소제곱콜로케이션법에 의한 왜곡량 모델링을 적용할 경우에는 다음식을 사용하도록 한다

[eqn: hat { delta } ~ = ~ C(S) C_D ^-1 l]

위 식에서 [eqn: hat { delta }]는 추정된 왜곡량이고, 벡터[eqn: C(S)]의 요소들은 데이터 점들과 보간점간의 거리를 사용하여 해석적 공분산 함수로부터 구한 값이다. 또한, 매트릭스 [eqn: C_D]의 요소들은 공통점들의 모든 조합에 의하여 거리를 사용하여 해석적 공분산 함수로부터 계산된다. 벡터 [eqn: l]은 각 공통점들에서 왜곡을 포함하고 있는 관측치이다

왜곡량 격자파일 작성

수치지형도 제작 지역의 범위를 포괄할 수 있도록 왜곡량의 범위(위도와 경도 범위)를 결정하여 왜곡량 격자 파일을 작성한다

X축, Y축의 잔여 왜곡량의 추정값에 X축, Y축의 경향방정식을 통해 계산된 X축, Y축의 왜곡량을 더한다

최종적으로 계산된 X축, Y축의 왜곡량을 0.005°(1/1,000 수치지형도의 도곽 크기)의 간격으로 격자 파일을 작성한다. 각 축별 격자 파일의 이름은 dx.txt와 dy.txt로 지정한다

<X(North)방향>                   <Y(East)방향>

MATLAB 이용한 왜곡량 모델링

왜곡량 산정의 절차

왜곡량의 계산

왜곡량은 공통점의 성과를 지역측지계 기준의 평면직각좌표를 세계측지계 기준의 평면직각좌표로 변환한 후, 공통점의 성과 중에서 세계측지계 기준의 평면직각좌표와 변환한 직각좌표간의 차이를 계산한 것이다

[eqn: TRIANGLE X=X _{GRS80} -X _{"GRS80"~"trans"}]	(1)

[eqn: TRIANGLE Y=Y _{GRS80} -Y _{"GRS80"~"trans"}]	(2)

여기서, [eqn: X_GRS80]과 [eqn: Y_GRS80]은 공통점의 성과중 세계측지계 기준의 평면직각좌표를 의미하며, [eqn: X_{"GRS80 trans"}]과 [eqn: Y_{"GRS80 trans"}]는 지역측지계 기준의 평면직각좌표를 세계측지계 기준의 평면직각좌표로 변환한 값을 의미한다

비상사성인 점 제거

비상사성인 점의 제거는 모든 공통점을 이용하여 X, Y방향의 왜곡량에 대한 표준편차를 계산한 후, 이 표준편차보다 □3□이상인 공통점을 제거한다

if [eqn: |TRIANGLE X| > 3sigma_X] or [eqn: |TRIANGLE Y| > 3sigma_Y] □ 비상사점	(3)

편중된 공통점 제거

왜곡량의 모델링을 통한 내삽을 수행할 때 공통점의 위치는 매우 중요하다. 만약 공통점이 편중되어 있다면, 이는 최소 제곱 콜로케이션 방법을 이용할 때 데이터의 편향을 일으킬 수 있다. 그렇기 때문에 편중되어 있는 공통점을 제거(data thinning)해야 한다

편중된 공통점의 제거하기 위하여서는 공통점의 최소거리는 100m로 한다

회귀식에 의한 지역적 변화제거

서로 다른 기준면을 모델화하기 위하여 좌표계의 차이에 의해서 발생되는 지역적인 왜곡은 다항식 표면으로 간주될 수 있다. 그렇기 때문에 회귀다항식을 이용하여 이러한 지역적 변화를 제거할 수 있다. 이러한 식은 아래와 같이 정의될 수 있다

[eqn: TRIANGLE X=a _{0} +b _{0} X _{"GRS80"~"trans"} +c _{0} Y _{"GRS80"~"trans"}]	(4)

[eqn: TRIANGLE Y=a _{1} +b _{1} X _{"GRS80"~"trans"} +c _{1} Y _{"GRS80"~"trans"}]	(5)

모든 공통점을 이용하여 다음과 같은 행렬식을 구성한다

[eqn: bmatrix{TRIANGLE X _{1}#TRIANGLE Y _{1}#TRIANGLE X _{2}#TRIANGLE Y _{2}#...#...#TRIANGLE X _{n}#TRIANGLE Y _{n}} = bmatrix{1&X _{1}&Y _{1}&0&0&0#0&0&0&1&X _{1}&Y _{1}#1&X _{2}&Y _{2}&0&0&0#0&0&0&1&X _{2}&Y _{2}#...&...&...&...&...&...#...&...&...&...&...&...#1&X _{n}&Y _{n}&0&0&0#0&0&0&1&X _{n}&Y _{n}} CDOT bmatrix{a _{0}#b _{0}#c _{0}#a _{1}#b _{1}#c _{1}}]□[eqn: bold y = bold B CDOT bold x]	(6)

여기서, [eqn: TRIANGLE X_n]과 [eqn: TRIANGLE Y_n]은 [eqn: n]번째 X와 Y방향으로의 왜곡량이며, [eqn: X _{n}]과 [eqn: Y _{n}]은 지역측지계 기준의 평면직각좌표를 세계측지계 기준의 평면직각좌표로 변환한 X와 Y방향의 값을 나타낸다

식 (6)에서부터 벡터 [eqn: bold x]는 정규방정식(normal equation)에 의해서 해를 구할 수 있다

[eqn: bold{x} = LEFT ( bold{B ^{T}} CDOT bold{B} RIGHT ) ^{-1} CDOT LEFT ( bold{B ^{T}} CDOT bold{y} RIGHT )]	(7)

벡터 [eqn: bold x]가 결정되면, 다음과 같은 방식으로 잔여왜곡량을 결정한다

[eqn: TRIANGLE X _{R} = TRIANGLE X- LEFT ( a _{0} +b _{0} X _{"GRS80"~"trans"} +c _{0} Y _{"GRS80"~"trans"} RIGHT )]	(8)

[eqn: TRIANGLE Y _{R} = TRIANGLE Y- LEFT ( a _{1} +b _{1} X _{"GRS80"~"trans"} +c _{1} Y _{"GRS80"~"trans"} RIGHT )]	(9)

여기서, [eqn: TRIANGLE X _{R}]과 [eqn: TRIANGLE Y _{R}]은 잔여왜곡량을 나타낸다

해석적 공분산함수 결정

데이터의 내삽을 위해서는 거리에 대한 가중치가 필요하다. 이러한 가중치는 경험적 공분산함수를 이용하여 해석적 공분산 함수로 최적화 시킬 수 있다. 잔여왜곡량은 랜덤변수로 가정될 수 있고, 랜덤변수는 가우시안 분포를 지니게 된다. 가우시안 분포는 각 데이터의 분산에 의해서 그 폭이 결정될 수 있게 되므로, 내삽을 위해 사용되는 거리의 가중치가 공분산함수를 이용할 때, 실제 데이터로부터 최적의 해석적 공분산 함수를 구하는 것이 매우 중요하다. 해석적 공분산 함수는 일반적으로 가우스함수 또는 Reilly의 공분산함수를 많이 사용하게 된다[1][3]

먼저, 해석적 공분산 함수를 결정하기 위해서는 실험적 공분산 함수를 결정해야 한다. 실험적 공분산 함수는 다음과 같이 결정될 수 있다

1) 모든 공통점에 대한 잔여왜곡량의 평균을 계산한다

[eqn: mu _{TRIANGLE X _{R}} = {1} over {N} sum _{k=1} ^{N} TRIANGLE X _{Rk}]	(10)

[eqn: mu _{TRIANGLE Y _{R}} = {1} over {N} sum _{k=1} ^{N} TRIANGLE Y _{Rk}]	(11)

2) 첫 번째 점에 대한 실험적 공분산 값을 계산한다. 이때 첫 번째 점과 다른 모든 점들간의 공분산([eqn: C _{X}], [eqn: C _{Y}])과 거리([eqn: r])를 계산한다

[eqn: C _{X1k} = LEFT ( TRIANGLE X _{R1} - mu _{TRIANGLE X _{R}} RIGHT ) CDOT LEFT ( TRIANGLE X _{Rk} - mu _{TRIANGLE X _{R}} RIGHT )]	(12)

[eqn: C _{Y1k} = LEFT ( TRIANGLE Y _{R1} - mu _{TRIANGLE Y _{R}} RIGHT ) CDOT LEFT ( TRIANGLE Y _{Rk} - mu _{TRIANGLE Y _{R}} RIGHT )]	(13)

[eqn: r _{1k} = sqrt {(X _{1} -X _{k} ) ^{2} +(Y _{1} -Y _{k} ) ^{2}}]	(14)

여기서, [eqn: k]는 1에서 [eqn: N]까지의 값이고, [eqn: C _{X1k}]과 [eqn: C _{Y1k}]는 공통점 1과 공통점 k사이의 X와 Y방향으로의 공분산이며, [eqn: r _{1k}]은 공통점 1과 공통점 k사이의 거리를 나타낸다. [eqn: C _{X1k}]과 [eqn: C _{Y1k}]은 [eqn: r _{1k}]의 함수이며, ([eqn: r _{1k}],[eqn: C _{X1k}]) 과 ([eqn: r _{1k}],[eqn: C _{Y1k}])로 표현할 수 있다

3) 두 번째 점부터 마지막 점까지 실험적 공분산 값을 계산한다. 여기서 주의할 점은 두 번째 점을 계산할 경우에는 첫 번째 점은 제외하여 계산한다. 즉, k번째 점을 계산할 경우, k-1번째 이하의 점들은 제외하여 계산한다. 이렇게 하여 모든 점에 대한 실험적 공분산 값을 결정한다

4) 2)와 3)에 의해서 결정된 실험적 공분산 값은 만약 공통점의 개수가 100개라면 X와 Y방향으로 각각 5050개가 결정되게 된다. 이를 그대로 이용하여 해석학적 공분산 함수를 결정하는 것은 무리가 있으며 바람직하지 못하다. 다음과 같이 거리([eqn: r])에 대하여 구간간격을 결정하고, 공분산 값을 평균한다. 예를 들면, 0∼500m, 500m∼1000m, 1000m∼1500m등과 같이 500m 간격내에 들어있는 모든 공분산값을 평균하여 X와 Y방향으로의 공분산을 계산하여 ([eqn: r],[eqn: C _{X}]) 과 ([eqn: r],[eqn: C _{Y}])를 결정한다

5) 1)에서 4)의 과정을 통해 실험적 공분산 값을 우리는 그대로 이용할 수 없다. 각점에는 위에서 언급한 오차들이 포함되어 있기 때문이다. 이들을 제거하기 위해서는 일반적으로 해석적 공분산함수를 이용한다. 해석적 공분산함수는 가우스함수 또는 Reilly의 공분산함수를 이용하게 된다(그림2)

해석적 공분산 함수의 비교

가우스 함수는 아래와 같이 정의된다

[eqn: C(s)`=`C _{0} e ^{-A ^{2} s ^{2}}]	(15)

여기서, [eqn: s]는 거리로 앞에서 설명한 [eqn: r]에 해당하며, [eqn: C_0]는 [eqn: s=0]에서의 분산이고, [eqn: A]는 상관거리로 가우스 함수의 폭을 결정하는 파라미터로 이 값이 크면 함수의 폭이 커져서 먼 거리에서도 많은 영향을 주는 것을 의미하며, 이 값이 작으면 짧은 거리에서만 영향을 주는 것을 의미한다

또한, Reilly의 공분산 함수는 아래와 같이 정의된다

[eqn: C(s)`=`C _{0}LEFT[1- 1over 2 LEFT( s over d RIGHT)^2 RIGHT ] CDOT e ^{- 1 over 2 LEFT( s over d RIGHT)^2}]	(16)

여기서, [eqn: d]는 실험데이터로부터 얻은 파라미터를 의미한다

위의 두 공분산 함수는 가까운 거리에서는 거의 비슷한 결과를 보이다가 거리가 멀어질수록 다른 값을 보이게 된다. 실제적으로는 거의 비슷한 값을 지니게 되지만, 많은 연구자들에 의해서 Reilly의 공분산함수를 사용하는 것이 더욱 나은 것으로 보고되고 있다[1],[2],[3]

위의 두 공분산 함수중 하나를 선택하여 4)에서 계산한 ([eqn: r],[eqn: C _{X}]) 과 ([eqn: r],[eqn: C _{Y}])에 대하여 최소제곱 곡선 접합법(Least Square Curve Fitting)을 수행한다

가우스 함수는 정규방정식을 이용하여 매우 쉽게 파라미터를 계산할 수 있지만, Reilly의 공분산 함수는 반복되는 과정을 통하여 얻어질 수 있다. 가우스 함수를 결정하는 과정을 알아보면, 먼저 식 (15)의 양변에 [eqn: log _{e}]를 취하면,

[eqn: log_e C(s)`=`log_e C _{0} -A ^{2} s ^{2}]	(17)

여기서, 미지수인 [eqn: log_e C]는 [eqn: a _{0}]로 [eqn: -A ^{2}]은 [eqn: b _{0}]로 놓으면, [eqn: s]의 이차함수로 표현가능하다

[eqn: log _{e} C(s)`=`a_0 +b_0 s ^{2}]	(18)

만약 점이 [eqn: (s _{k} ,C(s _{k} ))]가 주어졌을때, 아래의 행렬방정식을 구성한다

[eqn: BMATRIX {log_e C(s_1 ) #log_e C(s_1 ) #... #log_e C(s_N ) } = BMATRIX{1 & s_1^2 #1 & s_2^2 #...& ...#1 & s_N^2 } CDOT BMATRIX {a_0 # b_0}] □ [eqn: bold y = bold B CDOT bold x]	(19)

식 (7)에 의해서 미지수가 결정될 수 있다. 이렇게 미지수 [eqn: a _{0}]와 [eqn: b _{0}]가 결정되면, 아래의 수식에 의해서 [eqn: C_0]와 [eqn: A]를 계산할 수 있다

[eqn: C _{0} =e ^{a _{0}}], [eqn: A`=` sqrt {-b _{0}}]	(20)

위의 과정을 통하여 해석학적 공분산함수를 결정할 수 있다

최소 제곱 콜로케이션법에 의한 격자화

거리에 대한 가중치인 해석적 공분산함수가 결정되면, 이는 최소 제곱 콜로케이션 법에 의하여 격자화를 실시한다. 최소제곱 콜로케이션법은 아래와 같다

원하는 지점에 대한 잔여왜곡량의 계산은 근처에 인접하고 있는 점들을 이용하여 앞에서 계산한 해석적 공분산함수를 가중치로 하여 계산하는 방법으로 아래의 행렬방정식을 통해 계산할 수 있다

[eqn: TRIANGLE R _{P} = bold {C_I} CDOT bold C_D^-1 CDOT bold {TRIANGLE sigma}]	(21)

여기서,

[eqn: bold C_I = BMATRIX {C(d_P1 ) &C(d_P2 ) &...&C(d_Pn ) }]	(22)

[eqn: bold C_D = bmatrix{C(d _{11} )&C(d _{12} )&...&C(d _{1n} )#...&...&...&...#C(d _{n1} )&C(d _{n2} )&...&C(d _{nn} )}]	(23)

[eqn: bold{TRIANGLE sigma} ```=``` bmatrix{TRIANGLE sigma _{1}#TRIANGLE sigma _{2}#...#TRIANGLE sigma _{n}}]	(24)

여기서, [eqn: d _{pq}]는 [eqn: p]점과 [eqn: q]점 사이의 거리를 나타내고, [eqn: TRIANGLE R_P]는 점 [eqn: P]에서의 잔여왜곡량을 의미하며, [eqn: TRIANGLE sigma _{k}]는 k점에서의 X 또는 Y방향의 잔여왜곡량 [eqn: TRIANGLE X _{R}]과 [eqn: TRIANGLE Y _{R}]을 나타낸다

최소제곱 콜로케이션법은 근처의 점들을 이용하여 내삽하는 방식이므로, 근처의 점의 위치와 거리 점의 개수들이 매우 중요한 영향을 미치게 된다. 그러므로, 사용하는 점에 대한 제한 요소를 두어야 한다. 이러한 제한 요소에는 크게 두가지 파라미터를 사용하는 데 첫 번째 파라미터는 제한거리이며, 두 번째 파라미터는 최대사용점의 개수에 해당한다. 그림 3은 제한거리와 최대사용점의 개수에 대하여 잘 나타낸다. 이 파라미터에 해당하는 점들을 가지고 계산하기 위해서는 먼저 내삽하고자하는 점내의 제한거리안에 있는 점을 모두 찾는다. 찾아진 모든 점들의 개수가 최대사용점의 개수보다 적은지(이하인지)를 확인하고, 많다면(이상이면) 모든 점을 사용하는 것이 아니라 가장 가까운 점에서부터 최대사용점에 이를 때까지만 사용하거나, 거리와 방향을 모두 고려하여 결정하면 된다

최대사용점은 5∼10점이고, 제한거리는 4∼5Km로 한다

최소제곱 콜로케이션 파라미터

왜곡량 격자파일 제작

1/1,000수치지형의 도곽크기에 맞는 간격으로 최소제곱 콜로케이션 방법에 의하여 잔여왜곡량을 계산하고, 아래의 수식에 의해 왜곡량을 결정한다

[eqn: TRIANGLE X = TRIANGLE X _{R} + LEFT ( a _{0} +b _{0} X _{"GRS80"~"trans"} +c _{0} Y _{"GRS80"~"trans"} RIGHT )]	(25)

[eqn: TRIANGLE Y =TRIANGLE Y _{R}+ LEFT ( a _{1} +b _{1} X _{"GRS80"~"trans"} +c _{1} Y _{"GRS80"~"trans"} RIGHT )]	(26)

도곽 왜곡량을 이용한 내삽

1/1,000수치지형의 네 모서리 도곽의 왜곡량이 결정되면 그 왜곡량에 의해서 내부점 P는 내삽에 의해 결정될 수 있다. 내삽방법은 다음과 같다

[eqn: z _{P} =a _{0} +a _{1} X+a _{2} Y+a _{3} XY]	(27)

여기서,

[eqn: a _{0} `=`z _{1}]	(28)

[eqn: a _{1} `=`z _{2}-z_1]	(29)

[eqn: a _{2} =z _{4} -z _{1}]	(30)

[eqn: a _{3} =z _{1} +z _{3} -z _{2} -z _{4}]	(31)

[eqn: X`= {(lambda_P -lambda_1 )} over {(lambda_2 -lambda_1 )}]	(32)

[eqn: Y`= {(phi_P -phi_1 )} over {(phi_4 -phi_1 )}]	(33)

내부점 P에서의 왜곡량 계산법

MATLAB을 이용한 왜곡량 모델링 예시

테스트 지역 및 자료

시범지역별 공통점의 왜곡량 분석(의왕시)

측점명	Bessel 타원체(A)		GRS80타원체(B)		변환 GRS80타원체(C)		차이(B-C)
	X(North)	Y(East)	X(North)	Y(East)	X(North)	Y(East)	□X	□Y
L5	432268.79	202761.21	432574.59	202831.56	432573.89	202831.11	0.696	0.457
448	427722.12	205115.78	428028.05	205186.25	428027.30	205185.74	0.749	0.508
TL24	427541.69	200627.72	427847.49	200698.25	427846.79	200697.71	0.696	0.541
L65	430409.87	201037.40	430715.64	201107.84	430714.96	201107.34	0.683	0.496
TL29	427308.51	199360.33	427614.28	199430.89	427613.59	199430.33	0.682	0.553
TL28	426641.01	199208.91	426946.78	199279.49	426946.10	199278.92	0.684	0.564
L60	430484.73	199812.69	430790.47	199883.15	430789.80	199882.64	0.670	0.505
tl-16	430474.69	199588.40	430780.42	199658.86	430779.75	199658.35	0.665	0.506
T15	430470.25	201322.73	430776.03	201393.17	430775.34	201392.67	0.687	0.494
L62	429621.77	200879.64	429927.55	200950.11	429926.86	200949.60	0.689	0.510
L52	431218.55	200279.25	431524.29	200349.67	431523.62	200349.18	0.673	0.489
L54	431231.50	200912.65	431537.25	200983.06	431536.58	200982.58	0.675	0.483
TL14	431781.75	202270.23	432087.54	202340.61	432086.84	202340.14	0.693	0.467
TL9	432182.17	199817.82	432487.89	199888.22	432487.23	199887.74	0.660	0.478
L68	432034.09	197135.30	432339.73	197205.74	432339.10	197205.24	0.627	0.499
L69	431342.54	197349.38	431648.19	197419.84	431647.56	197419.34	0.631	0.506
TL46	431303.98	196849.33	431609.62	196919.80	431608.99	196919.29	0.629	0.512
TL44	430629.99	197781.61	430935.67	197852.09	430935.02	197851.57	0.643	0.516
TL33	422936.19	194457.29	423241.88	194528.05	423241.23	194527.40	0.647	0.651
BM32-10-1	421465.78	197866.29	421771.60	197937.05	421770.90	197936.40	0.697	0.649
5238	433456.49	187585.09	433761.81	187655.63	433761.31	187655.09	0.502	0.543

측점명	Bessel 타원체(A)		GRS80타원체(B)		변환 GRS80타원체(C)		차이(B-C)
	X(North)	Y(East)	X(North)	Y(East)	X(North)	Y(East)	□X	□Y
52381570	428511.56	198027.10	428817.27	198097.64	428816.61	198097.10	0.658	0.545
T325	422110.70	193938.89	422416.38	194009.68	422415.73	194009.01	0.645	0.666
TL1	432765.23	200740.20	433070.96	200810.58	433070.29	200810.11	0.668	0.464
TL8	432203.09	200719.14	432508.84	200789.52	432508.16	200789.05	0.671	0.472
L2-1	433046.10	201126.13	433351.84	201196.49	433351.17	201196.03	0.671	0.457
TL3	432852.81	201853.42	433158.58	201923.77	433157.89	201923.32	0.682	0.455
TL4	432974.24	202795.68	433280.02	202866.01	433279.33	202865.57	0.690	0.446
L50	432152.75	201114.79	432458.50	201185.17	432457.82	201184.70	0.676	0.470
L49	432248.07	201574.18	432553.84	201644.56	432553.16	201644.09	0.682	0.467
TL26	427024.70	200090.41	427330.49	200160.96	427329.80	200160.41	0.689	0.550
TL23	427930.76	200336.76	428236.51	200407.33	428235.85	200406.75	0.660	0.582
TL35	429710.60	196936.42	430016.27	197006.94	430015.63	197006.40	0.638	0.534
TL36	428936.51	196926.03	429242.19	196996.57	429241.54	196996.03	0.645	0.547
TL40	428865.99	196642.48	429171.66	196713.02	429171.02	196712.48	0.641	0.549
L67	429920.57	196554.49	430226.21	196625.00	430225.58	196624.47	0.628	0.533
TL38	427881.64	196273.17	428187.31	196343.76	428186.67	196343.19	0.638	0.565
BM26-2B	427264.36	200077.78	427570.15	200148.33	427569.46	200147.78	0.691	0.548
SW462	424920.57	202120.18	425226.45	202190.77	425225.73	202190.20	0.728	0.570
6220	434064.51	198768.61	434370.16	198838.97	434369.53	198838.51	0.637	0.462
310	434064.51	198768.62	434370.16	198838.98	434369.53	198838.52	0.635	0.458
AY310	434064.52	198768.62	434370.17	198838.98	434369.53	198838.52	0.635	0.458
SY325	422110.70	193938.89	422416.38	194009.68	422415.74	194009.02	0.645	0.666
311	428511.56	198027.08	428817.26	198097.62	428816.60	198097.08	0.658	0.544
L66	427285.72	199681.14	427591.50	199751.69	427590.81	199751.14	0.684	0.549
TL27	426905.62	199789.80	427211.41	199860.36	427210.72	199859.80	0.689	0.556
l-61	429726.19	200192.74	430031.95	200263.21	430031.27	200262.69	0.679	0.513
L15	430467.51	201355.05	430773.29	201425.49	430772.60	201424.99	0.687	0.494
TL19	429792.66	201235.23	430098.45	201305.69	430097.76	201305.18	0.691	0.504
L53	431328.25	200642.07	431634.00	200712.49	431633.33	200712.00	0.676	0.486
TL13	431543.29	201544.12	431849.06	201614.52	431848.38	201614.04	0.683	0.475
TL12	431380.01	201329.88	431685.78	201400.28	431685.10	201399.80	0.683	0.479
TL10	431175.49	199437.50	431481.21	199507.93	431480.54	199507.44	0.664	0.491
TL42	432116.68	197431.95	432422.32	197502.38	432421.69	197501.89	0.631	0.496
TL43	431275.51	197482.66	431581.16	197553.12	431580.53	197552.62	0.633	0.508
TL45	430475.39	196881.52	430781.04	196952.01	430780.40	196951.49	0.633	0.524
BM26-2A	427048.60	197198.97	427354.31	197269.57	427353.65	197269.00	0.657	0.572
BM26-3-1A	432092.88	194281.57	432398.42	194352.05	432397.83	194351.53	0.592	0.518
AY03	434862.05	197823.18	435167.66	197893.53	435167.04	197893.07	0.620	0.453
L64	427266.57	200068.61	427572.36	200139.16	427571.67	200138.61	0.690	0.548
TL25	427243.56	200340.26	427549.36	200410.81	427548.66	200410.26	0.697	0.548
TL17	429605.39	199804.49	429911.13	199874.97	429910.46	199874.46	0.674	0.517
L59	430488.35	200221.56	430794.10	200292.00	430793.43	200291.50	0.673	0.501
L58	430513.98	200592.01	430819.74	200662.46	430819.06	200661.96	0.678	0.499
TL18	429725.91	200489.77	430031.68	200560.23	430030.99	200559.72	0.684	0.511
TL11	431272.14	199849.01	431577.87	199919.43	431577.20	199918.94	0.666	0.491
L55	431365.73	201091.76	431671.56	201162.20	431670.81	201161.69	0.746	0.517
L56	431731.20	202036.52	432036.99	202106.90	432036.30	202106.43	0.690	0.470
L51	432201.47	200252.5	432507.20	200322.89	432506.53	200322.42	0.665	0.475
TL41	432223.01	196650.61	432528.63	196721.06	432528.01	196720.56	0.621	0.501
L70	431351.46	197192.06	431657.10	197262.52	431656.47	197262.01	0.630	0.507
L71	430462.49	197395.08	430768.16	197465.57	430767.52	197465.05	0.640	0.521
TL31	423879.30	194989.36	424184.99	195060.09	424184.34	195059.46	0.648	0.633
TL30	423833.81	194524.03	424139.48	194594.77	424138.84	194594.13	0.643	0.636
B302	437496.00	201094.88	437801.68	201165.10	437801.03	201164.71	0.644	0.392
L48	432816.54	202203.85	433122.32	202274.19	433121.63	202273.74	0.686	0.453
TL2	432751.47	201342.38	433057.22	201412.75	433056.55	201412.29	0.677	0.461
TL6	432257.88	202234.26	432563.66	202304.62	432562.97	202304.16	0.688	0.461
T5	432256.75	202763.34	432562.54	202833.70	432561.85	202833.24	0.696	0.458
TL21	427776.80	200172.91	428082.58	200243.44	428081.89	200242.90	0.691	0.540
TL20	428001.30	199458.02	428307.06	199528.55	428306.38	199528.01	0.680	0.543
L63	427915.63	199858.64	428221.40	199929.17	428220.72	199928.63	0.683	0.540
TL34	429835.16	196200.78	430140.79	196271.30	430140.17	196270.76	0.626	0.535
TL39	428934.72	196438.66	429240.38	196509.21	429239.74	196508.66	0.637	0.550
TL37	427881.26	196627.53	428186.94	196698.12	428186.29	196697.55	0.647	0.568
TL32	423073.66	194765.99	423379.35	194836.74	423378.70	194836.09	0.650	0.647
						평균	0.665	0.518
						표준편차	0.033	0.054

MATLAB 사용 S/W 명령어

1) 데이터 로딩 및 도시

>> load eui_exam.txt

>> data = eui_exam;

>> plot(data(:,7),data(:,6),'+');	% data를 화면에 디스플레이

>> num_data = size(data,1)

2) 데이터의 통계

>> mean_dx = mean(data(:,8))	%X방향의 왜곡량 평균

>> mean_dy = mean(data(:,9))	%Y방향의 왜곡량 평균

>> std_dx = std(data(:,8))	%X방향의 왜곡량 표준편차

>> std_dy = std(data(:,9))	%Y방향의 왜곡량 표준편차

>> max_dx = max(data(:,8))	%X방향의 왜곡량 최고값

>> max_dy = max(data(:,9))	%Y방향의 왜곡량 최고값

>> min_dx = min(data(:,8))	%X방향의 왜곡량 최저값

>> min_dy = min(data(:,9))	%Y방향의 왜곡량 최저값

	X	Y
평균	0.666	0.518
표준편차	0.033	0.054
최대치	0.749	0.666
최소치	0.502	0.392

3) 비상사점 찾기

>> tmp = find(abs(data(:,8)-mean_dx)>3*std_dx | abs(data(:,9)-mean_dy)>3*std_dy ); %비상사점 찾기

>> data(tmp,:) = [];	% 비상사인점 제거

>> hold on; plot(data(:,7),data(:,6),'og');	% data를 화면에 디스플레이

>> num_data = size(data,1)

□ 결과는 비상사인 점을 하나 발견. 측점명 5238

4) 편중된 공통점 제거

>> num_data = size(data,1)

>> k=1;	% 모든 점들 사이의 거리를 구한다.

for i=1:num_data,

for j=i+1:num_data,

h(k,1)=i;

h(k,2)=j;

h(k,3)=sqrt((data(i,6)-data(j,6))^2+(data(i,7)-data(j,7))^2);

k=k+1;

end

end

>> tmp = find(h(:,3)<100);	% 100m이하인 곳을 찾는다. m결정 필요

>> h_tmp = h(tmp,:);	% 100m이하인 h만 선택한다. 이것은 결정해야 할 것임

>> data_tmp = zeros(num_data,1);

>> for i=1:size(h_tmp,1),

data_tmp(h_tmp(i,1))=data_tmp(h_tmp(i,1))+1;

%data_tmp(h_tmp(i,2))=data_tmp(h_tmp(i,2))+1;

end

>> max_pt = max(data_tmp);	% max_pt가 1이면... 2칸 아래로..

>> tmp = find(data_tmp ==max_pt);

>> data(tmp,:) = [];	% 편중된 공통점중 가장 중복이 많은 점을 제거.. 그리고.. 다시 편중된 공통점 제거 처음으로..

>> num_data = size(data,1)	% 비상사점이 제거된 후의 점의 개수

>> data(h_tmp(:,2),:) = [];	% max_pt가 1인 경우에 이용됨. 마지막으로 1대1로 연결된 점중 아이디가 뒤에 있는 점 제거

>> num_data = size(data,1)	% 편중된 공통점이 제거된 후의 점의 개수

>> hold on; plot(data(:,7),data(:,6),'+r');	% data를 화면에 디스플레이

□ 결과적으로 78개의 점이 남았다.

5) 회귀식에 의한 지역적 변화제거

>> odd_line = 1:2:2*num_data-1;

>> even_line = 2:2:2*num_data;

>> mean_X = mean(data(:,6));	%

>> mean_Y = mean(data(:,7));

>> B = zeros(2*num_data,6);	% B행렬의 결정

>> B(odd_line,1)=1;

>> B(odd_line,2)=data(:,6)-mean_X;

>> B(odd_line,3)=data(:,7)-mean_Y;

>> B(even_line,4)=1;

>> B(even_line,5)=data(:,6)-mean_X;

>> B(even_line,6)=data(:,7)-mean_Y;

>> y = zeros(2*num_data,1);	% y벡터 결정

>> y(odd_line) = data(:,8);

>> y(even_line) = data(:,9);

>> x = inv(B'*B)*(B'*y)	% x벡터 결정

[eqn: a_0 = 0.666551282051284], [eqn: b_0 = -5.5362088287528e-006], [eqn: c_0 = 1.24362363073061e-005]

[eqn: a_1 = 0.517987179487181], [eqn: b_1 = -1.48187166441985e-005], [eqn: c_1 = -6.48804224419043e-006]

>> a0 = x(1); % 데이터 저장

>> b0 = x(2);

>> c0 = x(3);

>> a1 = x(4);

>> b1 = x(5);

>> c1 = x(6);

>> data(:,10) = data(:,8) - (a0+b0*(data(:,6)-mean_X)+c0*(data(:,7)-mean_Y));

>> data(:,11) = data(:,9) - (a1+b1*(data(:,6)-mean_X)+c1*(data(:,7)-mean_Y));

>> figure(2);

>> subplot(2,1,1), plot(data(:,10)), xlabel('position'), ylabel('X direction residual distortion(m)');

>> subplot(2,1,2), plot(data(:,11)), xlabel('position'), ylabel('Ydirection residual distortion(m)');

다음은 X방향과 Y방향으로의 잔여왜곡량을 나타낸다.

□ 여기서, 두 개의 점이 문제(비상사점)가 있다. 위치 31과 60번째에 있는 점이 다른 점들에 비해서 매우 높은 오차를 지니고 있다. 따라서, 오차가 큰 점에 대해서는 과대오차로 간주하여 제거한다

>> data([31 60],:) = [];

>> figure(3);

>> subplot(2,1,1), plot(data(:,10)), xlabel('position'), ylabel('X direction residual distortion(m)');

>> subplot(2,1,2), plot(data(:,11)), xlabel('position'), ylabel('Ydirection residual distortion(m)');

□ 제거된 데이터를 이용하여 다시 2차원 회귀식에 적용하면,

>> num_data = size(data,1)

>> odd_line = 1:2:2*num_data-1;

>> even_line = 2:2:2*num_data;

>> mean_X = mean(data(:,6));	%

>> mean_Y = mean(data(:,7));

>> B = zeros(2*num_data,6);	% B행렬의 결정

>> B(odd_line,1)=1;

>> B(odd_line,2)=data(:,6)-mean_X;

>> B(odd_line,3)=data(:,7)-mean_Y;

>> B(even_line,4)=1;

>> B(even_line,5)=data(:,6)-mean_X;

>> B(even_line,6)=data(:,7)-mean_Y;

>> y = zeros(2*num_data,1);	% y벡터 결정

>> y(odd_line) = data(:,8);

>> y(even_line) = data(:,9);

>> x = inv(B'*B)*(B'*y)	% x벡터 결정

>> a0 = x(1); % 데이터 저장

>> b0 = x(2);

>> c0 = x(3);

>> a1 = x(4);

>> b1 = x(5);

>> c1 = x(6);

>> data(:,10) = data(:,8) - (a0+b0*(data(:,6)-mean_X)+c0*(data(:,7)-mean_Y));

>> data(:,11) = data(:,9) - (a1+b1*(data(:,6)-mean_X)+c1*(data(:,7)-mean_Y));

>> figure(4);

>> subplot(2,1,1), plot(data(:,10)), xlabel('position'), ylabel('X direction residual distortion(m)');

>> subplot(2,1,2), plot(data(:,11)), xlabel('position'), ylabel('Ydirection residual distortion(m)');

6) 해석학적 공분산함수 결정

>> num_data = size(data,1)	% 편중된 공통점이 제거된 후의 점의 개수

>> mean_delX = mean(data(:,10));	% X방향의 잔여왜곡량의 평균 계산

>> mean_delY = mean(data(:,11));	% Y방향의 잔여왜곡량의 평균 계산

>> h = zeros(num_data,3);

>> k=1;	% 실험적 공분산을 구한다.

for i=1:num_data,

for j=i+1:num_data,

h(k,1)=sqrt((data(i,6)-data(j,6))^2+(data(i,7)-data(j,7))^2); % 거리 r의 계산

h(k,2)=(data(i,10)-mean_delX)*(data(j,10)-mean_delX);

h(k,3)=(data(i,11)-mean_delY)*(data(j,11)-mean_delY);

k=k+1;

end

end

>> dist = 200; num_dist = (2000/dist); % 간격설정 필요

>> histo = zeros(num_dist,3);

for i=1:num_dist,	% 간격설정 필요

tmp1 = find(h(:,1)>=(i-1)*dist & h(:,1) < i*dist);

histo(i,1) = dist*(i-1)+dist/2;

histo(i,2) = mean(h(tmp1,2));

histo(i,3) = mean(h(tmp1,3));

end

>> tmp = find(histo(:,1) < 100);

>> histo(tmp,:) = [];

>> figure(5);

>> subplot(2,1,1), plot(histo(:,1),histo(:,2),'o'), xlabel('distance (m)'), ylabel('X direction empirical covariance(m^2)');

>> subplot(2,1,2), plot(histo(:,1),histo(:,3),'o'), xlabel('distance (m)'), ylabel('Y direction empirical covariance(m^2)');

□ 위의 데이터를 fitting하여야 한다. fitting과정은 다음과 같은 과정을 통하여 계산하였다.

>> clear B x y

>> x = curvefit('gaussfn', [exp(-6) 0.003], histo(:,1), histo(:,2));

>> C0x = x(1);	% X방향의 잔류 왜곡량의 분산

>> Ax = x(2);	% X방향의 잔류 왜곡량의 상관거리

>> sx = 0:10:2000;

>> Csx = C0x*exp(-Ax^2*sx.^2);

>> figure(5); subplot(2,1,1), hold on, plot(sx,Csx,'r');

>> y = curvefit('gaussfn', [exp(-6) 0.003], histo(:,1), histo(:,3));

>> C0y = y(1);	% X방향의 잔류 왜곡량의 분산

>> Ay = y(2);	% X방향의 잔류 왜곡량의 상관거리

>> sy = 0:10:2000;

>> Csy = C0y*exp(-Ay^2*sy.^2);

>> figure(5); subplot(2,1,2), hold on, plot(sy,Csy,'r');

□ reilly의 공분산 함수, 여기서는 초기치가 매우 중요한 영향을 준다

>>  x = curvefit('reilly', [exp(-6) 0.003], histo(:,1), histo(:,2));

>> C0xr = x(1);	% X방향의 잔류 왜곡량의 분산

>> Axr = x(2);	% X방향의 잔류 왜곡량의 상관거리

>> sxr = 0:10:2000;

>> Csxr = C0xr*(1-Axr^2*sxr.^2).*exp(-Axr^2*sxr.^2);

>> figure(5); subplot(2,1,1), hold on, plot(sxr,Csxr,'g');

>> y = curvefit('reilly', [exp(-6) 0.003], histo(:,1), histo(:,3));

>> C0yr = y(1);	% X방향의 잔류 왜곡량의 분산

>> Ayr = y(2);	% X방향의 잔류 왜곡량의 상관거리

>> syr = 0:10:2000;

>> Csyr = C0yr*(1-Ayr^2*syr.^2).*exp(-Ayr^2*syr.^2);

>> figure(5); subplot(2,1,2), hold on, plot(syr,Csyr,'g');

7) 최소 제곱 콜로케이션에 의한 격자화

>> Rp(1) = mean(data(:,6));	% 내삽하고자 하는 X좌표

>> Rp(2) = mean(data(:,7));	% 내삽하고자 하는 Y좌표

>> clear h;

>> h=zeros(num_data,2);

>> h(:,1) = (1:1:num_data)';

>> h(:,2) = sqrt((data(:,6)-Rp(1)).^2+(data(:,7)-Rp(2)).^2);

>> h_tmp = find(h(:,2)<5000); % 5Km내의 점을 찾는다.

>> h(h_tmp,:) = [];

>> if size(h,1)>5,	% 5점 이상이라면,

h_tmp = sortrows(h,2);	% quicksort

h = h_tmp(1:5,:);

end

>> CIx = gaussfn([C0x Ax],h(:,2))';

>> clear tmp;

>> h_size = size(h,1);

tmp = zeros(h_size,h_size);

for i=1:h_size,

for j=1:h_size,

tmp(i,j)=sqrt((data(h(i,1),6)-data(h(j,1),6))^2+(data(h(i,1),7)-data(h(j,1),7))^2);

end

end

>> CDx = gaussfn([C0x Ax],tmp);

>> dsigx = data(h(:,1)',10);

>> Rp(3) = CIx*inv(CDx)*dsigx;

>> CIy = gaussfn([C0y Ay],h(:,2))';

>> CDy = gaussfn([C0y Ay],tmp);

>> dsigy = data(h(:,1)',11);

>> Rp(4) = CIy*inv(CDy)*dsigy;

>> Rp(5) = Rp(3) + (a0+b0*(Rp(1)-mean_X)+c0*(Rp(2)-mean_Y));

>> Rp(6) = Rp(4) + (a1+b1*(Rp(1)-mean_X)+c1*(Rp(2)-mean_Y));

□ 결국 우리가 원하는 점에서의 왜곡량은 위와 같이 결정된다.

□ 이러한 과정이 1/1,000수치지형도 도곽에 대하여 수행된다

<별첨 5>

세계측지계 변환 비용산정 예시

1. 설계제원

① 1/1,000 수치지형도 시가지 100도엽(약 25□)

② 정확도 평가를 위한 GPS 현지측량 50점

③ 1/1,000 수치지형도와 연계된 GIS 데이터 100도엽

2. 작업인원 수량산출

(단위 : 인)

구 분	특급 기술자	고급 기술자	중급 기술자	초급 기술자	정보 처리 기사	중급 기능사	초급 기능사	보통 인부	비 고
1/1,000수치지형도 좌표변환	12.0	20.0		8.0	4.0	40.0			( )은 내업
GPS 현지측량 (50점)		12.34 (4.56)	12.34 (6.44)	12.34 (7.51)			21.46	3.22
GIS 응용 데이터 좌표변환		8.0		16.0		48.0
계	12.0	44.90	18.78	43.85	4.0	88.0	21.46	3.22

※ 중급기능사는 지도제작분야, 초급기능사는 측량분야를 적용

① 1/1,000 수치지형도 좌표변환

- 특급기술자	: 0.15 □ 100도엽 □ 0.8 = 12.0인

- 고급기술자	: 0.25 □ 100도엽 □ 0.8 = 20.0인

- 초급기술자	:  0.1 □ 100도엽 □ 0.8 =  8.0인

- 정보처리기사	: 0.05 □ 100도엽 □ 0.8 =  4.0인

- 중급기능사(지도제작)	:  0.5 □ 100도엽 □ 0.8 = 40.0인

② GPS 현지측량 <건설표준품셈 21-3-4  4급기준점 측량 적용>

- 고급기술자(외업	) : 23.0 □ 1.15 □ (50□150) □ 1.4 = 12.34인

(내업	) :  8.5 □ 1.15 □ (50□150) □ 1.4 =  4.56인

- 중급기술자(외업	) : 23.0 □ 1.15 □ (50□150) □ 1.4 = 12.34인

(내업	) : 12.0 □ 1.15 □ (50□150) □ 1.4 =  6.44인

- 초급기술자(외업	) : 23.0 □ 1.15 □ (50□150) □ 1.4 = 12.34인

(내업	) : 14.0 □ 1.15 □ (50□150) □ 1.4 =  7.51인

- 초급기능사(외업	) : 40.0 □ 1.15 □ (50□150) □ 1.4 = 21.46인

- 보통인부  (외업	) :  6.0 □ 1.15 □ (50□150) □ 1.4 =  3.22인

③ GIS 응용데이터 좌표변환

- 고급기술자	: 0.1 □ 100도엽 □ 0.8 =  8.0인

- 초급기술자	: 0.2 □ 100도엽 □ 0.8 = 16.0인

- 중급기술자(지도제작)	: 0.6 □ 100도엽 □ 0.8 = 48.0인

3. 사업비 산출

(단위 : 원)

구 분	단 가	수 량	금 액	비 고
① 직접인건비 - 특급기술자 - 고급기술자 - 중급기술자 - 초급기술자 - 정보처리기사 - 중급기능사 - 초급기능사 - 보통인부	165,204 123,159 104,134 90,250 89,202 91,098 61,567 52,585	12.0 44.90 18.78 43.85 4.0 88.0 21.46 3.22	23,289,367 1,982,448 5,529,839 1,955,636 3,957,462 356,808 8,016,624 1,321,227 169,323
② 직접경비(여비) - 고급기술자 - 중급기술자 - 초급기술자 - 초급기능사	43,000 37,000 37,000 37,000	12.34 12.34 12.34 21.46	2,237,800 530,620 456,580 456,580 794,020
③ 제경비 (직접인건비□110%)			25,618,303
④ 기술료 (직접인건비+제경비)□20%			9,781,534
⑤ 소계			60,927,004
⑥ 부가가치세			6,092,700
⑦ 합 계			67,019,704

<첨부 1>

세계측지계 변환에 필요한 투영식

타원체면상의 경위도 좌표를 평면직각좌표로 변환

[eqn: Y(E)~=~□Y+k_0 ·N· left[A+ A^3 over6 (1-T+C)+ A^5 over 120 (5-18T + T^2 +72C - 58e'^2 ) right]]

[eqn: X(N)~=~□X+k_0 · left{&M-M_0 + N tan□· left( {A^2} over{2} + {A^4} over24 (5-T+9C+4C^2 ) ##&+ {A^6} over{720} (61-58T+T^2 +600C -330e'^2 ) right) right}]

여기서,

[eqn: &①~T=tan^2 □ ## &②~C= e^2 over{1-e^2 }cos^2 □ ## &③~A=(□-□_0 )cos□ (여기서, ~□와~ □_0 는 ~radian값임) ## &④~N(위도~□에서의 ~묘유선의~곡률~반경) = a over sqrt{1-e^2 sin^2 □} ## &⑤~M(자오선장)=a·left{left(1- {e^2} over {4} - {3e^4 }over {64} - {5e^6 }over {256} right)□ - left( {3e^2 } over{8} + {3e^4 }over {32} + {45e^6 }over {1024} right)sin2□ ##&~~~~~~~~~~~~~~+ left( {15e^4 } over{256} + {45e^6 }over{1024} right)sin4□ - {35e^6} over {3072} sin6□ right} ## &⑥~e^2 (제1이심률) = {a^2 - b^2} over a^2 ## &⑦~e'^2 (제2이심률) = {a^2 - b^2} over b^2]

그리고,

[eqn: &□ ~:~ 위도,~~□~:~경도,~~□_0 ~:~ 투영원점~위도,~~□_0 ~:~투영원점~경도, ##&a~:~타원체~장반경, ~~f ~:~ 편평률, ##&b~:~타원체~단반경(=a(1-f)), ##&k_0 ~:~원점축척계수, ##&rm□Y ~:~ Y축(East) ~원점~가산값, ##&rm□X~:~X축(North)~원점~가산값]

평면직각좌표를 타원체면상의 경위도 좌표로 변환

[eqn: □=□_1 - {N_1 tan□_1} overR_1 · left[&D^2 over2 - D^4 over24 (5+3T_1 +10C_1 -4C_1 ^2 -9e'^2 ) ##&+ D^6 over 720 (61+90T_1 + 298C_1 +45T_1 ^2 -252e'^2 - 3C_1 ^2 ) right]]

[eqn: □=□_0 + 1over cos□_1 &left[D- D^3 over6 (1+2T_1 +C_1 ) + D^5 over 120 (5-2C_1 + 28T_1 -3C_1 ^2 + 8e'^2 + 24T_1 ^2 ) right]]

여기서,

[eqn: &①~M = M_0 + (X-□X) over k_0 ##&②~□_1 = {M} over{a·left(1- e^2 over4 - 3e^4 over64 - 5e^6 over256 right) } ###&③~e_1 = {1-sqrt{1-e^2} }over {1+sqrt{1-e^2} } ###&④~□_1 = □_1 + ( 3e_1 over2 - 27e_1 ^3 over32 )sin2□_1 + ( 21e_1 ^2 over16 - 55e_1 ^4 over32 )sin4□_1 ##&~~~~+ ( 151e_1 ^3 over96 )sin6□_1 + ( 1097e_1 ^4 over512 )sin8□_1 ###&⑤~ R_1 (위도~□_1 에서의 ~자오선의~곡률~반경) ~=~ {a·(1-e^2 )}over {(1-e^2 sin^2 □_1 )^{3/2}} ## &⑥~ C_1 = e'^2 cos^2 □_1 ## &⑦~T_1 = tan^2 □_1 ##&⑧~N_1 (위도~□_1 에서의 ~묘유선의~곡률~반경) ~= ~a over sqrt{1-e^2 sin^2 □_1} ## & ⑨~D={Y-□Y}over{N_1 ·k_0 }]

그리고, [eqn: rmY ~:~ TM 좌표(East),~~ X~:~TM좌표(North)]

[eqn: &~ M_0 ~:~ 투영원점에~대한~자오선호장 , &k_0 ~:~원점축척계수]