水汽層析代數(shù)重構(gòu)算法

何 林，柳林濤，蘇曉慶，3，許超鈐，段鵬碩，3

1.貴州電力設(shè)計(jì)研究院，貴州貴陽(yáng)550000；2.中國(guó)科學(xué)院測(cè)量與地球物理研究所，湖北武漢430077；3.中國(guó)科學(xué)院大學(xué)，北京100049；4.武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢430079

1 引 言

自從文獻(xiàn)［1］提出地基GPS水汽探測(cè)技術(shù)后，其高精度、全天候、高時(shí)空分辨率等優(yōu)點(diǎn)在氣象學(xué)內(nèi)得到了廣泛應(yīng)用。地基GPS水汽探測(cè)所得的大氣可降水量（precipitable water vapor，PWV）是測(cè)站天頂方向總的水汽含量，無(wú)法反映出水汽的三維分布，因此限制了它在天氣分析和數(shù)值預(yù)報(bào)中的應(yīng)用［2］。文獻(xiàn)［3—7］研究表明，區(qū)域GPS網(wǎng)各測(cè)站斜路徑濕延遲含有一定的水汽三維分布信息，因此研究如何從斜路徑水汽含量反演水汽三維分布具有較大的科研價(jià)值和實(shí)用意義。

文獻(xiàn)［3］率先開(kāi)展水汽層析方面的研究，采用SVD方法求解層析方程，并用ECMWF預(yù)報(bào)值驗(yàn)證結(jié)果；文獻(xiàn)［4］分析了約束條件對(duì)GPS水汽層析結(jié)果的影響，并證明初值對(duì)層析結(jié)果的可靠性有重要影響；文獻(xiàn)［5］初步討論了 ART、SART、MART等在GNSS水汽層析中應(yīng)用，并給出了經(jīng)驗(yàn)松弛因子和迭代終止條件；文獻(xiàn)［6］驗(yàn)證了ART、MART、SIRT用于對(duì)流層水汽層析的可行性，也給出了參考松弛因子和迭代終止參數(shù)。

由于測(cè)站幾何分布、衛(wèi)星星座、網(wǎng)格劃分等因素的影響，層析方程系數(shù)矩陣往往稀疏且嚴(yán)重秩虧，導(dǎo)致其無(wú)法直接求解［5－8］。代數(shù)重構(gòu)算法是一種線性反演方法，采用迭代的形式避免了法方程的求逆。本文研究了代數(shù)重構(gòu)算法在水汽層析應(yīng)用中的各種問(wèn)題，如約束條件、松弛因子、終止條件等。本文首先介紹水汽層析原理，約束條件及層析方程構(gòu)建，然后分析代數(shù)重構(gòu)算法在層析方程求解中的各種問(wèn)題，包括初值的選擇、松弛因子的計(jì)算、迭代終止的確定等，并以香港SatRef網(wǎng)進(jìn)行驗(yàn)證分析。

2 層析原理

GPS信號(hào)在傳播過(guò)程中由于對(duì)流層的影響會(huì)造成一定的信號(hào)延遲（slant total delay，STD）［9－15］，該延遲可表示為

式中，N為大氣折射率，是靜力學(xué)折射率Nd和濕折射率Nw之和，其中Nd可由地面氣象觀測(cè)值計(jì)算。由此可得由水汽引起的濕延遲（slant wet delay，SWD）

通過(guò)高精度GPS數(shù)據(jù)處理軟件解算出STD、SWD等，然后由SWD計(jì)算斜路徑大氣可降水量（slant wet vapor，簡(jiǎn)稱SWV）。得到SWV后，根據(jù)實(shí)際情況對(duì)觀測(cè)方程進(jìn)行離散化，即將研究區(qū)域空間劃分網(wǎng)格。離散化后，觀測(cè)方程變?yōu)?/p>

式中，A為系數(shù)陣；x表示網(wǎng)格水汽密度；d為斜路徑水汽含量。

3 層析方程

3.1 約束條件

