Phương pháp collocation bình phương nhỏ nhất khác gì so với phương pháp bình sai bình phương nhỏ nhất truyền thống?

Bình sai bình phương nhỏ nhất cổ điển (Gauss-Markov) chỉ xử lý các tham số xác định có tính quy luật hình học và coi toàn bộ phần dư là sai số ngẫu nhiên không tương quan. Ngược lại, collocation bình phương nhỏ nhất phân rã phần dư thành hai thành phần riêng biệt gồm tín hiệu trường vật lý ngẫu nhiên có tương quan không gian và nhiễu đo ngẫu nhiên, cho phép vừa bình sai vừa nội suy trường thế tối ưu.

Tại sao hàm hiệp phương sai lại quan trọng hàng đầu trong collocation trắc địa vật lý?

Hàm hiệp phương sai xác định trực tiếp trọng số tương tác không gian và mối liên kết toán học giữa các điểm dữ liệu. Trong trắc địa vật lý, hàm hiệp phương sai bắt buộc phải thỏa mãn phương trình vi phân Laplace của trường thế, đảm bảo tính nhất quán vật lý giữa các đại lượng đo lường khác nhau như dị thường trọng lực, độ lệch dây dọi và độ cao geoid.

Phương pháp collocation bình phương nhỏ nhất được ứng dụng như thế nào tại Việt Nam?

Tại Việt Nam, phương pháp collocation bình phương nhỏ nhất được sử dụng chủ yếu trong quy trình loại bỏ và khôi phục (Remove-Compute-Restore) để xây dựng mô hình quasigeoid độ phân giải cao, tích hợp đồng thời dữ liệu trọng lực mặt đất, trọng lực hàng không và đo cao vệ tinh trên Biển Đông với độ chính xác đạt dưới 10 cm.

Phương pháp collocation bình phương nhỏ nhất là gì?

Phương pháp collocation bình phương nhỏ nhất là kỹ thuật ước lượng tuyến tính tối ưu trong trắc địa vật lý và địa vật lý, cho phép kết hợp đồng thời việc bình sai các tham số hệ thống, lọc tín hiệu ngẫu nhiên và khử nhiễu ngẫu nhiên từ các nguồn dữ liệu đo lường dị thường. Phương pháp này thiết lập lời giải giải tích thông qua không gian Hilbert tái tạo hạt nhân nhằm xác định các đại lượng của trường trọng lực Trái Đất như độ cao geoid, dị thường trọng lực hay độ lệch dây dọi. Bài viết này trình bày chi tiết về nguồn gốc lý thuyết, cấu trúc mô hình toán học, phương pháp xây dựng hàm hiệp phương sai, mối liên hệ tương đương với các kỹ thuật nội suy thống kê hiện đại và những ứng dụng thực tiễn nổi bật trong xây dựng mô hình quasigeoid độ phân giải cao tại Việt Nam.

Nguồn gốc lịch sử và bối cảnh hình thành lý thuyết

Lịch sử phát triển của phương pháp collocation bình phương nhỏ nhất bắt đầu từ thời điểm nhà toán học và trắc địa Torben Krarup đặt nền móng lý thuyết đầu tiên vào năm 1969, thông qua việc áp dụng không gian Hilbert với hạt nhân tái sinh vào việc biểu diễn thế trọng trường trong công trình nghiên cứu về nền tảng toán học của trắc địa vật lý. Tiếp nối công trình tiên phong của Krarup, Helmut Moritz đã hệ thống hóa và phát triển phương pháp này thành một lý thuyết tổng quát cho toàn ngành trắc địa.

Theo Moritz (1978), phương pháp collocation bình phương nhỏ nhất kết hợp đồng thời ba bài toán kinh điển của trắc địa: bình sai các tham số xu hướng có tính quy luật, lọc tín hiệu trường thực tế và dự báo nội suy các giá trị chưa biết tại các vị trí không có trạm đo. Sự kết hợp này mang lại một khung toán học thống nhất để xử lý các tập dữ liệu không đồng nhất về bản chất vật lý lẫn độ chính xác đo lường.

