利用NMF算法確定ZTD格網(wǎng)產(chǎn)品空間分辨率

劉志平，朱丹彤，王潛心，王 堅(jiān)

(中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院，江蘇 徐州 221116)

天頂對(duì)流層延遲(zenith troposphere delay，ZTD)包括對(duì)流層靜力延遲(zenith hydrostatic delay，ZHD)和對(duì)流層濕延遲(zenith wet delay，ZWD)，是GNSS高精度導(dǎo)航定位中的主要誤差源[1]。從2005年開(kāi)始，全球大地測(cè)量觀測(cè)系統(tǒng)(global geodetic observing system，GGOS)Atmosphere分析中心發(fā)布天頂對(duì)流層格網(wǎng)產(chǎn)品ZWD、ZHD[2-3](緯差2°×經(jīng)差2.5°)。很多學(xué)者利用GGOS對(duì)流層格網(wǎng)產(chǎn)品開(kāi)展了大量研究與應(yīng)用，包括對(duì)流層變化特征[4-5]、全球ZTD模型[6-7]、快速精密定位[8-9]、氣象研究[10]等。同時(shí)，也有文獻(xiàn)研究指出，利用0.5°分辨率的ECMWF氣象數(shù)據(jù)所計(jì)算出的ZTD結(jié)果更加接近真實(shí)情況[11]。由此可知，高空間分辨率的對(duì)流層格網(wǎng)產(chǎn)品自然有利于提高后續(xù)應(yīng)用與研究成果精度水平。然而，在實(shí)際應(yīng)用中也會(huì)加大數(shù)據(jù)存儲(chǔ)傳輸難度、降低計(jì)算實(shí)時(shí)性[12]。因此，應(yīng)綜合考慮產(chǎn)品精度、空間特征及數(shù)據(jù)量要求，定量探討對(duì)流層格網(wǎng)產(chǎn)品的最佳或合理空間分辨率。

綜上，為確定ZTD格網(wǎng)產(chǎn)品發(fā)布的空間分辨率，本文對(duì)GGOS天頂對(duì)流層和ECMWF氣象格網(wǎng)產(chǎn)品進(jìn)行研究。首先，顧及非負(fù)矩陣分解算法(non-negative matrix factorization,NMF)在數(shù)據(jù)壓縮[13]、特征提取[14-15]等方面的成功應(yīng)用，引入NMF算法提取對(duì)流層格網(wǎng)產(chǎn)品特征維數(shù)并針對(duì)該算法迭代效率低[16]等不足，提出一種基于隨機(jī)學(xué)習(xí)速率矩陣的NMF算法。其次，結(jié)合歐洲中尺度天氣預(yù)報(bào)中心(European Center for Medium Range Weather Forecasts，ECMWF)再分析資料(0.125°×0.125°)分析大氣可降水量數(shù)據(jù)(precipitable water vapor，PWV)、地表大氣壓數(shù)據(jù)(surface pressure，SP)與ZWD、ZHD相關(guān)系數(shù)，并通過(guò)ZTD經(jīng)緯向梯度值分組概率，初步確定格網(wǎng)產(chǎn)品發(fā)布的空間分辨率。再次，利用所提出的改進(jìn)NMF算法，提取不同空間分辨率PWV和SP產(chǎn)品的特征維數(shù)，并最終確定ZTD格網(wǎng)產(chǎn)品空間分辨率。最后，建議對(duì)ZTD格網(wǎng)產(chǎn)品進(jìn)行NMF分解并直接發(fā)布分解后的基矩陣W和系數(shù)矩陣H，并分析NMF發(fā)布方法的ZTD格網(wǎng)產(chǎn)品壓縮效果。

1 改進(jìn)的非負(fù)矩陣分解算法

E=A-WH

(1)

由于NMF算法是對(duì)待分解矩陣A的優(yōu)化逼近，故該算法需要在特定目標(biāo)函數(shù)的約束下進(jìn)行迭代，從而不斷逼近待分解矩陣A。目前應(yīng)用最廣泛的目標(biāo)函數(shù)是假設(shè)噪聲服從正態(tài)分布的歐幾里得距離[13,15]

(2)

式(2)分別對(duì)Wik、Hkj兩端求導(dǎo)，可得目標(biāo)函數(shù)在Wik、Hkj兩方向上的梯度

(3)

根據(jù)梯度下降法原理，設(shè)學(xué)習(xí)速率矩陣分別為λW、λH，則基矩陣W和系數(shù)矩陣H的迭代式可表示如下

