基于GPS-IR的美國中西部地區NDVI時間序列反演

吳繼忠，吳 瑋

?

基于GPS-IR的美國中西部地區NDVI時間序列反演

吳繼忠，吳 瑋

（南京工業大學測繪科學與技術學院，南京211816）

基于AVHRR（advanced very high resolution radiometer）、MODIS（moderate-resolution imaging spectroradiometer）等衛星遙感影像獲取的歸一化植被指數（normalized difference vegetation index，NDVI）存在大氣噪聲、土壤背景、飽和度等固有問題。GPS（global positioning system）衛星播發的L波段信號對土壤和植被水分含量變化較為敏感，GPS-IR（GPS-interferometric reflectometry）利用測地型接收機和天線記錄GPS反射信號的變化，進而反演測站環境參數。該文研究了利用GPS-IR反演區域NDVI時間序列的方法。采用4個GPS參考站2007－2015年近9 a的連續觀測數據，由偽距和相位觀測值計算了歸一化微波反射指數（normalized microwave reflection index，NMRI），傅立葉變換顯示NMRI具有明顯的周期特性，其中年周期和半年周期分量普遍較為突出。利用三角多項式擬合剔除NMRI中由積雪和降雨引起的粗差點后，其波動與同時間段內MODIS NDVI的趨勢一致。一元線性回歸結果顯示NMRI與NDVI之間存在顯著線性相關，相關系數在0.697～0.818（<0.001），NDVI反演誤差的均方根誤差在0.059～0.079，表明GPS-IR反演區域NDVI時間序列是可行的，該研究為獲取準實時、低成本和高時間分辨率的NDVI提供了新的思路。

模型；植被；遙感；GPS-IR；歸一化植被指數；歸一化微波反射指數；反演；相關分析

0 引 言

歸一化植被指數（normalized difference vegetation index，NDVI）可以反映植被生長情況、覆蓋情況、生物量情況和植被種類情況，已廣泛用于全球植被狀態的定量化研究和應用[1-7]。迄今為止，NDVI都是基于衛星遙感影像來獲取，常用的NDVI來源于AVHRR（advanced very high resolution radiometer）、MODIS（moderate-resolutionimaging spectroradiometer）等遙感數據產品。然而，NDVI也存在固有缺陷，特別是：1）受傳感器自身因素和大氣條件的影響，NDVI數據包含較大的噪聲影響[8-12]；2）NDVI數據易受土壤質地、土壤水分含量等土壤背景噪聲的影響[13-14]；3）在植被覆蓋較高的條件下，NDVI數值易迅速飽和[15-17]。

全球定位系統（global positioning system，GPS）是一種無線電導航定位系統，GPS衛星發射的載波信號位于微波的波段，能穿透植被，且對于植被和土壤含水量的變化較為敏感，經地表植被和土壤反射后的載波信號可被GPS接收機接收并記錄。GPS-IR（GPS -interferometric reflectometry）是近些年來發展的一種新方法，其主要思想是分析反射信號物理參數的變化，進而反演出環境變化信息，尤其是與地表土壤和植被相關的信息。Larson等研究利用GPS反射信號用于測量地表積雪厚度，提出了積雪厚度的正演模模型[18-19]。Chew等研究了GPS反射信號與土壤水分含量的回歸模型[20-22]。Small等最早研究了GPS反射信號（或多路徑信號）與植被狀態之間的關系，大量實驗結果表明，GPS偽距多路徑均方根誤差（root mean square error，RMS）與植被高度、植被含水量間都存在良好的負相關性[23]。Wan等研究建立了GPS信噪比數據振幅與草類植物含水量之間的線性回歸模型，利用這一模型反演植物含水量誤差小于1 kg/m2[24]。Larson等以GPS偽距和相位觀測值為基礎，提出了歸一化微波反射指數（normalized microwave reflection index，NMRI），發現NMRI與植被水分含量之間存在較好的相關性[25]。上述研究已發現GPS反射信號與植被水分含量之間存在相關性，能否利用GPS反射信號來反演NDVI數據則未見相關研究成果。本文研究目的在于分析NMRI的時頻分布特性，評價NMRI與NDVI的相關性，探索利用GPS-IR反演NDVI的可行性。論文利用4個不同區域的GPS參考站上長時間觀測數據進行了驗證和分析，其結論為區域NDVI數據獲取提供了新的思路。

