李 云 華 宇 燕保榮 郭 偉

（1.中國科學(xué)院國家授時(shí)中心，陜西西安710600；2.中國科學(xué)院大學(xué)，北京100049）

1 引 言

2018年，美國總統(tǒng)特朗普簽署了《國家安全與彈性授時(shí)法案》。該法案要求交通部在兩年內(nèi)建設(shè)針對(duì)GPS 的地基備用授時(shí)系統(tǒng)。這在一定程度上宣布了美國的重啟羅蘭計(jì)劃，同時(shí)提高羅蘭系統(tǒng)的精度使其與GPS 系統(tǒng)相當(dāng)。繼BPL 和長河二號(hào)系統(tǒng)的eLoran 升級(jí)改造之后，我國也開展了高精度地基授時(shí)方面的研究，計(jì)劃擴(kuò)建eLoran 發(fā)播臺(tái)滿足eLoran 信號(hào)的全國土覆蓋，同時(shí)要提高eLoran 系統(tǒng)的精度。

國內(nèi)外專家對(duì)ASF 網(wǎng)格進(jìn)行了一系列研究。2005年，Gregory William Johnson 等人通過大量測試數(shù)據(jù)分析二次時(shí)延的空間、時(shí)間和方向變化［3］。2012年，Hargreaves C 等人提出收集足夠測量數(shù)據(jù)做出高精度ASF 數(shù)據(jù)庫，能保證在海上定位結(jié)果優(yōu)于10m（95%）［4］。還有專家學(xué)者從節(jié)省ASF 網(wǎng)格數(shù)據(jù)庫的成本出發(fā)，分析用靜態(tài)數(shù)據(jù)和動(dòng)態(tài)數(shù)據(jù)建立ASF 網(wǎng)格兩種不同的方法［5～8］。

有了ASF 數(shù)據(jù)庫，用戶確立所在網(wǎng)格根據(jù)內(nèi)插算法解算出ASF 值，修正信號(hào)傳播時(shí)延，從而提高授時(shí)精度。在ASF 數(shù)據(jù)庫既定的情況下，內(nèi)插算法本身也會(huì)引入新的誤差，本文仿真分析不同空間插值算法帶來的誤差，并對(duì)算法進(jìn)行比較，得出雙線性插值算法誤差最小，反距離插值算法誤差最大。

2 eLoran 信號(hào)傳播時(shí)延

eLoran 信號(hào)從發(fā)射天線到接收天線所經(jīng)歷的時(shí)間稱為傳播時(shí)延Tp，如式（1）所示。假設(shè)信號(hào)在無限大空氣介質(zhì)中傳播時(shí)，由發(fā)射天線到接收天線所經(jīng)歷的時(shí)間稱為一次時(shí)延（Primary Factor，PF），又稱為基本時(shí)延；當(dāng)信號(hào)在海平面?zhèn)鞑r(shí)，傳播時(shí)延與一次時(shí)延之差，稱為二次時(shí)延（Secondary Factor，SF）反映了海平面對(duì)傳播時(shí)延的影響；當(dāng)信號(hào)在地平面?zhèn)鞑r(shí)，傳播時(shí)延與一次時(shí)延之差稱為附加二次時(shí)延（Added Secondary Factor，ASF）反映了地平面對(duì)傳播時(shí)延的影響［9］。

一次時(shí)延PF 計(jì)算公式如式（2）所示［10］

式中:C——真空中的光速，通常C=0.299 792 458km/μs；d——信號(hào)傳播的大地線距離；ns——地面的大氣折射指數(shù)，國際標(biāo)準(zhǔn)大氣值ns=1.000 315。

二次時(shí)延SF/附加二次時(shí)延ASF 的計(jì)算公式如式（3）所示

式中:ω——角速度，單位rad/s；W——衰減函數(shù)，與信號(hào)傳播距離、傳播路徑介質(zhì)的電導(dǎo)率、介電常數(shù)有關(guān)；argW——地波衰減函數(shù)W 的相位，單位rad；f——載波信號(hào)的頻率，100kHz；d——信號(hào)傳播的大圓距離，單位km；σ——信號(hào)傳播路徑上介質(zhì)的電導(dǎo)率；ε——信號(hào)傳播路徑上的介電常數(shù)。