(4)

在常規(guī)學(xué)習(xí)速率矩陣的基礎(chǔ)上，本文提出隨機(jī)學(xué)習(xí)速率矩陣

(5)

式中，α為隨機(jī)加速因子，α∈[0,1]。

將式(5)代入式(4)，可得W、H最終迭代估計(jì)式為

(6)

本文將式(6)稱(chēng)為基于隨機(jī)學(xué)習(xí)速率陣的NMF分解算法。其中，W、H迭代初值一般取0～1均勻分布矩陣，隨機(jī)加速因子α為[0,1]均勻分布的隨機(jī)數(shù)，且在每次迭代計(jì)算過(guò)程中獨(dú)立生成。顯見(jiàn)，當(dāng)隨機(jī)加速因子取恒定值α≡1時(shí)，式(6)即為標(biāo)準(zhǔn)的NMF分解算法[13,15]。

2 ZTD格網(wǎng)產(chǎn)品空間分辨率的確定

2.1 ZTD經(jīng)緯向梯度值概率統(tǒng)計(jì)

為初步研究天頂對(duì)流層數(shù)據(jù)的空間分布特征，本文對(duì)GGOS Atmosphere天頂對(duì)流層格網(wǎng)產(chǎn)品(緯差2°×經(jīng)差2.5°)和ECMWF氣象資料(緯差0.125°×經(jīng)差0.125°)計(jì)算等效天頂對(duì)流層數(shù)據(jù)隨緯度和經(jīng)度方向的梯度(緯向和經(jīng)向差分)，從而研究不同空間分辨率ZTD梯度絕對(duì)值的概率統(tǒng)計(jì)特征，以初步確定ZTD空間分辨率。試驗(yàn)數(shù)據(jù)選用2015年年積日1、UTC0時(shí)的GGOS Atmosphere天頂對(duì)流層格網(wǎng)數(shù)據(jù)和ECMWF提供的氣象格網(wǎng)數(shù)據(jù)PWV、SP，對(duì)上述不同分辨率的EZTD和GGOS提供的ZTD格網(wǎng)數(shù)據(jù)進(jìn)行梯度計(jì)算(方案1：EZTD緯向稀疏后計(jì)算梯度，3°×0.125°、2°×0.125°、1°×0.125°、0.5°×0.125°、0.25°×0.125°；方案2：EZTD經(jīng)向稀疏后計(jì)算梯度，0.125°×4°、0.125°×3°、0.125°×2°、0.125°×1°、0.125°×0.5°)，并統(tǒng)計(jì)梯度絕對(duì)值分組概率，所得結(jié)果見(jiàn)表1。從表1方案1的結(jié)果可看出，不同空間分辨率的EZTD格網(wǎng)產(chǎn)品的概率統(tǒng)計(jì)特征存在差異，各空間分辨率下小于3 cm的概率統(tǒng)計(jì)值分別為81.95%、80.34%、78.85%、78.29%、78.16%，可知緯差1°×經(jīng)差0.125°是緯向梯度絕對(duì)值概率統(tǒng)計(jì)特征的空間分辨率拐點(diǎn)。從方案2結(jié)果可看出，各空間分辨率下小于3 cm的概率統(tǒng)計(jì)值分別為92.20%、91.10%、90.40%、89.67%、89.50%，可得緯差0.125°×經(jīng)差2°是經(jīng)向梯度絕對(duì)值概率統(tǒng)計(jì)特征的空間分辨率拐點(diǎn)。同時(shí)，對(duì)比方案1和方案2可知，緯向梯度比徑向梯度復(fù)雜得多，這與文獻(xiàn)[7]中的結(jié)論一致。因此，初步確定緯差1°×經(jīng)差2°是保留天頂對(duì)流層格網(wǎng)產(chǎn)品特征信息的合理空間分辨率。此外，EZTD與GGOS Atmosphere提供的ZTD格網(wǎng)產(chǎn)品的分組概率統(tǒng)計(jì)值及變化趨勢(shì)比較吻合，進(jìn)一步驗(yàn)證了結(jié)合EZTD數(shù)據(jù)研究ZTD格網(wǎng)產(chǎn)品空間分辨率的合理性。

表1 梯度絕對(duì)值分組概率統(tǒng)計(jì) (%)

2.2 NMF分解與空間分辨率確定

