趙海生 許正文

（1.中國電波傳播研究所，電波環(huán)境特性及模化技術(shù)重點(diǎn)實(shí)驗(yàn)室，山東 青島266107；82.中國科學(xué)院空間天氣學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京100190）

1.引 言

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

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

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

2.多種數(shù)據(jù)資源

2.1 地基高分辨率電離層監(jiān)測網(wǎng)

研究發(fā)現(xiàn)，地震前兆的電離層異常特征尺度在百千米量級，具體尺度的大小因震級大小等因素而異。而我國現(xiàn)有的電波環(huán)境觀測網(wǎng)站之間相距至少600千米以上，難以有效開展電離層中小尺度不均勻性結(jié)構(gòu)的觀測與研究。因此，通過建設(shè)電離層地基高分辨率觀測網(wǎng)，實(shí)現(xiàn)百千米尺度的電離層監(jiān)測，對于開展電離層不均勻性研究，用于地震前兆研究，有著明顯的應(yīng)用背景。

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

圖1 電離層地基高分辨率觀測網(wǎng)分布圖

2.2 電子密度剖面反演方法

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

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

電離層圖上每個(gè)描跡對應(yīng)的高度并不是該點(diǎn)的實(shí)際高度，該高度稱為虛高，虛高和實(shí)高的對應(yīng)關(guān)系式為

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

2.3 三頻信標(biāo)TEC測量方法

三頻信標(biāo)與雙頻信標(biāo)相比多了一個(gè)載波ω3＝ω0，將三個(gè)頻率載波兩兩差分可以得到三組雙頻載波，三組載波利用差分多普勒技術(shù)可以得到三個(gè)差分相位表達(dá)式

式中：ΔP12、ΔP13、ΔP23為差分相位的絕對值，包含相對相位和相位積分常數(shù)。經(jīng)推導(dǎo)、整理可得（推導(dǎo)過程詳見文獻(xiàn) [19] ）

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

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

3.數(shù)據(jù)融合方法

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

表1 選取的垂測站點(diǎn)及反射點(diǎn)位置

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

2）由垂測儀得到底部電子密度剖面，由斜測儀得到反射點(diǎn)底部2～3個(gè)關(guān)鍵點(diǎn)的電子密度數(shù)據(jù)。

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

4）4個(gè)斜測反射點(diǎn)，底部數(shù)據(jù)由Chapman模型擬合得到，擬合過程中通過調(diào)整參數(shù)，使模型以最小二乘原理與斜測數(shù)據(jù)吻合；頂部數(shù)據(jù)由頂部探測儀數(shù)據(jù)得到。

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

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

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

4.數(shù)值模擬

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

時(shí)間：2009年11月17日世界時(shí)6時(shí)；

反演區(qū)域高度范圍：100～500 km；

衛(wèi)星緯度跨越：27°N～48°N；

反演區(qū)域緯度跨越：30°N～45°N；

臺(tái)站設(shè)置：沿120°E，在33°N～43°N之間，每2個(gè)緯度一個(gè)臺(tái)站；

經(jīng)度：120°E；

網(wǎng)格劃分：每個(gè)網(wǎng)格垂直10 km，緯度跨越0.5°；

地面站選?。喝绫?；

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

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

由于我國地震電磁試驗(yàn)衛(wèi)星尚未發(fā)射，目前無法得到實(shí)測的三頻信標(biāo)TEC數(shù)據(jù)和頂部探測數(shù)據(jù)，檢驗(yàn)該算法的數(shù)據(jù)只能通過模型得到。為了使模擬結(jié)果更符合實(shí)際，我們在數(shù)值模擬中加入了適當(dāng)?shù)碾S機(jī)誤差。反演結(jié)果和誤差比較如圖2和圖3所示。

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

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

5.結(jié) 論

垂直分辨率不高是電離層CIT技術(shù)面臨的一個(gè)重要問題，結(jié)合三頻信標(biāo)技術(shù)、地基高分辨率地震監(jiān)測網(wǎng)及頂部探測技術(shù)，提出了一種實(shí)用有效的聯(lián)合CIT算法，并利用IRI模型對算法進(jìn)行了檢驗(yàn)，結(jié)果表明：采用本文提出的算法，能夠提高CIT反演精度，有效彌補(bǔ)有限視角不足。本算法將頂部探測技術(shù)應(yīng)用于電離層CIT成像研究，進(jìn)一步增加了數(shù)據(jù)源的類型，提高了反演精度，是本算法創(chuàng)新之處，也為我國地震電磁衛(wèi)星電離層擾動(dòng)探測提供了一種新的算法。

[1] RAYMUND T D，AUSTEN J R，F(xiàn)RANKE 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）－h(huán)mF2 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] .武漢：武漢大學(xué)出版社，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.