3 空間插值算法

常用的空間插值算法，包括雙線性插值算法、反距離權(quán)值算法IDW（Inverse Distance Weighted，IDW）、Junkins 算法和自然鄰域插值算法等。這里重點(diǎn)對(duì)雙線性插值算法、反距離權(quán)值算法IDW、Junkins 算法進(jìn)行仿真分析。

3.1 雙線性插值算法

雙線性插值算法是在x，y兩個(gè)方向分別進(jìn)行一次線性插值，如圖1 所示。先x方向，后y方向插值，與先y方向，后在x方向插值，插值效果相同，其結(jié)果與插值順序無關(guān)［11］。

圖1 雙線性插值算法示意圖Fig.1 Bilinear interolation algorithm

四個(gè)頂點(diǎn)的坐標(biāo)分別為（x1，y1），（x2，y2），（x3，y3），（x4，y4），其ASF 分別為C1，C2，C3 和C4，待插值點(diǎn)（x，y）的ASF 如式（4）所示

其中，

3.2 反距離插值

反距離權(quán)重IDW 插值使用一組采樣點(diǎn)屬性的線性權(quán)重組合來確定待插值點(diǎn)的屬性，權(quán)重是一種反距離函數(shù)。這里反距離函數(shù)的冪值，取值為2。冪參數(shù)越高，可進(jìn)一步強(qiáng)調(diào)最近點(diǎn)對(duì)插值點(diǎn)的影響，因此，鄰近數(shù)據(jù)將受到更大影響，擬合的表面更加詳細(xì)（更不光滑），內(nèi)插值更接近相鄰采樣點(diǎn)的值。反之，冪參數(shù)越小，更遠(yuǎn)距離的采樣點(diǎn)對(duì)插值點(diǎn)的影響越大，擬合的表面更光滑。

IDW 插值方法假定每個(gè)輸入點(diǎn)都有著局部影響，這種影響隨著距離的增加而減弱。算法描述:現(xiàn)有n個(gè)離散點(diǎn)Z（x1，y1），Z（x2，y2）…Z（xn，yn），需要對(duì)插值點(diǎn)（x，y）進(jìn)行插值預(yù)測，求出Z（x，y），過程如下。

1）利用距離函數(shù)求出各離散點(diǎn)與插值點(diǎn)的距離，這里采用歐氏距離如式（5）所示

式中:（x，y）——插值點(diǎn)坐標(biāo)；（xi，yi）——已知點(diǎn)坐標(biāo)。

2）在冪參數(shù)確定的情況下，利用權(quán)重函數(shù)計(jì)算各個(gè)離散點(diǎn)的權(quán)重如式（6）所示

式中:n——離散點(diǎn)個(gè)數(shù)。

3）計(jì)算插值點(diǎn)的插值如式（7）所示

式中:Wi——第i采樣點(diǎn)的權(quán)重；Z（xi，yi）——第i點(diǎn)的屬性值。

3.3 Junkins 權(quán)值算法

Junkins 加權(quán)法是GPS 系統(tǒng)中用戶穿刺點(diǎn)電離層垂直延遲的計(jì)算方法。Junkins 加權(quán)法與反距離插值算法相比，計(jì)算量較小，與簡單的雙線性加權(quán)法相比，其空間相關(guān)性較強(qiáng)［12］，如圖2 所示。其中，1～4 分別指網(wǎng)格結(jié)點(diǎn)，其對(duì)應(yīng)的ASF 值分別為ASF1～ASF4；A 點(diǎn)為要擬合的點(diǎn)，其ASF 用ASF（φA，λA）表示為

圖2 Junkins 加權(quán)算法示意圖Fig.2 The Junkins weighting method

式中:W（x，y）——內(nèi)插的加權(quán)函數(shù)。

4 算法仿真分析

4.1 雙線性插值算法仿真

4.1.1 不同的網(wǎng)格大小

