MATLAB 고 스 크 뤼 그 투영 반산 실현

가우스 투영 (가우스 - 크 레 그 투영) 의 반산
2020 - 06 을 업데이트 하고 좌표 시스템 을 WGS - 84 좌표계 로 통일 하여 스 크 립 트 함 수 를 정리 합 니 다.
가우스 투영 의 반산 은 현지의 국부 좌표계 (x, y) 에서 현지의 지리 좌표계 (B: 위도, L: 경도) 로 전환 하 는 것 을 말한다.기 존의 보 문 MATLAB 가 고 스 크 뤼 그 투영 을 실현 한 것 은 고 스 투영 에 대해 간략하게 설명 한 것 이기 때문에 본 보 문 은 고 스 크 뤼 그 투영 의 원 리 를 소개 하지 않 고 고 고 스 투영 반산 의 알고리즘 절차 와 실 현 된 MATLAB 스 크 립 트 만 제시 했다.본 박문 참고 문헌 자 료 는 다음 과 같다.
[1] 공상 원, 곽 국제 명, 유종 천 등. 대지 측량 학 기초 (제2판) [M]. 무한: 무한대학 출판사, 2009. [2] 정 붕 비, 성 영연, 문 한강 등. 2000 국가 대지 좌표계 실 용 보전 [M]. 북경: 측량 출판사, 2008. [3] 우 리 연. 측량 좌표 전환 모델 연구 와 전환 시스템 실현 [D]. 장안대학, 2010. [4] 이 후 박,변 소 봉 고 스 투영 의 복 변 함수 표시 [J]. 측량 학보, 2008 (01): 5 - 9
가우스 투영 반산 공식
현지 위도 B B B 의 계산 공식 은 다음 과 같다.
B = B f − ρ t f 2 M f y ( y N f ) [ 1 − 1 12 ( 5 + 3 t f 2 + η f 2 − 9 η f 2 t f 2 ) ( y N f ) 2 + 1 360 ( 61 + 90 t f 2 + 45 t f 4 ) ( y N f ) 4 ] \begin{aligned} B=& B_{f}-\frac{\rho t_{f}}{2 M_{f}} y\left(\frac{y}{N_{f}}\right)\left[1-\frac{1}{12}\left(5+3 t_{f}^{2}+\eta_{f}^{2}-9 \eta_{f}^{2} t_{f}^{2}\right)\left(\frac{y}{N_{f}}\right)^{2}\right.\\ &+\left.\frac{1}{360}\left(61+90 t_{f}^{2}+45 t_{f}^{4}\right)\left(\frac{y}{N_{f}}\right)^{4}\right] \end{aligned} B=​Bf​−2Mf​ρtf​​y(Nf​y​)[1−121​(5+3tf2​+ηf2​−9ηf2​tf2​)(Nf​y​)2+3601​(61+90tf2​+45tf4​)(Nf​y​)4]​
y = 0 y = 0 y = 0 시의 위 도 를 저 점 위도 B f B 라 고 한다f. BF, 밑 점 위도 의 교체 계산 공식 은 다음 과 같다.
X = a 0 B − a 2 2 sin ⁡ 2 B + a 4 4 sin ⁡ 4 B − a 6 6 sin ⁡ 6 B + a 8 8 sin ⁡ 8 B X=a_{0} B - \ frac {a {2}} {2} \ sin 2 B + \ frac {a {4}} {4} \ sin 4 B - \ frac {a {6}} {6} \ sin 6 B + \ frac {a {8} {8} \ sin 8 B X = a0 B − 2a2 sin2B + 4a4 sin4B − 6a6 sin6B + 8a8 sin8B 중,{ a 0 = m 0 + 1 2 m 2 + 3 8 m 4 + 5 16 m 6 + 35 128 m 8 a 2 = 1 2 m 2 + 1 2 m 4 + 15 32 m 6 + 7 16 m 8 a 4 = 1 8 m 4 + 3 16 m 6 + 7 32 m 8 a 6 = 1 32 m 6 + 1 16 m 8 a 8 = 1 128 m 8 ; { m 0 = a ( 1 − e 2 ) m 2 = 3 2 e 2 m 0 m 4 = 5 4 e 2 m 2 m 6 = 7 6 e 2 m 4 m 8 = 9 8 e 2 m 6 \left\{\begin{aligned} a_{0}=m_{0}+\frac{1}{2} m_{2}+\frac{3}{8} m_{4}+\frac{5}{16} m_{6}+\frac{35}{128} m_{8} \\ a_{2}=\frac{1}{2} m_{2}+\frac{1}{2} m_{4}+\frac{15}{32} m_{6}+\frac{7}{16} m_{8} \\ a_{4}=\frac{1}{8} m_{4}+\frac{3}{16} m_{6}+\frac{7}{32} m_{8} \\ a_{6}=\frac{1}{32} m_{6}+\frac{1}{16} m_{8} \\ a_{8}=\frac{1}{128} m_{8} \end{aligned};\quad \right. \left\{\begin{aligned} m_{0}=a(1-e^2) \\ m_{2}=\frac{3}{2} e^2 m_{0} \\ m_{4}=\frac{5}{4} e^2 m_{2}\\ m_{6}=\frac{7}{6} e^2 m_{4} \\ m_{8}=\frac{9}{8} e^2 m_{6} \end{aligned}\right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧​a0​=m0​+21​m2​+83​m4​+165​m6​+12835​m8​a2​=21​m2​+21​m4​+3215​m6​+167​m8​a4​=81​m4​+163​m6​+327​m8​a6​=321​m6​+161​m8​a8​=1281​m8​​;⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎧ m0 = a (1 − e2) m2 = 23 e2m0 m4 = 45 e2m2 m6 = 67 e2m4 m8 = 89 e2m6 B f B B f Bf 교체 계산 과정 은: 교체 초기: 령 B f 0 = x a 0 B f ^ 0 = \ \ frac {x}B f0 = a 0 x 의 교체 과정: F i = = − a 2 의 죄 가 있 었 고, 2 B 의 f i + a 4 의 죄 가 있 었 고, 4 B 의 f i 와 6 의 죄 가 있 었 고, 6 의 죄 가 있 었 고, 6 B 의 f i + a 8 의 죄 가 있 었 고, 8 B f f f i ^ {i} - \ frac {a {2} {2}} {2} \ \ sin 2 B {f}}} \ \ \ sin 2 B {f}}}}} {{{i}}}} \ \ \ \ \ \ \ \ frac {a {a}}}} {{4} \ \ \ \ \ \ \ \ \ \ \ \ frac {a 4 B {f}}}}}}} \ \ \ \ \ \ \ \ \ \ \ \ sin 6 B {f} ^ {i} + \ frac {a {8}} {8} \ sin 8 B {f} ^ {i}Fi = − 2a2 sin2Bfi + 4a4 sin4Bfi − 6a6 sin6Bfi + 8a8 sin8Bfi; 령 B f i + 1 = X − F i a 0 B f ^ {i + 1} = \ \ frac {X - F ^ i} {a 0} Bfi + 1 = a0 X − Fi 교체 정지:ϵ |B_f^{i+1} - B_f^i|∣Bfi+1​−Bfi​∣<ϵ시 교체 중지. 주:ϵ \epsilon ϵ주어진 한도 값 을 위해 한도 값 은 보통 1E - 8 보다 작 습 니 다.
현지 경도 l l 의 계산 공식 은 다음 과 같다.
l = ρ cos ⁡ B f ( y N f ) [ 1 − 1 6 ( 1 + 2 t f 2 + η f 2 ) ( y N f ) 2 + 1 120 ( 5 + 28 t f 2 + 24 t f 4 + 6 η f 2 + 8 η f 2 t f 2 ) ( y N f ) 4 ] \begin{aligned} l=& \frac{\rho}{\cos B_{f}}\left(\frac{y}{N_{f}}\right)\left[1-\frac{1}{6}\left(1+2 t_{f}^{2}+\eta_{f}^{2}\right)\left(\frac{y}{N_{f}}\right)^{2}+\frac{1}{120}\left(5+28 t_{f}^{2}\right.\right.\\ &+\left.\left.24 t_{f}^{4}+6 \eta_{f}^{2}+8 \eta_{f}^{2} t_{f}^{2}\right)\left(\frac{y}{N_{f}}\right)^{4}\right] \end{aligned} l=​cosBf​ρ​(Nf​y​)[1−61​(1+2tf2​+ηf2​)(Nf​y​)2+1201​(5+28tf2​+24tf4​+6ηf2​+8ηf2​tf2​)(Nf​y​)4]​
상기 식 중:

