基于動態(tài)閾值的HJ－1B圖像云檢測算法研究

韓 杰，楊磊庫，李慧芳，梁洪有，馬曉紅，謝玉娟

(河南理工大學(xué)礦山空間信息技術(shù)國家測繪局重點實驗室，焦作 454000)

基于動態(tài)閾值的HJ－1B圖像云檢測算法研究

韓 杰，楊磊庫，李慧芳，梁洪有，馬曉紅，謝玉娟

(河南理工大學(xué)礦山空間信息技術(shù)國家測繪局重點實驗室，焦作 454000)

針對環(huán)境與災(zāi)害監(jiān)測預(yù)報小衛(wèi)星星座1B星(HJ－1B)的傳感器波段特征，通過對云檢測閾值的年度波動性分析，提出了一種基于波譜標(biāo)準(zhǔn)差異常的“動態(tài)閾值云檢測算法”。該算法的實現(xiàn)主要是采用圖像配準(zhǔn)、波段運算、線性回歸和誤差分析等方法，獲取云異常區(qū)域，進而剔除圖像中的云像元。研究結(jié)果表明，該算法可以很好地檢測出不同時相、不同類型下墊面上空的云像元，為有效利用HJ－1B數(shù)據(jù)并提高其圖像分類精度發(fā)揮作用。

標(biāo)準(zhǔn)差異常;動態(tài)閾值;云檢測;環(huán)境與災(zāi)害監(jiān)測預(yù)報小衛(wèi)星星座1B星(HJ－1B)

0 引言

環(huán)境與災(zāi)害監(jiān)測預(yù)報小衛(wèi)星星座1B星(以下簡稱HJ－1B)的圖像憑借其獲取方便、監(jiān)測區(qū)域大、分辨率較高等特點，已在生態(tài)破壞、環(huán)境污染檢測和災(zāi)后評估等諸多領(lǐng)域得到廣泛應(yīng)用［1－2］。然而因受云層遮蓋或云陰影遮擋的干擾，HJ－1B圖像中常常會出現(xiàn)許多“盲區(qū)”，使該區(qū)域內(nèi)的地物信息丟失或變得模糊，嚴(yán)重影響到地表和大氣各種參數(shù)反演的準(zhǔn)確性和可靠性。因此，在遙感定量化應(yīng)用前，采用快速有效的云檢測方法剔除圖像中云覆蓋區(qū)域，是圖像預(yù)處理中的一個重要環(huán)節(jié)［3］。

目前，關(guān)于云檢測的算法很多，例如 ISCCP(international satellite cloud climatology project)方法［4］、APPOLLO(AVHRR processing scheme over cloud land and ocean)方法［5］、CLAVR(NOAA cloud advanced very high resolution radiometer) 方法［6］等，這些算法主要是利用預(yù)先設(shè)定的經(jīng)驗閾值進行云檢測。然而由于不同時相、不同高度云層的反射率與亮溫有所不同，利用固定閾值判識每景圖像中的云像元是不可靠的，會導(dǎo)致檢測的精度好壞不一［7］。

為解決上述問題，本文針對HJ－1B CCD/IRS的波段設(shè)置，提出了一種基于波譜標(biāo)準(zhǔn)差異常的“動態(tài)閾值云檢測算法”。該算法將云覆蓋區(qū)域作為超出晴空輻射值正常變化范圍的異常區(qū)域，利用該區(qū)域波譜標(biāo)準(zhǔn)差異常特征獲取動態(tài)閾值，針對不同類型下墊面逐步建立云掩模。研究結(jié)果表明，該算法能夠準(zhǔn)確地剔除HJ－1B圖像中的云像元，具有一定的通用性。

1 研究區(qū)與數(shù)據(jù)源

1.1 研究區(qū)概況

研究區(qū)位于河南省焦作市的西北部，E111°07'～113°44'，N34°08'～ 35°04'之間，海拔在85～940 m之間;屬暖溫帶大陸性季風(fēng)氣候，四季分明。該區(qū)域的下墊面類型較為復(fù)雜，包括植被、土壤、城市用地、水體等，具有較好的代表性，適用于檢驗本文算法的通用性。

1.2 HJ－1B 數(shù)據(jù)介紹