1 原理與方法

1.1 反射區與歸一化微波反射指數計算

GPS-IR使用的觀測設備是常規的測量型接收機及天線，不需要兩幅天線，也無需更改天線的朝向。若GPS接收機天線高為，反射信號覆蓋區域可用第一菲涅耳區進行描述，其形狀是由短半軸和長半軸定義的橢圓

圖1 第一菲涅耳區

Fig.1 First Fresnel zones

NMRI是評價反射信號振幅變化的一個綜合性指標，其核心是計算1載波上偽距多路徑指標MP1的RMS（root mean square）值，MP1的表示為[26]

式中1是1載波上偽距觀測值，m；11575.42 MHz，21227.60 MHz；10.19 m，20.24 m；1和2是1、2載波相位觀測值，周。與導航定位計算不同，計算MP1的RMS值不需要使用對應整周模糊度的真實值，在無周跳的情況下整周模糊度保持常數不變，其真值的絕對大小不影響RMS的計算，因此周跳的探測處理非常關鍵。采用經過改進的TurboEdit方法[27]將每顆單日觀測數據周跳探測處理，得到重新劃分的若干個“干凈”弧段，將每個弧段的MP1做去均值化處理后進行弧段合并，計算每顆衛星單日MP1的RMS值并作加權平均，最終得到單日MP1的RMS值。NMRI的計算以MP1的RMS值為基礎，其計算方法如下

（3）

式中RMSMP1是單日MP1的RMS值，max(RMSMP1)是RMSMP1序列數值由大到小排列前5%的RMSMP1的平均值，因此NMRI絕大部分值在0～1之間，有少數值為負數。

1.2 歸一化植被指數計算

歸一化植被指數由可見光波段和近紅外波段二者反射率的反差來表征植被生物量的測度，其計算公式為

式中ch1、ch2分別為紅光波段和近紅外波段經過大氣校正的地面反射率。紅光波段（波長620～670 nm）處于入射輻射的光譜吸收區，近紅外波段（波長841～876 nm）處于入射輻射的光譜反射區。NDVI的數值范圍是[-1,1]，數值越大則綠色植被越密集。

2 試 驗

2.1 數據來源與計算

GPS數據使用了美國板塊邊緣觀測計劃（Plate Boundary Observatory，PBO）中的4個參考站自2007－2015年近9 a的觀測數據，數據采樣間隔為30 s，衛星截止高度角為5°。由于PBO參考站用途是用于地殼運動監測，其觀測環境較好，所用測站周邊都是草地覆蓋且沒有高大障礙物的影響，部分參考站的觀測環境視圖如圖2所示。

4個參考站的概略情況如表1所示。P042和P048在2007－2015年間使用的是Trimble NETRS接收機。P041在2007－2012年使用的是Trimble NETRS接收機，2013年2月更換成Trimble NETR9接收機，由于不同型號接收機對多路徑信號的處理性能不同，為保持數據一致性，P041站僅使用2007－2012年的數據，類似的原因，P054使用2012－2015年的數據。將4個參考站在上述時間段內的觀測數據下載后，按照式（3）計算單日RMSMP1。

表 1 GPS 參考站概況

目前MODIS NDVI被公認為數據質量較高的植被指數產品之一，本文使用的NDVI數據來自于美國地質調查局陸地過程分布式數據中心的MODIS植被指數產品MOD13Q1，其空間分辨率為250 m，時間分辨率為16 d。利用ENVI軟件處理2007－2015年間的影像，各個參考站NDVI值取以測站為中心的3′3像素NDVI的均值。

2.2 數據分析

