国产日韩欧美一区二区三区三州_亚洲少妇熟女av_久久久久亚洲av国产精品_波多野结衣网站一区二区_亚洲欧美色片在线91_国产亚洲精品精品国产优播av_日本一区二区三区波多野结衣 _久久国产av不卡

?

造紙廢水處理過程微小故障檢測方法研究

2017-09-15 08:56王齡松馬璞璠葉鳳英熊智新趙小燕劉鴻斌
中國造紙 2017年8期
關(guān)鍵詞:廢水處理矩陣變量

王齡松 馬璞璠 葉鳳英 熊智新 趙小燕 劉鴻斌,2,*

(1.南京林業(yè)大學(xué)江蘇省制漿造紙科學(xué)與技術(shù)重點(diǎn)實驗室,江蘇南京,210037;2.華南理工大學(xué)制漿造紙工程國家重點(diǎn)實驗室,廣東廣州,510640)

·廢水處理過程檢測·

造紙廢水處理過程微小故障檢測方法研究

王齡松1馬璞璠1葉鳳英1熊智新1趙小燕1劉鴻斌1,2,*

(1.南京林業(yè)大學(xué)江蘇省制漿造紙科學(xué)與技術(shù)重點(diǎn)實驗室,江蘇南京,210037;2.華南理工大學(xué)制漿造紙工程國家重點(diǎn)實驗室,廣東廣州,510640)

微小故障的檢測是過程監(jiān)測領(lǐng)域的一個重要研究方向。傳統(tǒng)的多元統(tǒng)計過程監(jiān)測方法無法對過程微小故障進(jìn)行有效監(jiān)控。本課題將多元累積和控制(CUSUM)方法及多元指數(shù)加權(quán)移動平均(EWMA)方法分別與主成分分析(PCA)相結(jié)合用于造紙廢水處理過程中微小故障的過程監(jiān)測。研究結(jié)果證實了多元累積和控制方法和多元指數(shù)加權(quán)移動平均方法的有效性。

造紙廢水處理過程;主成分分析;累積和控制;指數(shù)加權(quán)移動平均;故障檢測

(*E-mail: hbinjm@163.com)

隨著科技的發(fā)展,現(xiàn)代控制系統(tǒng)和設(shè)備的復(fù)雜性不斷增加、規(guī)模不斷擴(kuò)大,生產(chǎn)過程中的控制系統(tǒng)一旦發(fā)生故障將會造成巨大損失。在過程監(jiān)測領(lǐng)域,故障的定義是系統(tǒng)的某個或某些特征或者變量發(fā)生了超過允許范圍的偏差[1]。過程監(jiān)測和故障診斷技術(shù)可以提高系統(tǒng)的可靠性,降低事故的風(fēng)險,對于工業(yè)生產(chǎn)過程的正常運(yùn)行起到非常重要的作用[2]。

流程工業(yè)生產(chǎn)過程中的測量變量往往相互關(guān)聯(lián),傳統(tǒng)的單變量統(tǒng)計過程控制技術(shù)無法滿足工業(yè)上的監(jiān)控要求?;诙嘣队袄碚摰亩嘣y(tǒng)計過程監(jiān)控技術(shù)逐漸受到學(xué)術(shù)界和工業(yè)界的廣泛重視,并且在化工生產(chǎn)過程中得到了成功的應(yīng)用[3-5]。化工生產(chǎn)過程都存在著緩慢變化,比如設(shè)備的老化、催化劑的性能變化以及化學(xué)藥品濃度的變化等,這些緩慢變化的微小故障如果不能及時地檢測出來,輕則會導(dǎo)致產(chǎn)品質(zhì)量的下降,重則甚至可能導(dǎo)致生產(chǎn)事故的發(fā)生。

累積和控制(Cumulative sum,CUSUM)與指數(shù)加權(quán)移動平均(Exponent weighted moving average,EWMA)是兩種經(jīng)典的單變量統(tǒng)計過程控制方法,廣泛應(yīng)用于單變量過程系統(tǒng)中微小信號的故障檢測[6]。然而,單變量CUSUM和EWMA并沒有考慮生產(chǎn)過程中多種變量之間的相關(guān)性,對于復(fù)雜的工業(yè)生產(chǎn)過程來講,僅使用單變量監(jiān)控不僅操作麻煩,圖形眾多,而且會帶來頻繁的誤報和漏報問題[7]。為此,Chen等人[8]提出將CUSUM、EWMA與主成分分析(Principal component analysis,PCA)相結(jié)合來對微小故障進(jìn)行檢測。