Trước khi lý thuyết collocation ra đời, các bài toán trắc địa thường được giải quyết tách biệt thông qua phương pháp bình sai gián tiếp Gauss-Markov cho các tham số hình học hoặc thông qua các công thức tích phân biên kinh điển của Stokes và Vening Meinesz cho bài toán thế. Tuy nhiên, các công thức tích phân biên đòi hỏi phải phủ kín dữ liệu trọng lực trên toàn cầu với mật độ dày đặc. Collocation bình phương nhỏ nhất đã giải quyết triệt để nút thắt này bằng cách cho phép tính toán trực tiếp từ các tập điểm đo phân tán không đồng đều.

Mô hình toán học tổng quát và phương trình quan trắc

Theo chuyên khảo của Hofmann-Wellenhof và Moritz (2006), mô hình toán học tổng quát của collocation bình phương nhỏ nhất có dạng phương trình tuyến tính trong không gian định chuẩn:

\mathbf{x} = \mathbf{A}\mathbf{X} + \mathbf{s} + \mathbf{n}

Trong cấu trúc phương trình trên, các đại lượng được định nghĩa chặt chẽ như sau:

$\mathbf{x}$ là vector chứa các giá trị quan trắc thực tế thu thập từ thực địa, ví dụ như các giá trị đo dị thường trọng lực, thế trọng trường hay độ lệch dây dọi.
$\mathbf{A}$ là ma trận thiết kế (design matrix) biểu diễn mối quan hệ hàm số giữa các đại lượng quan trắc và các tham số xu hướng hệ thống.
$\mathbf{X}$ là vector các tham số chưa biết cần bình sai, đại diện cho phần xu hướng xác định (deterministic trend).
$\mathbf{s}$ là vector tín hiệu ngẫu nhiên (signal) tại chính các điểm quan trắc, mang bản chất là trường ngẫu nhiên liên tục có kỳ vọng toán bằng 0.
$\mathbf{n}$ là vector nhiễu ngẫu nhiên của phép đo (noise), mang kỳ vọng toán bằng 0 và không tương quan với tín hiệu trường.

Mục tiêu cốt lõi của bài toán là tìm kiếm vector tham số ước lượng $\hat{\mathbf{X}}$ và dự báo tín hiệu $\hat{s}_P$ tại một điểm tính toán $P$ bất kỳ trong không gian, sao cho phương sai của sai số dự báo đạt giá trị cực tiểu. Theo Hofmann-Wellenhof và Moritz (2006), ước lượng tối ưu tuyến tính không chệch của vector tham số hệ thống được xác định theo nguyên lý bình phương nhỏ nhất tổng quát:

\hat{\mathbf{X}} = (\mathbf{A}^T \bar{\mathbf{C}}^{-1} \mathbf{A})^{-1} \mathbf{A}^T \bar{\mathbf{C}}^{-1} \mathbf{x}

Trong đó ma trận $\bar{\mathbf{C}}$ biểu diễn tổng ma trận hiệp phương sai của toàn bộ vector quan trắc, được phân rã thành hai thành phần độc lập:

\bar{\mathbf{C}} = \mathbf{C}_{ss} + \mathbf{C}_{nn}

Tại đây, $\mathbf{C}_{ss}$ là ma trận tự hiệp phương sai giữa các điểm tín hiệu quan trắc, còn $\mathbf{C}_{nn}$ là ma trận phương sai - hiệp phương sai của sai số đo lường. Sau khi đã tách biệt phần xu hướng hệ thống, vector tín hiệu dự báo $\hat{s}_P$ tại một vị trí bất kỳ $P$ được tính toán thông qua mối quan hệ:

\hat{s}_P = \mathbf{C}_{Ps} \bar{\mathbf{C}}^{-1} (\mathbf{x} - \mathbf{A}\hat{\mathbf{X}})

Trong biểu thức này, $\mathbf{C}_{Ps}$ là vector hiệp phương sai chéo giữa tín hiệu tại điểm cần dự báo $P$ và các tín hiệu tại những điểm quan trắc đã có. Ưu điểm nổi bật của công thức này là không chỉ cho ra giá trị dự báo tối ưu mà còn cung cấp kèm theo phương sai sai số hậu nghiệm của giá trị tính toán, giúp đánh giá độ tin cậy định lượng của kết quả tại từng vị trí không gian.

Hàm hiệp phương sai và tính chất toán học trong trường thế