HJ－1B衛(wèi)星平臺上裝載了一臺CCD相機和一臺紅外(IRS)相機，其波段設(shè)置及空間分辨率等參數(shù)［8］如表1所示。本研究中所使用算法測試數(shù)據(jù)的成像時間為2010年11月12日。

表1 HJ－1B衛(wèi)星相機主要參數(shù)Tab.1 Principal parameters of HJ－1B camera

2 特征波段閾值波動性分析

云檢測的基本原理是利用云與植被、土壤以及水域等下墊面介質(zhì)的光譜差異進行云檢測的。云在可見光范圍有較高的反射率，但在熱紅外波段的亮溫較低［9］。針對HJ－1B波段設(shè)置，結(jié)合以往相關(guān)研究結(jié)論［10－12］，本文選取第 3 波段(紅光波段)和第8波段(熱紅外波段)作為特征波段進行云檢測。

以往大多數(shù)研究是利用預(yù)先設(shè)定的閾值判識云像元的，即將圖像中可見光波段反射率大于某一閾值或者熱紅外波段亮溫小于某一閾值的像元作為云像元［13］。為判斷特征波段閾值是否穩(wěn)定，本文選取研究區(qū)不同季度獲取的HJ－1B圖像進行云檢測閾值的波動性分析。從每季度圖像中挑選10 d數(shù)據(jù)，利用ENVI中區(qū)域直方圖統(tǒng)計及密度分割技術(shù)，結(jié)合目視解譯，獲取各特征波段的最佳閾值。其中，紅光波段反射率閾值及熱紅外波段亮溫閾值的波動曲線如圖1所示。

圖1 紅波段反射率(左)和熱紅外波段亮溫(右)閾值波動曲線Fig.1 Threshold volatility curves of reflectance of red band(left)and brightness temperature of infrared band(right)

從圖1可以看出，特征波段閾值變化主要有2個特征:①不同時相的云檢測閾值是不穩(wěn)定的(這主要是因為不同季節(jié)、不同高度云層的反射率及亮溫值有所差異，進而導(dǎo)致閾值的波動);②亮溫曲線表現(xiàn)出明顯升溫降溫現(xiàn)象(主要與太陽輻射時間和強度有關(guān))。根據(jù)斯特藩－玻爾茲曼(Stefan－Boltzmann)定律可知:在比輻射率一定的前提下，輻射總量與物體溫度的4次方成正比。同種物體的比輻射率雖然會隨著波長的變化而變化，但變化很小。因此，云在不同光譜波段向外輻射的能量都會增加，且在熱紅外波段(11 μm處)大氣水汽吸收的能量比較少，導(dǎo)致傳感器接收到的熱輻射增多，使云層表現(xiàn)出亮溫增高的現(xiàn)象。冬季該過程減弱，則表現(xiàn)出降溫的趨勢。

從特征波段閾值年度波動性分析結(jié)果可以看出，幾乎每一景圖像的云檢測閾值都是不同的，因此要想利用統(tǒng)一、固定的閾值即可很好地識別任意一景圖像中的云像元是不可能的。

3 基于波譜標(biāo)準(zhǔn)差異常的動態(tài)閾值云檢測方法

本文通過對特征波段閾值波動性分析，提出了一種基于波譜標(biāo)準(zhǔn)差異常的“動態(tài)閾值云檢測方法”。該方法的基本思路是:首先對同一區(qū)域多日晴空圖像進行輻射定標(biāo)、重采樣、裁剪、圖像配準(zhǔn)等預(yù)處理;然后將不同時相特征波段疊加在一起，形成紅波段反射率和熱紅外波段亮溫數(shù)據(jù)集，利用ENVI獲取2個數(shù)據(jù)集的均值圖像(晴空背景場)和標(biāo)準(zhǔn)差圖像;分析待檢測圖像與背景場的差異，確定因云遮蓋或云陰影遮擋而產(chǎn)生的異常特征，進而利用標(biāo)準(zhǔn)差圖像及圖像誤差分析方法獲取動態(tài)閾值;最后針對不同下墊面類型，采用分類別逐像元搜索算法實現(xiàn)云掩模的自動生成。具體流程見圖2。

圖2 基于波譜標(biāo)準(zhǔn)差異常的動態(tài)閾值云檢測流程Fig.2 Flow chart of dynamic threshold cloud detection based on spectral standard deviation anomaly

3.1 數(shù)據(jù)預(yù)處理