針對造紙廢水處理系統(tǒng)的多變量特性,本課題首先將單變量CUSUM和單變量EWMA進(jìn)行多變量擴(kuò)展,形成多變量CUSUM(Multivariate CUSUM, MCUSUM)和多變量EWMA(Multivariate EWMA, MEWMA),然后將MCUSUM和MEWMA作為數(shù)據(jù)預(yù)處理方法,最后將MCUSUM和MEWMA與PCA方法分別結(jié)合,構(gòu)成MCUSUM-PCA和MEWMA-PCA方法。對造紙廢水處理過程控制系統(tǒng)加入偏移和漂移兩類故障信號,采用統(tǒng)計指標(biāo)對系統(tǒng)進(jìn)行過程監(jiān)測和故障診斷。

1 故障檢測方法簡介

1.1 多變量統(tǒng)計過程監(jiān)測

對于一組正常工況的數(shù)據(jù)Xm×n(式中,m為樣本個數(shù);n為變量個數(shù)),先將原始數(shù)據(jù)矩陣進(jìn)行標(biāo)準(zhǔn)化處理,即將矩陣各列轉(zhuǎn)化為均值為0,方差為1的向量。然后運(yùn)用PCA,將X分解成等式(1):

(1)

式中,T∈Rm×r、P∈Rn×r分別為主元得分矩陣和載荷矩陣;E∈Rm×n為預(yù)測殘差矩陣;r為主元個數(shù),其值可用交叉檢驗或觀察累積方差貢獻(xiàn)率來確定。

原始數(shù)據(jù)矩陣經(jīng)過降維處理后,變量之間相關(guān)性得到消除。通過對主元得分矩陣和預(yù)測殘差矩陣進(jìn)行變換,構(gòu)建出兩個統(tǒng)計量T2和Q,以檢測造紙廢水處理過程的微小故障。

對于測試數(shù)據(jù)xk=[xk,1,xk,2,…,xk,n]T,按照公式(2)計算統(tǒng)計量T2:

(2)

式中,Λ=diag(λ1,λ2,…,λr),λ為X的協(xié)方差矩陣的特征值。T2的統(tǒng)計控制限CL服從F分布,由式(3)計算得到:

(3)

式中,F(xiàn)r,m-r,α表示自由度為r和m-r的F分布的上側(cè)α分位數(shù)。當(dāng)統(tǒng)計量超出控制限的時候,視為過程失控,系統(tǒng)發(fā)出報警。

統(tǒng)計量Q,又稱平方預(yù)測誤差統(tǒng)計量,其計算如公式(4):

(4)

式中,I為對角線全為1的單位矩陣。

統(tǒng)計量Q的統(tǒng)計控制限CL按公式(5)計算:

(5)

1.2 多元累積和控制圖(MCUSUM)

工業(yè)生產(chǎn)過程中的緩慢變化和微小變化,如催化劑的緩慢失效,生產(chǎn)設(shè)備的逐年老化,以及測量設(shè)備的失靈等,會對生產(chǎn)過程造成不容忽視的不良影響。這些細(xì)微的變化如果不能被及時發(fā)現(xiàn)和處理,很有可能會影響到生產(chǎn)的進(jìn)行以及產(chǎn)品的質(zhì)量,甚至出現(xiàn)不可挽回的安全事故。對于微小故障的監(jiān)控,多元累積和控制圖(CUSUM控制圖)和指數(shù)加權(quán)移動平均控制圖(EWMA控制圖)是單變量統(tǒng)計過程監(jiān)控領(lǐng)域中的有效方法。該類方法只能對單變量過程進(jìn)行監(jiān)控,應(yīng)用范圍受到限制,無法滿足當(dāng)代復(fù)雜大型生產(chǎn)設(shè)備以及生產(chǎn)過程的復(fù)雜性要求。本課題采用如下多元累積和控制圖(MCUSUM)方法對造紙廢水處理系統(tǒng)進(jìn)行微小故障檢測。通過對多元數(shù)據(jù)中各個變量在d個時刻內(nèi)的數(shù)據(jù)進(jìn)行加和,一方面,可以放大數(shù)據(jù)中的微小變化,提高后續(xù)PCA的過程監(jiān)測能力;另一方面,也能夠降低異常數(shù)據(jù)對過程監(jiān)測統(tǒng)計量的影響,增強(qiáng)監(jiān)測系統(tǒng)的魯棒性。