為比較原NMF分解算法和本文所提的基于隨機(jī)學(xué)習(xí)速率矩陣的NMF分解算法的優(yōu)劣，采用2015年年積日1、UTC0時(shí)的大氣可降水量格網(wǎng)產(chǎn)品PWV(緯差2°×經(jīng)差2.5°，數(shù)據(jù)矩陣維數(shù)為91×145)，進(jìn)行NMF分解算法改進(jìn)前后的效果對(duì)比分析。其中，特征維數(shù)設(shè)置為9～91、步長(zhǎng)取2，均方根誤差和誤差下降率閾值分別為0.5 cm、1×10-7。此外，鑒于NMF分解結(jié)果的不唯一性，為使對(duì)比結(jié)果具有統(tǒng)計(jì)意義，每個(gè)特征維數(shù)進(jìn)行100次重復(fù)NMF分解計(jì)算。方案1：利用原NMF算法進(jìn)行迭代分解，統(tǒng)計(jì)迭代終止時(shí)的迭代次數(shù)和均方根誤差；方案2：利用改進(jìn)的NMF算法進(jìn)行迭代分解，統(tǒng)計(jì)迭代終止時(shí)的迭代次數(shù)和均方根誤差。計(jì)算結(jié)果如圖1所示。

圖1 兩種方法迭代次數(shù)和均方根誤差

從圖1可以看出，改進(jìn)后NMF算法迭代次數(shù)得到了極大減少、迭代終止時(shí)的均方根誤差也得到了降低，表明基于隨機(jī)學(xué)習(xí)速率矩陣的NMF算法優(yōu)于原NMF算法。其次，原NMF算法在特征維數(shù)范圍內(nèi)無(wú)法保證迭代收斂，說(shuō)明該方法不能成功提取PWV的空間特征；而改進(jìn)的NMF算法在特征維數(shù)大于某一值后，能夠保證迭代收斂。而特征維數(shù)是空間分布特征的直觀表達(dá)，通過(guò)對(duì)比不同空間分辨率下的特征維數(shù)，能夠完整表達(dá)PWV空間特征的最佳空間分辨率。因此，下面采用改進(jìn)的NMF算法進(jìn)行對(duì)流層延遲特征維數(shù)的計(jì)算分析。

為對(duì)比不同分辨率對(duì)流層延遲的特征維數(shù)，選擇2015年4個(gè)年積日(冬季DOY1、春季DOY90、夏季DOY182和秋季DOY274)、每日UTC0與UTC12時(shí)共8個(gè)時(shí)刻的可降水量產(chǎn)品PWV、地表大氣壓產(chǎn)品SP(緯差0.125°×經(jīng)差0.125°)，進(jìn)而對(duì)PWV和SP數(shù)據(jù)進(jìn)行稀疏化，得到緯差×經(jīng)差分別為2.5°×5°、2°×2.5°、2°×2°、1°×2°、1°×1°、0.5°×1°、0.5°×0.5°共7種空間分辨率的全球格網(wǎng)產(chǎn)品，分別對(duì)PWV、SP數(shù)據(jù)進(jìn)行NMF分解(均方根誤差和誤差下降率閾值分別為0.5 cm、1×10-7)并研究特征維數(shù)與空間分辨率的關(guān)系。需說(shuō)明的是，PWV與SP數(shù)據(jù)已轉(zhuǎn)換為等效ZWD與等效ZHD(單位：厘米)。表2為相應(yīng)的PWV固有特征維數(shù)統(tǒng)計(jì)結(jié)果。表3為相應(yīng)的PWV固有特征新維數(shù)統(tǒng)計(jì)結(jié)果。

表2統(tǒng)計(jì)了不同時(shí)刻、不同分辨率PWV固有特征維數(shù)。對(duì)比同一年積日不同UTC下的特征維數(shù)，可以發(fā)現(xiàn)不同UTC的空間特征基本保持不變，表現(xiàn)為全球范圍內(nèi)大氣可降水量的短時(shí)穩(wěn)定性。同時(shí)，對(duì)比同一UTC不同年積日下的空間特征，不同季節(jié)的特征維數(shù)存在細(xì)小差異：夏季的空間特征最為復(fù)雜，其固有特征維數(shù)為42，春季、冬季的固有特征維數(shù)為36，空間特征最簡(jiǎn)單，說(shuō)明不同季節(jié)大氣可降水量的空間分布特征存在差異，從側(cè)面反映了全球范圍的季節(jié)性差異。