雙線性插值算法在不同大小網(wǎng)格內(nèi)的插值誤差分析如圖3 所示，從上到下依次是網(wǎng)格大小為0.5°，0.2°和0.1°，信號(hào)傳播路路徑是平均陸地電導(dǎo)率σ=0.003，ε=22 誤差范圍分別為（0～6）ns，（0～1）ns 和（0～0.2）ns，顯然網(wǎng)格越小，ASF 網(wǎng)格的精度越高。誤差最小的區(qū)域分布在網(wǎng)格的四個(gè)頂點(diǎn)附近，誤差最大的區(qū)域集中在網(wǎng)格的中心部位。因此相對(duì)來說離網(wǎng)格頂點(diǎn)越近，誤差越小，距離網(wǎng)格頂點(diǎn)越遠(yuǎn)誤差越大；距離網(wǎng)格中心越近誤差越大，距離網(wǎng)格中心越遠(yuǎn)，誤差越小。

4.1.2 不同的傳輸路徑

圖3 雙線性插值算法（不同的網(wǎng)格大?。〧ig.3 The bilinear interpolation algorithm （different mesh sizes）

雙線性插值算法應(yīng)用在海水路徑、平均陸地和甚干燥地路徑三種路徑上的0.1°網(wǎng)格中的插值結(jié)果的誤差如圖4 所示。海水介質(zhì)的誤差范圍在（0～0.1）ns 左右，平均陸地和甚干燥地在（0～0.2）ns 左右。隨著路徑的電導(dǎo)率變小，路徑的傳輸介質(zhì)對(duì)信號(hào)的阻力增強(qiáng)，同時(shí)ASF 網(wǎng)格的誤差也變大。因此海水路徑對(duì)信號(hào)的阻力最小，ASF 值最小，相應(yīng)的ASF 網(wǎng)格計(jì)算的誤差最??；隨著信號(hào)傳輸路徑上介質(zhì)對(duì)信號(hào)的阻力逐漸變大，ASF 值變大，相應(yīng)的ASF 網(wǎng)格計(jì)算結(jié)果的誤差逐漸變大。從網(wǎng)格中心處向外部誤差逐漸變小，呈現(xiàn)處出一定的層次感。

圖4 雙線性插值算法（不同傳輸路徑）Fig.4 The bilinear interpolation algorithm （different path）

4.2 反距離插值算法仿真

4.2.1 不同的網(wǎng)格大小

對(duì)反距離權(quán)值算法的仿真如圖5 所示，從上到下依次是網(wǎng)格大小為0.5°，0.2°和0.1°，誤差范圍分別為（0～40）ns，（0～15）ns 和（0～8）ns。顯然網(wǎng)格越小，ASF 網(wǎng)格能實(shí)現(xiàn)的精度越高。IDW 算法網(wǎng)格誤差分布較分散，網(wǎng)格中靠近頂點(diǎn)的分散區(qū)域誤差最小，同時(shí)網(wǎng)格的中心部分誤差最小，誤差最大的區(qū)域包括網(wǎng)格的側(cè)邊的幾個(gè)分散區(qū)域。

圖5 反距離加權(quán)插值算法（不同網(wǎng)格大?。〧ig.5 IDW interpolation algorithm （different mesh sizes）

4.2.2 不同的傳輸路徑

針對(duì)海水路徑、平均陸地和甚干燥地三種路徑上，網(wǎng)格大小0.1°時(shí)的誤差情況如圖6 所示。不同的信號(hào)傳輸路徑ASF 網(wǎng)格內(nèi)插結(jié)果的誤差:海水介質(zhì)的誤差范圍在（0～2.5）ns，平均陸地的誤差范圍在（0～6）ns 和甚干燥地的誤差范圍在（0～8）ns。從甚干燥地到海水隨著路徑的電導(dǎo)率變小，路徑的傳輸介質(zhì)對(duì)信號(hào)的阻力減小，同時(shí)ASF 網(wǎng)格的誤差也變小。因此甚干燥地路徑對(duì)信號(hào)的阻力最大，相應(yīng)的ASF網(wǎng)格計(jì)算的誤差最大，隨著信號(hào)傳輸路徑上介質(zhì)對(duì)信號(hào)的阻力逐漸變小，ASF 網(wǎng)格計(jì)算的誤差逐漸變小。

4.3 Junkins 插值算法仿真

4.3.1 不同的網(wǎng)格大小

