王 維 王解先

(同濟(jì)大學(xué)測(cè)量與國(guó)土信息工程系，上海 200092)

聯(lián)合迭代重構(gòu)算法在對(duì)流層水汽三維重構(gòu)中的應(yīng)用研究*

王 維 王解先

(同濟(jì)大學(xué)測(cè)量與國(guó)土信息工程系，上海 200092)

通過數(shù)值模擬仿真將SIRT算法應(yīng)用于水汽層析過程，結(jié)果顯示該算法可以較好地重構(gòu)水汽的空間分布，并且不受觀測(cè)方程性質(zhì)的影響，收斂較快，較之求逆的層析算法更易于實(shí)現(xiàn)。

聯(lián)合迭代重構(gòu)算法;圖像重構(gòu);水汽;層析;迭代

1 引言

全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS)的觀測(cè)信號(hào)穿過大氣時(shí)由于受到水汽折射影響發(fā)生信號(hào)偏折，從而產(chǎn)生延遲。目前已經(jīng)有許多研究通過利用信號(hào)延遲來(lái)獲取對(duì)流層水汽信息，并反演出水汽的空間分布。

目前，借助GNSS層析技術(shù)反演水汽密度的方法得到發(fā)展，其求解法方程有以下幾種方法：

1)在理想情況下，每個(gè)網(wǎng)格范圍內(nèi)存在至少一個(gè)測(cè)站，法方程系數(shù)陣滿秩可以直接求逆，但在目前GNSS監(jiān)測(cè)網(wǎng)條件下較難滿足。

2)在測(cè)站分布不均勻的時(shí)候，觀測(cè)值在網(wǎng)格中的分布不均勻，會(huì)出現(xiàn)覆蓋過剩與覆蓋不足，用奇異值分解法［1］(SVD)將法方程系數(shù)陣對(duì)角化，分離系數(shù)陣欠定部分，超定部分可以直接求逆［2］。

3)附加約束條件方程可以提供一部分信息同時(shí)使法方程系數(shù)陣穩(wěn)定，但是隨著不同類型觀測(cè)值的增加，使得最小二乘求逆和先驗(yàn)權(quán)的確定難度加大，采用Kalman濾波算法［3］對(duì)層析網(wǎng)格參數(shù)進(jìn)行求解可以避開約束條件的限制［4］，但是狀態(tài)向量初值及其方差協(xié)方差矩陣難以確定。

近年來(lái)，代數(shù)重構(gòu)算法不斷用于層析反演中，該方法穩(wěn)定并且迭代次數(shù)較少，可以進(jìn)行快速求解［5］。聯(lián)合迭代重構(gòu)算法是代數(shù)重構(gòu)算法中的一種，由于其計(jì)算的簡(jiǎn)便性，越來(lái)越多地用于圖像重構(gòu)，如反演電離層電子密度分布［6，7］。該算法在圖像重構(gòu)過程中不再像代數(shù)重構(gòu)算法逐個(gè)網(wǎng)格參數(shù)依次進(jìn)行改正，而是在某一輪迭代中，所有網(wǎng)格參數(shù)的函數(shù)值均用上一次迭代的改正數(shù)進(jìn)行整體修正，因而提高了計(jì)算效率和反演結(jié)果的精度。本文將聯(lián)合迭代重構(gòu)算法用于對(duì)流層水汽重構(gòu)以期縮短計(jì)算時(shí)間，提高近實(shí)時(shí)短期氣象預(yù)報(bào)與災(zāi)害監(jiān)測(cè)的效率。

2 水汽層析原理

假設(shè)在指定觀測(cè)時(shí)間段內(nèi)每個(gè)網(wǎng)格包含的水汽密度是均勻的，根據(jù)GPS信號(hào)在所經(jīng)過網(wǎng)格內(nèi)穿過的長(zhǎng)度，列出觀測(cè)方程［2］為

用矩陣的形式表示為