假設(shè)從過程中采集到的數(shù)據(jù)為Xm×n,對該數(shù)據(jù)進(jìn)行處理,定義t時刻的累積和變量值為:

(6)

式中,d表示各樣本累計和的步數(shù),即從當(dāng)前時刻往前推d個時刻;xi表示第i個采樣;m為總的樣本個數(shù)。為保證數(shù)據(jù)時效性及減少實際應(yīng)用中的計算量,不對由yt構(gòu)成的數(shù)據(jù)集進(jìn)行歸一化處理,直接將t≤d時的數(shù)據(jù)舍棄,然后對其進(jìn)行傳統(tǒng)PCA過程監(jiān)測分析,從而實現(xiàn)與PCA的結(jié)合。

1.3 多元指數(shù)加權(quán)移動平均控制圖(MEWMA)

當(dāng)故障的幅值小于其所能夠檢測到的臨界故障幅值時,傳統(tǒng)的PCA故障檢測方法無法有效地檢測出該故障,而指數(shù)加權(quán)移動平均控制圖(EWMA控制圖)則適用于緩慢漂移和微小故障的檢測。將MEWMA與PCA相結(jié)合可以提高主元分析的故障檢測性能[8]。

作為單變量EWMA的擴(kuò)展,MEWMA定義見公式(7):

zk=ωxk+(1-ω)zk-1

(7)

式中,0<ω≤1且z0=0。

X∈Rm×n經(jīng)過MEWMA處理后變?yōu)閆∈Rm×n,Z的協(xié)方差矩陣為SZ[ω/(2-ω)]S,S為X的協(xié)方差矩陣。對經(jīng)過MEWMA預(yù)處理過的數(shù)據(jù)矩陣進(jìn)行傳統(tǒng)PCA過程監(jiān)測分析,從而實現(xiàn)與PCA的結(jié)合。

統(tǒng)計量T2采用公式(8)來計算[8]:

(8)

統(tǒng)計量Q采用公式(9)來計算[8]:

(9)

2 造紙廢水處理過程故障檢測

2.1 造紙廢水處理過程數(shù)據(jù)

本課題所用的數(shù)據(jù)取自于廣東東莞的一家造紙廠好氧段廢水的監(jiān)測數(shù)據(jù)。該數(shù)據(jù)含有170組測量值,8個過程及水質(zhì)變量,如圖1所示。圖1中左邊縱坐標(biāo)對應(yīng)于進(jìn)水化學(xué)需氧量(CODinf)、出水化學(xué)需氧量(CODeff)、進(jìn)水懸浮固形物(SSinf)、出水懸浮固形物(SSeff)值;右邊縱坐標(biāo)表示流量Q(104m3/d)、pH值、溫度T(℃)和溶解氧濃度DO(mg/L)[9]。

圖1 造紙廢水處理過程數(shù)據(jù)

2.2 故障數(shù)據(jù)的構(gòu)造

從故障診斷的角度分析,故障可以歸納為三類:傳感器故障、執(zhí)行器故障以及過程參數(shù)故障[1]。本課題以傳感器故障檢測為例來驗證所提出微小故障檢測方法的有效性。常見的傳感器故障類型有偏移故障、漂移故障、完全故障以及精度下降故障[10]。本課題對傳感器變量CODinf和pH值分別構(gòu)建數(shù)據(jù)故障。對傳感器變量CODinf加入均值的20%數(shù)據(jù)故障,對傳感器變量pH值加入0.05系數(shù)的漂移故障,得到偏移和漂移兩種故障類型(表1,圖2)。在MATLAB軟件中對170個造紙廢水樣本數(shù)據(jù)進(jìn)行分析處理,前100個樣本作為訓(xùn)練集,后70個數(shù)據(jù)作為測試集,測試集中取一部分?jǐn)?shù)據(jù)加入故障信息(本課題從第121個采樣點(diǎn)開始加入故障信息)。

表1 偏移和漂移故障引入公式

圖2 偏移故障和漂移故障

2.3 MCUSUM和MEWMA中關(guān)鍵參數(shù)的確定

2.3.1 MCUSUM預(yù)處理中d值的確定