前面已經(jīng)提到層析方程是秩虧的，即大量網(wǎng)格沒(méi)有GPS信號(hào)穿過(guò)。對(duì)于那些沒(méi)有信號(hào)穿過(guò)的網(wǎng)格，方程解算的結(jié)果與初值一樣。為了解決該問(wèn)題，需引入平滑約束條件。

3.1.1 水平平滑約束

水汽在空間分布上具有連續(xù)性，根據(jù)地理學(xué)第一定律，距離越近的網(wǎng)格，相關(guān)性越大，故可定義平滑約束方程［16］

式中，H為水平平滑約束矩陣；l為約束條件數(shù)?？紤]到水汽緩慢均勻的變化特性，可采用均值濾波器進(jìn)行平滑處理［7］。對(duì)于那些處于研究區(qū)域中間的網(wǎng)格來(lái)說(shuō)，約束條件由該網(wǎng)格所在層面周?chē)?個(gè)網(wǎng)格的二階Laplacian算子給出，即

對(duì)于處于探測(cè)平面邊緣的網(wǎng)格來(lái)說(shuō)，該算子必須要經(jīng)過(guò)適當(dāng)調(diào)整才能使用。這里給出每個(gè)平面頂點(diǎn)處的算子，分別為左上角H1，右上角H2，左下角H3，右下角H4

對(duì)于那些處于邊線而非頂點(diǎn)處的網(wǎng)格來(lái)說(shuō)，該算子調(diào)整為上邊線H5，下邊線H6，左邊線H7，右邊線H8

在附加該水平平滑約束后，有多個(gè)信號(hào)穿過(guò)的網(wǎng)格將給臨近網(wǎng)格提供改正值，沒(méi)有信號(hào)穿過(guò)的網(wǎng)格將從其臨近網(wǎng)格獲取改正信息，從而使整體結(jié)果得到修正。

3.1.2 垂直約束

探空資料提供的剖面信息顯示水汽主要集中在6km以下，10km處水汽含量接近于零，因此將上邊界的水汽密度強(qiáng)制約束為0.01g／m3。根據(jù)上述說(shuō)明，構(gòu)建垂直約束方程

3.2 層析方程構(gòu)建

根據(jù)以上約束條件和觀測(cè)方程，構(gòu)造出層析方程

式中，As為斜方向投影函數(shù)；SWV為斜路徑水汽含量；Ad為天頂方向投影函數(shù)；PWV為天頂方向水汽含量；H表示水平約束條件；V表示垂直約束條件。

4 代數(shù)重構(gòu)算法

這里簡(jiǎn)單介紹代數(shù)重構(gòu)算法中常用的8種算法，它 們 分 別 為 Kaczmarz、Randkaczmarz、Symkaczmarz、SART、Landweber、Cimmino、CAV、DROP［5，17－18］。 其 中 Kaczmarz、Randkaczmarz、Symkaczmarz算法都屬于 ART（algebraic reconstruction techniques）算法，區(qū)別僅在于內(nèi)循環(huán)順序不一樣，Kaczmarz為依次循環(huán)，Randkaczmarz為隨機(jī)循環(huán)，而Symkaczmarz包括一次順序 循 環(huán) 和 一 次 逆 序 循 環(huán)［10］。 而 SART、Landweber、Cimmino、CAV 和DROP算法都屬于SIRT（simultaneous iterative reconstruction techniques）算法，與ART算法的不同在于它不再單獨(dú)對(duì)每條觀測(cè)路徑上的網(wǎng)格逐個(gè)進(jìn)行修正，故沒(méi)有內(nèi)循環(huán)，而是一次性對(duì)所有觀察路徑進(jìn)行迭代。各算法的具體計(jì)算公式請(qǐng)查閱相關(guān)參考文獻(xiàn)，限于篇幅，這里不做詳細(xì)介紹。

5 層析方程求解

5.1 迭代初值選取