x x x x 는 북방 방향의 좌표 이 고

y y y 는 동방 방향의 좌표 (중국 지역 내수 에서 500000 m 감소),

ρ = 180 × 3600 / π \rho=180×3600/\pi ρ=180×3600 / pi 는 라디안 초,

a a a 는 지구 타원체 의 긴 반지름 이다.

e e e 는 지구 타원 공의 첫 번 째 편심 율 이다.

e ′ e ′ e 는 지구 타원체 의 두 번 째 편 심 률 이다.

N f N f Nf 는 저 점 위도 B f B f Bf 에서 지구의 묘 유 권 곡률 반지름 이다.

M f M f Mf 는 저 점 위도 B f B f Bf 에서 지구의 자 오 권 곡률 반지름 이다.

η f 2 = e ′ 2 c o s 2 B f \eta_f^2 = e'^2cos^2{B_f} ηf2​=e′2cos2Bf​，

t f = t a n B f t_f=tan{B_f} tf​=tanBf​。

가우스 투영 반산 의 MATLAB 함수

function Coord = GaussProInverse(Pos)
% INPUT // Units of longitude and latitude is RAD (°)
% REF[1]//    , .2000           [M].  :     ,2008.
% REF[2]//    ,   .       (   )[M].  :       ,2010.
% Longitude of the Earth central meridian
MerLon = 114;          % Wuhan is 114 deg
% Extract the local projected coordinate (x & y)
Coord = Pos;
Eth.D2R = 0.0174532925199433;  % pi/180
x = Pos(1);
y = Pos(2) - 500000;
%% Earth orientation parameters of WGS 84 Coordinate System
Eth.R0 = 6378137.0;
Eth.f = 1/298.257223563;
Eth.Rp = 6356752.3142452;      % R0*(1 - f);
% First Eccentricity and its Squared
Eth.e12 = 0.006694379990141;   % (2f - f*f);
Eth.e11 = 0.081819190842622;   % sqrt(2f - f*f);
% Second Eccentricity and its Squared
Eth.e22 = 0.00673949674227643; % (2f - f*f)/(1 + f*f - 2*f);
Eth.e21 = 0.08209443794969570; % sqrt((2f - f*f)/(1 + f*f - 2*f));
Bf = Meridian2Latitude(x, Eth.R0, Eth.e12);
tnBf = tan(Bf); tn2Bf = tnBf*tnBf; tn4Bf = tn2Bf*tn2Bf;
csBf = cos(Bf); cs2Bf = csBf*csBf; Eta2 = Eth.e22*cs2Bf;
COE = sqrt(1 - Eth.e12*sin(Bf)*sin(Bf));
Nf = Eth.R0/COE;
Mf = Eth.R0*(1 - Eth.e12)/COE^3;
%% Calculate Latitude(B)
YNf = y/Nf;
YNf2 = YNf*YNf;
YNf4 = YNf2*YNf2;
TYMN = 0.5*tnBf*y*YNf/Mf;
COE1 = (5 + 3*tn2Bf + Eta2 -9*Eta2*tn2Bf)*YNf2/12;
COE2 = (61 + 90*tn2Bf + 45*tn4Bf)*YNf4/360;
Lat = Bf - TYMN*(1 - COE1 + COE2);
%% Calculate Longitude(L)
YBNf = YNf/csBf;
COE3 = (1 + 2*tn2Bf + Eta2)*YNf2/6;
COE4 = (5 + 28*tn2Bf + 24*tn4Bf + 6*Eta2 + 8*Eta2*tn2Bf)*YNf4/120;
Lon = MerLon*Eth.D2R + YBNf*(1 - COE3 + COE4);
Coord(1) = Lat/Eth.D2R;
Coord(2) = Lon/Eth.D2R;
end

function Bf = Meridian2Latitude(x,a,e12)
m0 = a*(1 - e12);    m2 = 3*e12*m0/2;
m4 = 5*e12*m2/4;     m6 = 7*e12*m4/6;
m8 = 9*e12*m6/8;     a8 = m8/128;
a6 = m6/32 + m8/16;  a4 = m4/8 + 3*m6/16 + 7*m8/32;
a0 = m0 + m2/2 + 3*m4/8 + 5*m6/16 + 35*m8/128;
a2 = m2/2 + m4/2 + 15*m6/32 + 7*m8/16;
B0 = x/a0;
while 1
    F = -a2*sin(2*B0)/2 + a4*sin(4*B0)/4 - a6*sin(6*B0)/6 + a8*sin(8*B0)/8;
    Bf = (x - F)/a0;
    if abs(B0 - Bf)<1e-10
        break;
    end
    B0 = Bf;
end
end

이 방법 은 전형 적 인 고 스 반산 알고리즘 으로 고 스 반산 에 대한 상세 한 추 도 는 대지 측량 교재 나 최신 문헌 을 참고 할 수 있다.