球坐標多面函數對大區域數據插值的應用研究

孟 濤,孟祥恩,曾現敏,馬治忠*,吳嘯龍

（1. 國家海洋局北海海洋技術保障中心，山東 青島 266033；2. 西安科技大學，陜西 西安 710054）

目前大地測量學領域主要采用的觀測手段是空間大地測量觀測技術，該技術打破了以往經典大地測量的時空局限性，從而實現精確、直接、全面的在超大尺度乃至全球尺度上的大地測量監測數據采集，而研究基于球坐標多面函數的空間數據插值方法也是為了適應這一技術帶來的變化[1-4]。

1 球坐標多面函數觀測方程

多面函數在笛卡爾坐標系中的數學形式可表示為[5]：

式中，n為參與多面函數建模的已知節點數；Q(X,Y,Xj,Yj)為多面函數模型的核函數；αj為多面函數待定系數。

一種常用多面函數的雙曲型核函數形式為：

上式中光滑因子δ在控制著曲面平緩程度的同時，在其幾何意義上還表示這曲面上最小值點到X-Y坐標平面的距離，即δ2=(z-zj)2。在幾何上該參數也可以視為不同半徑的兩個同心球層在同一點處的半徑差異[6]。因此，在笛卡爾坐標系下，結合笛卡爾平面坐標和球坐標之間的轉換關系，式（2）可以寫為如下形式：

其中，

式中，R為地球球狀模型的平均曲率半徑；r 為最優內部球層半徑常數；θ為余緯度；λ為經度。

式（3）中的r 為針對某一研究區域（范圍）最優的內部球層常半徑，Hardy 等給出了最優半徑迭代公式為：

式（6），（7）分別對應為全球尺度下、局部地區的第一輔助方程；式（8）為第二輔助方程；n為研究區域內用于多面函數建模的已知觀測數據個數；A為用于建立多面函數模型的已知數據所覆蓋區域的球面面積。由式（5）～（8）可知，當給定分布在球面上的建模數據點的個數n以及地球平均曲率半徑后，可以得到針對當前已知數據分布情況的最優r 解析式。然而，式（5）為無顯式解的隱函數，需要根據二分法求解當前已知結點分布水平對應的內球層最優半徑。

經過坐標變換后，基本的多面函數模型仍可以表達為線性結構。在實驗范圍內，對n個實測點選取m個點數據（m≤n），構建出多面函數模型式（9），其在最小二乘算法下的最優解如式（10）～（13）所示[7～9]：

假定已知n個采樣點上的地殼垂直運動速率和位置，并確定了核函數的具體形式，則速率面擬合就是依據n個采樣點的速率觀測值求解一組待定參數。

當m=n時，

2 算例與分析

2.1 算例1

本算例旨在驗證球坐標系下多面函數模型對于大、中尺度空間離散數據插值應用的可行性及其精確程度。研究范圍選取為120°～135°E，18°～28°N 的843 個已知點大地水準面高程數據。選取其中分布均勻的94個點作為外部檢查點，用于評定所建立模型的擬合精度。研究分別采用不同數目的已知數據構建數學模型，計算并評價球坐標系和笛卡爾坐標系下多面函數的擬合精度；球坐標多面函數與其他空間數據插值算法的擬合精度。

外部檢查點均方根誤差計算公式[10-11]為：

式中，n為檢核數據個數；XK為待檢核點位置的已知量；XK為待檢核點位置的模型擬合結果。

本文分別給出了起算結點個數為843,281,211,169,141,121,106時，在球坐標系與笛卡爾坐標系的多面函數中，對外部檢核點的均方根誤差分別進行計算，結果走勢如圖1 所示。三角形符號表示球坐標多面函數的建模結果，星形符號表示笛卡爾平面坐標多面函數的建模結果。從圖中可以看出，隨著用于建立多面函數模型的已知結點數目的減少，多面函數計算外部檢核點的精度逐漸變差，其均方根誤差在球坐標多面函數中建模結果的變化幅度要顯著小于笛卡爾坐標系的結果。

圖1 不同坐標系統的多面函數建模均方根誤差變化趨勢對比

不同方案建模結果精度評定的詳細統計結果見表1，結合圖1 所示趨勢與表1 中的擬合殘差統計結果可知，多面函數在球坐標系下相較于笛卡爾坐標系的建模結果，對于離散點的穩定性以及擬合精度更優。

表1 不同方案擬合殘差統計/m

