多數據源電離層層析成像方法

趙海生 許正文

（1.中國電波傳播研究所，電波環境特性及模化技術重點實驗室，山東 青島266107；82.中國科學院空間天氣學國家重點實驗室，北京100190）

1.引 言

電離層層析成像（簡稱CIT）已成為一種重要的電離層探測手段，地面接收站通過接收衛星信標信號，根據信號的差分相位、時延、法拉第旋轉等，提取從衛星信標到地面接收站之間的電子密度總含量（TEC）。根據地面站測得的不同星地連線的斜向TEC，進行反演區域的電離層CIT成像。電離層CIT成像技術從提出至今已有20多年的歷史，眾多學者提出各種各樣的CIT算法[1－4]，但直到現在仍沒有徹底解決電離層CIT技術的有限視角、稀疏布站等不適定性問題。因此，CIT成像結果的精度仍然不高，從而影響了電離層CIT技術的廣泛應用。

大地震前存在電離層異常擾動的觀點已逐漸被認可，從1964年發現地震與電離層存在某些聯系，至今已有許多設備觀測到了震前電離層異常擾動[5－11]。但是地震與電離層擾動存在不對應性，很多時候電離層出現異常擾動并沒有地震發生。地震和電離層擾動的不對應性，是試圖通過觀測電離層異常變化來預報地震的學者首先要解決的難題。為了減輕地震對我國人民生命財產造成的重大損失，也為了解決這一難題，我國啟動了電磁衛星計劃，其目標是力爭到2020年，初步形成高低軌道優化配置、穩定在軌運行的地震電磁衛星星座，從而使空間技術與地面觀測臺網結合起來，實施地震立體觀測，獲取更加全面系統的觀測信息，推動防震減災能力的提高。

要從眾多的電離層擾動中，精確提取地震信息，就要深入研究地震引起電離層擾動的特點。比如，地震擾動的強度、高度、形狀、持續時間等可能與其他原因引起的電離層擾動存在差異，因此，只有深入了解了地震擾動的特征才能將它精確識別。我國決定在地震電磁試驗衛星上攜帶三頻信標發射機，試圖通過電離層CIT成像技術，進行地震擾動特征的提取。然而目前CIT技術的現狀，并不具備高精度CIT成像的能力，在地震擾動特征提取上存在很大難度。多年的實踐經驗表明：單靠CIT算法的改進，無法徹底解決有限視角、稀疏布站等問題，無法從根本上提高CIT成像精度。許多學者將目光轉向多數據源融合方法，試圖從另一途徑解決CIT技術面臨的問題，其中有學者提出全球定位系統（GPS）數據和地基垂測儀數據聯合CIT成像的方法[12]，還有人提出三頻信標數據和我國華北地震監測網垂測、斜測數據聯合反演的方法[13]，這些方法一定程度上提高了反演精度，但沒有徹底擺脫模型對CIT技術的制約，初值往往依靠模型生成，初值若與實際吻合得好反演結果精度高，反之精度下降。本文提出的融合地面垂測、斜測數據、頂部探測儀數據和三頻衛星信標數據的多數據源電離層聯合CIT方法，從根本上擺脫了模型的限制，初值由地面垂測、斜測數據和頂部探測儀數據生成，因此，該方法能夠獲得穩定的高精度反演結果。模擬結果證明了該方法的可行性。

2.多種數據資源

2.1 地基高分辨率電離層監測網

研究發現，地震前兆的電離層異常特征尺度在百千米量級，具體尺度的大小因震級大小等因素而異。而我國現有的電波環境觀測網站之間相距至少600千米以上，難以有效開展電離層中小尺度不均勻性結構的觀測與研究。因此，通過建設電離層地基高分辨率觀測網，實現百千米尺度的電離層監測，對于開展電離層不均勻性研究，用于地震前兆研究，有著明顯的應用背景。

我國的首個電離層地基高分辨率觀測網于2009年4月建成并投入觀測，采用電離層斜向探測技術，依托中國電波傳播研究所已建的北京、長春、青島、新鄉、蘇州5個垂直探測站，在中國地震局的首都圈地震監測臺站布設20臺電離層斜向探測儀，每個斜測站同時接收5個垂測站的信號，鏈路反射點總共可形成100個電離層監測點，鏈路大圓距離均在2000 km以內，其中600 km以內的斜測鏈路達50條。經計算，按照100 km的空間分辨率，各鏈路的反射點基本可以覆蓋整個華北地區。垂測站、斜測站分布圖如圖1所示。

圖1 電離層地基高分辨率觀測網分布圖

2.2 電子密度剖面反演方法

電子密度剖面就是電子密度隨高度的變化，它是電離層探測所獲得的最有意義的結果之一。由某高度的反射頻率可以精確計算該高度的電子密度

式中：N為電子密度；m為電子質量；ε0為自由空間介電常數；e為電子電荷量；fp為等離子體頻率。

電離層圖上每個描跡對應的高度并不是該點的實際高度，該高度稱為虛高，虛高和實高的對應關系式為

