利用NMF算法確定ZTD格網產品空間分辨率

2018-06-04 03:12:29劉志平朱丹彤王潛心

測繪通報 2018年5期

劉志平，朱丹彤，王潛心，王堅

(中國礦業大學環境與測繪學院，江蘇徐州 221116)

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

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

1 改進的非負矩陣分解算法

E=A-WH

(1)

由于NMF算法是對待分解矩陣A的優化逼近，故該算法需要在特定目標函數的約束下進行迭代，從而不斷逼近待分解矩陣A。目前應用最廣泛的目標函數是假設噪聲服從正態分布的歐幾里得距離[13,15]

(2)

式(2)分別對Wik、Hkj兩端求導，可得目標函數在Wik、Hkj兩方向上的梯度

(3)

根據梯度下降法原理，設學習速率矩陣分別為λW、λH，則基矩陣W和系數矩陣H的迭代式可表示如下

(4)

在常規學習速率矩陣的基礎上，本文提出隨機學習速率矩陣

(5)

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

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

(6)

本文將式(6)稱為基于隨機學習速率陣的NMF分解算法。其中，W、H迭代初值一般取0～1均勻分布矩陣，隨機加速因子α為[0,1]均勻分布的隨機數，且在每次迭代計算過程中獨立生成。顯見，當隨機加速因子取恒定值α≡1時，式(6)即為標準的NMF分解算法[13,15]。

2 ZTD格網產品空間分辨率的確定

2.1 ZTD經緯向梯度值概率統計

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

表1 梯度絕對值分組概率統計 (%)

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

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

圖1 兩種方法迭代次數和均方根誤差

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

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

表2統計了不同時刻、不同分辨率PWV固有特征維數。對比同一年積日不同UTC下的特征維數，可以發現不同UTC的空間特征基本保持不變，表現為全球范圍內大氣可降水量的短時穩定性。同時，對比同一UTC不同年積日下的空間特征，不同季節的特征維數存在細小差異：夏季的空間特征最為復雜，其固有特征維數為42，春季、冬季的固有特征維數為36，空間特征最簡單，說明不同季節大氣可降水量的空間分布特征存在差異，從側面反映了全球范圍的季節性差異。

表2 PWV固有特征維數統計

為了分析SP特征維數穩定性與格網間距的關系，表3統計了不同時刻、不同空間分辨率SP固有特征維數。從該表可以得出，SP特征維數穩定性高于PWV，不同UTC、不同季節的特征維數基本保持不變，而且隨著格網間距降至1°以后特征維數基本穩定。此外，同一空間分辨率下、不同時刻之間的特征維數變現出較高的一致性，這也從側面說明了地表大氣壓在全球范圍內具有較高的穩定性。

表3 SP固有特征維數統計

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

為評價NMF的數據壓縮效果，表4統計了原發布方法和NMF發布新方法的數據量以及數據壓縮比。其中，原壓縮比為分辨率為1°×2°與2°×2.5°原發布方法數據量的比值，NMF壓縮比為1°×2°分辨率NMF發布方法與2°×2.5°分辨率原發布方法數據量的比值。不難看出，對流層格網產品采用1°×2°分辨率會導致數據量變為原來的2.5倍，而采用NMF發布形式的NMF壓縮比分別為1.1、1.9，NMF發布方法數據壓縮效率達到56%、24%。

表4 原數據與NMF發布方法比較

3 結論

(1) 針對常規非負矩陣分解算法計算效率低、難以提取固有特征維數的問題，基于隨機學習速率矩陣，本文提出了一種改進的非負矩陣分解算法。計算結果表明所提出的改進算法迭代效率高、迭代收斂誤差小。

(2) 分析了ECMWF PWV、SP與GGOS ZWD、ZHD產品之間的強相關性，結合GGOS ZTD和ECMWF EZTD計算了經向、緯向梯度值分組概率并初步確定了ZTD格網產品的空間分辨率，在此基礎上利用NMF算法提取了7種空間分辨率、8個時刻PWV、SP的特征維數，從而確定了ZWD和ZHD格網產品的固有特征維數和最佳空間分辨率。

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

參考文獻：

[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] 陳俊勇，黨亞民.完善大地坐標框架和地球重力場時變測量的進展——2005年國際大地測量協會(IAG)科學大會札記[J].測繪通報，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] 姚宜斌，何暢勇，張豹,等.一種新的全球對流層天頂延遲模型GZTD[J].地球物理學報，2013，56(7):2218-2227.

[7] 姚宜斌，胡羽豐，張豹.利用多源數據構建全球天頂對流層延遲模型[J].科學通報，2016,61(24):2730-2741.

[8] 高旺，高成發，潘樹國,等.北斗三頻寬巷組合網絡 RTK 單歷元定位方法[J].測繪學報，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] 陳欽明，宋淑麗，朱文耀.亞洲地區ECMWF/NCEP資料計算ZTD的精度分析[J].地球物理學報，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] 杜海順，張旭東,等.圖嵌入正則化投影非負矩陣分解人臉圖像特征提取[J].中國圖象圖形學報，2014，19(8):1176-1184.

[15] 施蓓琦，劉春，孫偉偉,等.應用稀疏非負矩陣分解聚類實現高光譜影像波段的優化選擇[J].測繪學報，2013，42(3):351-358.

[16] 李樂，章毓晉.非負矩陣分解算法綜述[J].電子學報，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.