為了進一步研究大、中空間尺度下球坐標系多面函數的插值精度，本文對同一組已知數據，分別采用局部多項式法、最鄰近點法、克里格法、改進謝別德法等算法對研究區域內的250 個外部檢核點進行建模擬合[12-14]，其擬合殘差分布如圖2所示，紅色曲線為球坐標多面函數對外部檢核點的擬合殘差分布。由此可知，在研究區域內部，球多面函數對于空間數據的插值精度達到甚至優于一些常規算法。

圖2 不同算法的外部檢核點擬合殘差對比

2.2 算例2

水準數據是典型空間大地測量離散的觀測結果，采用動態平差方法獲得的區域地表垂直運動信息在時間域、空間域方面是不連續的[15-16]。實際上，地面垂直運動表現為地表上連續形狀變化的曲面，各水準點則是這個曲面上的離散采樣結果。為實現在時間域、空間域中的連續性垂直變化信息，多面函數被廣泛應用于給予離散數據獲取連續形變速率面的研究當中。

然而，兩點高差為水準測量的直接觀測量，利用水準觀測量平差得到的水準速率僅有相對意義，所以必須確定參考基準。不同參考基準的選擇可能獲得完全不同趨勢的形變速率場。因此，在根據水準數據判斷區域地面變化的穩定性時，常常引入垂直速率梯度。垂直運動高梯度帶及速率等值線的急劇轉折部位往往是地面構造薄弱、易失穩區域。

梯度作為一個向量，依據定義，地殼垂直運動速率梯度應為相應坐標分量的偏導數。黃立人等曾基于笛卡爾坐標系下多面函數模型給出了垂直速率在N方向和E方向的偏導數分量計算模型，并針對地震高危區域進行了深入研究，取得了大量有益結果。然而，地面垂直形變速率場的徑向梯度對于區域地面構造穩定性的判定、指向特性還有待進一步研究。

根據Hardy與Gopfert[4]，任何由式（3）、（4）進一步推導得到的函數關系仍然為線性模型。在利用球多面函數理論構建區域垂直速率面模型的基礎上，可以進一步求解其速率在徑向的梯度分量。假設式（3）表示垂直形變速率面，則其對徑向坐標分量R的偏導數可以按照下式進行計算：

式中，先驗參數r和ψj可以由式（4）、（5）進行計算。

為了便于驗證該模型的客觀性與正確性性，本文選取分布在福建省的水準垂直速率場數據資料，利用式（16）計算出研究區域地面垂直形變速率徑向梯度分布，其結果如圖3所示。

圖3 徑向梯度分布

從圖3 中可以看出，地面垂直形變速率場的徑向梯度分布規律與區域內的重大活動構造（地殼的薄弱區域）具有良好的相關性。福建省東南沿海地殼地震活動主要受控于北東向的斷裂系，及與其正交的北西斷裂系的共同作用[17]。從徑向梯度帶的密集變化區域分布可以看出，高梯度帶地區與歷史地震分布具有良好的空間相關性。也就是說，斷裂帶的交匯部位，徑向梯度分布密集變化劇烈，對應的地震時間分布頻度高。

3 結 語

根據本文研究，可以得到以下幾點結論：

1）從正雙曲型核函數的幾何意義出發，同時控制著正雙曲面的平緩程度以及曲面上所有點到x-y坐標平面的最小距離。因此，可以將視為兩點z坐標分量的差異，從而最終獲得球坐標系下多面函數的基本模型。

2）為了進一步分析球坐標多面函數對大、中尺度上的球多面函數的擬合精度，本文選取分布在我國東南部范圍為120°～135°E，18°～28°N的大地水準面離散數據進行插值精度分析。通過計算球坐標多面函數與笛卡爾平面坐標多面函數在不同建模方案下的整體擬合殘差以及均方根誤差，本文研究認為針對較大矩形區域內的離散數據，采用球坐標系多面函數的得到的數據插值精度與穩定性均優于相比于笛卡爾坐標系的計算結果。另一方面，本文進一步對比了球多面函數與幾種常見空間數據插值模型的插值精度，結果顯示本文模型相對于其他算法，精確度與穩定性都較好。

3）由于水準測量結果僅具有相對意義，難以單純地依靠速率值大小來衡量地殼運動的穩定性。因此常將地殼運動的垂直運動高梯度帶，及其等值線的急劇轉折部位作為未來區域地面差異性變化或構造應力顯著作用區的判斷依據。本文在球坐標系多面函數擬合垂直速率面的基礎上，進一步推導了垂直速率的徑向梯度分量。并選取福建省東南沿海地區的垂直速率觀測數據進行徑向梯度分布的解算。結果顯示判定垂直速率的徑向梯度分布與地殼構造薄弱、易失穩的區域就有良好的空間對應關系。因此，該研究方法在理論上同樣適用于礦山采空區地面沉降監測及相關地質災害、危險區區劃等工作的研究。