圖3給出了4個參考站的RMSMP1時間序列，可以看出RMSMP1序列總體上具有明顯的周期特性，這是植被生長隨季節變化的體現；同時還可以看到圖上也有明顯偏離趨勢線的粗差點，這些粗差點主要是受積雪、降雨等天氣的影響而產生的[25]。由圖3中P042和P054的變化還可以發現，2012年內有部分RMSMP1數值相對于其他年份明顯偏大，且持續時間較長，這一現象極有可能是2012年美國中西部地區遭遇極端嚴重的干旱氣候所造成的[28]，圖4給出了各個參考站2007－2015年間各年份內累計降雨量，其中2012年的降雨量分別是前5年平均降雨量的68%、46%、103%、72%，可見P041、P042和P054 3個站在2012年的降雨量均低于往年的平均值，相應的RMSMP1數值也明顯高于往年，這一特點在圖上有清晰地體現，其直接原因是在干旱條件下，地表土壤和植被水分含量很低，偽距多路徑效應增大，從而造成RMSMP1值增加[23]。

為揭示RMSMP1序列在頻率域的分布規律，利用快速傅立葉變換將離散時間序列轉換為離散頻譜，圖5是快速傅立葉變換后得到的頻譜圖，由頻譜圖可以看出，RMSMP1序列中普遍存在明顯的周期信號，P041、P048和P054三個站上有年周期信號和半年周期信號，P041和P054兩個站的年周期信號比半年周期信號突出，而P048站的半年周期信號比年周期信號突出，振幅高出20%；P042站上僅有明顯的年周期信號，其最大振幅對應的信號頻率為1.02 a-1，周期約為358 d。

由于RMSMP1容易受到積雪、降雨等影響，在計算NMRI之前需要預處理。最直接的方法是參照氣象觀測記錄，將降雨降雪時間段的數據直接剔除，這種方法只能人工進行，工作量較大。考慮到RMSMP1序列顯著的周期性，用三角多項式對RMSMP1序列進行擬合，將擬合殘差大于2倍中誤差的觀測值剔除，再次進行擬合并進行反復迭代，直到擬合殘差無超限為止，圖6顯示了經過誤差點檢測后的RMSMP1序列，對比圖3可看出粗差點有明顯的減少。圖7顯示了P042站上的NMRI，由式（3）可知NMRI與RMSMP1是線性相關，因此NMRI和RMSMP1具有相同的周期特性。

NMRI的時間分辨率為1 d，高于NDVI的時間分辨率。為便于對比，需將二者的采樣率進行統一。為此將采樣率高的NMRI進行插值處理，用三角多項式插值的方法獲得與NDVI采樣完全同步的NMRI數據。圖8a－8d分別顯示了4個參考站上NMRI和NDVI的散點分布圖。

從圖8可以看出，二者在時域內的總體變化趨勢基本一致，NMRI和NDVI的峰值、谷值出現的時間吻合較好，可初步判定兩者間存在相關性，這一特征說明由NMRI來反演NDVI是可行的。

2.3 反演模型的建立和驗證

為便于模型驗證，以時間順序將試驗數據的前60%用于建模，后40%用于模型檢驗。根據上述思路，分別將4個參考站的NMRI為自變量，NDVI為因變量，進行一元線性回歸，表2給出了回歸分析的結果。

表2 回歸分析結果

由表2可以看出不同測站上的回歸系數是不一樣的，這與數據量有關，還與參考站硬件處理反射信號的方式有關，硬件的類型乃至固件的版本都會引起RMSMP1數值的尺度變化[29-30]，圖3中不同測站上RMSMP1振幅的差異也能印證這一結論，但不同參考站上相關系數沒有明顯差異，介于0.697～0.818之間，在顯著性水平選擇為0.001的條件下，檢驗全部通過，表明NMRI和NDVI之間存在顯著相關。考慮到一元線性回歸的檢驗、檢驗和檢驗的等價性，檢驗和檢驗也必定通過。

將NMRI數據的后40%代入表2中建立的回歸模型，得到計算的NDVI，將從遙感影像獲取的NDVI值作為真實值，計算每個站上NDVI的反演誤差的均方根誤差，其大小分別為0.059、0.061、0.079、0.069，其數值與回歸分析殘差的均方根誤差基本接近。為分析誤差的區間分布規律，計算NDVI反演相對誤差（即反演誤差與真實值的比值），以20%的間隔統計每個區間內相對誤差所占的比例，其統計結果見圖9。

從圖9可以看出，NDVI相對誤差較小的所占比例較大，以±20%以內的相對誤差為例，4個站上所占比例分別為60%、75%、70%、81%，隨著相對誤差的增大，其比例越來越小，誤差分布總體上接近于正態分布，說明反演模型是有效的。