Theo Darbeheshti và Featherstone (2008), hàm hiệp phương sai đóng vai trò cốt lõi trong collocation bình phương nhỏ nhất vì nó quy định toàn bộ trọng số tương tác không gian giữa các điểm dữ liệu. Nếu hàm hiệp phương sai không phản ánh đúng đặc tính vật lý của trường thế, kết quả nội suy sẽ bị chệch hoặc mất tính tối ưu thống kê.

Trong trắc địa vật lý, thế trọng trường nhiễu loạn $T$ thỏa mãn phương trình điều hòa vi phân Laplace tại vùng không gian bên ngoài vỏ Trái Đất:

\nabla^2 T = \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2} = 0

Do thế nhiễu loạn là một hàm điều hòa, hàm tự hiệp phương sai của nó cũng phải tuân thủ điều kiện hàm điều hòa không gian đối với cả hai vị trí điểm tính toán. Dựa trên tính chất đối xứng cầu toàn cầu, hàm hiệp phương sai của thế trọng trường thường được khai triển theo chuỗi đa thức Legendre thông qua các phương sai bậc:

K(P, Q) = \sum_{n=2}^{\infty} \sigma_n \left( \frac{R_B^2}{r_P r_Q} \right)^{n+1} P_n(\cos \psi)

Trong đó $R_B$ là bán kính mặt cầu Bjerhammar nằm hoàn toàn bên trong Trái Đất, $r_P$ và $r_Q$ là khoảng cách từ tâm Trái Đất tới hai điểm $P$ và $Q$ , $\psi$ là góc cầu giữa hai điểm, $P_n(\cos \psi)$ là đa thức Legendre bậc $n$ và $\sigma_n$ là phương sai bậc của thế nhiễu loạn. Các mô hình giải tích nổi tiếng như mô hình Tscherning-Rapp hay mô hình giải tích của Moritz cho phép tính toán giải tích các đạo hàm riêng của thế, tạo điều kiện tính chuyển hiệp phương sai trực tiếp giữa các đại lượng vật lý khác nhau.

Tuy nhiên, trên phạm vi khu vực cục bộ có địa hình phức tạp, giả thiết về tính đẳng hướng và tính dừng của trường trọng lực thường không hoàn toàn thỏa mãn. Theo nghiên cứu của Darbeheshti và Featherstone (2008), việc áp dụng mô hình hàm hiệp phương sai phi dừng trong collocation 2 chiều cho phép thích ứng linh hoạt với sự thay đổi địa chất và độ dốc địa hình, giúp giảm thiểu đáng kể sai số nội suy ở các vùng núi cao.

Mối tương quan với Kriging và các mô hình học máy hiện đại

Một góc nhìn lý thuyết thú vị trong toán học ứng dụng là sự đồng quy của nhiều phương pháp xuất phát từ các ngành khoa học độc lập. Collocation bình phương nhỏ nhất trong trắc địa vật lý có cấu trúc toán học hoàn toàn tương đồng với phương pháp Kriging vạn năng (Universal Kriging) do Georges Matheron phát triển trong địa thống kê mỏ khoáng, cũng như phương pháp dự báo tuyến tính không chệch tốt nhất (BLUP) trong thống kê ứng dụng.

Xét trong bối cảnh trí tuệ nhân tạo và khoa học dữ liệu đương đại, collocation bình phương nhỏ nhất chính là một trường hợp chuyên biệt của hồi quy quá trình Gaussian (Gaussian Process Regression). Cả hai phương pháp đều đặt giả thiết trường dữ liệu là một quá trình ngẫu nhiên dừng với phân phối chuẩn nhiều chiều, sử dụng hàm nhân hiệp phương sai (covariance kernel) để ánh xạ không gian và tối ưu hóa hàm hợp lý cực đại để xác định các siêu tham số.

Điểm khác biệt căn bản làm nên sức mạnh độc tôn của collocation trong khoa học Trái Đất nằm ở nguyên lý tương thích vật lý trường thế. Trong khi Kriging và học máy thông thường sử dụng các hàm nhân mang tính kinh nghiệm thống kê thuần túy (như hàm nhân Gaussian, Matérn hay đa thức bậc hai), collocation bình phương nhỏ nhất bắt buộc hàm nhân phải thỏa mãn các định luật vật lý của cơ học thiên thể và phương trình Laplace, nhờ đó bảo toàn tính nhất quán giữa thế năng, gia tốc trọng trường và hướng của dây dọi.