式中：h′（f）為虛高；μ′為群折射指數；h0為電離層底高；fp為等離子體頻率；h為實高。推算電子密度剖面關鍵是計算每個反射頻率對應的反射高度。反射高度的精確計算至今還是一個難題[14－17]，目前廣泛應用的解算電子密度剖面的方法是 “分片法”，算法的具體推導過程參考文獻[18] 。

2.3 三頻信標TEC測量方法

三頻信標與雙頻信標相比多了一個載波ω3＝ω0，將三個頻率載波兩兩差分可以得到三組雙頻載波，三組載波利用差分多普勒技術可以得到三個差分相位表達式

式中：ΔP12、ΔP13、ΔP23為差分相位的絕對值，包含相對相位和相位積分常數。經推導、整理可得（推導過程詳見文獻 [19] ）

式中：Δφ12、Δφ13分別為f1與f2，f1與f3的差分相位的小數部分。在三頻衛星信標中n1＝3，n2＝8，f0＝16.668 MHz，將其帶入式（4）得

式中：8.3165×1016（單位為el/m2）稱為三頻信標測量TEC的相位模糊系數；k2為一正整數，它與雙頻信標的相位積分常數一樣，是需要確定的，把k2稱為三頻相位積分常數。文獻[19] 給出了k2的計算方法。

3.數據融合方法

從我國首都圈地震監測網選取8個站點（如表1），其中包括4個垂測站和4個斜測站，各個站反射點分布如表1所示。攜帶頂部探測儀的地震電磁衛星沿120°E經度飛行，在每個反射點及27°N和48°N上空完成一次掃頻，分別得到一張頂部電離圖，反演得到頂部電子密度剖面。同一時刻的頂部探測數據和底部探測數據聯合，經過數據融合處理可得到整個反演區域二維電子密度剖面，我們以此作為迭代初值進行了CIT成像研究。具體數據融合算法如下：

表1 選取的垂測站點及反射點位置

1）由頂部探測儀得到8個反射點和27°N、48°N上空頂部電子密度剖面。

2）由垂測儀得到底部電子密度剖面，由斜測儀得到反射點底部2～3個關鍵點的電子密度數據。

3）由4個垂測儀得到的底部電子密度剖面與同一時刻同一緯度位置的頂部電子密度剖面結合，得到4個完整電子密度剖面。

4）4個斜測反射點，底部數據由Chapman模型擬合得到，擬合過程中通過調整參數，使模型以最小二乘原理與斜測數據吻合；頂部數據由頂部探測儀數據得到。

5）在27°N、48°N處底部數據完全由Chapman模型得到，頂部數據由頂部探測儀數據得到。

6）對以上方法得到的10個電子密度剖面分別進行濾波平滑處理，然后用二維插值函數擬合出0.5°分辨率的完整的二維電子密度剖面。

由于站點分布稀疏、數據擬合誤差等原因，經過數據融合處理的二維電子密度剖面存在不可避免的誤差。但是該電子密度剖面以實測數據為基礎，完全克服了模型的限制，以此作為迭代初值與以模型作為初值相比，更具有穩定性和可靠性，從而提高了CIT成像結果的精度。

4.數值模擬

為了驗證算法的有效性，我們采用乘法代數算法（MART）進行電離層CIT反演，根據我國首都圈地震監測網的布站情況，采用沿120°E經線，30°N～45°N區域作為反演區域，選取8個地面站點（表1）。為了檢驗算法對小尺度擾動的監測能力，我們在34°N～38°N和40°N～42°N附近分別加入一個水平尺度約400 km和200 km的電離層擾動。再結合衛星信標的具體情況設置試驗參數如下：

時間：2009年11月17日世界時6時；

反演區域高度范圍：100～500 km；

衛星緯度跨越：27°N～48°N；

反演區域緯度跨越：30°N～45°N；

臺站設置：沿120°E，在33°N～43°N之間，每2個緯度一個臺站；

經度：120°E；

網格劃分：每個網格垂直10 km，緯度跨越0.5°；

地面站選取：如表1；

TEC采樣間隔：約0.064°（每隔0.064°一次TEC采樣）；

迭代初值：20091017UT0600和20090517UT0600時刻IRI2007模型算得的電子密度分布。

由于我國地震電磁試驗衛星尚未發射，目前無法得到實測的三頻信標TEC數據和頂部探測數據，檢驗該算法的數據只能通過模型得到。為了使模擬結果更符合實際，我們在數值模擬中加入了適當的隨機誤差。反演結果和誤差比較如圖2和圖3所示。

從圖2、圖3可以看出采用本文算法，基于多數據融合，利用實測數據生成迭代初值的反演結果明顯優于由模型生成的迭代初值的反演結果，20091017UT0600時刻的IRI模型作為初值的反演誤差為：1.2367×1011el/m3，20090517UT0600時刻的IRI模型作為初值的反演誤差為：1.4689×1011el/m3，而由實測數據生成的迭代初值反演誤差僅1.3809×1010el/m3.另外，在水平尺度為200 km的小尺度擾動的監測方面，本文算法對小尺度擾動形狀、尺度的反演更接近模型，明顯優于另外兩個反演結果。

