Nghiên cứu ứng dụng thuật toán Gauss-Jourdal trong xử lý số liệu trắc địa công trình - Nguyễn Việt Hà

Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số Đặc san CBES, 04 - 2018 225 NGHIÊN CỨU ỨNG DỤNG THUẬT TOÁN GAUSS-JOURDAL TRONG XỬ LÝ SỐ LIỆU TRẮC ĐỊA CÔNG TRÌNH Nguyễn Việt Hà1,*, Phạm Tuấn Cường2 Tóm tắt: Bài báo này đã ứng dụng thuật toán Gauss_Jordan để xử lý số liệu đo trong trắc địa công trình. Cụ thể là ứng dụng thuật toán Gauss_Jordan để tiến hành đồng thời giải bài toán tuyến tính và nghịch đảo ma trận trong bài toán bình sai và đánh giá độ chính xác trị đo. Kết quả việc ứng dụng là giải thuật, modul phần mềm tính toán bình sai và đánh giá độ chính xác trị đo, làm cho bài toán đơn giản hơn và tính toán nhanh hơn trên máy tính, ngoài ra còn dễ dàng hơn cho việc quản lý, lập chương trình và cập nhật chương trình trên máy tính. Từ khóa: Xử lý số liệu trắc địa; Đánh giá độ chính xác; Nghịch đảo ma trận. 1. MỞ ĐẦU Đối với công tác xử lý số liệu trắc địa công trình, tính đúng đắn của số liệu sau xử lý không những phụ thuộc vào độ chính xác đo đạc thực địa, mà còn chịu ảnh hưởng rất lớn bởi phương pháp xử lý số liệu. Từ trước đến nay đã có nhiều chương trình tính toán xử lý số liệu trắc địa công trình như chương trình bình sai mặt bằng, độ cao; chương trình đánh giá độ chính xác mặt bằng, độ cao; chương trình tính chuyển tọa độ,... các chương trình trên thường có một nội dung tính toán rất quan trọng là giải bài toán tuyến tính và nghịch đảo ma trận cho việc tính toán các hàm trọng số phục vụ đánh giá độ chính xác của các hàm ẩn số. Trước đây, việc giải bài toán tuyến tính và nghịch đảo ma trận thường được tiến hành theo thuật toán Gauss, thuật toán khai căn Choleski hoặc thuật toán truy hồi [1, 2, 4]. Với các thuật toán trên việc giải bài toán tuyến tính và nghịch đảo ma trận được tiến hành độc lập và theo một thuật toán phức tạp, điều này sẽ làm cho bài toán dài hơn và chậm hơn, ngoài ra còn gây khó khắn cho việc quản lý và lập chương trình và cập nhật chương trình trên máy tính. Trong khuôn khổ của bài báo khoa học này, nhóm nghiên cứu chúng tôi sẽ đề xuất phương pháp ứng dụng thuật toán Gauss_Jordal [2] để giải quyết bài toán giải phương trình tuyến tính và nghịch đảo ma trận trên cùng một phép tính. Điều này mang lại nhiều ưu điểm như tiết kiệm thời gian tính toán, thuật toán đơn giản, dễ quản lý, dễ hiểu. 2. BÀI TOÁN TUYẾN TÍNH VÀ CÁC PHƯƠNG PHÁP GIẢI 2.1. Bài toán tuyến tính trong trắc địa công trình Trong toán học một hệ phương trình gồm n phương trình tuyến tính với n ẩn số x 1 , x 2 ,...,x n có dạng như sau: (1) Hệ phương trình tuyến tính trên có thể viết dưới dạng phương trình ma trận là AX = b, trong đó: 11 12 1 21 22 2 1 2 ... ... ... .... ... ... ... n n n n nn a a a a a a A a a a             ; 1 2 ... n x x X x             ; 1 2 ... n b b b b                       nn nn 2n2 1n1 2n 2n 222 121 1n 1n 212 111 b xa . . . xa xa ............................ . . . . . . . . . . . . . . . . b xa . . . xa xa b xa . . . xa xa Toán học, Cơ học & Ứng dụng N. V. Hà, P. T. Cường, “Nghiên cứu ứng dụng thuật toán số liệu trắc địa công trình.” 226 Nếu det(A) ≠ 0 thì nghiệm của hệ (1) có thể tính theo công thức 1X A b . Áp dụng công thức tính ma trận nghịch đảo ta có thể biến đổi và dẫn đến lời giải được diễn tả bằng định lý Cramer như sau: Định lý Cramer. Gọi A j là ma trận nhận được từ ma trận A bằng cách thay cột thứ j bằng cột b, khi đó hệ (1) có nghiệm duy nhất và x j được tính bởi công thức det(A ) det(A) j jx  (2) Tuy nhiên, trong thực hành người ta không dùng công thức này để tính nghiệm vì số phép tính quá lớn. Người ta dùng những phương pháp hữu hiệu hơn mà chúng tôi sẽ giới thiệu sau đây. 2.2. Một vài phương pháp giải bài toán tuyến tính trong trắc địa công trình a. Phương pháp khử Gauss Phương pháp khử Gauss dùng cách khử dần các ẩn số để đưa hệ phương trình tuyến tính về một dạng tam giác trên rồi giải hệ tam giác này từ dưới lên trên, không phải tính một định thức nào. Phương pháp này được thực hiện qua các bước sau: Quá trình xuôi: - Bước 1: Dùng phương trình trên cùng để khử x 1 trong n-1 phương trình còn lại. Giả sử a 11 ≠0. (Để cho đơn giản, trước khi khử ta có thể chia phương trình thứ nhất cho a 11 ). Cụ thể để khử x 1 ở các hàng tiếp theo thứ k( k=2,3,n) ta phải tính lại các hệ số a kj ở hàng thứ k (với j=1,2,..n+1) như sau: a kj =a kj -a 1j *a k1 /a 11 - Bước 2: Dùng phương trình thứ 2 để khử x 2 trong n-2 phương trình còn lại tiếp theo. Giả sử a 22 ≠0. (Để cho đơn giản, trước khi khử ta có thể chia phương trình thứ hai cho a 22 ). Cụ thể để khử x 2 ở hàng tiếp theo thứ k (k=3,4,n) ta phải tính lại các hệ số a kj ở hàng thứ k (j=2,..n+1) như sau: a kj =a kj -a 2j *a k2 /a 22 - Bước i: Dùng phương trình i để khử x i trong các phương trình tiếp theo thứ i+1,i+2, ..., n. Giả sử a ii ≠0. Để cho công thức đơn giản, trước khi khử ta có thể chia phương trình thứ i cho a ii . Cụ thể để khử x i ở hàng thứ k (k=i+1,n) ta phải tính lại các hệ số a kj ở hàng thứ k (j=i,..n+1) như sau: a kj =a kj -a ij *a ki /a ii - Bước n-1: Dùng phương trình thứ n-1 để khử x n-1 trong phương trình thứ n. Giả sử a n-1 n-1 ≠0. (Để cho đơn giản, trước khi khử x n-1 ta có thể chia phương trình thứ n-1 cho a n-1 n-1 ) Cụ thể để khử x n-1 ở hàng thứ n ta phải tính lại các hệ số a nj ở hàng thứ n (j=n- 1,n,n+1) như sau: a nj =a nj -a n-1j *a n-1i /a n-1n-1 Kết thúc quá trình khử. Quá trình giải ngược x n = b n /a nn hoặc ( x n =b n ) . . . x i = (b i -(aΣ+=nij1 ij x j ) )/a ii ) hoặc (b i -(aΣ+=nij1 ij x j ) ), i =n-1, n-2, ..., 1 Ta áp dụng các phép biến đổi sơ cấp như vừa trình bày để biến đổi ma trận chữ nhật cấp nx(n+1) trên đây về dạng: Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số Đặc san CBES, 04 - 2018 227 b. Phương pháp khử Gauss-Jordan Phương pháp khử Gauss-Jordan dùng cách khử dần các ẩn để đưa hệ phương trình đã cho về một dạng ma trận đường chéo rồi giải hệ phương trình này, không phải tính một định thức nào Phương pháp này được thực hiện qua các bước sau: - Bước 1: Dùng phương trình đầu tiên để khử x 1 trong n-1 phương trình còn lại, cách làm tương tự như phương pháp khử để tính định thức... (Để cho công thức đơn giản, trước khi khử ta có thể chia phương trình thứ nhất cho a 11 ). Cụ thể để khử x 1 ở hàng thứ k( k=2,3,n) ta phải tính lại các hệ số a kj ở hàng thứ k (j=1,2,..n+1) như sau: a kj =a kj -a 1j *a k1 /a 11 - Bước i: Dùng phương trình i để khử x i trong các phương trình thứ 1,2, i- 1,i+1,i+2,...,n.. (Để cho công thức đơn giản, trước khi khử ta có thể chia phương trình thứ i cho a ii ). Cụ thể để khử x i ở hàng thứ k (k=1,2, i-1,i+1,i+2,...,n.) ta phải tính lại các hệ số a kj ở hàng thứ k (j=i,..n+1) như sau: a kj =a kj -a ij *a ki /a ii . - Bước n: Dùng phương trình thứ n để khử x n trong phương trình thứ 1,2, ..., n-1.. (Để cho công thức đơn giản, trước khi khử ta có thể chia phương trình thứ n cho a nn ) Cụ thể để khử x n ở hàng thứ k( k=1,2, ..,n-1.) ta phải tính lại các hệ số a kj ở hàng thứ k (j=n,n+1) như sau: a kj =a kj -a nj *a kn /a nn Tương tự phép khử Gauss tại mỗi bước, trước khi khử ta phải chọn trụ tối đại. Cụ thể tại bước i ta luôn chọn hàng có phần tử a ri có giá trị tuyệt đối lớn nhất rồi đổi cho hàng thứ i cho hàng thứ r. Hệ phương trình sau khi khử có dạng: Hoặc nếu tại các bước (bước i) ta chia cho hệ số a ii : Tức là ta đã có các nghiệm mà không cần phải tính toán thêm. Cũng như trong phương pháp khử Gauss, khi cài đặt trên máy tính ta dùng một mảng a thay cho cả ma trận A và vec tơ b. Tức là:                    1 ' 12 ' 2 ' 11 ' 1 ' 12 1......00 .................. .................. ......10 ......1 nn nn nn a aa aaa           nn nn 2 2 22 1 1 11 b xa ....................................... bxa b xa 1 1 2 2 x b x b .......................................   n n x b        Toán học, Cơ học & Ứng dụng N. V. Hà, P. T. Cường, “Nghiên cứu ứng dụng thuật toán số liệu trắc địa công trình.” 228 Ta áp dụng các phép biến đổi sơ cấp như vừa trình bày để biến đổi ma trận chữ nhật cấp nx(n+1) trên đây về dạng Vậy ta có x i = a' i,(n+1). 3. ÁP DỤNG PHƯƠNG PHÁP KHỬ GAUSS - JORDAN ĐỂ TÍNH MA TRẬN NGHỊCH ĐẢO Để giải hệ n phương trình n ẩn Ax = b, trong phương pháp khử Gauss-Jordan ta đã dùng các phép biến đổi sơ cấp để đưa phương trình này về dạng: Ex = b' (3) Vì Ex = x, do đó, ta có x=b'. Nếu B là một ma trận chữ nhật cấp n x k tùy ý, ta có thể áp dụng phương pháp khử Gauss-Jordan để giải đồng thời k hệ n phương trình n ẩn: AX = B (4) trong đó: Ta viết ma trận B bên phải ma trận A như sau: Nếu ma trận A không suy biến, ta có thể áp dụng các phép biến đổi sơ cấp để đưa ma trận này về dạng:                                     nnnnn n n nnnnnn nn nn baaa baaa baaa aaaa aaaa aaaa ...... .................. .................. ...... ...... ...... .................. .................. ...... ...... 21 222221 111211 121 1222221 1111211                    1 ' 12 ' 11 ' 1......00 .................. .................. 0......10 0......01 nn n n a a a              nknn k k xxx xxx xxx X ... ............ ... ... 21 22221 11211              nknn k k bbb bbb bbb B ... ............ ... ... 21 22221 11211             nknnnnnn kn kn bbbaaa bbbaaa bbbaaa ...... ........................ ...... ...... 2121 2222122221 1121111211 11 12 1 21 22 2 1 2 1 0 ... 0 ' ' ... ' 0 1 ... 0 ' ' ... ' ... ... ... ... ... ... ... ... 0 0 ... 1 ' ' ... ' k k n n nk b b b b b b b b b             Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số Đặc san CBES, 04 - 2018 229 Khi đó, ta có: Đặt: Xét trường hợp đặc biệt B = E, ta có ma trận B' chính là ma trận nghịch đảo của ma trận A. Thật vậy, nếu X là nghiệm của (4) thì: X = A -1 B (5) Nếu B = E thì X = A -1 . Do đó, việc tìm ma trận nghịch đảo của ma trận A tương đương với việc giải phương trình AX = E (6) Vậy để tính ma trận nghịch đảo ta cần thực hiện các bước sau: • Viết thêm ma trận đơn vị E bên cạnh ma trận A • Áp dụng phép biến đổi sơ cấp trên hàng cho đến khi ma trận có dạng Khi đó ta có Chú ý. Trong quá trình biến đổi ta có thể đổi các hàng của ma trận. Điều này không ảnh hưởng đến kết quả thu được: Ma trận C vẫn là ma trận nghịch đảo của ma trận A ban đầu. Lý do là vì để tìm ma trận nghịch đảo ta chỉ cần xác định ma trận nghiệm. Ma trận nghiệm không bị thay đổi nếu ta đổi chỗ các hàng. 4. LẬP CHƯƠNG TRÌNH TÍNH TOÁN ÁP DỤNG THUẬT TOÁN GAUSS- JORDAN ; '...'' ............ '...'' '...'' ... ............ ... ... 21 22221 11211 21 22221 11211                           nknn k k nknn k k bbb bbb bbb xxx xxx xxx X              nknn k k bbb bbb bbb B '...'' ............ '...'' '...'' ' 21 22221 11211             1...00... ........................ 0...10... 0...01... 21 22221 11211 nnnn n n aaa aaa aaa             nnnn n n ccc ccc ccc ...1...00 ........................ ...0...10 ...0...01 21 22221 11211              nnnn n n ccc ccc ccc A ... ............ ... ... 21 22221 11211 1 Toán học, Cơ học & Ứng dụng N. V. Hà, P. T. Cường, “Nghiên cứu ứng dụng thuật toán số liệu trắc địa công trình.” 230 a. Môđun lập hệ phương trình số hiệu chỉnh và lập hệ phương trình chuẩn void Pthc(int i,int sdm) { int m,i1,i2; for( m=1;m<=sdm+1;m++) a[m]=0.0; i1=trido[i].ds; i2=trido[i].dt; a[i1]=-1.0; a[i2]=1.0; a[sdm+1]=trido[i].hh-(diem[i2].h-diem[i1].h); } //-------------------------------------------------------} void Ptc(int i,int sdm) { double p; int j,i1,i2; p=1.0/trido[i].tram; for(i1=1;i1<=sdm+1;i1++) for(i2=i1;i2<=sdm+1;i2++) rnm[i2*(i2-1)/2+i1]=rnm[i2*(i2-1)/2+i1]+p*a[i1]*a[i2]; } b. Môđun giải ẩn và nghịch đảo ma trận bằng thuật toán Gauss-Jordal Đây là môđun để giải ẩn và nghịch đảo ma trận trên một thuật toán, với đoạn chương trình rất ngắn này có thể thay thế được những môđun rất dài khi viết bằng thuật toán Gauss. Theo nhóm nghiên cứu thì cùng một bài toán mà số dòng viết bằng thuật toán Gauss-Jordal/ số dòng viết bằng thuật toán Gauss là 18/80. Điều này cho thấy thuật toán mà tác giải lựa chọn là tối ưu hơn về mặt quản lý và tốc độ. void Gauss_Jordan_Invers(int sdm) { int i,j,k,m = 2*sdm+1; for (i = 0; i < sdm; i++) for (j = 0; j < sdm; j++) //tao ma tran E phia sau if (j == i) GA[i][sdm+j+1] = 1; else GA[i][sdm+j+1] = 0; for (k = -1; k < sdm - 1; k++) { i = k + 1; for (j = m; j >= i; j--) GA[i][j] = GA[i][j]/GA[i][i]; for (i = 0; i < sdm; i++) { if (i == k + 1) continue; for (j = m; j > k; j--) GA[i][j] = GA[i][j] - GA[i][k+1]*GA[k+1][j]; } } } 5. TÍNH TOÁN THỰC NGHIỆM VÀ KIỂM TRA ĐỘ CHÍNH XÁC CỦA CHƯƠNG TRÌNH Sau khi hoàn thành phần mềm do nhóm nghiên cứu thành lập, để thực nghiệm kiểm chứng độ chính xác và tin cậy của phần mềm, nhóm đã chạy chương trình lấy số liệu quan Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số Đặc san CBES, 04 - 2018 231 trắc lún công trình Tòa nhà hỗn hợp HH3 khu đô thị mới Mỹ Đình-Mễ Trì-Từ Liêm-Hà Nội làm số liệu thực nghiệm, vẫn số liệu đo trên được nhóm nghiên cứu tiến hành bình sai bằng hai phần mềm khác là phần mềm PICKNET 3.0 và phần mềm BSDC của TS Nguyễn Việt Hà để đánh giá độ chính xác và độ tin cậy của chương trình, kết quả được thể hiện ở bảng dưới. Chỉ tiêu đánh giá Kết quả do nhóm nghiên cứu Kết quả phần mềm BSDC (Nguyễn Việt Hà) Kết quả phần mềm Picknet 3.0 Sai số trung phương trọng số đơn vị 0.21 (mm) 0.21 (mm) 0.21 (mm) Sai số điểm yếu nhất 0.64 (mm) 0.64 (mm) 0.64 (mm) Chênh cao có số hiệu chỉnh lớn nhất 1.47 (mm) 1.5 (mm) 1.5 (mm) 6. KẾT LUẬN Thuật toán Gauss_Jordan mà nhóm tác giả nghiên cứu cho thấy có thể ứng dụng trong công tác xử lý số liệu trắc địa công trình, tính đúng đắn của số liệu sau xử lý đã được tính toán thực nghiệm và so sánh với một số phần mềm tin cậy hiện nay. Với thuật toán Gauss_Jordan cho phép giải bài toán tuyến tính và nghịch đảo ma trận được tiến hành đồng thời bằng một thuật toán đơn giản, điều này sẽ làm cho bài toán ngắn hơn và tính toán nhanh hơn trên máy tính, ngoài ra còn dễ dàng hơn cho việc quản lý, lập chương trình và cập nhật chương trình trên máy tính. TÀI LIỆU THAM KHẢO [1]. Hoàng Ngọc Hà “Tính toán trắc địa và cơ sở dữ liệu”. NXB Giáo dục -2002 [2]. Nguyễn Trọng San, Đào Quang Hiếu , Nguyễn Tiến Năng “Giáo trình trắc địa phổ thông”. Đại học Mỏ - Địa Chất -1992 [3]. Nguyễn Trọng San, Đào Quang Hiếu, Đinh Công Hòa “Trắc địa cơ sở'' Tập 1, Tập 2. Hà Nội – 2002. [4]. Trần Khánh, Nguyễn Việt Hà “Bài giảng ứng dụng tin học trong trắc địa công trình” Hà Nội 2012. [5]. Phạm Hồng Thái “Bài giảng ngôn ngữ lập trình C/C++ ’’. Hà Nội 2003. ABSTRACT ON THE APPLICATION OF GAUSS-JOURDAL FOR PROCESSING DATA IN ENGINEERING SURVEYING In this paper, an application of Gauss-Jordan for processing surveyed measurements in the field of engineering surveying is presented. This algorithm is applied for determining response variables and inverse covariance matries in purpose of the adjustment and evaluation of accuracy. A modul using this algorithm is developed for surveying applications which has many benefits such as faster for running, easier for control of software and update on data. Keywords: Surveying processing; Evaluation of accuracy; Inverse matrix. Nhận bài ngày 25 tháng 02 năm 2018 Hoàn thiện ngày 16 tháng 3 năm 2018 Chấp nhận đăng ngày 20 tháng 3 năm 2018 Địa chỉ: 1 Khoa Trắc địa, Bản đồ và Quản lý đất đai, Trường Đại học Mỏ - Địa chất; 2 Khoa Khoa học cơ bản, Trường Đại học Mỏ - Địa chất ; * Email: nguyenvietha@humg.edu.vn, phamtuancuong@humg.edu.vn