在MCUSUM算法中,參數(shù)d對故障檢測結(jié)果有最直接的影響。通過對比偏移和漂移兩組故障檢測結(jié)果發(fā)現(xiàn),MCUSUM-PCA的監(jiān)控效果隨著d值的增大而增強(qiáng)。較大的d值會使監(jiān)控過程更加穩(wěn)定,不會使控制圖因為某一樣本值的突變而發(fā)生劇烈變化。另一方面,d值的增大也意味著需要積累的數(shù)據(jù)量步數(shù)增大,故障檢測的滯后也隨之變大。綜合考慮,本文d值選定為12。

圖3 偏移故障的PCA檢測結(jié)果

圖4 偏移故障的MCUSUM-PCA檢測結(jié)果(d=12)

圖5 偏移故障的MEWMA-PCA檢測結(jié)果(ω=0.2)

2.3.2 MEWMA預(yù)處理中λ值的確定

在MEWMA算法中,權(quán)值參數(shù)ω對MEWMA-PCA統(tǒng)計量期望值的變化影響很大。當(dāng)ω取較大值時(接近1),統(tǒng)計量能夠及時地反映過程中的突然變化,從而對突變故障具有良好的檢測性能,然而對于在檢測過程中出現(xiàn)的緩變故障,卻會出現(xiàn)漏報情況;而當(dāng)ω取較小值時(接近0),MEWMA能夠起到很好地平滑作用,能夠過濾過程中的噪聲,但同時也會引起一定的檢測滯后。因此,ω參數(shù)的選取要綜合考慮噪聲和檢測滯后兩個方面,二者權(quán)衡,以求達(dá)到最好的檢測效果。經(jīng)過對造紙廢水處理數(shù)據(jù)的多次模擬試驗,最終確定ω=0.2。

2.4 故障檢測結(jié)果與分析

2.4.1 偏移故障檢測結(jié)果

使用PCA、MCUSUM-PCA和MEWMA-PCA分別對偏移故障數(shù)據(jù)進(jìn)行檢測,結(jié)果如圖3~圖5所示。對于偏移故障,計算結(jié)果表明,MCUSUM-PCA與MEWMA-PCA兩種方法都能夠在故障(21~70樣本點(diǎn))出現(xiàn)后的3~4個樣本點(diǎn)內(nèi)檢測出故障,監(jiān)測過程超過控制線而發(fā)出故障報警(圖4和圖5)。相比較而言,PCA方法對偏移故障的檢測效果較差(圖3)。

圖6 漂移故障的PCA檢測結(jié)果

圖7 漂移故障的MCUSUM-PCA檢測結(jié)果(d=12)

圖8 漂移故障的MEWMA-PCA檢測結(jié)果(ω=0.2)

PCA、MCUSUM-PCA和MEWMA-PCA的故障正確檢測率如表2所示。MCUSUM-PCA和MEWMA-PCA的Q統(tǒng)計指標(biāo)的正確檢測率分別為0.96和0.94,MEWMA-PCA檢測出的異常數(shù)據(jù)比MCUSUM-PCA滯后一個樣本點(diǎn)。從圖4和圖5中可以看出,兩種方法對出現(xiàn)故障之后的持續(xù)異常狀態(tài)都取得了有效的監(jiān)控。然而,PCA方法的Q統(tǒng)計指標(biāo)的正確檢測率只有0.16,這是因為偏移故障的幅度變化比較微小,常規(guī)的PCA監(jiān)測方法無法對故障造成的微小變化進(jìn)行積累,導(dǎo)致無法對故障進(jìn)行準(zhǔn)確的檢測。

表2 PCA、MCUSUM-PCA和MEWMA-PCA對偏移故障

2.4.2 漂移故障檢測結(jié)果

使用PCA、MCUSUM-PCA和MEWMA-PCA分別對漂移故障數(shù)據(jù)進(jìn)行檢測,結(jié)果如圖6~圖8所示。對于漂移故障,MCUSUM-PCA與MEWMA-PCA兩種方法都能夠在故障(21~70樣本點(diǎn))出現(xiàn)后的10~15個樣本點(diǎn)內(nèi)檢測出故障,監(jiān)測過程超過控制線而發(fā)出故障報警(圖7和圖8)。相比較而言,PCA方法對偏移故障的檢測效果較差(圖6)。