其中，A為系數(shù)陣，由所有觀測(cè)值穿過網(wǎng)格的長(zhǎng)度組成。X為參數(shù)陣，表示了網(wǎng)格的水汽密度。SWV為常數(shù)陣，每個(gè)分量表示一條斜路徑方向上的水汽含量。由于網(wǎng)格數(shù)量較大，并且斜路徑觀測(cè)值在分布上的不均勻以及時(shí)段內(nèi)所積累的觀測(cè)值數(shù)量等均不能保證每個(gè)網(wǎng)格都有觀測(cè)值穿過。為了解決某些網(wǎng)格沒有觀測(cè)值的問題，增加水平及垂直約束方程，并與觀測(cè)方程一起解算，模型為：

式中H、V為水平約束與垂直約束方程系數(shù)陣。

3 聯(lián)合迭代重構(gòu)算法

聯(lián)合迭代算法是代數(shù)重構(gòu)算法的一種改進(jìn)算法，與代數(shù)重構(gòu)算法不同的是它不再單獨(dú)對(duì)每一條觀測(cè)路徑上的水汽密度逐個(gè)進(jìn)行修正，而是一次性對(duì)所有觀測(cè)路徑進(jìn)行迭代［10］，根據(jù)每一次迭代的修正量再對(duì)對(duì)流層水汽密度分布圖像做整體修正。迭代式如下：

式中，x表示網(wǎng)格水汽密度，m表示觀測(cè)方程的行數(shù)(包括約束方程)，n表示觀測(cè)方程的列數(shù)，即參數(shù)的個(gè)數(shù)，aij表示方程系數(shù)陣第i行第j列個(gè)元素，λ為松弛因子，SWV表示常數(shù)項(xiàng)。

平均偏差MDIF：

中誤差RMS：

對(duì)于水汽的三維重構(gòu)問題，初值可以通過以下方式獲得：1)根據(jù)測(cè)站的氣象文件，內(nèi)插各個(gè)網(wǎng)格內(nèi)的氣象元素，進(jìn)而獲得網(wǎng)格初始水汽值;2)利用標(biāo)準(zhǔn)大氣直接求得初值，該方法雖然簡(jiǎn)單，但精度較差;3)依靠數(shù)值預(yù)報(bào)模型的計(jì)算結(jié)果，如NCEP/ECMWF/MM5提供的氣壓、溫度、濕度等網(wǎng)格數(shù)據(jù)，可以較準(zhǔn)確地內(nèi)插出研究區(qū)域內(nèi)網(wǎng)格的水汽密度。

4 數(shù)值模擬實(shí)驗(yàn)

4.1 實(shí)驗(yàn)概況

利用上海地區(qū)GPS監(jiān)測(cè)網(wǎng)進(jìn)行數(shù)值模擬實(shí)驗(yàn)。研究范圍水平跨度0.6°，垂直方向取10 km?？臻g網(wǎng)格水平分辨率為0.15°，垂直分辨率為500 m，東西和南北方向各分為4格，垂直方向分為20層，空間區(qū)域內(nèi)網(wǎng)格總數(shù)為320個(gè)。該區(qū)域內(nèi)包含測(cè)站9個(gè)，從平面來(lái)看平均分布在各個(gè)網(wǎng)格內(nèi)(圖1)。實(shí)驗(yàn)取2007年9月13日UTC 0：00—0：30時(shí)段內(nèi)的數(shù)據(jù)，具體模擬過程如下：

1)觀測(cè)方程系數(shù)陣。衛(wèi)星發(fā)射信號(hào)到達(dá)地面的接收機(jī)，信號(hào)穿過大氣層時(shí)在研究區(qū)域內(nèi)穿過網(wǎng)格產(chǎn)生了截距，這些距離構(gòu)成系數(shù)陣A。

2)模擬觀測(cè)值。歐洲中尺度天氣預(yù)報(bào)(ECMWF)資料以網(wǎng)格形式提供了全球大氣層的氣壓、溫度、比濕等信息。在指定的研究區(qū)域內(nèi)，沿觀測(cè)路徑穿過的網(wǎng)格內(nèi)插ECMWF氣象資料，并在該路徑上積分［6］出觀測(cè)值SWVs。在實(shí)際觀測(cè)中，觀測(cè)值受白噪聲影響存在一定誤差，而模擬觀測(cè)值是一種理想無(wú)噪聲狀態(tài)的觀測(cè)量，為了與實(shí)際情況相仿，給模擬觀測(cè)值增加一定量的隨機(jī)誤差：