在綜合考慮HJ－1B的重訪周期及地物覆蓋類型變化的基礎(chǔ)上，本文將季度作為時間跨度，建立季度背景場。由于積雪在云檢測中是一個較難判識的地物，為檢驗本文算法的有效性，選用研究區(qū)冬季的晴空圖像建立冬季晴空背景場。最終獲取的特征波段均值、標(biāo)準(zhǔn)差圖像如圖3所示。

圖3 特征波段均值及標(biāo)準(zhǔn)差圖像Fig.3 Mean and standard deviation images of the feature bands

本文所選測試圖像的特征波段如圖4所示。

圖4 待檢測圖像紅波段反射率(左)和熱紅外波段亮溫(右)Fig.4 Reflectance of red band(left)and brightness temperature of infrared band(right)of the image to be detected

3.2 線性回歸分析

在利用背景場進行去云處理前，分析背景場與待檢測圖像特征波段在無云區(qū)域的相關(guān)性是十分必要的［14］。特征波段之間的線性回歸公式為

式中:Xi為均值背景場第i個特征波段;Yi為待檢測圖像第i個特征波段;i為特征波段數(shù)(i=1，2);ai和bi為回歸系數(shù)。線性回歸分析的結(jié)果見表2。

表2 線性回歸分析結(jié)果Tab.2 Result of HJ －1B Camera linear regression analysis

從表2可以看出，特征波段在晴空區(qū)域的相關(guān)性較高，說明該背景場可以作為云檢測的參考底圖。

3.3 差值圖像分析

將待檢測圖像與背景場對應(yīng)波段做差值運算，獲取HJ－1B各特征波段的差值圖像(圖5)。

圖5 紅波段反射率(左)和熱紅外波段亮溫(右)差值圖像Fig.5 D －value images of reflectance of red band(left)and brightness temperature of infrared band(right)

通過對圖5目視判讀可以得出:①在紅波段反射率差值圖像的左上角出現(xiàn)了明顯且集中的“正值區(qū)”，在熱紅外波段亮溫差值圖像的左上角則出現(xiàn)位置大致相同的“負(fù)值區(qū)”;②其他區(qū)域差值較為接近。

為了定量分析上述現(xiàn)象，本文從特征波段差值圖像中隨機挑選1行(第204行)數(shù)據(jù)進行分析，其剖面圖如圖6所示。

圖6 紅波段反射率(左)和熱紅外波段亮溫(右)差值圖像(第204行)剖面圖Fig.6 Profiles of D －value images(line 204)of reflectance of red band(left)and brightness temperature of infrared band(right)

從圖6可以看出，在特征波段差值圖像的剖面曲線上出現(xiàn)了明顯的波峰與波谷;其他區(qū)域差值雖有一定的波動，但波動不大。此現(xiàn)象與目視解譯的結(jié)論相吻合。

3.4 云區(qū)波譜異常特征分析

現(xiàn)實中云層遮蓋和云陰影遮擋都會產(chǎn)生異常信息，因此在提取云異常區(qū)域前，需要分析和判斷因云層遮蓋或云陰影遮擋而產(chǎn)生的異常特征。

1)對于待檢測圖像中的云陰影而言，該區(qū)域受云陰影的遮擋，使得在成像時刻其反射率比晴空區(qū)域低，陰影異常會出現(xiàn)在紅波段反射率差值圖像中的“負(fù)值區(qū)”。但陰影區(qū)的溫度不會因陰影的瞬間遮擋而有大幅度降低，同樣由斯特藩－玻爾茲曼定律可知其亮溫也不會有較大幅度下降，因此該區(qū)域在熱紅外波段亮溫差值圖像中不會表現(xiàn)出明顯的異常。

2)當(dāng)待檢測圖像中存在云層遮蓋時，該區(qū)域的紅波段反射率較高、熱紅外波段亮溫較低，而背景場對應(yīng)區(qū)域的紅波段反射率較低、熱紅外波段亮溫較高，因此因云遮蓋引起的異常在紅波段反射率和熱紅外波段亮溫差值圖像中應(yīng)分別出現(xiàn)在“正值區(qū)”和“負(fù)值區(qū)”。

通過對不同因素引起的異常特征分析，確定了云異常在特征波段差值圖像中出現(xiàn)的區(qū)域。然而如何選取合適的閾值提取該異常區(qū)域，還需要進一步分析。