文獻(xiàn)［4］的研究表明，水汽層析結(jié)果的可靠性在很大程度上依賴于先驗(yàn)信息的精度，因此初值的選擇也是代數(shù)重構(gòu)算法的重要部分。對(duì)于水汽的層析解算，初值的獲取一般有以下3種方式：①利用無(wú)線電探空儀提供的大氣廓線信息，該方法精度最高；②利用數(shù)值預(yù)報(bào)模型的計(jì)算結(jié)果，該方法也較準(zhǔn)確；③直接利用標(biāo)準(zhǔn)大氣分布模型，該方法精度較差。

5.2 松弛因子選擇

松弛因子在迭代中起到調(diào)節(jié)改正量的作用，影響著收斂速度和迭代結(jié)果的質(zhì)量。文獻(xiàn)［5］證明，松弛因子的選擇甚至比算法還重要。這里介紹 一 種 計(jì) 算 最 優(yōu) 松 弛 因 子 的 方 法［5，17－18］，松 弛 因子能以最快的速度收斂到最小方差，計(jì)算松弛因子需要一定的先驗(yàn)知識(shí)。

該算法［18］由兩步組成，首先確定能獲取整體最小方差所有可行的λ，然后確定能通過(guò)最少的迭代次數(shù)就能達(dá)到該最小方差的λ。如果用ηλ表示對(duì)應(yīng)給定λ的最小方差，則由于迭代的離散本質(zhì)，ηλ會(huì)隨著λ的變化而緩慢變化，但是對(duì)于不同的λ，ηλ的偏差很小。以SIRT算法為例，試驗(yàn)發(fā)現(xiàn)，若要達(dá)到整體方差最小，λSIRT＝1／ρ（ATMA）比較合適，則可定義最小方差的上限＝1.01ηλSIRT。

假設(shè)kλ是給定的松弛因子λ為達(dá)到所需要的迭代次數(shù)，接下來(lái)就是在最小化kλ中計(jì)算出最優(yōu)的λ。試驗(yàn)表明kλ是λ的單峰函數(shù)，故可用以下經(jīng)過(guò)改進(jìn)的黃金分割搜索法來(lái)計(jì)算最優(yōu)的λ。該方法首先給定λ一個(gè)初始搜索域（α，β），即SIRT的收斂區(qū)間，然后計(jì)算出兩個(gè)初始內(nèi)部點(diǎn)α′＝α＋r（β－α）和β′＝α＋（1－r）（β－α），其中r＝（3 －）／2。假定此時(shí)的迭代次數(shù)為κα′和κβ′，相關(guān)的最小方差為ηα′和ηβ′，然后根據(jù)以下判斷來(lái)縮小搜索區(qū)間。

ηα′＞：最小方差λ＝α′沒(méi)有達(dá)到整體最小方差水平，將區(qū)間縮減到（α′，β）。

ηβ′＞：最小方差λ＝ηβ′沒(méi)有達(dá)到整體方差最小值，將區(qū)間縮減到（α，β′）。

κα′≥κβ′：α′和β′對(duì)于λ都是可行的值，根據(jù)單峰性，將區(qū)間縮減到（α′，β）。

κα′＜κβ′：α′和β′對(duì)于λ都是可行的值，根據(jù)單峰性，將區(qū)間縮減到（α，β′）。

按照上述步驟循環(huán)搜索松弛因子，直到其所在的區(qū)間寬度足夠?。ㄈ?0－6），然后以區(qū)間中點(diǎn)為最優(yōu)的λ。該訓(xùn)練算法對(duì)于ART算法也類似，但是對(duì)于Randkaczmarz方法，由于內(nèi)循環(huán)指數(shù)i的前后不一致性，故需要通過(guò)試驗(yàn)來(lái)獲取不同的λ。此外，在實(shí)現(xiàn)該訓(xùn)練算法中需設(shè)定最大循環(huán)數(shù)kmax，這是由于對(duì)于某些問(wèn)題，可能沒(méi)有在指定的循環(huán)數(shù)內(nèi)達(dá)到最小方差，此時(shí)就需要在k＝kmax下確定能夠達(dá)到最小方差的最優(yōu)松弛因子。