表2 PWV固有特征維數(shù)統(tǒng)計(jì)

為了分析SP特征維數(shù)穩(wěn)定性與格網(wǎng)間距的關(guān)系，表3統(tǒng)計(jì)了不同時(shí)刻、不同空間分辨率SP固有特征維數(shù)。從該表可以得出，SP特征維數(shù)穩(wěn)定性高于PWV，不同UTC、不同季節(jié)的特征維數(shù)基本保持不變，而且隨著格網(wǎng)間距降至1°以后特征維數(shù)基本穩(wěn)定。此外，同一空間分辨率下、不同時(shí)刻之間的特征維數(shù)變現(xiàn)出較高的一致性，這也從側(cè)面說(shuō)明了地表大氣壓在全球范圍內(nèi)具有較高的穩(wěn)定性。

表3 SP固有特征維數(shù)統(tǒng)計(jì)

綜合上述分析，PWV和SP數(shù)據(jù)的最佳分辨率均為1°×2°，固有特征維數(shù)分別為42、70。需指出的是，產(chǎn)品發(fā)布若采用原有的發(fā)布形式會(huì)導(dǎo)致數(shù)據(jù)量增大(數(shù)據(jù)矩陣由91×145增大為181×181，后者為前者的2.48倍)，而數(shù)據(jù)量激增不利于存儲(chǔ)效率和傳輸速度。顧及ZWD固有特征維數(shù)為42、ZHD有特征維數(shù)為70，故經(jīng)NMF分解后的基矩陣與系數(shù)矩陣可獲得較好數(shù)據(jù)壓縮效果(詳細(xì)分析見(jiàn)表4)。因此，本文建議GGOS Atmosphere對(duì)流層格網(wǎng)產(chǎn)品的空間分辨率應(yīng)采用1°×2°，發(fā)布數(shù)據(jù)內(nèi)容為格網(wǎng)產(chǎn)品經(jīng)NMF算法分解后的基矩陣和系數(shù)矩陣。

為評(píng)價(jià)NMF的數(shù)據(jù)壓縮效果，表4統(tǒng)計(jì)了原發(fā)布方法和NMF發(fā)布新方法的數(shù)據(jù)量以及數(shù)據(jù)壓縮比。其中，原壓縮比為分辨率為1°×2°與2°×2.5°原發(fā)布方法數(shù)據(jù)量的比值，NMF壓縮比為1°×2°分辨率NMF發(fā)布方法與2°×2.5°分辨率原發(fā)布方法數(shù)據(jù)量的比值。不難看出，對(duì)流層格網(wǎng)產(chǎn)品采用1°×2°分辨率會(huì)導(dǎo)致數(shù)據(jù)量變?yōu)樵瓉?lái)的2.5倍，而采用NMF發(fā)布形式的NMF壓縮比分別為1.1、1.9，NMF發(fā)布方法數(shù)據(jù)壓縮效率達(dá)到56%、24%。

表4 原數(shù)據(jù)與NMF發(fā)布方法比較

3 結(jié) 論

(1) 針對(duì)常規(guī)非負(fù)矩陣分解算法計(jì)算效率低、難以提取固有特征維數(shù)的問(wèn)題，基于隨機(jī)學(xué)習(xí)速率矩陣，本文提出了一種改進(jìn)的非負(fù)矩陣分解算法。計(jì)算結(jié)果表明所提出的改進(jìn)算法迭代效率高、迭代收斂誤差小。

(2) 分析了ECMWF PWV、SP與GGOS ZWD、ZHD產(chǎn)品之間的強(qiáng)相關(guān)性，結(jié)合GGOS ZTD和ECMWF EZTD計(jì)算了經(jīng)向、緯向梯度值分組概率并初步確定了ZTD格網(wǎng)產(chǎn)品的空間分辨率，在此基礎(chǔ)上利用NMF算法提取了7種空間分辨率、8個(gè)時(shí)刻PWV、SP的特征維數(shù)，從而確定了ZWD和ZHD格網(wǎng)產(chǎn)品的固有特征維數(shù)和最佳空間分辨率。

(3) 建議ZTD格網(wǎng)產(chǎn)品發(fā)布的空間分辨率調(diào)整為1°×2°，以更好地描述ZTD空間分布特征。同時(shí)，針對(duì)空間分辨率提高導(dǎo)致的數(shù)據(jù)量增大問(wèn)題，推薦對(duì)產(chǎn)品進(jìn)行NMF分解并直接發(fā)布分解后的基矩陣W和系數(shù)矩陣H，可使ZWD和ZHD壓縮效率分別達(dá)到56%、24%。