3 結果與討論

NDVI作為一種重要的遙感參數，迄今為止均來源于AVHRR、MODIS等遙感數據產品，難以克服大氣噪聲、土壤背景和飽和度等固有缺陷。本文研究了利用GPS-IR（GPS- interferometric reflectometry）反演區域NDVI數據，并通過長時間的實測數據進行了檢驗和分析，主要結論為：

1）基于GPS-IR生成的NMRI具有明顯的周期性，其中年周期和半年周期較為突出；

2）在干旱氣候條件下土壤和植被水分含量低，造成反射信號振幅增大，NMRI值相應變小；

3）NMRI和NDVI在時域內波動趨勢趨于一致，峰值、谷值出現的時間吻合較好，依此建立了NDVI反演的一元線性回歸模型，相關系數在0.728～0.776（<0.05），NDVI擬合殘差的均方根誤差在0.056～0.091；

由于NMRI的獲取更為便捷，利用GPS-IR反演區域NDVI時間序列具備可行性。與遙感衛星大尺度大范圍的影像獲取不同，GPS-IR記錄的是測站周圍地表反射信號的特征參數，其作用范圍是以測站為中心的圓形區域，有效面積為數千平方米。目前，以GPS為代表的GNSS連續運行參考站已經成為一種空間信息基礎設施，參考站的數量和密度在不斷增加，GPS-IR可獲取的觀測范圍也在不斷增大，全球范圍內大量分布的參考站可望成為潛在的NDVI傳感器，且具有低成本、準實時、高時間分辨率的優點，不受傳感器自身因素和大氣條件的影響，但GPS-IR受雨雪、地形條件的影響，GPS-IR與土壤、植被間相互作用的物理機制尚不完全明確，上述問題還需要在后續工作中進一步深入探索。

[1] 邊金虎，李愛農，宋孟強，等. MODIS 植被指數時間序列Savitzky-Golay濾波算法重構[J]. 遙感學報，2010，14(4)：725－741.

Bian Jinhu, Li Ainong, Song Mengqiang, et al. Reconstruction of NDVI time-series datasets of MODIS based on Savitzky-Golay filter[J]. Journal of Remote Sensing, 2010, 14(4): 725－741. (in Chinese with English abstract)

[2] Steven M, Malthus T, Baret F, et al. Intercalibration of vegetation indices from different sensor systems[J]. Remote Sensing of Environment, 2003, 88(4): 412－422.

[3] Lunetta R, Knight J, Edirwickrema J, et al. Land-cover change detection using multi-temporal MODIS NDVI data[J]. Remote Sensing of Environment, 2006, 105(2): 142－154.

[4] Montandon L, Small, E. The impact of soil reflectance on the quantification of the green vegetation fraction from NDVI[J]. Remote Sensing of Environment, 2008, 112(4): 1835－1845.

[5] Leon J, van L, Casady G. Using MODIS-NDVI for the modeling of post-wildfire vegetation response as a function of environmental conditions and pre-fire restoration treatments[J]. Remote Sensing, 2012, 4(3): 598－621.

[6] 杜加強，趙晨曦，賈爾恒·阿哈提，等. 近30 a新疆月NDVI動態變化及其驅動因子分析[J]. 農業工程學報，2016，32(5)：172－181.