5.3 迭代終止條件

代數(shù)重構(gòu)算法需要通過(guò)迭代來(lái)達(dá)到試驗(yàn)數(shù)據(jù)的響應(yīng)，考慮到有可能迭代到某一個(gè)點(diǎn)后結(jié)果變差，或出現(xiàn)參數(shù)震蕩，因此必須確定達(dá)到試驗(yàn)數(shù)據(jù)的最佳響應(yīng)時(shí)刻，從而最優(yōu)的結(jié)束循環(huán)［5，18－20］。一種方式是當(dāng)結(jié)果趨于穩(wěn)定，即當(dāng)xk＋1和xk兩者之差小于某個(gè)閾值時(shí)結(jié)束；另一種方式是當(dāng)｜Axk－mk｜＝min時(shí)結(jié)束循環(huán)［5］。但這些標(biāo)準(zhǔn)并不適用于該問(wèn)題，為此，文獻(xiàn)［5—6］定義了一系列參數(shù)作為判斷迭代終止的依據(jù)。

本文采用正規(guī)化累積周期圖（normalized cumulative periodogram，NCP）作為迭代終止條件［18－20］，接下來(lái)詳細(xì)介紹該方法。在 NCP中，將每次迭代的剩余量rk＝b－Axk看作時(shí)間序列，由于其表示的信號(hào)與白噪聲有著明顯不同，因此可以判斷出最優(yōu)迭代數(shù)k。定義＾r(nóng)k∈Cm為rk的離散傅里葉變換，q為滿足q≤m／2的最大整數(shù)，則rk的NCP元素（Ck∈Rq）為

如果rk由白噪聲構(gòu)成，則其期望功率譜為平的，期望 NCP元素位于從（0，0）到（q，1）的直線上。實(shí)際噪聲的功率譜沒(méi)有理想中的平坦，但NCP元素分布仍然接近于直線。該問(wèn)題的關(guān)鍵在于確定剩余量從信號(hào)主導(dǎo)到白噪聲主導(dǎo)的迭代數(shù)目k［19－20］，一旦rk呈現(xiàn)出白噪聲的分布特點(diǎn)，即NCP分布接近于直線，迭代就可以終止。決定rk是否為白噪聲的計(jì)算是實(shí)際NCP元素的Kolmogorov－Smirnoff檢驗(yàn)，若實(shí)際的NCP元素位于理論NCP元素的Kolmogorov－Smirnoff區(qū)間內(nèi)，則表示可接受rk為白噪聲這個(gè)假設(shè)，表明此時(shí)迭代可終止，從而可以確定出最佳迭代數(shù)kopt。若假設(shè)檢驗(yàn)的顯著性水平為5%，則Kolmogorov－Smirnoff區(qū)間的界線為±1.36q－1／2。

6 試驗(yàn)分析

試驗(yàn)采用香港衛(wèi)星定位參考站網(wǎng)（SatRef）提供的2012年148～155d的數(shù)據(jù)，用GAMIT進(jìn)行解算，引入BJFS、SHAO和KUNM 3個(gè)IGS站，衛(wèi)星高度截止角取為10°，天頂靜力學(xué)延遲采用Saastamoinen模型，對(duì)流層延遲采用線性估計(jì)模型，映射函數(shù)采用VMF1模型。層析試驗(yàn)中，將研究區(qū)域劃分為6×6×10＝360個(gè)網(wǎng)格，如圖1所示，采用京士柏氣象站（HKKP）的無(wú)線電探空儀資料作為檢驗(yàn)層析結(jié)果質(zhì)量的標(biāo)準(zhǔn)。按第3節(jié)的方法構(gòu)建層析方程，層析初值選擇為層析時(shí)刻前3天的無(wú)線電探空儀數(shù)據(jù)的均值，接下來(lái)分別分析松弛因子、迭代終止規(guī)則和不同的代數(shù)重構(gòu)算法對(duì)層析結(jié)果的影響。

