Nâng cao độ chính xác định vị GPS động bằng lọc Kalman

Nghiên cứu Tạp chí Khoa học Tài nguyên và Môi trường - Số 19 - năm 201820 NÂNG CAO ĐỘ CHÍNH XÁC ĐỊNH VỊ GPS ĐỘNG BẰNG LỌC KALMAN Đinh Xuân Vinh, Cao Minh Thủy Trườ ng Đạ i họ c Tà i nguyên và Môi trườ ng Hà Nộ i Tó m tắ t Nội dung thảo luận về những phương trình vi phân trong bộ lọc Kalman phù hợp với quá trình định vị GPS động. Quá trình định vị GPS động được mô tả như những chuyển động ngẫu nhiên theo thời gian. Các trị đo GPS được biểu diễn trong một phương trình vi phân kèm theo nhiễu thực. Một mô hình biểu diễn thời gian của chuyển động được xây dựng kèm theo nhiễu trắng. Sự kết hợp các phương trình vi phân của thế giới thực và mô hình được xây dựng dựa trên các phương trình biểu diễn sự phù hợp với chuyển động của anten máy thu GPS động. Khảo sát này có thể góp phần nâng cao độ chính xác thành lập bản đồ tỷ lệ lớn, giảm chi phí nhân lực, tăng tiến độ thi công và bổ sung ứng dụng cho bộ lọc Kalman. Từ khóa: Lọc Kalman; GPS động. Abstract Enhance the accuracy of kinematic GPS positioning with Kalman fi ltering This paper discusses the differential equations in the Kalman fi lter that are suitable for kinematic GPS positioning. Kinematic GPS positioning is described as random motion over time. GPS measurements are expressed in a differential equation with real noise. A time model of motion is constructed with white noise. The combination of real world and model is based on equations, that are representing the motion of the kinematic GPS receiver antenna. This paper can contribute to improving the accuracy of large scale mapping, reducing employee costs, increasing schedules and adding Kalman fi lter applications. Keywords: Kalman fi lter; Kinematic GPS. A. Giớ i thiệ u lọ c Kalman Phé p lọ c tuyế n tí nh cá c giá trị đo đượ c củ a mộ t tậ p hợ p cá c biế n ngẫ u nhiên để ướ c lượ ng, hay nó i chí nh xá c hơn là để dự bá o cá c giá trị không đo đượ c củ a mộ t tậ p hợ p khá c, đã đượ c cá c nhà khoa họ c thế giớ i quan tâm từ rấ t sớ m. Phương phá p đầ u tiên đị nh hì nh ướ c lượ ng tố i ưu từ dữ liệ u có nhiễ u là phương phá p Bì nh phương nhỏ nhấ t củ a Gauss (1777 - 1855). Tí nh chấ t chắ c chắ n củ a trị đo có chứ a sai số (nhiễ u) đượ c xá c nhậ n bở i Galileo (1564 - 1642). Đầ u thế kỉ 20, Kolmogorov (1903 - 1987) và Wiener (1894 - 1964) đã sá ng tạ o ra lý thuyế t dự bá o, là m mề m và lọ c theo quy trì nh ngẫ u nhiên Markov. Kalman (1960) đề xuấ t phương phá p lọ c tuyế n tí nh mớ i [1], giả i quyế t được bà i toá n vi phân bậ c hai tuyế n tí nh. Trạ ng thá i tứ c thờ i củ a mô hì nh tuyế n tí nh độ ng vớ i sự tham gia củ a nhiễ u trắ ng được ước lượng bằ ng sử dụ ng cá c trị đo ở trạ ng thá i tương quan tuyế n tí nh xen lẫ n nhiễ u trắ ng. Phương trình hệ thống lọc Kalman rời rạc là một ước lượng trạng thái xRn theo một quy trình bị chi phối bởi phương trình vi phân tuyến tính ngẫu nhiên sau: 1 1 1k k k kx Fx Gu w     (1) Nghiên cứu Tạp chí Khoa học Tài nguyên và Môi trường - Số 19 - năm 2018 21 Với trị đo zRm tuân theo phương trình sau k k kz Hx v  (2) trong đó: xk là vector chỉ trạng thái hệ thống; ma trận F kích thước (n x n) trong phương trình vi phân là ma trận hệ số của ẩn tại trạng thái trước đó (k-1) so với trạng thái hiện thời k. Ma trận G là ma trận hệ số đầu vào điều chỉnh tùy ý của ẩn uRl liên hệ với trạng thái của ẩn x, trong trắc địa thì nó biểu thị các nguyên nhân gây nên biến đổi hệ thống, ảnh hưởng tới quy trình ngẫu nhiên của hệ thống. Ma trận H kích thước (m x n) trong phương trình trị đo là ma trận hệ số của trị đo zk, wk-1 là nhiễu trắng hệ thống và nó được biểu diễn như một vector; vk là nhiễu trắng trị đo được biểu diễn dưới dạng vector. Chỉ số k chỉ thời điểm của hệ thống và k-1 là thời điểm trước đó. Phương trình (1) phù hợp với mô hình vận động (là mô hình có ngoại lực tác động gây biến đổi vận tốc và gia tốc) và không thể tìm thấy trong mô hình động (là mô hình không có ngoại lực tác động) thành phần Guk-1 vì không có nguyên nhân gây biến dạng nào được tính đến trong mô hình. Cũng không thể tìm thấy trong mô hình tĩnh thành phần Fxk-1 vì vật thể phản ứng ngay tức thì với những thay đổi đầu vào. Trong mô hình đồng nhất không có nguyên nhân gây biến dạng, nên ma trận hệ thống được xác định là ma trận đơn vị. Vector trạng thái tự nhiên của xk lẽ dĩ nhiên là biến không đo được, còn zk là giá trị đo được. Biến ngẫu nhiên wk-1 và vk biểu diễn nhiễu hệ thống và nhiễu trị đo, chúng được giả thiết là độc lập với nhau, là nhiễu trắng và tuân theo phân phối chuẩn, nghĩa là p(w)~N(0,Q) (3) p(v)~N(0,Q) (4) Ta có ma trận nhiễu hệ thống Q liên quan tới vector nhiễu hệ thống theo: Q=E[wwT] (5) Ma trận nhiễu trị đo R có liên hệ với vector nhiễu trị đo v theo: R=E[vvT] (6) Nếu chúng ta mang những trị đo với chu kỳ T s để đưa vào phép lọc, thì việc đầu tiên là ta phải tìm được ma trận cơ sở Φ. Ma trận cơ sở của hệ thời gian bất biến có thể tìm được từ ma trận hệ thống động [3] như sau: (t)=ଂ-1[(sl-F)-1] (7) Ở đây, I là ma trận đơn vị, ଂ-1 là biến đổi Laplace nghịch đảo, F là ma trận hệ thống động. Có thể chứ ng minh đượ c cá c phương trì nh Riccati biể u diễ n Hiệ p phương sai tiên nghiệ m, Hiệ p phương sai hậ u nghiệ m và giá trị Hiệ u í ch củ a cá c bướ c lọ c Kalman. Phương trì nh Riccati như sau: 1 k T k k k k'M P Q   (8) 1T T k k k kK M H ( HM H R )   (9) k k kP ( I K H )M  (10) Ở đây, Pk là ma trận hiệp phương sai mô tả sai số trong ước lượng trạng thái sau khi cập nhật; Mk là ma trận hiệp phương sai mô tả sai số trong ước lượng trạng thái trước khi cập nhật. Ma trận nhiễu rời rạc Qk có thể tìm được từ ma trận nhiễu liên tục Q và ma trận cơ sở theo 0 sT TQk ( )Q ( )d      (10) Để bắt đầu phương trình Riccati, ta cần ma trận hiệp phương sai ban đầu P 0 . Nghiên cứu Tạp chí Khoa học Tài nguyên và Môi trường - Số 19 - năm 201822 B. Thự c nghiệ m quan trắ c GPS độ ng Mục tiêu áp dụng phương pháp xử lý sau trong công tác thu tín hiệu GPS động với loại máy thu GPS thông dụng, rẻ tiền, nhưng có thể cho ta chất lượng vị trí điểm đạt độ chính xác cỡ xen ti met, phù hợp yêu cầu xây dựng lưới khống chế đo vẽ tỷ lệ lớn, hoặc đo vẽ chi tiết thành lập bản đồ tỷ lệ 1:500. Chúng tôi tiến hành một cách cẩn thận công tác thu tín hiệu GPS tại khu vực trường Đại học Tài nguyên và Môi trường Hà Nội, ngày 07 tháng 9 năm 2017. Thiết bị thu tín hiệu GPS gồm 3 máy thu loại tín hiệu một tần số X20 của hãng Huace - Trung Quốc, số hiệu các máy là: 100957, 100961 và 101533. Loại máy X20 tương đối cũ, chỉ thu được duy nhất tín hiệu GPS khoảng cách giả và GPS pha sóng tải, không thu được tín hiệu Glonass và Beidou. Đặt góc chân trời 150, tần số lấy mẫu 5 giây. Thời gian đo bắt đầu lúc 10 giờ 20 phút, giờ Hà Nội. Kết thúc đo lúc 11 giờ 38 phút. Trong 50 phút đầu tiên, ba máy thu ở chế độ tĩnh. Phương pháp đo tương đối cho phép xác định chính xác vị trí điểm máy thu. Khoảng 25 phút cuối, máy số 100961 di động theo hai hướng gần vuông góc nhau. Đó là do máy 100961 đặt ở ngã ba đường. Quá trình di động máy theo hai con phố khoảng 18 phút, sau đó đặt trở lại máy vào chân ba chạc vẫn giữ nguyên trên mốc khoảng 7 phút, sau đó kết thúc ca đo. Đặc thù phố nhỏ, bề ngang phố khoảng 5 mét. Một con phố có nhiều cây to hai vệ đường, dẫn tới tín hiệu GPS bị mất trong vài phút. Điều kiện vệ tinh và chất lượng máy thu rất khiêm tốn. Hầu hết thời gian đo chỉ thu được tín hiệu 5 vệ tinh. Đây gần như là giới hạn cuối của chất lượng ca đo. Đối với ca đo tĩnh thì không vấn đề gì, nhưng với ca đo động thì chất lượng tín hiệu rất tồi. Khoảng cách giữa các điểm trắc địa khá gần nhau, từ 229 mét đến 280 mét và phân bố như hình 1. Cấu tạo mốc trắc địa như hình 2. Hình 3 là máy thu 100961 và sinh viên đo đạc. Hình 4 là tập trị đo và hướng di chuyển của máy 100961. Xử lý dữ liệu đo bằng phần mềm Compas kèm theo máy. Số liệu đo tĩnh theo phương pháp tương đối đạt kết quả tốt. Sai số vị trí điểm thu được có độ chính xác ±1 mm. Lý do là khoảng cách giữa các điểm khá gần nhau (chưa đến 300 mét). Về trị đo động, lựa chọn 249 trị đo động có thời gian từ 4:15:00 GPST đến 4:38:15 GPST. Khoảng thời gian từ 4:28:20 GPST đến 4:30:55 GPST không thu được tín hiệu. Lý do, máy động di chuyển trên con phố nhỏ, có nhiều cây xanh ven đường, máy X20 chỉ nhận được tín hiệu GPS của 5 vệ tinh, các vệ tinh Glonass và Beidou không cấu trúc trong máy thu X20. Sau khi di chuyển theo hai con phố và quay trở lại mốc trắc địa cũ, máy 100961 đặt trở lại trên chân máy vẫn giữa nguyên trên mốc trắc địa và thu tín hiệu GPS với thời gian khoảng 7 phút. Sau đó tắt máy và kết thúc thực nghiệm. Nghiên cứu Tạp chí Khoa học Tài nguyên và Môi trường - Số 19 - năm 2018 23 Hình 1: Sơ đồ điểm thu GPS Hình 2: Mốc khống chế gần cổng trường Hình 3: Máy thu 100961 Hình 4: Tập trị đo và hướng di chuyển của 100961 Bảng 1. Một đoạn trị đo GPS động (WGS84) Kinematic GPST x-ecef(m) y-ecef(m) z-ecef(m) 1 4:15:00 AM -1617957.304 5731130.610 2276349.631 2 4:15:05 AM -1617958.055 5731134.222 2276350.955 3 4:15:10 AM -1617957.632 5731132.958 2276350.265 4 4:15:15 AM -1617958.075 5731134.736 2276350.733 5 4:15:20 AM -1617957.721 5731134.335 2276350.712 6 4:15:25 AM -1617957.813 5731135.291 2276351.113 7 4:15:30 AM -1617957.576 5731134.737 2276350.986 8 4:15:35 AM -1617957.569 5731134.896 2276351.047 9 4:15:40 AM -1617957.597 5731135.530 2276351.200 Nghiên cứu Tạp chí Khoa học Tài nguyên và Môi trường - Số 19 - năm 201824 Hình 5: Số lượng vệ tinh trong ca đo thực nghiệm C. Xử lý số liệu Lọc Kalman được triển khai theo phương pháp sử dụng tích phân Euler và phương pháp Bình phương nhỏ nhất. Các phương trình Riccati nhằm nâng cao hiệu suất của bộ lọc và tối ưu hóa sau mỗi bước lọc được triển khai dưới dạng đa thức. Do vậy, lọc Kalman cũng được triển khai ở dạng đa thức. Thự c hiệ n lọ c Kalman dạ ng đa thứ c vớ i cá c bậ c lầ n lượ t là 0, 1, 2. Đố i vớ i bậ c 0, ta có phương trì nh ma trậ n cơ sở ban đầ u như sau: 1 1k k k k kˆ ˆ ˆx x K ( z x )    (12) Ký hiệu: kxˆ ước lượng Kalman tại thời đoạn k; 1kxˆ  ước lượng Kalman thời đoạn k-1; K1k hiệu ích của ước lượng Kalman bậc 0 (trạng thái 1) tại thời đoạn k; zk trị đo tại thời đoạn k. Độ lệch của lọc bậc 0 được định nghĩa: 1k k kˆRe s z x   Lọ c bậ c 1 kèm theo vận tốc có dạ ng như sau: (13) Ký hiệu: tần số đo (khoảng cách giữa các thời đoạn); hiệu ích Kalman đối với vị trí điểm; hiệu ích Kalman đối với vận tốc chuyển động của điểm; ước lượng vận tốc điểm Kalman tại thời đoạn k; ước lượng gia tốc điểm Kalman tại thời đoạn k; ước lượng vận tốc của điểm tại thời đoạn k-1. Độ lệch của lọc bậc 1 được định nghĩa: Hiệu ích của lọc Kalman bậc 1 được tính theo phương pháp bình phương nhỏ nhất đệ quy: k=1, 2, 3, ..., n. Lọ c bậ c 2 kèm theo gia tốc và vận tốc có dạ ng như sau: Nghiên cứu Tạp chí Khoa học Tài nguyên và Môi trường - Số 19 - năm 2018 25 (14) Ký hiệu: ước lượng gia tốc điểm Kalman tại thời đoạn k; ước lượng gia tốc của điểm tại thời đoạn k-1; K3k hiệu ích Kalman đối với gia tốc chuyển động của điểm. Độ lệch của lọc bậc 2 được định nghĩa: (15) Hiệu ích của lọc Kalman bậc 2 được tính:   2 1 3 3 3 2 1 2 1 2k ( k k )K , k , , ,n, k k ( k )      (16) 2 18 2 1 1 2k s ( k )K , k( k )( k )T    (17)   3 2 60 1 2k s K , k k k T    (18) Hình 6: Giá trị tọa độ Y (nét liền) và lọc Kalman (nét đứt) theo thời gian Hình 7: Giá trị tọa độ X (nét liền) và lọc Kalman (nét đứt) theo thời gian Lưu ý, khoảng 7 phút cuối máy 100961 ở trạng thái tĩnh. Theo đồ thị 6 và 7, trị đo X và Y của máy biến động rất lớn, thậm chí hơn 10 mét. Sử dụng bộ lọc Kalman cho ta một ước lượng tốt hơn. Vị trí điểm 100961 có tọa độ (VN2000) xác định trong ca đo tĩnh và đo động là: Bảng 2. So sánh tọa độ hai ca đo tĩnh và động của máy 100961 STT Tọa độ X (mét) Tọa độ Y (mét) 1 Ca đo tĩnh 2328352.7827 579494.7087 2 Ca đo động 2328352.7540 579494.6767 3 Độ lệch 0.0287 0.0320 Nghiên cứu Tạp chí Khoa học Tài nguyên và Môi trường - Số 19 - năm 201826 C. KẾ T LUẬ N Lọc Kalman thực sự rất hiệu quả khi phân tích dữ liệu vận động theo thời gian. Những nghiên cứu bước đầu nhằm sử dụng máy thu GPS rẻ tiền, thông dụng thay thế những máy thu GPS RTK đắt tiền. Kết quả có thể chấp nhận được trong điều kiện số lượng vệ tinh rất ít (5 vệ tinh), chất lượng tín hiệu yếu và đôi khi mất tín hiệu. Với độ chính xác ±3 cm có thể ứng dụng trong đo vẽ chi tiết hoặc xây dựng lưới khống chế đo vẽ. Nếu có thể sử dụng loại máy thu được các vệ tinh GPS, Glonass, Beidou thì chất lượng sẽ nâng lên đáng kể, thời gian đo cũng vì thế tăng nhanh hơn. Lời cảm ơn: Nghiên cứu này được thực hiện với sự hỗ trợ của đề tài nghiên cứu khoa học cấp cơ sở “Nghiên cứu nâng cao độ chính xác xác định vị trí điểm máy thu trong đo GPS động xử lý sau phục vụ thành lập bản đồ tỷ lệ lớn ở Việt Nam”, mã số 13.01.17.O.03. Tác giả chân thành cảm ơn nhóm sinh viên ĐH5QĐ9 đã nhiệt tình tham gia thực nghiệm. TÀ I LIỆ U THAM KHẢ O [1]. R.E. Kalman (1960). A New Approach to Linear Filtering and Prediction Problems, Journal of Basic Engineering, 82 (series D):34-45. Copyright @ 1960 by ASME. [2]. Norman Morrison (1969). Intro to Sequential Smoothing and Prediction. McGraw-Hil Book Company, New York. [3]. Đinh Xuân Vinh, Phan Văn Hiến, Nguyễn Bá Dũng (2016). Lý thuyết và phương pháp phân tích biến dạng. Nhà xuất bản Tài nguyên Môi trường và Bản đồ Việt Nam. ISBN: 978-604-904-875-3. [4]. Phan Văn Hiến, Đinh Xuân Vinh (2010). Ứng dụng lọc Kalman trong phân tích biến dạng nhà cao tầng do bức xạ nhiệt mặt trời. Tạp chí Xây dựng, số 5-2010. ISSN 0866-0762. [5]. Đề tài nghiên cứu khoa học cấp cơ sở “Nghiên cứu nâng cao độ chính xác xác định vị trí điểm máy thu trong đo GPS động xử lý sau phục vụ thành lập bản đồ tỷ lệ lớn ở Việt Nam”, mã số 13.01.17.O.03. Chủ nhiệm đề tài: TS. Đinh Xuân Vinh. BBT nhận bài: 22/01/2018; Phản biện xong: 27/2/2018 [5]. Gilad Even-Tzur, Lior Shahar (2015). Application of extended free net adjustment constraints in two-step analysis of deformation network. Acta Geod Geophys DOI 10.1007/s40328-015-0119-3, Springer. [6]. Dagogo M.J. FUBARA (1973). Geodetic numerical and statistical analysis of data. Space Program Offi ce, Battelle Columbus Laboratories, Columbus, Ohio 43201. https://link.springer.com/ article/10.1007%2FBF02522077. [7]. Haim B. Papo (1986). Extended free net adjustment constraints. US Department of commerce, National Oceanic and Atmospheric Administration, National Ocean Service. [8]. Hoàng Ngọc Hà, Trương Quang Hiếu (2003). Cơ sở toán học xử lý số liệu trắc địa. NXB Giao thông vận tải. [9]. Hoàng Ngọc Hà (2006). Bình sai tính toán lưới trắc địa và GPS. NXB Khoa học kỹ thuật. [10]. Huang Shengxiang, Yin Hui, Jiang Zheng (2004). Xử lý số liệu trong quan trắc biến dạng. NXB Đại học Vũ Hán (Bản dịch của PGS. TS Phan Văn Hiến. NXB Khoa học và kỹ thuật, 2010). [11]. Trần Khánh, Nguyễn Quang Phúc (2010). Quan trắc chuyển dịch và biến dạng công trình. NXB Giao thông vận tải. [12]. Tao Benzao (2017). Bình sai lưới tự do và phân tích biến dạng. Phan Văn Hiến và Phạm Quốc Khánh dịch. NXB Tài nguyên - Môi trường và Bản đồ Việt Nam. BBT nhận bài: 11/01/2018; Phản biện xong: 20/02/2018 KHẢO SÁT ỨNG DỤNG CỦA THUẬT TOÁN BJERHAMMAR... (tiếp theo trang 81)