孟 濤,孟祥恩,曾現(xiàn)敏,馬治忠*,吳嘯龍

（1. 國(guó)家海洋局北海海洋技術(shù)保障中心，山東 青島 266033；2. 西安科技大學(xué)，陜西 西安 710054）

目前大地測(cè)量學(xué)領(lǐng)域主要采用的觀測(cè)手段是空間大地測(cè)量觀測(cè)技術(shù)，該技術(shù)打破了以往經(jīng)典大地測(cè)量的時(shí)空局限性，從而實(shí)現(xiàn)精確、直接、全面的在超大尺度乃至全球尺度上的大地測(cè)量監(jiān)測(cè)數(shù)據(jù)采集，而研究基于球坐標(biāo)多面函數(shù)的空間數(shù)據(jù)插值方法也是為了適應(yīng)這一技術(shù)帶來的變化[1-4]。

1 球坐標(biāo)多面函數(shù)觀測(cè)方程

多面函數(shù)在笛卡爾坐標(biāo)系中的數(shù)學(xué)形式可表示為[5]：

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

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

上式中光滑因子δ在控制著曲面平緩程度的同時(shí)，在其幾何意義上還表示這曲面上最小值點(diǎn)到X-Y坐標(biāo)平面的距離，即δ2=(z-zj)2。在幾何上該參數(shù)也可以視為不同半徑的兩個(gè)同心球?qū)釉谕稽c(diǎn)處的半徑差異[6]。因此，在笛卡爾坐標(biāo)系下，結(jié)合笛卡爾平面坐標(biāo)和球坐標(biāo)之間的轉(zhuǎn)換關(guān)系，式（2）可以寫為如下形式：

其中，

式中，R為地球球狀模型的平均曲率半徑；r 為最優(yōu)內(nèi)部球?qū)影霃匠?shù)；θ為余緯度；λ為經(jīng)度。

式（3）中的r 為針對(duì)某一研究區(qū)域（范圍）最優(yōu)的內(nèi)部球?qū)映０霃?，Hardy 等給出了最優(yōu)半徑迭代公式為：

式（6），（7）分別對(duì)應(yīng)為全球尺度下、局部地區(qū)的第一輔助方程；式（8）為第二輔助方程；n為研究區(qū)域內(nèi)用于多面函數(shù)建模的已知觀測(cè)數(shù)據(jù)個(gè)數(shù)；A為用于建立多面函數(shù)模型的已知數(shù)據(jù)所覆蓋區(qū)域的球面面積。由式（5）～（8）可知，當(dāng)給定分布在球面上的建模數(shù)據(jù)點(diǎn)的個(gè)數(shù)n以及地球平均曲率半徑后，可以得到針對(duì)當(dāng)前已知數(shù)據(jù)分布情況的最優(yōu)r 解析式。然而，式（5）為無顯式解的隱函數(shù)，需要根據(jù)二分法求解當(dāng)前已知結(jié)點(diǎn)分布水平對(duì)應(yīng)的內(nèi)球?qū)幼顑?yōu)半徑。

經(jīng)過坐標(biāo)變換后，基本的多面函數(shù)模型仍可以表達(dá)為線性結(jié)構(gòu)。在實(shí)驗(yàn)范圍內(nèi)，對(duì)n個(gè)實(shí)測(cè)點(diǎn)選取m個(gè)點(diǎn)數(shù)據(jù)（m≤n），構(gòu)建出多面函數(shù)模型式（9），其在最小二乘算法下的最優(yōu)解如式（10）～（13）所示[7～9]：

假定已知n個(gè)采樣點(diǎn)上的地殼垂直運(yùn)動(dòng)速率和位置，并確定了核函數(shù)的具體形式，則速率面擬合就是依據(jù)n個(gè)采樣點(diǎn)的速率觀測(cè)值求解一組待定參數(shù)。

當(dāng)m=n時(shí)，

2 算例與分析

2.1 算例1

本算例旨在驗(yàn)證球坐標(biāo)系下多面函數(shù)模型對(duì)于大、中尺度空間離散數(shù)據(jù)插值應(yīng)用的可行性及其精確程度。研究范圍選取為120°～135°E，18°～28°N 的843 個(gè)已知點(diǎn)大地水準(zhǔn)面高程數(shù)據(jù)。選取其中分布均勻的94個(gè)點(diǎn)作為外部檢查點(diǎn)，用于評(píng)定所建立模型的擬合精度。研究分別采用不同數(shù)目的已知數(shù)據(jù)構(gòu)建數(shù)學(xué)模型，計(jì)算并評(píng)價(jià)球坐標(biāo)系和笛卡爾坐標(biāo)系下多面函數(shù)的擬合精度；球坐標(biāo)多面函數(shù)與其他空間數(shù)據(jù)插值算法的擬合精度。

外部檢查點(diǎn)均方根誤差計(jì)算公式[10-11]為：

式中，n為檢核數(shù)據(jù)個(gè)數(shù)；XK為待檢核點(diǎn)位置的已知量；XK為待檢核點(diǎn)位置的模型擬合結(jié)果。

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