3.5 云異常信息提取

提取云異常信息的過程主要考慮晴空輻射值的正常變化及HJ－1B圖像自身誤差的允許范圍兩方面。對于晴空輻射值而言，由于太陽照射強度、大氣中水汽含量、氣溶膠濃度等外界因素的變化，會導(dǎo)致晴空反射率、亮溫發(fā)生正常的改變。本文通過對研究區(qū)不同時相的HJ－1B差值圖像和標(biāo)準(zhǔn)差圖像進行比較，利用ENVI軟件中區(qū)域直方圖統(tǒng)計、剖面圖分析和波段運算等方法分析兩者在云異常區(qū)域的光譜特征，最終選取2倍標(biāo)準(zhǔn)差圖像作為界定晴空正常變化范圍的閾值(即超出該閾值的區(qū)域被判識為異常區(qū)域)。實際上，該閾值是動態(tài)的，因為隨著檢測季節(jié)和檢測區(qū)域的變化，HJ－1B差值圖像和標(biāo)準(zhǔn)差圖像都會有相應(yīng)的變化。

然而，在衛(wèi)星發(fā)射過程中及發(fā)射升空后，受周圍環(huán)境、運動狀態(tài)、光譜響應(yīng)等變化的影響，傳感器元件的性能也會改變，其靈敏度會下降，致使測量值與真實值存在一定的誤差［15］。因此，在獲取超出晴空正常變化范圍的異常圖像后，還要考慮該圖像自身誤差的允許范圍。本文參考測量平差中的誤差分析方法［16］，取置信度在95.4%(即超出該圖像自身2倍標(biāo)準(zhǔn)差)的區(qū)域作為最終的云異常區(qū)域(如圖7所示，圖7中μ與δ分別代表HJ－1B圖像均值與標(biāo)準(zhǔn)差)。由于各景異常圖像的均值、標(biāo)準(zhǔn)差不同，因此該閾值也是動態(tài)變化的。

圖7 HJ－1B圖像誤差允許范圍示意圖Fig.7 Sketch map of acceptable range of HJ－1B image error

根據(jù)上述獲取云異常信息的過程，提出了基于波譜標(biāo)準(zhǔn)差異常的“動態(tài)閾值云檢測模型”，即

式中:DIr，SDIr和 CMIr分別代表 HJ－1B紅波段反射率的差值圖像、標(biāo)準(zhǔn)差圖像和云異常圖像;DIbt，SDIbt和CMIbt分別代表熱紅外波段亮溫的差值圖像、標(biāo)準(zhǔn)差圖像和云異常圖像;mean()與std()函數(shù)分別用來獲取HJ－1B圖像的均值與方差;GE()為紅波段云異常判識函數(shù)，它將圖像中大于或等于某一數(shù)值的區(qū)域判識為1，其他區(qū)域判識為0;LE()為熱紅外波段亮溫云異常判識函數(shù)，它將圖像中小于或等于某一數(shù)值的區(qū)域判識為1，其他區(qū)域判識為0。

利用上述云檢測模型，最終獲得2個特征波段的云異常圖像(圖8)，其中白色為云異常區(qū)域(值為1)，黑色為晴空區(qū)(值為0)。

圖8 紅波段反射率(左)和熱紅外波段亮溫(右)云異常區(qū)域Fig.8 Cloud abnormal regions of reflectance of red band(left)and brightness temperature of infrared band(right)

3.6 云掩模生成

雖然2個特征波段的云異常區(qū)域已被提取出來，但不能直接取兩者的交集作為最終的云掩模，因為不同類型下墊面會對云異常提取的結(jié)果產(chǎn)生較大影響。

1)當(dāng)云像元位于水體、植被等暗區(qū)域中時，該像元在HJ－1B紅波段反射率差值圖像中會出現(xiàn)明顯的“正值區(qū)”，從而會被準(zhǔn)確地判識為云像元。

2)當(dāng)云像元位于積雪、城市用地等亮區(qū)域中時，由于積雪及云在HJ－1B可見光波段的反射率都較高，在待檢測圖像與背景場相減后，云像元與背景場中亮地物的反射率相抵消，因而會被誤判為晴空區(qū)。在這種情況下，要利用云在亮溫差值圖像中表現(xiàn)出的明顯“負(fù)值區(qū)”進行分析。