junkins 算法的結(jié)果誤差分析如圖7 所示，從上到下依次是網(wǎng)格大小為0.5°，0.2°和0.1°，誤差范圍分別為（ -20～20）ns，（ -10～10）ns 和（ -5～5）ns，網(wǎng)格越小，ASF 網(wǎng)格能實(shí)現(xiàn)的精度越高。誤差分布的規(guī)律性介于雙線性插值差算法和反距離插值算法之間，網(wǎng)格內(nèi)左下角部分網(wǎng)格誤差小，右上角區(qū)域網(wǎng)格誤差大，誤差的分布較集中。

圖6 反距離加權(quán)插值算法（不同傳輸路徑）Fig.6 IDW interpolation algorithm （different path）

圖7 Junkins 算法（不同網(wǎng)格大小）Fig.7 The junkins weighting method（different mesh size）

4.3.2 不同的傳輸路徑

按照網(wǎng)格0.1°時(shí)計(jì)算的海水路徑、平均陸地和甚干燥地路徑的誤差情況如圖8 所示，海水介質(zhì)的誤差范圍在（ -2 ～2）ns，平均陸地的誤差范圍在（ -5～5）ns和甚干燥地的誤差范圍在（ -5～5）ns。隨著路徑的電導(dǎo)率變小，路徑的傳輸介質(zhì)對(duì)信號(hào)的阻力增強(qiáng)，同時(shí)ASF 網(wǎng)格的誤差也變大。因此，海水路徑對(duì)信號(hào)的阻力最小，相應(yīng)的ASF 網(wǎng)格計(jì)算誤差最小，隨著信號(hào)傳輸路徑上介質(zhì)對(duì)信號(hào)的阻力逐漸變大，ASF 網(wǎng)格的計(jì)算誤差逐漸變大。

圖8 Junkins 算法（不同傳輸路徑）Fig.8 The Junkins weighting method（different path）

5 算法比較

首先在待測試區(qū)域建立4 個(gè)網(wǎng)格，每個(gè)網(wǎng)格經(jīng)緯度方向的大小0.2°，待測試區(qū)域到發(fā)播臺(tái)的信號(hào)傳輸路徑屬于平均陸地。比較三種網(wǎng)格應(yīng)用算法如圖9 所示，從上到小依次是IDW 插值算法、Junkins 插值算法和Bilinear 插值算法，從最左邊一列看出三種算法誤差分布IDW 插值算法誤差分布的規(guī)律性最差，誤差最小的區(qū)域集中在網(wǎng)格中心部位，誤差大的區(qū)域分散在網(wǎng)格的邊緣區(qū)域；Junkins 算法的規(guī)律性介于二者中間，誤差最小的區(qū)域集中在網(wǎng)格的左下部分，誤差大的區(qū)域集中在網(wǎng)格的右上部分；雙線性算法誤差分布的規(guī)律性最好:離網(wǎng)格頂點(diǎn)越近，誤差越小，離網(wǎng)格頂點(diǎn)越遠(yuǎn)，誤差越大，從網(wǎng)格中心向外部擴(kuò)散開來誤差逐漸變小。

圖9 中間一列是色柱，反距離插值算法的誤差范圍是（0～15）ns；Junkins 算法的誤差范圍是（ -10 ～10）ns；雙線性插值算法的誤差范圍是（0～1）ns，雙線性插值算法引入的誤差最大。

圖9 中，右側(cè)3 個(gè)圖分別是三種算法的誤差概率分布函數(shù)曲線圖。IDW 插值算法誤差不大于20ns，95%的誤差小于12ns，Junkins 插值算法誤差不大于15ns；95%的誤差小于10.7ns，Bilinear 插值算法誤差不大于1.2ns，95%的誤差小于0.9943ns。三種算法分析比較得出IDW 插值算法誤差最大，Junkins 算法誤差次之，雙線性插值算法的誤差最小。

圖9 三種插值算法的比較示意圖Fig.9 The comparison of three interpolation algorithms

6 結(jié)束語

通過以上仿真分析比較，雙線性插值算法的誤差最小，而且呈現(xiàn)一定的規(guī)律型從網(wǎng)格中心到網(wǎng)格邊緣誤差逐漸減??；反距離插值算法優(yōu)于要計(jì)算距離算法計(jì)算量最大，誤差最大，而且誤差的規(guī)律性差；Junins 算法計(jì)算量、誤差大小以及誤差的規(guī)律性在三種算法中居中。