圖4、圖5進一步給出了反演結果的峰值高度比較和在32°N、37°N和42°N三個緯度點電子密度剖面比較（圖5橫坐標的單位為1×1011el/m3）。從兩圖能夠清晰看出本文算法在提高反演結果垂直分辨率方面的效果，這一方法有效彌補了CIT成像有限視角的不足。

5.結 論

垂直分辨率不高是電離層CIT技術面臨的一個重要問題，結合三頻信標技術、地基高分辨率地震監測網及頂部探測技術，提出了一種實用有效的聯合CIT算法，并利用IRI模型對算法進行了檢驗，結果表明：采用本文提出的算法，能夠提高CIT反演精度，有效彌補有限視角不足。本算法將頂部探測技術應用于電離層CIT成像研究，進一步增加了數據源的類型，提高了反演精度，是本算法創新之處，也為我國地震電磁衛星電離層擾動探測提供了一種新的算法。

[1] RAYMUND T D，AUSTEN J R，FRANKE S J，et al.Application of computerized tomography to the investigation of ionospheric structures[J] .Radio Sci.，1990，25（5）：771－789.

[2] FREMOUW F J，SECAN J A，HOWE B M.Application of stochastic inverse theory to ionospheric tomography[J] .Radio Sci.，1992，27（5）：721－732.

[3] KUNITAKE M，OHTAKA K，MARUYAMA T，et al.Tomographic imaging of the ionosphere over Japan by the modified truncated SVD method[J] .Ann.Geophys.，1995，13（12）：1303－1310.

[4] FRIDMAN S V，L NICKISCH J，HAUSMAN M.Personal－computer－based system for real－time reconstruction of the three－dimensional ionosphere using data from diverse sources[J] .Radio Sci.，2009，44（3）：RS3008.doi：10.1029/2008RS004040.

[5] PULINETS S A.Strong earthquakes prediction possibility with the help of topside sounding from satellites[J] .Adv.Space Res.，1988，21（3），455－458.

[6] LIU J Y，CHEN Y I，PULINETS S A，et al.Seismoionospheric signatures prior to M 6 Taiwan earthquakes[J] .Geophys.Res.Lett.，2000，27（19）：3113－3116.

[7] SILINA A S，LIPEROVSKAYA E V，LIPEROVSKY V A，et al.Ionospheric phenomena before strong earthquakes[J] .Nat.Hazards Earth Syst.Sci.，2001，1（3）：113－118.

[8] PULINETS S A，BOYARCHUK K.Ionospheric Precursors of Earthquakes [M] . Berlin Heidelbert：Springer－Verlag，2004.

[9] PARROT M，BERTHELIER J J，LEBRETON J P，et al.Examples of unusual ionospheric observations made by the DEMETER satellite over seismic regions[J] .Phys.Chem.Earth，2006，31（4－9）：486－495.

[10] ZHAO B，WANG M，YU T，et al.Is an unusual large enhancement of ionospheric electron density linked with the 2008 great Wenchuan earthquake？[J] .J.Geophys.Res.，2008，113：A11304.doi：10.1029/2008JA013613.

[11] LIU J Y，CHEN Y I，CHEN C H，et al.Seismoionospheric GPS total electron content anomalies observed before the 12 May 2008 Mw7.9 Wenchuan earthquake[J] .J.Geophys.Res.，2009，114：A04320.doi：10.1029/2008JA013698.

[12] MA X F，MARUYAMA T.Three－dimensional ionospheric tomography using observation data of GPS ground receivers and ionosonde by neural network[J] .J.Geophys.Res.，2005，110：A05308.doi：10.1029/2004JA010797.

[13] ZHAO H S，XU Z W，WU J，et al.Ionospheric tomography by combining vertical and oblique ionograms with TEC retrieved from a tri－band beacon[J] .J.Geophys.Res.，2010，115：A10303.doi：10.1029/2010JA015285.

[14] DAVIES K.Ionospheric Radio[M] .London：Peter Peregrinus，1990.

[15] OBROU O K，BILITZA D，ADENIYI J O，et al.E－quatorial F2 layer peak height and correlation with vertical ion drift and M（3000）F2[J] .Adv.Space Res.，2003，31（3）：513－520.

[16] RAWER K and EYFRIG R.Improving the M（3000）－hmF2 relation[J] .Adv.Space Res.，2004，33（6）：878－879.

[17] XU Z W，WU J，IGARASHI K，et al.Long－term ionospheric trends based on ground－based ionosonde observations at Kokubunji，Japan[J] .J.Geophys.Res.， 2004， 109： A09307. doi：10.1029/2004JA010572.

[18] 熊年祿，唐存琛，李行健.電離層物理概論[M] .武漢：武漢大學出版社，1999.

[19] ZHAO H S，XU Z W，J WU，et al.Ionospheric tomography of small－scale disturbances with a triband beacon：a numerical study[J] .Radio Sci.，2010，45：RS3008.doi：10.1029/2009RS004247.