3)網(wǎng)格參數(shù)。網(wǎng)格內(nèi)的水汽密度作為待估參數(shù)代入式(4)進(jìn)行解算。根據(jù)測(cè)站的氣象文件，內(nèi)插各個(gè)網(wǎng)格內(nèi)的氣壓、溫度及水汽壓，獲取網(wǎng)格水汽的初值。

4)考慮沒有觀測(cè)值通過的網(wǎng)格，給原觀測(cè)方程加上水平約束和垂直約束。

4.2 算法實(shí)驗(yàn)與結(jié)果分析

作為迭代收斂與否的關(guān)鍵因素之一，松弛因子的選取尤為重要。本文在一定范圍內(nèi)，以平均偏差MDIF和平均中誤差RMS作為迭代質(zhì)量?jī)?yōu)劣的標(biāo)準(zhǔn)對(duì)松弛因子進(jìn)行最優(yōu)搜索。如圖2所示，松弛因子λ的搜索范圍為0～0.1，當(dāng)λ從0.01增加到0.12時(shí)，迭代處于收斂狀態(tài)，同時(shí)MDIF和RMS不斷減小，當(dāng)λ=0.12時(shí)達(dá)到最小;當(dāng)松弛因子大于0.12時(shí)，迭代發(fā)散，因此λ最佳值為0.12。圖3為迭代收斂速度。從圖3可以看出，隨著迭代次數(shù)的增加，SIRT算法迅速收斂。平均偏差收斂較快，在迭代10次以后趨于穩(wěn)定，而RMS則收斂較慢，迭代100次以后趨于平緩。迭代次數(shù)為100時(shí)，MDIF達(dá)到0.09 mm，RMS為1.57 mm，之后精度變化緩慢趨于穩(wěn)定。

圖2 松弛因子對(duì)迭代質(zhì)量的影響Fig.2 Influence of the relaxation parameter on the quality of reconstruction

圖3 SIRT算法的收斂情況Fig.3 Convergence of SIRT algorithm

圖4 層析水汽值與ECMWF結(jié)果的比較Fig.4 Comparison between the vertical profile of water vapor obtained from tomography and from ECMWF

為了檢驗(yàn)SIRT算法層析水汽分布的質(zhì)量，將迭代后的網(wǎng)格水汽值與ECMWF資料進(jìn)行比較。從圖4可以看出聯(lián)合迭代重構(gòu)算法較好地給出了水汽垂直輪廓線，正確地反映出水汽在垂直方向上的變化趨勢(shì)，網(wǎng)格內(nèi)水汽值與ECMWF資料變化趨勢(shì)一致，在數(shù)值上略有上下波動(dòng)。在大氣5.5 km范圍內(nèi)層析結(jié)果小于ECMWF的結(jié)果，差值較小;5.5 km以上的網(wǎng)格層析結(jié)果則大于ECMWF的結(jié)果，差值有所增大。圖5是整個(gè)空間區(qū)域?qū)游銎拭鎴D，每幅圖均自西向東方向給出某一緯度的垂直剖面。圖上正確地反映了水汽在空間垂直分布的特性，即大部分水汽聚集在近地面3 km區(qū)域內(nèi)。隨高度的增加，水汽值下降，當(dāng)高度達(dá)到10 km時(shí)，水汽值接近于0。另外，圖5中緯度B=31.18°對(duì)應(yīng)的剖面圖出現(xiàn)了一些異常，第4層出現(xiàn)水汽突然變大的網(wǎng)格，該網(wǎng)格水汽高于周圍網(wǎng)格，不符合水汽在垂直方向的變化規(guī)律，反映出SIRT算法迭代存在不足之處。

圖5 區(qū)域水汽垂直剖面圖Fig.5 Tomographic reconstruction

5 結(jié)論

1)松弛因子起到調(diào)節(jié)修正的作用，由于它影響著迭代速度和迭代結(jié)果的質(zhì)量，因此松弛因子的正確選取尤為重要。本文給出的松弛因子選取范圍在0.01～0.12，最佳值為0.12;

2)SIRT算法比較依賴初值，合理地選取初值不僅可以提高計(jì)算速度加快算法收斂，還可以在觀測(cè)值對(duì)網(wǎng)格覆蓋不佳的情況下，獲得更為準(zhǔn)確的結(jié)果;