So sánh phương pháp collocation bình phương nhỏ nhất với các kỹ thuật khác

Để thấy rõ ưu thế và phạm vi thích hợp của phương pháp collocation bình phương nhỏ nhất, bảng dưới đây so sánh phương pháp này với các kỹ thuật xấp xỉ trường trọng lực và nội suy không gian phổ biến:

Tiêu chí so sánh	Collocation bình phương nhỏ nhất (LSC)	Tích phân Stokes biến đổi (FFT)	Kriging địa thống kê
Bản chất toán học	Ước lượng tuyến tính tối ưu trong không gian Hilbert	Tích phân số giải tích bài toán biên	Nội suy tối ưu dựa trên biến thiên hàm variogram
Khả năng xử lý dữ liệu dị thường	Tích hợp đồng thời nhiều loại dữ liệu khác nhau	Chỉ xử lý trực tiếp một loại dữ liệu (dị thường trọng lực)	Chủ yếu xử lý dữ liệu cùng đại lượng vật lý
Phân bố không gian của dữ liệu	Hoạt động trên lưới điểm tùy ý, thưa hoặc dày	Bắt buộc dữ liệu phải trải đều trên lưới tọa độ đều	Hoạt động trên tập điểm rời rạc tùy ý
Mô hình hóa sai số đo lường	Khử nhiễu tự động thông qua ma trận hiệp phương sai sai số	Không có cơ chế tách nhiễu trực tiếp trong tích phân	Khử nhiễu đo thông qua hiệu ứng nugget
Ràng buộc cơ sở vật lý	Thỏa mãn nghiêm ngặt phương trình Laplace của trường thế	Thỏa mãn bài toán biên vật lý toàn cầu	Chỉ mang tính chất kinh nghiệm thống kê hình học
Chi phí tính toán ma trận	Rất lớn do phải đảo ma trận toàn phương kích thước N x N	Rất nhanh nhờ thuật toán biến đổi Fourier nhanh (FFT)	Trung bình đến lớn tùy thuộc số lượng điểm lân cận

Ứng dụng thực tiễn trong xác định mô hình geoid tại Việt Nam

Trong thực tiễn hiện đại hóa hệ thống độ cao quốc gia, phương pháp collocation bình phương nhỏ nhất đóng vai trò nền tảng để xây dựng mô hình bề mặt quasigeoid địa phương. Bề mặt quasigeoid là bề mặt chuẩn phục vụ việc chuyển đổi trực tiếp độ cao trắc địa đo bằng vệ tinh GNSS sang độ cao thủy chuẩn phục vụ các công trình xây dựng hạ tầng, giao thông và thủy lợi.

Kỹ thuật tính toán chuẩn mực thường áp dụng quy trình loại bỏ và khôi phục (Remove-Compute-Restore). Trong quy trình này, phần bước sóng dài của trường trọng lực được loại bỏ bằng một mô hình trọng trường toàn cầu, phần hiệu ứng địa hình cục bộ được loại bỏ bằng mô hình độ cao số. Thành phần trường dư còn lại có tính dừng cao, rất thuận lợi cho việc áp dụng phương pháp collocation bình phương nhỏ nhất để nội suy lên lưới tọa độ chi tiết trước khi khôi phục lại các thành phần đã loại bỏ.

Theo nghiên cứu của Vu và cộng sự (2019), mô hình quasigeoid trọng lực phân giải cao của Việt Nam xây dựng bằng phương pháp collocation bình phương nhỏ nhất đạt độ lệch chuẩn 9,7 cm khi kiểm chứng độc lập trên 812 điểm độ cao GNSS thủy chuẩn phân bố trên toàn lãnh thổ quốc gia. Độ chính xác vượt trội này khẳng định năng lực xử lý vượt trội của phương pháp collocation khi kết hợp các nguồn dữ liệu trọng lực đo mặt đất, đo hàng không và dữ liệu đo cao vệ tinh trên vùng biển Đông.

Những ưu thế này mở ra tiềm năng to lớn trong việc ứng dụng collocation bình phương nhỏ nhất để theo dõi biến dạng vỏ Trái Đất, phân tích sụt lún mặt đất tại vùng đồng bằng sông Cửu Long và xây dựng cơ sở dữ liệu trọng lực biển phục vụ bảo đảm an ninh hàng hải và nghiên cứu địa chất dầu khí trong tương lai.