PCA、MCUSUM-PCA和MEWMA-PCA的故障正確檢測率如表2所示。MCUSUM-PCA方法的Q統(tǒng)計指標(biāo)的正確檢測率為0.76,MEWMA-PCA方法的Q統(tǒng)計指標(biāo)的正確檢測率為0.72。兩種方法對漂移故障的正確檢測率均小于其對偏移故障的正確檢測率,原因在于漂移故障是一個緩慢變化的微小故障,需要一定的時間才能表現(xiàn)出故障的存在,在故障發(fā)生的初始時刻,統(tǒng)計指標(biāo)并沒有超出控制限,監(jiān)測系統(tǒng)無法檢測到此時的故障,直到故障信號積累到了一定程度,監(jiān)測系統(tǒng)才能監(jiān)測到故障發(fā)出報警。之后,兩種方法對出現(xiàn)故障之后的持續(xù)異常狀態(tài)都取得了有效的監(jiān)控。而PCA方法的Q統(tǒng)計指標(biāo)的正確檢測率為0,這是因為漂移故障在用PCA方法檢測時,殘差空間并不能有效地表達(dá)過程的變化趨勢,Q統(tǒng)計指標(biāo)全部位于統(tǒng)計限之下,故沒有進(jìn)行故障報警。而對于T2統(tǒng)計指標(biāo),PCA方法對漂移故障的正確檢測率也低于MCUSUM-PCA和MEWMA-PCA的T2和Q統(tǒng)計指標(biāo)。

3 結(jié) 論

本課題采用多變量累積和控制(MCUSUM)、多變量指數(shù)加權(quán)移動平均(MEWMA)分別與PAC結(jié)合構(gòu)成的MCUSUM-PCA、MEWMA-PCA及PCA三種多元統(tǒng)計方法對造紙廢水處理過程進(jìn)行監(jiān)測研究。對于偏移和漂移的微小故障,MCUSUM-PCA和MEWMA-PCA的故障正確檢測率均高于傳統(tǒng)PCA方法的故障正確檢測率。在故障檢測的靈敏度方面,MCUSUM-PCA和MEWMA-PCA的表現(xiàn)也優(yōu)于傳統(tǒng)PCA監(jiān)測方法。MCUSUM-PCA和MEWMA-PCA都能在故障出現(xiàn)2~3個樣本點(diǎn)內(nèi)就可檢測出偏移故障,在10~15個樣本點(diǎn)內(nèi)就可檢測出漂移故障;而傳統(tǒng)PCA方法卻因?qū)ξ⑿」收系臋z測效果不明顯而出現(xiàn)大量的漏檢現(xiàn)象。

[1] HUANG Dao-ping, QIU Yu, LIU Yi-qi, et al. Review of data-driven fault diagnosis and prognosis for wastewater treatment[J]. Journal of South China University of Technology (Natural Science Edition), 2015, 43(3): 111. 黃道平, 邱 禹, 劉乙奇, 等. 面向污水處理的數(shù)據(jù)驅(qū)動故障診斷及預(yù)測方法綜述[J]. 華南理工大學(xué)學(xué)報 (自然科學(xué)版), 2015, 43 (3): 111.

[2] ZHOU Dong-hua, HU Yan-yan. Fault diagnosis techniques for dynamic systems[J]. Acta Automatica Sinica, 2009, 35(6): 748. 周東華, 胡艷艷. 動態(tài)系統(tǒng)的故障診斷技術(shù)[J]. 自動化學(xué)報, 2009, 35(6): 748.

[3] Liu Tian-long, Shen Wen-hao. A review of applications of fault diagnostic expertsystem in wastewater treatment[J]. Paper Science & Technology, 2011, 30(2): 75. 劉天龍, 沈文浩. 污水處理過程中故障診斷專家系統(tǒng)的應(yīng)用[J]. 造紙科學(xué)與技術(shù), 2011, 30(2): 75.

[4] QIN S. Joe. Survey on data-driven industrial process monitoring and diagnosis[J]. Annual Reviews in Control, 2012, 36(2): 220.

[5] JI Hong-quan, HE Xiao, ZHOU Dong-hua. Fault detection techniques based on multivariates statistical analysis[J]. Journal of Shanghai Jiaotong University, 2015, 49 (6): 842. 紀(jì)洪泉, 何 瀟, 周東華. 基于多元統(tǒng)計分析的故障檢測方法[J]. 上海交通大學(xué)學(xué)報, 2015, 49 (6): 842.