參考文獻(xiàn)：

[1] LI W，YUAN Y，OU J，et al.New Versions of the BDS/GNSS Zenith Tropospheric Delay Model IGGtrop[J].Journal of Geodesy，2015，89(1):73-80.

[2] 陳俊勇，黨亞民.完善大地坐標(biāo)框架和地球重力場(chǎng)時(shí)變測(cè)量的進(jìn)展——2005年國(guó)際大地測(cè)量協(xié)會(huì)(IAG)科學(xué)大會(huì)札記[J].測(cè)繪通報(bào)，2005(12):1-4.

[3] CHEN Q，SONG S，HEISE S，et al.Assessment of ZTD Derived from ECMWF/NCEP Data with GPS ZTD over China[J].GPS Solutions，2011，15(4):415-425.

[4] ZHANG J Y，WANG L，YANG S，et al.Decadal Changes of the Wintertime Tropical Tropospheric Temperature and Their Influences on the Extratropical Climate[J].Science Bulletin，2016，61(9):737-744.

[5] YUCHECHEN A E，LAKKIS S G，LAVORATO M B.On the Stability of the Troposphere/Lower Stratosphere and Its Relationships with Cirrus Clouds and Three Mandatory Levels over Buenos Aires[J].International Journal of Remote Sensing，2016，37(7):1541-1552.

[6] 姚宜斌，何暢勇，張豹,等.一種新的全球?qū)α鲗犹祉斞舆t模型GZTD[J].地球物理學(xué)報(bào)，2013，56(7):2218-2227.

[7] 姚宜斌，胡羽豐，張豹.利用多源數(shù)據(jù)構(gòu)建全球天頂對(duì)流層延遲模型[J].科學(xué)通報(bào)，2016,61(24):2730-2741.

[8] 高旺，高成發(fā)，潘樹(shù)國(guó),等.北斗三頻寬巷組合網(wǎng)絡(luò) RTK 單歷元定位方法[J].測(cè)繪學(xué)報(bào)，2015，44(6):641-648.

[10] DOUSA J，VACLAVOVIC P.Real-time Zenith Tropospheric Delays in Support of Numerical Weather Prediction Applications[J].Advances in Space Research，2014，53(53):1347-1358.

[11] JIN S，LUO O F，GLEASON S.Characterization of Diurnal Cycles in ZTD from a Decade of Global GPS Observations[J].Journal of Geodesy，2009，83(6):537-545.

[12] 陳欽明，宋淑麗，朱文耀.亞洲地區(qū)ECMWF/NCEP資料計(jì)算ZTD的精度分析[J].地球物理學(xué)報(bào)，2012，55(5):1541-1548.

[13] LEE D D，SEUNG H S.Learning the Parts of Objects by Non-negative Matrix Factorization[J].Nature，1999，401(6755):788-91.

[14] 杜海順，張旭東,等.圖嵌入正則化投影非負(fù)矩陣分解人臉圖像特征提取[J].中國(guó)圖象圖形學(xué)報(bào)，2014，19(8):1176-1184.

[15] 施蓓琦，劉春，孫偉偉,等.應(yīng)用稀疏非負(fù)矩陣分解聚類(lèi)實(shí)現(xiàn)高光譜影像波段的優(yōu)化選擇[J].測(cè)繪學(xué)報(bào)，2013，42(3):351-358.

[16] 李樂(lè)，章毓晉.非負(fù)矩陣分解算法綜述[J].電子學(xué)報(bào)，2008，36(4):737-743.

[17] ZHAN C，LI W，OGUNBONA P.Local Representation of Faces through Extended NMF[J].Electronics Letters，2012，48(7):373-375.

[18] YAO Y，XU C，ZHANG B，et al.GTM-III：a New Global Empirical Model for Mapping Zenith Wet Delays onto Precipitable Water Vapour[J].Geophysical Journal International，2014，197(1):202-212.

[19] ZHANG D，GUO J，CHEN M，et al.Quantitative Assessment of Meteorological and Tropospheric Zenith Hydrostatic Delay Models[J].Advances in Space Research，2016，58(6):1033-1043.

[20] 苗啟廣，王寶樹(shù).圖像融合的非負(fù)矩陣分解算法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2005，17(9):2029-2032.