通過上述分析，結(jié)合各種典型地物的光譜特性及紅波段閾值的波動性分析結(jié)果，采用分類別逐像元搜索算法生成最終云掩模。當(dāng)背景場像元的紅波段反射率ρ≤0.18時，該像元被劃定為暗地物，這時若紅波段異常結(jié)果中對應(yīng)像元為異常，則該像元被判識為云像元;當(dāng)ρ＞0.18時，該像元被判識為積雪、沙土等亮地物，這時若亮溫結(jié)果中為異常像元，該像元也被認(rèn)定為云像元。利用該算法對整景圖像進行逐像元搜索檢測，采用IDL生成最終云掩模(圖9)。

圖9 云掩模生成流程圖Fig.9 Flow chart of cloud mask formation

4 檢測結(jié)果及精度分析

根據(jù)云像元判識流程獲取待檢測HJ－1B圖像的云掩模，利用該掩模對待檢測圖像進行去云處理，結(jié)果如圖10所示。

圖10 HJ－1B的B3(R)，B2(G)，B1(B)模擬真彩色合成圖像(左)和去云結(jié)果圖像(右)Fig.10 HJ －1B simulated true color image composed of B3(R)，B2(G)，B1(B)(left)and corresponding cloud removed result(right)

經(jīng)過目視解譯可以清晰地看出，該方法能將較亮的城市用地與云像元分開，可有效地去除厚云覆蓋區(qū)域，但在薄云處會出現(xiàn)局部漏判現(xiàn)象。為了更好地評價該算法的適應(yīng)性，本文又在冬季HJ－1B圖像中隨機選取其他幾天獲取的圖像進行去云處理，得出了相似的結(jié)論。與此同時，本文利用ENVI中評價分類精度的方法，對所有的檢測結(jié)果進行定量分析，在原始圖像和去云圖像中對每類地物隨機選取135個采樣點，將相同點進行比較，其統(tǒng)計結(jié)果如表3所示。

表3 云檢測結(jié)果精度分析Tab.3 Accuracy assessment of cloud detection result

從表3可以看出，上述方法在薄云區(qū)域出現(xiàn)有較多的誤判，這主要是由于傳感器所記錄的薄云信息在HJ－1B圖像的任何波段表現(xiàn)出來的響應(yīng)程度均不明顯［17］。但從總體上看，云像元檢測精度和無云像元檢測精度都在90%以上，可以滿足一般應(yīng)用的需要。

5 結(jié)論

1)基于波譜標(biāo)準(zhǔn)差異常的云檢測特征波段選擇、動態(tài)閾值獲取和云檢測的方法是可靠的，可以準(zhǔn)確地剔除HJ－1B遙感圖像中水體、城市用地及土壤等下墊面上空云像元。

2)本文提出的基于波譜標(biāo)準(zhǔn)差異常的“動態(tài)閾值云檢測算法”對于厚云檢測效果較好，但對薄云的檢測效果仍不夠理想。在以后的研究中，建議嘗試?yán)?S輻射傳輸模型對剔除厚云區(qū)域后的圖像進行大氣校正處理，以降低薄云對遙感定量化應(yīng)用的負(fù)面影響。

3)在評價云檢測精度方面，現(xiàn)存的大多數(shù)方法是通過在去云圖像中選取隨機點與原始圖像進行對比分析、經(jīng)目視解譯和精度統(tǒng)計后評價檢測結(jié)果的。但如何判斷檢測結(jié)果的真實性仍是現(xiàn)階段研究的一個難點，今后應(yīng)在這方面開展相關(guān)的研究工作。

［1］易 玲，汪 瀟，劉 斌.HJ－1衛(wèi)星數(shù)據(jù)質(zhì)量及其在土地利用中的應(yīng)用研究［J］.國土資源遙感，2009(3):74 －77.

［2］孫中平，熊文成，魏 斌，等.環(huán)境一號衛(wèi)星CCD影像質(zhì)量評價研究［J］.紅外，2010，31(9):30 －36.

［3］單 娜，鄭天垚，王貞松.快速高準(zhǔn)確度云檢測算法及其應(yīng)用［J］.遙感學(xué)報，2009，13(6):1147 －1155.