[6] LI Juan, ZHOU Dong-hua, SI Xiao-sheng, et al. Review of incipient fault diagnosis methods[J]. Control Theory & Applications, 2012, 29(12): 1517. 李 娟, 周東華, 司小勝, 等. 微小故障診斷方法綜述[J]. 控制理論與應(yīng)用, 2012, 29(12): 1517.

[7] GE Zhi-qiang, YANG Chun-jie, SONG Zhi-huan. Research and application of small shifts detection method based on MEWMA-PCA[J]. Information and Control, 2007, 36(5): 650. 葛志強(qiáng), 楊春節(jié), 宋執(zhí)環(huán). 基于MEWMA-PCA的微小故障檢測方法研究及其應(yīng)用[J]. 信息與控制, 2007, 36(5): 650.

[8] CHEN Junghui, LIAO Chien-Mao, LIN Franz Ren Jen, et al. Principle component analysis based control charts with memory effect for process monitoring[J]. Industrial & Engineering Chemistry Research, 2001, 40(6): 1516.

[9] YANG Hao, MO Wei-lin, XIONG Zhi-xin, et al. Soft Sensor Modeling of Papermaking Effluent Treatment Processes Using RPLS[J]. China Pulp & Paper, 2016, 35(10): 31. 楊 浩, 莫衛(wèi)林, 熊智新, 等. 基于RPLS的造紙廢水處理過程軟測量建模[J]. 中國造紙, 2016, 35(10): 31.

(責(zé)任編輯:馬 忻)

Incipient Fault Detection in Papermaking Wastewater Treatment Processes

WANG Ling-song1MA Pu-fan1YE Feng-ying1XIONG Zhi-xin1ZHAO Xiao-yan1LIU Hong-bin1,2,*

(1.JiangsuProvincialKeyLabofPulpandPaperScienceandTechnology,NanjingForestryUniversity,Nanjing,JiangsuProvince, 210037; 2.StateKeyLabofPulpandPaperEngineering,SouthChinaUniversityofTechnology,Guangzhou,GuangdongProvince, 510640)

Incipient fault detection has been an important research topic in the field of process monitoring. However, traditional multivariate statistical process monitoring methods fail to detect incipient faults. In this paper, two methods were developed for the incipient fault detection of a papermaking wastewater treatment process: multivariate cumulative sum combined with principal component analysis (PCA) and multivariate exponent weighted moving average combined with PCA. The results proved the effectiveness of the proposed fault monitoring methods.

papermaking wastewater treatment process; principal component analysis; cumulative sum (CUSUM); exponent weighted moving average (EWMA); fault detection

王齡松先生,在讀碩士研究生;研究方向:制漿造紙過程控制。

2017- 04-25(修改稿)

制漿造紙工程國家重點(diǎn)實驗室開放基金資助項目(201610);南京林業(yè)大學(xué)高層次人才科研啟動基金(163105996);江蘇省制漿造紙科學(xué)與技術(shù)重點(diǎn)實驗室開放基金項目(201530)。

TS79

A

10.11980/j.issn.0254- 508X.2017.08.004

*通信作者:劉鴻斌,副教授;研究方向:制漿造紙過程建模、監(jiān)測和控制。

猜你喜歡
廢水處理矩陣變量
抓住不變量解題
也談分離變量
一種O-A-A-O工藝在焦化廢水處理中的應(yīng)用
汽車涂裝廢水處理技術(shù)的研究
廢水處理技術(shù)與水環(huán)境保護(hù)的研究
上旋流厭氧反應(yīng)器在造紙廢水處理中的應(yīng)用
初等行變換與初等列變換并用求逆矩陣
矩陣
矩陣
矩陣
澄迈县| 德昌县| 新竹县| 阿坝县| 尼木县| 沁阳市| 满城县| 理塘县| 太仓市| 弥渡县| 许昌县| 雷波县| 阳泉市| 佳木斯市| 台东市| 衡水市| 巩留县| 祁连县| 皋兰县| 习水县| 重庆市| 常熟市| 巴里| 遂宁市| 江达县| 祥云县| 肇东市| 留坝县| 鹤岗市| 徐闻县| 福海县| 凤台县| 肇庆市| 普洱| 定兴县| 宁国市| 呈贡县| 明水县| 静安区| 吴桥县| 弥渡县|