圖1 不同坐標(biāo)系統(tǒng)的多面函數(shù)建模均方根誤差變化趨勢(shì)對(duì)比

不同方案建模結(jié)果精度評(píng)定的詳細(xì)統(tǒng)計(jì)結(jié)果見表1，結(jié)合圖1 所示趨勢(shì)與表1 中的擬合殘差統(tǒng)計(jì)結(jié)果可知，多面函數(shù)在球坐標(biāo)系下相較于笛卡爾坐標(biāo)系的建模結(jié)果，對(duì)于離散點(diǎn)的穩(wěn)定性以及擬合精度更優(yōu)。

表1 不同方案擬合殘差統(tǒng)計(jì)/m

為了進(jìn)一步研究大、中空間尺度下球坐標(biāo)系多面函數(shù)的插值精度，本文對(duì)同一組已知數(shù)據(jù)，分別采用局部多項(xiàng)式法、最鄰近點(diǎn)法、克里格法、改進(jìn)謝別德法等算法對(duì)研究區(qū)域內(nèi)的250 個(gè)外部檢核點(diǎn)進(jìn)行建模擬合[12-14]，其擬合殘差分布如圖2所示，紅色曲線為球坐標(biāo)多面函數(shù)對(duì)外部檢核點(diǎn)的擬合殘差分布。由此可知，在研究區(qū)域內(nèi)部，球多面函數(shù)對(duì)于空間數(shù)據(jù)的插值精度達(dá)到甚至優(yōu)于一些常規(guī)算法。

圖2 不同算法的外部檢核點(diǎn)擬合殘差對(duì)比

2.2 算例2

水準(zhǔn)數(shù)據(jù)是典型空間大地測(cè)量離散的觀測(cè)結(jié)果，采用動(dòng)態(tài)平差方法獲得的區(qū)域地表垂直運(yùn)動(dòng)信息在時(shí)間域、空間域方面是不連續(xù)的[15-16]。實(shí)際上，地面垂直運(yùn)動(dòng)表現(xiàn)為地表上連續(xù)形狀變化的曲面，各水準(zhǔn)點(diǎn)則是這個(gè)曲面上的離散采樣結(jié)果。為實(shí)現(xiàn)在時(shí)間域、空間域中的連續(xù)性垂直變化信息，多面函數(shù)被廣泛應(yīng)用于給予離散數(shù)據(jù)獲取連續(xù)形變速率面的研究當(dāng)中。

然而，兩點(diǎn)高差為水準(zhǔn)測(cè)量的直接觀測(cè)量，利用水準(zhǔn)觀測(cè)量平差得到的水準(zhǔn)速率僅有相對(duì)意義，所以必須確定參考基準(zhǔn)。不同參考基準(zhǔn)的選擇可能獲得完全不同趨勢(shì)的形變速率場(chǎng)。因此，在根據(jù)水準(zhǔn)數(shù)據(jù)判斷區(qū)域地面變化的穩(wěn)定性時(shí)，常常引入垂直速率梯度。垂直運(yùn)動(dòng)高梯度帶及速率等值線的急劇轉(zhuǎn)折部位往往是地面構(gòu)造薄弱、易失穩(wěn)區(qū)域。

梯度作為一個(gè)向量，依據(jù)定義，地殼垂直運(yùn)動(dòng)速率梯度應(yīng)為相應(yīng)坐標(biāo)分量的偏導(dǎo)數(shù)。黃立人等曾基于笛卡爾坐標(biāo)系下多面函數(shù)模型給出了垂直速率在N方向和E方向的偏導(dǎo)數(shù)分量計(jì)算模型，并針對(duì)地震高危區(qū)域進(jìn)行了深入研究，取得了大量有益結(jié)果。然而，地面垂直形變速率場(chǎng)的徑向梯度對(duì)于區(qū)域地面構(gòu)造穩(wěn)定性的判定、指向特性還有待進(jìn)一步研究。

根據(jù)Hardy與Gopfert[4]，任何由式（3）、（4）進(jìn)一步推導(dǎo)得到的函數(shù)關(guān)系仍然為線性模型。在利用球多面函數(shù)理論構(gòu)建區(qū)域垂直速率面模型的基礎(chǔ)上，可以進(jìn)一步求解其速率在徑向的梯度分量。假設(shè)式（3）表示垂直形變速率面，則其對(duì)徑向坐標(biāo)分量R的偏導(dǎo)數(shù)可以按照下式進(jìn)行計(jì)算：

式中，先驗(yàn)參數(shù)r和ψj可以由式（4）、（5）進(jìn)行計(jì)算。

為了便于驗(yàn)證該模型的客觀性與正確性性，本文選取分布在福建省的水準(zhǔn)垂直速率場(chǎng)數(shù)據(jù)資料，利用式（16）計(jì)算出研究區(qū)域地面垂直形變速率徑向梯度分布，其結(jié)果如圖3所示。

圖3 徑向梯度分布

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

3 結(jié) 語

根據(jù)本文研究，可以得到以下幾點(diǎn)結(jié)論：

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

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

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