［4］Rossow W B，Mosher F，Kinsella E，et al.ISCCP Cloud Algorithm Intercomparison［J］.Journal of Climate and Applied Meteorology，1985，24(9):877 －903.

［5］Kriebel K T，Gesell G，K?stner M，et al.The Cloud Analysis Tool APOLLO:Improvements and Validations［J］.International Journal of Remote Sensing，2003，24(12):2389 －2048.

［6］Stowe L L，McClain E P，Carey R，et al.Global Distribution of Cloud Cover Derived from NOAA/AVHRR Operational Satellite Data［J］.Adv Space Res，1991，11(3):51 －54.

［7］劉 健.FY－2云檢測中動態(tài)閾值提取技術(shù)改進方法研究［J］.紅外與毫米波學(xué)報，2010，29(4):288 －292.

［8］任 平，楊存健，周介銘.HJ－1A/B星CCD多光譜遙感數(shù)據(jù)特征評價及應(yīng)用研究［J］.遙感技術(shù)與應(yīng)用，2010，25(1):138－142.

［9］何全軍，曹 靜，黃 江，等.基于多光譜綜合的MODIS數(shù)據(jù)云檢測研究［J］.國土資源遙感，2006(3):19－22.

［10］楊昌軍，許健民，趙鳳生.時間序列在FY2C云檢測中的應(yīng)用［J］.大氣與環(huán)境光學(xué)學(xué)報，2008，3(5):377 －391.

［11］郭洪濤，王 毅，劉向培，等.衛(wèi)星云圖云檢測的一種綜合優(yōu)化方法［J］.解放軍理工大學(xué)學(xué)報:自然科學(xué)版，2010，11(2):221－227.

［12］劉顯通，劉 奇，傅云飛，等.基于TRMM VIRS可見光和紅外五通道的白天云檢測方案［J］.大氣與環(huán)境光學(xué)學(xué)報，2010，5(2):128－140.

［13］宋小寧，趙英時.MODIS圖象的云檢測及分析［J］.中國圖象圖形學(xué)報(A 輯)，2003，8(9):1079 －1082.

［14］李炳燮，馬張寶，齊清文，等.Landsat TM遙感影像中厚云和陰影去除［J］.遙感學(xué)報，2010，14(3):1 －6.

［15］鞏 慧.HJ－1星CCD相機在軌輻射定標(biāo)與真實性檢驗研究［D］.北京:中國科學(xué)院遙感應(yīng)用研究所，2010.

［16］葛永慧.測量平差［M］.北京:中國礦業(yè)大學(xué)出版社，2005.

［17］張從容，張 正，張為良，等.利用衛(wèi)星遙感進行海上透明云薄云的檢測［J］.海洋預(yù)報，2005，22(z1):87 －93.

Research on Algorithm of Cloud Detection for HJ－1B Image Based on Dynamical Thresholding

HAN Jie，YANG Lei－ku，LI Hui－fang，LIANG Hong－you，MA Xiao－h(huán)ong，XIE Yu－juan
(State Bureau of Surveying and Mapping Key Laboratory of Mine Spatial Information，Henan Polytechnic University，Jiaozuo 454000，China)

Through analyzing the annual volatility of cloud detection thresholds and integrating the band characteristics of HJ－1B，the authors proposed a new algorithm of cloud detection for HJ－1B image based on dynamical thresholding according to the spectral standard deviation anomaly.Using image registration，band math，linear regression and error analysis，the authors acquired the cloud abnormal regions which could be used to remove cloud pixels from the image.The results show that this approach can detect cloud pixels over different periods and in different scenes successfully，thus promoting the use of HJ－1B data and improving the precision of image classification.

standard deviation anomaly;dynamic threshold;cloud detection;HJ－1B

TP 751.1

A

1001－070X(2012)02－0012－07

2011－08－01;

2011－08－27

河南理工大學(xué)礦山空間信息技術(shù)國家測繪局重點實驗室開放基金資助項目(編號:KLM201112)和中國礦業(yè)大學(xué)江蘇省資源環(huán)境信息工程重點實驗室開放基金資助項目(編號:JS200907)共同資助。

10.6046/gtzyyg.2012.02.03

韓 杰(1987－)，男，碩士研究生，主要從事遙感圖像處理方面的研究。E－mail:rainbow871027@163.com。

(責(zé)任編輯:劉心季)