圖1 SatRef網(wǎng)格劃分圖Fig.1 SatRef schematic meshing

6.1 松弛因子的試驗(yàn)分析

采用文中5.2節(jié)的方法計(jì)算出各算法的最優(yōu)松弛因子（其中Randkaczmarz算法的最優(yōu)松弛因子由試驗(yàn)得到），并給出各算法松弛因子的有效取值范圍、默認(rèn)值、文獻(xiàn)［5—6］提供的經(jīng)驗(yàn)參考值，如表1所示。

表1 各算法的松弛因子Tab.1 The relaxation factors for algorithms

表1中，為A陣最大奇異值的估值為D1／2A最大奇異值的估值為D1／2SA最大奇異值的估值為S－1ATMA譜半徑的估值，各參數(shù)定義請(qǐng)參考文獻(xiàn)［18］。

首先分析不同松弛因子對(duì)解算結(jié)果的影響，分別為默認(rèn)值、文獻(xiàn)［5—6］分別給出的經(jīng)驗(yàn)參考值，以及由5.2節(jié)方法計(jì)算的最優(yōu)松弛因子。由于文獻(xiàn)［5—6］只給出了Kaczmarz和SART兩種算法的參考松弛因子，故此處僅對(duì)這兩種算法進(jìn)行分析（表2、表3）。

表2 不同松弛因子達(dá)到NCP迭代結(jié)果所需迭代次數(shù)對(duì)比表Tab.2 The comparison of required iterations number to achieve NCP result for different relaxation factor

表3 不同松弛因子的NCP迭代結(jié)果精度對(duì)比表Tab.3 The comparison of NCP results accuracy for different iterative relaxation factor g／m3

從表2、表3可以看出，在采用同樣的迭代初值和終止規(guī)則下，不同的松弛因子會(huì)導(dǎo)致收斂速度和結(jié)果精度不同。在Kaczmarz算法中，系統(tǒng)默認(rèn)松弛因子最大，收斂最快，但結(jié)果精度最低，文獻(xiàn)［6］的參考迭代因子最小，結(jié)果精度最高，但收斂速度最慢，而黃金分割算法給出的松弛因子精度與文獻(xiàn)［6］的參考迭代因子相當(dāng)，但收斂速度遠(yuǎn)大于它，綜合表現(xiàn)最好；而在SART算法中，4個(gè)松弛因子的解算結(jié)果精度相當(dāng)，但最優(yōu)松弛因子的收斂速度遠(yuǎn)高于另外3種算法，綜合表現(xiàn)也最好。

6.2 迭代終止規(guī)則的驗(yàn)證分析

接下來(lái)分析不同的迭代終止規(guī)則對(duì)水汽層析結(jié)果的影響。試驗(yàn)選取了最大迭代數(shù)（即迭代達(dá)到最大迭代數(shù)后終止，本試驗(yàn)中設(shè)為10 000），文獻(xiàn)［5］和文獻(xiàn)［6］分別給出的參考終止條件，NCP迭代終止條件4種不同迭代終止規(guī)則，采用最優(yōu)松弛因子進(jìn)行解算（表4）。

表4 不同迭代終止規(guī)則的水汽層析結(jié)果精度比較Tab.4 The comparison of accuracy of the results for different iterations termination rules g／m3

從表4結(jié)果可以看出，采用不同迭代終止規(guī)則，水汽層析結(jié)果有著較大的差異。在這兩種算法中，NCP迭代終止規(guī)則的層析結(jié)果精度最高；其次為文獻(xiàn)［5］提出的參考終止條件，文獻(xiàn)［6］與文獻(xiàn)［5］精度比較接近，精度最差的是最大迭代數(shù)。這說(shuō)明在代數(shù)重構(gòu)算法中，并非迭代次數(shù)越多結(jié)果越好，可能到某一點(diǎn)后迭代會(huì)發(fā)散，導(dǎo)致結(jié)果變差，因此選擇合理的迭代終止條件很重要，本文提出的NCP終止規(guī)則較為理想。

