TRƯỜNG ĐẠI HỌC BÁCH KHOA
THÀNH PHỐ HỒ CHÍ MINH
KHOA KỸ THUẬT XÂY DỰNG
BÁO CÁO
BÁO CÁO BÀI TẬP LỚN TRẮC ĐỊA CAO CẤP 1
GVHD: PGS. TS NGUYỄN NGỌC LÂU
SVTH: TRẦN ĐÔNG GIANG
MSSV: 81100916
Lớp: XD11TD
Tháng 05/2014
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
LỜI CẢM ƠN
Để hoàn thành bài báo cáo này, em xin chân thành cảm ơn Thầy Nguyễn Ngọc
Lâu đã tận tình giảng dạy và hướng dẫn em trong suốt quá trình hoàn thành bài
báo cáo này. Mặc dù rất cố gắng nhưng do kiến thức chuyên môn bản thân còn
hạn chế nên bài báo cáo này không tránh khỏi những thiếu sót. Kính xin Thầy
góp ý để em rút ra những kinh nghiệm quý báu, em xin chân thành cảm ơn./.
Trần Đông Giang.
2
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
BÀI TẬP LỚN SỐ 1
1. YÊU CẦU
Xác định kích thước của mặt Spheroid dựa vào chiều dài 2 cung kinh tuyến và vĩ độ
tương ứng với chiều dài 2 cung kinh tuyến ấy.
Số liệu: bộ số liệu có số thứ tự thứ 14
STT
14
Vĩ độ (0)
1.5 – 2.5
86.5 -87.5
Độ dài KT (m)
110576.904
111693.063
Xác định kích thước của mặt Spheroid tương ứng với bộ số liệu trên:
- a=?
- e2 = ?
2. GIẢI QUYẾT BÀI TOÁN ĐẶT RA:
Để xác định kích thước của mặt Spheroid trái đất ta có 4 phương pháp:
- Phương pháp Eratosthenes.
- Phương pháp Viện Hàn Lâm Khoa Học Pháp.
- Phương pháp sử dụng mạng lưới tam giác quốc gia.
- Phương pháp đo trọng lực.
Khi biết vĩ độ Ф và độ dài cung kinh tuyến S, để xác định kích thước mặt Spheroid
ta có thể tìm hiểu 2 cách giải: cách giải gần đúng và cách giải chính xác.
CÁCH GIẢI GẦN ĐÚNG: dựa trên giá trị xấp xỉ của tích phân.
- Chiều dài cung kinh tuyến quan hệ với vĩ độ:
2
S1 M d ~ M m1 2 (2 1)
1
4
S2 M d ~ M m 43 (4 3)
3
Với bán kính cung kinh tuyến tại vĩ độ trung bình
3
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
M
1 2
m
M m1 2 =
a. 1 e2
2
2 1 2
)
1 e sin (
2
a. 1 e2
2
2 3 4
)
1 e sin (
2
;
3/2
3/2
;
GVHD: PGS.TS NGUYỄN NGỌC LÂU
12
m
1 2
2
34
m
3 4
2
Viết lại:
S1 = a(1-e2).
S2 = a(1-e2).
Ta lấy:
( 2 1 )
1 e sin ( )
2
2
12
m
3/2
( 2 1 )
1 e
2
sin 2 (34
m )
(1)
3/2
3/2
S1 (1 e2 .sin 2 (34
m ))
3/2
S2 (1 e2 .sin 2 (12
m ))
→ giải tìm e2
Thay vào hệ (1) giải ra a.
THỰC HIỆN:
+ Ta có:
12
m
1 2
200 '0 ;
2
12 2 1 10
+ Tương tự ta có:
34
m
4 3
8700 '0 ;
2
34 4 3 10
3/2
S1 (1 e2 .sin 2 (34
m ))
3/2
S2 (1 e2 .sin 2 (12
m ))
110576.904 (1 e .sin (87 ))
=
2
2
0 3/2
111693.063 (1 e .sin (2 ))
2
2
0
3/2
Khai triển phương trình, tính toán ta tìm được: e2 0.006699662637 (kết quả
tính bằng máy tính tay Casio).
Thay e2 vừa tính được vào phương trình thứ nhất của hệ (1), tìm được a:
4
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
( 2 1 )
S1 = a(1-e2).
a=
1 e sin ( )
S . 1 e sin ( )
; Với (
a
6378244.44919 m.
2
1
GVHD: PGS.TS NGUYỄN NGỌC LÂU
2
2
12
m
2
12
m
3/2
3/2
(1 e2 ).( 2 1 )
2
1 ) ở đơn vị rad.
5
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
CÁCH GIẢI CHÍNH XÁC (tính lặp): có thể chính xác đến mm.
Công thức tính chiều dài cung kinh tuyến từ xích đạo đến điểm có vĩ độ Ф là:
S A0 . A2 .sin 2 A4 .sin 4 A6 .sin 6 A8.sin8
Suy ra:
S1 S 2 S1 A0 (2 1 ) A2 (sin 22 sin 21 ) A4 (sin 42 sin 41 ) ... (1)
S2 S 4 S3 A0 (4 3 ) A2 (sin 24 sin 23 ) A4 (sin 44 sin 43 ) ... (2)
Với:
e 2 3e 4 5e6 175e8
A0 a 1
...
4 64 256 16384
3a
e 4 15e6 35e8
A2 e 2
...
8
4 128 512
15a 4 3e6 35e8
...
e
256
4
64
8
35a 6 5e
A6
...
e
3072
4
A4
A8
315a 8
(e ...)
131072
(Ta chỉ lấy đến bậc 8).
Lấy
A' ( 1 ) A2' (sin2 2 sin21 ) ...
(1)
0' 2
(2)
A0 ( 4 3 ) A2' (sin2 4 sin23 ) ...
Ta thấy: Ai a. fi (e2 )
f0 (e2 )( 2 1 ) f 2 (e2 )(sin2 2 sin21 ) ...
g (e 2 )
2
2
f0 (e )( 4 3 ) f 2 (e )(sin2 4 sin23 ) ...
Biến đổi phương trình về dạng:
e2 g ' e4 , e6 , e8
(0)
Cho e = 0.0067 tính lặp:
2
e (0) = 0.0067 e (1) → e ( 2) →…e2 chính xác.
2
2
2
( i 1)
(i )
Nếu thỏa mãn điều kiện | e
- e | < 10-10 thì dừng lại, không thỏa lặp tiếp.
Khi có e2 phương trình (1) hoặc (2) tìm ra a.
Kết thúc bài toán.
2
2
6
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
THỰC HIỆN:
Ta có:
S1 S2 S1 A0 (2 1 ) A2 (sin 22 sin 21 ) A4 (sin 42 sin 41) ...
= 0.01745329252 – 0.01742074307e2 – 0.00003249480611e4
– 0.000000054553258e6 – 9.1169355 10-11e8.
S2 S 4 S3 A0 (4 3 ) A2 (sin 24 sin 23 ) A4 (sin 44 sin 43 ) ...
= 0.01745329252 + 0.008654277117e2 + 0.006436690135e4 +
0.005319108761e6 + 0.004615207846e8.
Lập tỷ số
S1 0.01745329252 - 0.01742074307e2 - 0.00003249480611e4 0.000000054553258e6 - 9.11693551011 e8
S2
0.01745329252+0.008654277117e2 0.006436690135e4 0.005319108761e6 0.004615207846e8
= 0.9900069085.
e2
0.00001744123494 0.006404862507432e4 0.0052660089737e6 0.00456908774254e8
0.025988537202616
Tiến hành lặp tìm ra e2 và a.
7
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
-
GVHD: PGS.TS NGUYỄN NGỌC LÂU
Chương trình Matlab phục vụ công việc lặp tìm ra e2 và a:
Tên chương trình : tinh_kthuoc.m
Tham số nhập: không có.
Tham số xuất: a, e2
% CTC TINH KICH THUOC MAT SPHEROID:
function [a,e2]=tinh_kthuoc
format long
fid=fopen('kichthuoc.txt','w');
% Cho gia tri khoi dau:
e20=0.0067;
% Cho e2 gia tri chenh lech voi e20 de khoi dong vong lap:
e2=e20+1;
% Tien hanh lap de tim e2:
while abs(e2-e20)>1e-10
% Han sai cho phep la: epsilon=10^-10
e20=e2;
% Bieu thuc lap cuoi cung:
e2=(-1.744123493666698e004+0.006404862507432*e2.^2+0.005266008973700*e2.^3+0.00456908774
2540*e2.^4)/-0.025988537202616;
end
% Tim a:
a = 110576.904./(0.017453292519943-0.017420743069254*e23.249480610668151e-005*e2.^2-5.455325872419023e-008*e2.^39.116933862582735e-011*e2.^4);
% Ghi ket qua:
fprintf(fid,'Kich thuoc spheroid la:\n');
% Ghi gia tri ban truc a:
fprintf(fid,'a= ');
fprintf(fid,'%f\n',a);
% Ghi gia tri tam sai e2:
fprintf(fid,'e2= ');
fprintf(fid,'%f',e2);
fclose(fid);
disp('Gia tri e2:'); e2
disp('Gia tri ban truc a: '); a
Sau khi chạy chương trình trên ta tìm được a và e2 như sau:
a = 6378244.998786067 m.
e2 = 0.006700001641419.
8
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
ĐÁNH GIÁ VÀ KẾT LUẬN:
Kích thước
a (m)
e2
Tính gần đúng
6378244.44919
0.006699662637
Tính lặp
6378244.998786067
0.006700001641419
NHẬN XÉT:
- Từ chiều dài của 2 cung kinh tuyến và vĩ độ tương ứng ta có thể tìm được kích
thước Spheroid tương ứng.
- 2 công thức cho kết quả e2 lệch nhau không nhiều nhưng a lại lêch nhau
khoảng 0.5496 m.
- Do công thức tính gần đúng tính dựa xấp xỉ trên chiều dài cung kinh tuyến.
→ chiều dài càng lớn thì cho sai số càng lớn.
- Công thức tính lặp cho kết quả chính xác hơn và có thể áp dụng cho chiều dài
cung kinh tuyến lớn hơn.
→ Công thức tính lặp thì phù hợp hơn.
9
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
BÀI TẬP LỚN SỐ 2
1. YÊU CẦU:
Khảo sát độ chính xác của các công thức Schreiber, Gauss, Robbins, Vincenty cho
các cự ly từ 10km đến 20000km và phương vị trắc địa từ 0-1800. Từ đó hãy cho biết
rằng với độ chính xác 1 mm ở chiều dài và 0.001” ở góc. Cự ly giới hạn của các
công thức trên là bao nhiêu?
2. PHƯƠNG PHÁP GIẢI QUYẾT:
GIỚI THIỆU CÁC CÔNG THỨC
Công thức
Schreiber (1878)
Gauss (1846)
Robbins (1962)
Phương pháp tính
Trực tiếp
Lặp
Trực tiếp
Vincenty (1975)
Lặp
2.1.
Cự ly/ độ chính xác
S < 50Km
S < 100Km
Độ chính xác 20mm n
cho cự ly 1500 km, 16 m
cho đường đáy 4500 km
Độ chính xác vài mm
cho cự ly 20000km
Kiểm tra các Chương trình con (so với ví dụ trong Giáo trình
TĐCC1)
Ở bước này, các công thức bài toán thuận gồm: Schreiber, Gauss1, Robbins1
và Vincenty1 sử dụng cùng một bộ số liệu :
1 100 2839.61, 1 10604816.25, S 50000(m), A12 450
Nên ta tạo file data.txt dùng chung cho cả 4 chương trình con (CTC).
Với thứ tự số liệu như sau:
10
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
data.txt
1
1
2
3
S
a
1
A12
e2
data.txt
11
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
Ta viết thêm 1 CTC Matlab docfile.m để nhập dữ liệu từ file data.txt
function
[phi1,lam1,S,A12,a,e2]=docfile
% Mo flie:
fid=fopen('data.txt','r');
% Doc dong 1:
phi1=fscanf(fid,'%f',1);
phi1=dms2rad(phi1);
lam1=fscanf(fid,' %f',1);
lam1=dms2rad(lam1);
% Doc dong 2:
S=fscanf(fid,'%f',1);
A12=fscanf(fid,'%f',1);
A12=dms2rad(A12);
% Doc dong 3:
a=fscanf(fid,'%f',1);
e2=fscanf(fid,'%f',1);
% Dong file:
fclose(fid);
12
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
TIẾN HÀNH KIỂM TRA CÁC CÔNG THỨC
CÔNG THỨC SCHREIBER: (chỉ có Công thức bài toán thuận)
Áp dụng cho các công thức cự ly < 50 km, sử dụng các công thức tam giác cầu để
tính toán. Xấp xỉ tam giác Spheroid thành tam giác cầu để tính toán.
Thuật toán:
Chương trình : Schreiber.m
+ Tham số nhập: phi1,lam1,S,A12.
+ Tham số xuất: phi2,lam2,A21.
Tính toán:
Tiến hành chạy chương trình: do lấy dữ liệu cố định từ file data.txt nên ta chỉ việc
chạy code, không cần nhập tham số đầu vào (do code chỉ để kiểm tra nên ta cho dữ
liệu đầu vào cố định, giống với giáo trình).
function [phi2,lam2,A21]=Schreiber(phi1,lam1,S,A12)
Kết quả nhận được sau khi chạy chương trình:
phi2 = 10.474970056734399
lam2 = 1.070740082370997e+002
A12 = 2.250334836719714e+002
Vậy:
2 100 47 '49.7006 ''
2 107007 '40.0824 ''
A21 225003'34.8367"
13
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
CÔNG THỨC GAUSS: gồm CT cho bài toán Thuận và Nghịch.
- Áp dụng cho các công thức cự ly < 100 km, dựa vào một điểm nằm giữa
đường trắc địa để tính toán.
- Là công thức lặp có công thức thuận và công thức nghịch để tính toán kiểm
tra.
Thuật toán BT thuận:
Chương trình: Gauss1.m
+ Tham số nhập: phi1,lam1,S,A12.
+ Tham số xuất: phi2,lam2,A21.
Tính toán:
Tiến hành chạy chương trình: tương tự ta cũng chạy cho dữ liệu vào cố định
từ file dữ liệu data.txt.
function [phi2,lam2,A21]=Gauss1(phi1,lam1,S,A12)
Kết quả nhận được sau khi chạy chương trình:
phi2 = 10.474970045877257
lam2 = 1.070740082370979e+002
A12 = 2.250334836699251e+002
Vậy:
2 100 47 '49.70046 ''
2 107007 '40.08237 ''
A21 225003'34.8367"
Thuật toán BT nghịch:
Chương trình: Gauss2.m.
+ Tham số nhập: phi1, lam1, phi2, lam2.
+ Tham số xuất: A12_1,A21_1,S_1.
14
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
Tính toán:
Tiến hành chạy chương trình: ta lấy phi1, lam1 từ file data.txt; lấy phi2, lam2
từ kết quả chạy bài toán thuận của CT Gauss.
function [A12_1,A21_1,S_1]=Gauss2(phi1,lam1,phi2,lam2)
Kết quả nhận được sau khi chạy chương trình:
A12_1 = 45.000000038177028
A21_1 = 2.250334837080817e+002
S_1 = 4.999999986901564e+004
Vậy:
A12 45000 '00.00038"
A21 225003'34.83708"
S 49999.9999m
NHẬN XÉT:
So sánh 2 công thức thuận và nghịch ta thấy: so với số liệu khởi tính ban đầu,
chiều dài đường trắc địa sai 0.1 mm, góc phương vị thuận sai 0.00038’’ và phương
vị nghịch sai 0.00008’’ đều xấp xỉ với sai số tính toán đã cho là 0.0001’’.
15
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
CÔNG THỨC ROBBINS:
- Áp dụng cho các công thức cự ly trung bình.
- Là công thức tính trực tiếp, có công thức thuận và công thức nghịch để tính
toán kiểm tra.
Thuật toán BT thuận:
Chương trình: Robbins1.m
+ Tham số nhập: phi1, lam1, S, A12.
+ Tham số xuất: phi2, lam2, A21.
Tính toán:
Tiến hành chạy chương trình: tương tự ta cũng chạy cho dữ liệu vào cố định
từ file dữ liệu data.txt.
function [phi2,lam2,A21]=Robbins1(phi1,lam1,S,A12)
Kết quả nhận được sau khi chạy chương trình:
phi2 = 10.474970049371874
lam2 = 1.070740082331966e+002
A12 = 2.250334836684963e+002
Vậy:
2 100 47 '49.70049 ''
2 107007 '40.08233''
A21 225003'34.8367"
Thuật toán BT nghịch:
Chương trình: Robbins2.m.
+ Tham số nhập: phi1, lam1, phi2, lam2.
+ Tham số xuất: A12_1,A21_1,S_1.
16
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
Tính toán:
Tiến hành chạy chương trình: ta lấy phi1, lam1 từ file data.txt; lấy phi2,
lam2 từ kết quả chạy bài toán thuận của CT Robbins.
function [S_1,A12_1,A21_1]=Robbins2(phi1,lam1,phi2,lam2)
Kết quả nhận được sau khi chạy chương trình:
A12_1 = 45.000000015962030
A21_1 = 2.250334836844141e+002
S_1 = 4.999999987186867e+004
Vậy:
A12 450 00 '00.0002"
A21 225003'34.8368"
S 49999.9999m
NHẬN XÉT:
So sánh 2 công thức thuận và nghịch ta thấy: so với số liệu khởi tính ban đầu,
chiều dài đường trắc địa sai 0.1 mm, góc phương vị thuận sai 0.0002’’ và phương
vị nghịch sai 0.0001’’ đều xấp xỉ với sai số tính toán đã cho là 0.0001’’.
17
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
CÔNG THỨC VINCENTY:
- Áp dụng hầu hết cho các cự ly , có độ chính xác cao nhất.
- Là công thức tính trực lặp, có công thức thuận và công thức nghịch để tính
toán kiểm tra.
Thuật toán BT thuận:
Chương trình: Vincenty1.m
+ Tham số nhập: phi1, lam1, S, A12.
+ Tham số xuất: phi2, lam2, A21.
Tính toán:
Tiến hành chạy chương trình: tương tự ta cũng chạy cho dữ liệu vào cố định
từ file dữ liệu data.txt.
function [phi2,lam2,A21]=Vincenty1(phi1,lam1,S,A12)
Kết quả nhận được sau khi chạy chương trình:
phi2 = 10.474970045514842
lam2 = 1.070740082370895e+002
A12 = 2.250334836699234e+002
Vậy:
2 100 47 '49.70046 ''
2 107007 '40.08237 ''
A21 225003'34.8367"
Thuật toán BT nghịch:
Chương trình: Vincenty2.m.
+ Tham số nhập: phi1, lam1, phi2, lam2.
+ Tham số xuất: A12_1, A21_1, S_1.
18
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
GVHD: PGS.TS NGUYỄN NGỌC LÂU
Tính toán:
Tiến hành chạy chương trình: ta lấy phi1, lam1 từ file data.txt; lấy phi2,
lam2 từ kết quả chạy bài toán thuận của CT Vincenty.
function [S_1,A12_1,A21_1]=Vincenty2(phi1,lam1,phi2,lam2)
Kết quả nhận được sau khi chạy chương trình:
A12_1 = 45.000000000117247
A21_1 = 2.250334836700407e+002
S_1 = 4.999999999714783e+004
Vậy:
A12 45000 '00.00000"
A21 225003'34.83670"
S 49999.999997m
NHẬN XÉT:
So sánh 2 công thức thuận và nghịch ta thấy: so với số liệu khởi tính ban đầu,
chiều dài đường trắc địa sai 0.003 mm, góc phương vị thuận và nghịch hầu như
không sai lệch (lấy đến phần trăm ngàn của giây) → độ chính xác cao nhất.
19
BÁO CÁO BTL TRẮC ĐỊA CAO CẤP
2.2.
GVHD: PGS.TS NGUYỄN NGỌC LÂU
Khảo sát độ chính xác (cự ly giới hạn) của 4 công thức:
KHẢO SÁT CÔNG THỨC VÀ KẾT QUẢ
Với độ chính xác ở góc giới hạn là 0.01’’, độ chính xác 1mm ở chiều dài. Ta
tiến hành khảo sát lần lượt các công thức.
Chọn:
1 100
1 1060
Ta thực hiện khảo sát các cự ly từ: 10 km đến 20 000 km; góc phương vị từ: 00 đến
1800 .
Để thực hiện khảo sát, ta sử dụng các chương trình khảo sát sau:
% KHAO SAT CONG THUC ABC:
format long
% Mo file de ghi du lieu
fid=fopen('ABC.txt','w');
a=6378137; e2=0.0067;
phi1=dms2rad(10); lam1=dms2rad(106);
for S=10000:10000:20000000
for A12=0:5:180
A12=dms2rad(A12);
% Giai BT Thuan
[phi2,lam2,A21]=CT_thuan(phi1,lam1,S,A12);
% Giai BT Nghich
[S_1,A12_1,A21_1]=CT_nghich(phi1,lam1,phi2,lam2);
% Kiem tra
if abs(S_1-S)>0.001 || abs(A12_1-A12)>(0.01/206265)
|| abs(A21_1-A21)>(0.01/206265)
break;
end
end
if abs(S_1-S)>0.001
break;
end
% Ghi ket qua
fprintf(fid,'%f\t%f\t%f\t%f\n',S,S_1,S_1-S,A12);
end
disp('Cu ly gioi han cua CT ABC:'); S
disp('Goc phuong vi:'); A12
% Dong file
fclose(fid);
20
- Xem thêm -
.PDF 
30 
83 
89 