3)根據(jù)平均誤差MDIF和RMS作為迭代收斂條件，迭代100次時(shí)MDIF和RMS均趨于穩(wěn)定。由于SIRT是一種整體修正的方法，雖然迭代次數(shù)越多可以使結(jié)果越精確，但同時(shí)也容易造成對(duì)個(gè)別網(wǎng)格的過度修正而偏離正確的結(jié)果。

層析過程受到多方面因素影響，并且SIRT算法尚有需要改進(jìn)的空間，因此聯(lián)合迭代重構(gòu)算法在水汽重構(gòu)中的應(yīng)用還需要進(jìn)一步的研究。

致謝 感謝朱文耀和宋淑麗研究員的指導(dǎo)和提供氣象資料!

1 Flores A，A Rius and G Ruffini.4D tropospheric tomography using GPS slant wet delays，Ann［J］.Geophys.，2000，18：223-224.

2 Miidla P，Rannat K and Uba P.Tomographic approach for troposheric water vapor detection［J］.Computational Methods in Applied Mathematics，2008，(8)：263-278.

3 Lubomir P，Grandinarsky and Per Jarlemark.Ground-based GPS tomography of water vapor：analysis of simulated and real data［J］.Journal of the Meteorological Society of Japan，2004，82：551-560.

4 張雙成，等.基于Kalman濾波的斷層掃描初步層析水汽折射率分布［J］.武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2008.8，(8)：798-802.(Zhang Shuangcheng，et al.Preliminary tomography of special wet refractivity distribution based on kalman filter［J］.Geomatics and Information Science of Wuhan University，2008.8，(8)：798-802)

5 Stolle C，et al.A GPS based three-dimensional ionospheric imaging tool：Process and assessment［J］.Adv Space Res.，2006，38(11)：2 313-2 317.

6 Debao Wen，Sanzhi Liu and Pingying Tang.Tomographic reconstruction of ionospheric electron density based on constrained algebraic reconstruction technique［J］.J GPS Solut.，2010，doi：10.1007/s10291-010-0161-0.

7 Liu Sanzhi，et al.Inversion of ionosphere electron density based on a constrained simultaneous iteration reconstruction technique［J］.IEEE Transactions on Geoscience and Remote Sensing，2010，48：2 455-2 459.

8 宋淑麗，等.上海GPS網(wǎng)層析水汽三維分布改善數(shù)值預(yù)報(bào)濕度場(chǎng)［J］.科學(xué)通報(bào)，2005，(20)：2 271-22 77.(Song Shuli，et al.Water vapor tomography from GPS data and it’s application in improving numberical weather forecast［J］.Science Bulletin，2005，(20)：2 271-2 277)

9 Bust G S and Mitchell C N.History，current state，and future directions of ionospheric imaging［J］.Rev Geophys.，2008，46，RG1003，doi：10.1029/2006RG000212.

10 Heise S，et al.Integrated water vapor from IGS groundbased GPS observations：initial results from a global 5-min data set［J］.Ann Geophys.，2009，27：2 851-2 859.

APPLICATION OF SIMULTANEOUS ITERATIONS RECONSTRUCTION TECHNIQUE FOR 3D WATER VAPOR TOMOGRAPHY SYSTEM

Wang Wei and Wang Jiexian
(Department of Surveying and Geo-informatics，Tongji University，Shanghai 200092)

Simultaneous iterations reconstruction technique(SIRT)is one of the several reconstruction algorithms of the ART family.It is used widely for tomography technique because of its convenience in dealing with large sparse matrix.Thus SIRT is carried out to reconstruct the three dimension distribution of water vapor by simulation.The results show that SIRT can work effectively in water vapor tomography，get rapid convergence，and implement more easily than inversion.

SIRT;reconstruction;water vapor;tomography;iteration

1671-5942(2011)06-0100-05

2011-03-26

國(guó)家自然科學(xué)基金(40974018);國(guó)家863計(jì)劃項(xiàng)目(2009AA12Z307)

王維，女，1982年生，博士，主要從事GPS數(shù)據(jù)處理與GPS氣象學(xué)研究.E-mail：wangwei_nj@126.com

P207

A