Du Jiaqiang, Zhao Chenxi, Jiaerheng Ahati, et al. Analysis on spatio-temporal trends and drivers in monthly NDVI during recent decades in Xinjiang, China based two datasets [J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(5): 172－181. (in Chinese with English abstract)

[7] Otto M, H?pfner C, Curio J, et al. Assessing vegetation response to precipitation in northwest Morocco during the last decade: an application of MODIS NDVI and high resolution reanalysis data[J]. Theoretical and Applied Climatology, 2016, 123(1): 23－41.

[8] Roerink G J, Menenti M,Verhoef W. Reconstructing cloud free NDVI composites using Fourier analysis of time series[J]. International Journal of Remote Sensing, 2000, 21(9): 1911－1917.

[9] Atzberger C, Eilers P. Evaluating the effectiveness of smoothing algorithms in the absence of ground reference measurements [J]. International Journal of Remote Sensing, 2011, 32(13): 3689－3709.

[10] Chen J, J?nssonc P, Tamurab M, et al. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golayfilter[J]. Remote Sensing of Environment, 2004, 91: 332－344.

[11] Jennifer N, Gregory J. Noise reduction of NDVI time series: An empirical comparison of selected techniques [J]. Remote Sensing of Environment, 2009, 113(1): 248－258.

[12] 吳文斌，楊鵬，唐華俊，等. 兩種NDVI 時間序列數據擬合方法比較[J]. 農業工程學報，2009，25(11)：183－188.

Wu Wenbin, Yang Peng, Tang Huajun, et al.Comparison of two fitting methods of NDVI time series datasets [J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2009, 25(11): 183－188. (in Chinese with English abstract)

[13] Cheng Tao, Rivard B, Sanchez-Azofeifa A. Spectroscopic determination of leaf water content using continuous wavelet analysis [J]. Remote Sensing of Environment, 2011, 115(2): 659－670.

[14] Libonati R, Dacamara C, Pereira C, et al. On a new coordinate system for improved discrimination of vegetation and burned areas using MIR/NIR information [J]. Remote Sensing of Environment, 2011, 115(6): 1464－1477.

[15] 齊述華，王長耀，牛錚. 利用溫度植被旱情指數（TVDI）進行全國旱情監測研究[J]. 遙感學報，2003，7(5)：420－427.

Qi Shuhua, Wang Changyao, Niu Zheng. Evaluating soil moisture status in china using the temperature/vegetation dryness index (TVDI)[J]. Journal of Remote Sensing, 2003, 7(5): 420－427. (in Chinese with English abstract)

[16] Sandholt I, Rasmussen K, Andersen J. A simple interpretation of the surface temperaturevegetation index space for assessment of surface moisture status[J]. Remote Sensing of Environment, 2002, 79(2): 213－224.

[17] Liu Fang, Qin Qiming, Zhan Zhiming. A novel dynamic stretching solution to eliminate saturation effect in NDVI and its application in drought monitoring[J]. Chinese Geographical Science, 2012, 22(6): 683－694.

[18] Larson K, Nievinski F. GPS snow sensing: Results from the earthscope plate boundary observatory[J]. GPS Solutions, 2013, 17(3): 41－52.

[19] Nievinski F, Larson K. Forward modeling for near-surface reflectometry and positioning applications[J]. GPS Solution, 2014, 18(6): 309－322.

[20] Chew C, Small E, Larson K, et al. Effects of near-surface soil moisture on GPS snr data: Development of a retrieval algorithm for volumetric soil moisture[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(1): 537－543.

[21] Larson K, Small E, Gutmann E, et al. Using GPS multipath to measure soil moisture fluctuations: initial results[J]. GPS Solutions, 2008, 12(3): 173－177.

[22] 敖敏思，朱建軍，胡友健，等. 利用SNR觀測值進行GPS 土壤濕度監測[J]. 武漢大學學報：信息科學版，2015，40(1)：117－120.

Ao Minsi, Zhu Jianjun, Hu Youjian, et al. Comparative experiments on soil moisture monitoring with GPS SNR observations[J]. Geomatics and Information Science of Wuhan University, 2015, 40(1): 117－120. (in Chinese with English abstract)

[23] Small E, Larson K, Braun J. Sensing vegetation growth with reflected GPS signals[J]. Geophysical Research Letters, 2010, 37(12): 1－5.

[24] Wan Wei, Larson K. Using geodetic GPS receivers to measure vegetation water content[J]. GPS Solution, 2015, 19(2): 237－248.

[25] Larson K, Small E. Normalized microwave reflection index, part 1: A vegetation measurement derived from GPS networks[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(5): 1501－1511.

[26] Estey L, Meertens C. TEQC: The multi-purpose toolkit for GPS/GLONASS data[J]. GPS Solutions, 1999, 3(1): 42－49.

[27] 吳繼忠，施闖，方榮新. TurboEdit單站GPS數據周跳探測方法的改進[J]. 武漢大學學報：信息科學版，2011，36(1)：29－33.

Wu Jizhong, Shi Chuang, Fang Rongxin. Improving the single station data cycle slip detection approach turboedit[J]. Geomatics and Information Science of Wuhan University, 2011, 36(1): 29－33. (in Chinese with English abstract)

[28] Mallya G, Zhao L, Song X, et al. 2012 Midwest drought in the United States[J]. Journal of Hydrologic Engineering, 2013, 18(7): 737－745.

[29] Bilich A, Larson K. Mapping the GPS multipath environment using the signal-to-noise ratio (SNR)[J]. Radio Science, 2007, 42(6): 1－16.

[30] Vázquez G, Bennett R, Spinler J. Assessment of pseudorange multipath at continuous GPS stations in Mexico[J]. Positioning, 2013, 4(3): 253－265.

Retrieving NDVI in midwestern America using GPS-interferometric reflectometry

Wu Jizhong, Wu Wei

(,,211816)

the NDVI (normalized difference vegetation index) data, routinely derived from the AVHRR (advanced very high resolution radiometer) or MODIS (moderate resolution imaging spectroradiometer) imagery, is a key indicator of vegetation status and a useful parameter in studies of terrestrial vegetation cover, it has been widely used in remote sensing studies to reflect regional and global vegetation dynamics. However, the inherent defects of NDVI, including the atmospheric noise, soil effects and saturation problems are unavoidable, and thus impede further analysis and have a risk to generating erroneous results. Global Positioning System-Interferometric Reflectometry (GPS-IR) is a bistatic radar remote sensing technique that relates temporal changes in reflected GPS signals to changes in environmental parameters surrounding a ground-based GPS site. All GPS satellites transmit signals at L-band, which is similar to those used in active microwave radar applications. L-band signals have a higher correlation with vegetation water content, therefore GPS reflections will be sensitive to water within and on the surface of vegetation, as well as water in soil and snow. The sensing footprint of GPS-IR is on the order of a thousand square meters, which depends on the antenna height and satellite elevation angle. Other than specially-designed antenna or receiver in order to estimate environmental parameters, GPS-IR utilizes geodetic-quality GPS receivers and antennas, which are currently used at many of the already-existing GPS stations. This article presents a new method to retrieve regional NDVI data using NMRI (normalized microwave reflection index), which is an index derived from GPS observations. An experiment was conducted to evaluate the feasibility of the NDVI retrieval using NMRI. In the experiment, continuous GPS observations of four plate boundary observatory GPS reference stations in midwestern America during the interval 2008-2012 and MOD13Q1 product within the same time from MODIS were used. In the first step, the NMRI time series were calculated with the GPS pseudoranges and carrier phase observations preprocessed with an improved Turboedit method, and then NDVI time series were extracted from MOD13Q1 product. In the second step, NMRI and NDVI were compared and analyzed. The temporal fluctuations of NMRI showed a clear periodicity as well as sudden drops, which were not compatible with the gradual process of vegetation change. Fast Fourier transform revealed that the annual and semi-annual periodicities exhibited dominant amplitude. To obtain cleaned NMRI data, trigonometric polynomial fitting method was adopted to remove outliers. A relatively high correlation coefficient between NMRI and NDVI was found, the coefficients of determination varied from 0.697 to 0.818 (with a significance level of<0.001), showing a near linear relationship involving these variables. With regression analysis, a linear retrieve model for NDVI could be established on each reference station, the root mean square of NDVI retrieve errors varied from 0.059 to 0.079. The outcomes of this study suggested that GPS-IR would be almost equally capable of retrieving regional NDVI data, in contrast, GPS-IR had the potential to be in near real time, with low price and high temporal resolution, and what’s more, existing GPS networks around the world had the potential to be the NDVI sensors, which could be regarded as a new opportunity to obtain NDVI data.

models; vegetation; remote sensing; GPS-interferometric reflectometry; normalized microwave reflection index; normalized difference vegetation index; retrieve; correlation analysis

10.11975/j.issn.1002-6819.2016.24.024

P228.4; P237.9

A

1002-6819(2016)-24-0183-06

2016-05-01

2016-11-29

國家自然科學基金資助項目（41504024）；江蘇省測繪地理信息科研項目（JSCHKY201413）

吳繼忠，男，湖北紅安人，博士，副教授，碩士生導師，主要從事衛星導航定位應用研究。南京 南京工業大學測繪科學與技術學院，211816。Email：jzwumail@163.com