6.3 不同解算方案的結(jié)果分析

根據(jù)松弛因子和終止規(guī)則的不同設(shè)計(jì)了4種層析方程解算方案，分別為默認(rèn)松弛因子和最大迭代數(shù)（即方案1）；最優(yōu)松弛因子和最大迭代數(shù)目（即方案2）；默認(rèn)松弛因子和NCP終止規(guī)則（方案3）；最優(yōu)松弛因子和NCP終止規(guī)則（方案4），并用8種常見(jiàn)的代數(shù)重構(gòu)算法分別按這4種方案進(jìn)行解算，將結(jié)果整理成表5。

表5 不同解算方案的結(jié)果精度比較Tab.5 The comparison of accuracy of results for different solutions g／m3

由表5可以看出，對(duì)于同樣的層析方程，采用不同的方案進(jìn)行解算，得到的結(jié)果精度有較大差異。其中方案1解算結(jié)果精度最差，方案4精度最高，后兩種方案的結(jié)果精度都遠(yuǎn)高于前兩種解算方案。這說(shuō)明在代數(shù)重構(gòu)算法中，迭代終止規(guī)則的確定比松弛因子的選擇更為重要。在方案4中，精度最高的為CAV算法，其次為Cimmino算法。

為了便于更加直觀地分析解算結(jié)果，這里給出2012年155d0h京士柏氣象站處的水汽層析結(jié)果垂直廓線，如圖2所示。

圖2 采用方案4得到的水汽廓線Fig.2 Vapor profile got by strategy 4

從圖2可以看出，這8種代數(shù)重構(gòu)算法的結(jié)果都能較好地反映出水汽的分布趨勢(shì)，與無(wú)線電探空儀的結(jié)果接近，這表明代數(shù)重構(gòu)算法能夠較好地實(shí)現(xiàn)水汽層析。

7 結(jié) 論

本文討論了代數(shù)重構(gòu)算法在水汽層析應(yīng)用中的各種問(wèn)題，包括約束條件，層析方程，松弛因子，終止條件等，并提出了相應(yīng)解決方案，比較分析了8種常見(jiàn)的代數(shù)重構(gòu)算法。

松弛因子影響著迭代收斂速度和結(jié)果精度，在代數(shù)重構(gòu)算法中占有重要地位，本文給出了一種計(jì)算最優(yōu)松弛因子的黃金分割搜索算法，試驗(yàn)表明該算法能以較快的速度收斂到最優(yōu)結(jié)果。迭代終止條件在代數(shù)重構(gòu)算法中也很重要，本文給出了一種NCP迭代終止規(guī)則，試驗(yàn)結(jié)果表明該終止條件結(jié)果穩(wěn)定可靠，精度較高。試驗(yàn)還證明，迭代終止規(guī)則的確定比松弛因子的選擇更加重要。在各種層析方程解算方案中，采用文中計(jì)算最優(yōu)松弛因子的算法和NCP迭代終止規(guī)則結(jié)果精度最高。各種常見(jiàn)的代數(shù)重構(gòu)算法在采用該解算方案時(shí)，精度最高的為CAV算法，其次為Cimmino算法。

［1］ BEVIS M，BUSINGER S，HERRING T，et al.GPS Meteorology：Remote Sensing of Atmospheric Water Vapor Using the Global Positioning System［J］.Journal of Geophysical Research：Atmospheres，1992，97（D14）：15787－15801.

［2］ DING Jincai.GPS Meteorology and Its Applications［M］.Beijing：Meteorological Press，2009：148－161.（丁金才.GPS氣象學(xué)及其應(yīng)用［M］.北京：氣象出版社，2009：148－161.）

［3］ FLORES A，RUFFINIG G，RIUS A.4DTropospheric Tomography Using GPS Slant Wet Delays［J］.Annales Geophysicae，2000，18：223－224.

［4］ YU Shengjie，LIU Lintao，LIANG Xinghui.Influence Analysis of Conditions on GPS Water Vapor Tomography［J］.Acta Geodaetica et Cartographica Sinica，2010，39（5）：491－496.（于勝杰，柳林濤，梁星輝.約束條件對(duì)GPS水汽層析解算的影響分析［J］.測(cè)繪學(xué)報(bào)，2010，39（5）：491－496.）

［5］ BENDER M，DICK G，GE M R，et al.Development of a GNSS Water Vapour Tomography System Using Algebraic Reconstruction Techniques［J］.Advances in Space Research，2011，47（10）：1704－1720.

［6］ WANG Wei，WANG Jiexian.Ground－based GPS Water Vapor Tomography Based on Algebraic Reconstruction Technique［J］.Journal of Computer Applications，2011，31（11）：3149－3151.（王維，王解先.基于代數(shù)重構(gòu)技術(shù)的對(duì)流層水汽層析［J］.計(jì)算機(jī)應(yīng)用，2011，31（11）：3149－3151.）

［7］ ROHM W.The Ground GNSS Tomography：Unconstrained Approach［J］.Advances in Space Research，2013，51（9）：501－513.

［8］ WEN D B，LIU S Z，TANG P Y.Tomographic Reconstruction of Ionospheric Electron Density Based on Constrained Algebraic Reconstruction Technique［J］.GPS Solution，2012，14（4）：375－380.

［9］ TROLLER M，GEIGER A，BROCKMANN E，et al.Tomographic Determination of the Spatial Distribution of Water Vapor Using GPS Observations［J］.Advances in Space Research，2006，37：2211－2217.

［10］ ROHM W，BOSY J.The Verification of GNSS Tropospheric Tomography Model in a Mountainous Area［J］Advances in Space Research，2008，47：1721－1730.

［11］ NOTARPIETRO R，CUCCA M，GABELLA M，et al.Tomographic Reconstruction of Wet and Total Refractivity Fields from GNSS Receiver Networks［J］.Advances in Space Research，2011，47（1）：898－912.

［12］ GRADINARSKY L P，JARLEMARK P.Ground Based GPS Tomography of Water Vapor：Analysis of Simulated and Real Data［J］.Journal of the Meteorological Society of Japan，2004，82（18）：551－560.

［13］ HIRAHARA K.Local GPS Tropospheric Tomography［J］.Earth Planets and Space，2000，52：935－939.

［14］ TROLLER M，GEIGER A，BROCKMANN E，et al.Tomographic Determination of the Spatial Distribution of Water Vapour Using GPS Observations［J］.Advances in Space Research，2006，37：2211－2217.

［15］ DENG Z，BENDER M，DICK G，et al.Retrieving Tropospheric Delays from GPS Networks Densified with Single Frequency Receivers［J］.Geophysical Research Letters，2009，36（19）：802－806.

［16］ HOBIGER T，KONDO T，KOYAMA Y.Constrained Simultaneous Algebraic Reconstruction Technique（C－SART）：A New and Simple Algorithm Applied to Ionospheric Tomography［J］.Earth Planets Space，2008，60：727－735.

［17］ MING Y.Convergence Analysis of SART：Optimization and Statistics［J］.International Journal of Computer Mathematics，2013，90（1）：30－47.

［18］ HANSEN P C，SAXILD－HANSE M.AIR Tools：A MATLAB Package of Algebraic Iterative Reconstruction Methods［J］.Journal of Computational and Applied Mathematics，2012，236（8）：2167－2178.

［19］ HANSEN P C，KILMER M E，KJELDSEN R H.Exploiting Residual Information in the Parameter Choice for Discrete Ill－posed Problems［J］.BIT Numerical Mathematics，2006，46：41－59.

［20］ RUST B W，O’LEARY D P.Residual Periodograms for Choosing Regularization Parameters for Ill－posed Problems［J］.Inverse Problems，2008，24（3）：34005－34034.