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

?

EEMD 修正爆破加速度零漂信號中的最優(yōu)白噪聲系數(shù)*

2019-09-25 03:24:06王志亮陳貴豪黃佑鵬
爆炸與沖擊 2019年8期
關(guān)鍵詞:噪聲系數(shù)頻段修正

王志亮,陳貴豪,黃佑鵬

(合肥工業(yè)大學(xué)土木與水利工程學(xué)院,安徽 合肥 230009)

在巖石鉆孔爆破過程中,爆源近區(qū)巖石的沖擊振動信號具有持時短、頻帶寬的特點(diǎn)[1]。目前廣泛使用的壓電式加速度傳感器在采集此類信號時,大多會出現(xiàn)“零漂”現(xiàn)象[2-6]。造成壓電傳感器出現(xiàn)零漂誤差的原因主要有[2-3]:(1)受到高頻沖擊波作用時,傳感器由于自身諧振反應(yīng)導(dǎo)致內(nèi)部壓電晶體磁疇變化,這種變化反映到輸出信號上,出現(xiàn)零漂現(xiàn)象;(2)高沖擊作用下,傳感器壓電元件受到的作用力超過其在裝配時受到的預(yù)壓力,進(jìn)而造成壓電晶體與慣性質(zhì)量塊之間出現(xiàn)微量運(yùn)動,造成零漂現(xiàn)象;(3)傳感器的質(zhì)量-彈簧系統(tǒng)受到高沖擊作用時的動剛度不足會導(dǎo)致其幅值線性度變差,可能造成零漂;(4)受限于傳感器壓電材料和制造工藝,信號放大器電路易過早進(jìn)入飽和堵塞狀態(tài),RC 電路放電時間延長到?jīng)_擊結(jié)束后,表現(xiàn)為零漂;(5)傳感器基座與被測物體表面連接剛度不足造成零漂誤差。由于爆破過程中巖石結(jié)構(gòu)會產(chǎn)生不可逆的損傷,因而每次爆破信號的采集是單次性的。所以,如果爆破振動信號采集過程中出現(xiàn)零漂現(xiàn)象,采用合理方法對誤差信號進(jìn)行修正非常有必要[4-6]。

近年來,不少學(xué)者對于零漂信號中產(chǎn)生的趨勢項(xiàng)去除方法進(jìn)行研究:龍源等[6]基于經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, EMD)、小波法、最小二乘法分別研究了爆破震動測試信號中趨勢項(xiàng)的去除算法,并對三種方法去除趨勢項(xiàng)效果及重構(gòu)信號的時頻特征進(jìn)行對比分析,指出多項(xiàng)式趨勢項(xiàng)中的多項(xiàng)式階數(shù)選取對于信號的修正效果有較大影響;王燕等[7]采用EMD 法和最小二乘法相結(jié)合,以侵徹結(jié)束后的零漂趨勢作為侵徹過程中的零漂趨勢,對侵徹不同介質(zhì)的實(shí)測過載數(shù)據(jù)進(jìn)行去除零漂的處理,通過對修正后的過載數(shù)據(jù)進(jìn)行一次積分和二次積分,得到的侵徹速度和位移與實(shí)測結(jié)果吻合較好;王濟(jì)[8]研究發(fā)現(xiàn)滑動平均法應(yīng)用較為簡便,但其本質(zhì)僅是為了曲線基線回歸零點(diǎn),在處理信號時有一定的盲目性;謝全民等[9]認(rèn)為小波法需要針對振動信號特點(diǎn)事先選取合適的小波基函數(shù)和分解層數(shù);吳祖堂等[10]通過理論分析和實(shí)驗(yàn)驗(yàn)證得出沖擊結(jié)束后的壓電式傳感器信號呈指數(shù)形式衰減的結(jié)論,但是信號的零漂往往在沖擊信號達(dá)到峰值前就已經(jīng)形成,沖擊結(jié)束后信號的指數(shù)衰減規(guī)律和沖擊過程中的趨勢項(xiàng)類型是否一致,目前尚無定論;Huang 等人[11]揭示EMD 方法及其改進(jìn)算法—集合經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition, EEMD)方法的基函數(shù)都來自于其信號本身,對于信號的先驗(yàn)性信息要求較低,在處理爆破沖擊這類非平穩(wěn)信號時優(yōu)勢明顯。實(shí)際上,EEMD 方法在EMD 方法基礎(chǔ)上,通過添加在頻譜上均勻分布的白噪聲有效地解決了EMD 方法分解得到的固有模態(tài)分量(intrinsic mode function, IMF) 中存在的模態(tài)混疊和端點(diǎn)震蕩問題[12],故認(rèn)為EEMD 方法處理零漂信號的效果優(yōu)于EMD 方法。但EEMD 方法中需要給出添加白噪聲的系數(shù),目前沒有確定的量化指標(biāo)去評價零漂信號的修正效果,從而無法確定白噪聲系數(shù)的合理取值范圍。

基于此,本文先給出不同白噪聲系數(shù)下EEMD 結(jié)合高低頻處理方法修正巖石爆破中加速度零漂信號的過程;再輔以沖擊響應(yīng)譜分析,提出表征修正前后頻域平均偏差幅度的修正指數(shù);在修正指數(shù)分析的基礎(chǔ)上確定EEMD 修正爆破零漂信號中白噪聲系數(shù)的最優(yōu)取值范圍,再用加速度信號積分后的速度信號對比討論最優(yōu)噪聲系數(shù)范圍的合理性,最后對比EEMD 法與滑動平均法,討論二者的優(yōu)點(diǎn)與不足。

1 算法原理

EEMD 是針對EMD 算法中所存在的不足進(jìn)行改進(jìn)的算法,EMD 算法實(shí)現(xiàn)過程如圖1 所示[11]。EEMD 主要利用高斯白噪聲在頻譜上均勻分布的特性來彌補(bǔ)原始信號中一些缺失的尺度,使得不同時間尺度的信號自動分解到其相適應(yīng)的參考尺度上去,從而有效地抑制EMD 分解信號時存在的模態(tài)混疊問題。由于添加的白噪聲總體均值為零,所以將多次添加白噪聲后所求得的IMF 分量的集合平均值作為最終結(jié)果,以抵消白噪聲對于真實(shí)IMF 分量的影響。

圖 1 EMD 算法流程圖Fig. 1 Flow chart of EMD algorithm

EEMD 算法的實(shí)現(xiàn)步驟如下[12]。

(1) 將隨機(jī)高斯白噪聲信號nj(t)添加到原始信號x(t):

式中:xj(t)是第j 次加噪信號,k、m 分別是添加白噪聲的系數(shù)和總次數(shù),σ 是原始信號的標(biāo)準(zhǔn)差。

(2) 用EMD 方法將xj(t)分解成一系列IMF 分量ci,j和余項(xiàng)δj:

式中:ci,j和δj分別為第j 次加入白噪聲后信號分解得到的第i 個IMF 分量和余項(xiàng),gj為第j 次加入白噪聲后信號的分解尺度。

(3) 如果j<m,則重復(fù)步驟(1)和(2)。

(4) 獲取I=min(g1, g2,···, gm),計算I 相應(yīng)的分解的IMF 的集合平均值作為最終結(jié)果:

2 信號處理

2.1 信號介紹

本文所處理的零漂信號來源于花崗巖爆破試驗(yàn),其數(shù)據(jù)采集過程和時程曲線分別如圖2 和圖3 所示。爆源為1.5 g 太安炸藥,爆速約為6 000 m·s?1。圖3 中加速度信號A 測點(diǎn)的爆心距約為60 cm。信號放大器的采樣頻率為5×105Hz,加速度傳感器的頻響范圍是40~4×104Hz,設(shè)置低通3×104Hz。由圖3可以看出峰值加速度超過了4×104m·s?2,同時加速度信號在沖擊結(jié)束后沒有迅速歸零,而是持續(xù)一段時間后才緩慢歸零,這是高速沖擊環(huán)境下的典型零漂信號。

圖 2 試驗(yàn)與數(shù)據(jù)采集Fig. 2 Test and data collection

圖 3 爆破加速度信號AFig. 3 Blast acceleration signal A

2.2 白噪聲系數(shù)范圍

Wu 等[12]在分析EEMD 理論時,建議k 近似取0.20;Jiang 等[13]用EEMD 方法對振動信號分析時指出,如果信號主頻段由高頻分量組成,則k 的取值范圍建議為0~0.20。

巖石的爆破試驗(yàn)結(jié)果表明[14-15],巖石質(zhì)點(diǎn)的振動主頻率與炸藥裝藥量成反比關(guān)系,與炸藥爆速成正比關(guān)系;裝藥深度對振動主頻率也有很大影響,深孔爆破產(chǎn)生的地震波主頻較低,淺孔爆破產(chǎn)生的地震波主頻較高。本文爆破試驗(yàn)的裝藥量小、深度淺、炸藥爆速大且傳播距離短,故認(rèn)為其主頻段主要由高頻分量組成。因此,白噪聲系數(shù)k 選在0~0.20 范圍內(nèi)討論。

2.3 修正過程

2.3.1 EEMD 分解

為分析白噪聲系數(shù)對于EEMD 方法修正零漂信號的影響,在0~0.20 范圍內(nèi)選取0.05、0.10、0.15 和0.20 四個k 值對應(yīng)的白噪聲添加到原始信號中進(jìn)行EEMD 分解,添加次數(shù)均為100 次。圖4 給出k=0.10 時分解得到的IMF 分量和余項(xiàng)δ。

2.3.2 IMF 分量頻譜分析

將分解得到的IMF 分量分別進(jìn)行FFT 變換,可得各IMF 分量的頻譜信息。不同白噪聲系數(shù)下分解得到的各IMF 分量的主頻如表1 所示。限于篇幅,僅給出部分IMF 分量的頻譜圖,如圖5 所示。

從表1 可以看出,白噪聲系數(shù)對前兩階IMF 分量的主頻影響很小,前兩階IMF 分量主頻均在18 800 Hz 左右。結(jié)合圖5(a) 可知,白噪聲系數(shù)主要影響其主頻幅值,系數(shù)越大,主頻幅值越小;從IMF2 至IMF11 分量,隨著分解階數(shù)的增加,IMF 分量主頻迅速衰減,且不同白系數(shù)噪聲下分解的同階IMF 分量主頻與白噪聲系數(shù)呈現(xiàn)正相關(guān)關(guān)系;IMF11 之后的分量的主頻降至低頻段,且基本保持不變。同時,從圖5 可看出,相同白噪聲系數(shù)下分解得到的IMF 分量隨著分解階數(shù)增加,其頻段降低,頻帶漸窄。

圖 4 EEMD 分解信號(k=0.10)Fig. 4 Signals decomposed by EEMD method (k=0.10)

2.3.3 IMF 分量處理

(1)低頻處理。由表1 知,IMF11 至IMF13 分量和余項(xiàng)δ 的主頻集中在40 Hz 以下,低于試驗(yàn)中加速度傳感器的頻響下限 (40 Hz),故認(rèn)為其可信度很低。因此,將IMF11 至IMF13 分量和余項(xiàng)δ 消除。

(2)高頻處理。由于原始信號就含有大量高頻噪聲,故對不同白噪聲系數(shù)下的IMF1 至IMF10 分量分別進(jìn)行小波閾值去噪處理,閾值由heursure 函數(shù)自適應(yīng)獲取,去噪方式采用軟閾值處理[16-17]。

表 1 不同白噪聲系數(shù)下IMF 分量主頻Table 1 Dominant frequencies of IMF components under different white noise coefficients

圖 5 部分IMF 分量頻譜圖Fig. 5 Spectrograms of partial IMF components

分別對處理過的分量IMF1 至IMF10 進(jìn)行重構(gòu),即可得到不同白噪聲系數(shù)下EEMD 修正后的爆破沖擊信號。信號A 修正后的四個加速度信號都很好地消除了零漂趨勢,但由于彼此在加速度時程圖上的波形差異很小導(dǎo)致難以區(qū)分,故僅給出白噪聲系數(shù)為0.10 時的修正信號,如圖6 所示。

3 修正信號分析

3.1 頻域分析

圖 6 EEMD 修正前后的信號Fig. 6 Signals before and after EEMD correction

圖7 給出了修正前后信號的沖擊響應(yīng)譜。從沖擊響應(yīng)譜上來看,四個修正的信號與原始信號A 的頻域“差異”主要集中在4 000 Hz 以下。在4 000 Hz 以上的頻段,修正前后的頻域“差異”不明顯,很好的保留了原始信號的峰值和變化規(guī)律。需要指出,這里的“差異”是指處理后的信號與原始信號在沖擊響應(yīng)譜幅值上的相對偏差度,而非絕對偏差值。由于40 Hz 以下頻段已經(jīng)超出加速度傳感器的頻率響應(yīng)范圍,這部分原始信號認(rèn)為是不可信的,采用的修正方法應(yīng)將這部分的信號盡可能削弱,同時應(yīng)該盡量保留中高頻段可信信號的頻譜特征。因此,為兼顧防止信號中高頻段過度修正和盡量削弱低頻失信段的目的,需要定義一個具體指標(biāo)來評價信號的修正效果。

圖 7 不同白噪聲系數(shù)下EEMD 修正信號的沖擊響應(yīng)譜Fig. 7 Shock response spectra of signals modified by EEMD under different white noise coefficients

3.1.1 修正指數(shù)定義

設(shè)所選頻段的頻率序列為{f1, f2,···, fl},對應(yīng)的沖擊響應(yīng)譜序列為{P1, P2, ···, Pl},其中f 為頻率,P 為每個頻率對應(yīng)的沖擊響應(yīng)譜幅值,l 為頻段的頻率點(diǎn)總數(shù)。則反映信號修正前后沖擊響應(yīng)譜平均偏差幅度的修正指數(shù)定義如下:

3.1.2 指數(shù)應(yīng)用

為了評價不同白噪聲系數(shù)下信號在不同頻段上的修正效果,在(0, 0.20]中以0.001 為間隔取200 個k 值進(jìn)行信號處理。原始信號A 在這些白噪聲系數(shù)下處理后,可得200 個修正信號。然后,將這些修正信號與原始信號的沖擊響應(yīng)譜代入式(5)計算,即可得到不同白系數(shù)噪聲下相應(yīng)頻段的修正指數(shù),如圖8 所示。在給出的五個頻段上,修正指數(shù)與白噪聲系數(shù)都呈現(xiàn)近似冪指數(shù)的關(guān)系,可假定修正指數(shù)與白噪聲系數(shù)的函數(shù)關(guān)系為:

式中:a0、a1、和b 為常數(shù)。

結(jié)合圖8 與表2,對不同的噪聲系數(shù)k 取值范圍作如下討論:

(1) 0<k<k1。在此區(qū)間內(nèi),失信頻段上的修正指數(shù)隨白噪聲系數(shù)變化的幅度很大,擬合函數(shù)最大值小于相應(yīng)的修正指數(shù)上限值的88%。故認(rèn)為在此區(qū)間失信頻段修正效果不夠穩(wěn)定,且對信號失信頻段的削弱效果不足,不滿足處理方法對于較好去除其低頻趨勢項(xiàng)的要求。

(2) k1≤k≤k2。在此區(qū)間內(nèi),失信頻段上的修正指數(shù)隨白噪聲系數(shù)變化的幅度相比于前一區(qū)間減小很多。修正指數(shù)平均值達(dá)到80.63%,比失信頻段的修正指數(shù)下限8.86%提高了810%,已經(jīng)達(dá)到修正指數(shù)上限值的92.98%。雖然在置信頻段及其三個分段上,修正指數(shù)會隨著噪聲系數(shù)取值區(qū)間的上升而增大,可能導(dǎo)致在這些頻段上出現(xiàn)信號修正過度的問題,但在這四個頻段上前后兩個區(qū)間的修正指數(shù)平均值分別只相差0.69%、3.73%、0.30%和0.13%,修正指數(shù)上升的幅度不大。故認(rèn)為這個區(qū)間內(nèi)修正后信號較好地去除了低頻信號,且在置信頻段尤其是高頻段的頻域失真很小,是可以接受的。

(3) k2<k ≤0.20。這一區(qū)間與[k1, k2]區(qū)間相比,噪聲系數(shù)的增大會使置信頻段尤其是高頻段的頻域失真增大,但與此同時信號失信頻段的削弱效果不能得到有效提升,白噪聲系數(shù)增大的意義很小。

綜上分析,可以確定信號A 對應(yīng)的最優(yōu)的白噪聲系數(shù)k 取值區(qū)間為[0.256,0.066]。

圖 8 不同白系數(shù)噪聲下不同頻段上的沖擊響應(yīng)譜修正指數(shù)Fig. 8 Correction indices of shock response spectra on different frequency bands under different white noise coefficients

表 2 不同頻段修正指數(shù)范圍Table 2 Ranges of correction index on different frequency bands

3.1.3 修正指數(shù)驗(yàn)證

通過圖9 看出,隨著噪聲系數(shù)的提高,EEMD 方法修正后信號與原始信號B 在沖擊響應(yīng)譜上的差異呈上升趨勢,修正指數(shù)與噪聲系數(shù)依然呈現(xiàn)如式(5) 的冪指數(shù)函數(shù)關(guān)系。但是當(dāng)噪聲系數(shù)k 上升到0.10 之后,EEMD 方法對失信頻段0~40 Hz 上信號的頻域削弱幾乎不再增加,但是對于置信頻段上的三個分段上信號的頻域削弱依然在增加,這也就意味著在(0.10, 0.20)范圍內(nèi),增大噪聲系數(shù)不能進(jìn)一步削弱失信頻段信號,反而會導(dǎo)致修正后信號在中高頻失真加大。如采用3.1.2 節(jié)中的方法,同樣可以確定一個大致的最優(yōu)噪聲系數(shù)取值范圍,即為[0.098,0.114]。

圖 9 不同白系數(shù)噪聲下不同頻段上的沖擊響應(yīng)譜修正指數(shù)Fig. 9 Correction indices of shock response spectrum on different frequency bands under different white noise coefficients

3.2 時域分析

從圖3、圖6 和圖10 中可以看出,僅從加速度零漂信號本身來說,EEMD 處理零漂信號A 和B 后得到的信號均很好地去除了原始信號中的零漂趨勢項(xiàng),其波形在爆破沖擊結(jié)束后都很快回歸基線,保留了原始爆破信號的衰減特征。

事實(shí)上,零漂偏差在加速度信號上表現(xiàn)得不如加速度時域積分后得到的速度信號顯著,這是因?yàn)椴粌H加速度信號基線漂移會導(dǎo)致積分后速度信號基線漂移嚴(yán)重,而且原始信號上下振幅的不對稱也會加劇積分曲線漂移。因此,對于加速度信號A 及其在不同白噪聲系數(shù)下處理得到的修正信號進(jìn)行積分,得到的速度信號如圖11 所示。

圖 10 爆破加速度信號B 及其修正信號Fig. 10 Blast acceleration signal B and its corrected signal

圖 11 EEMD 修正前后速度曲線Fig. 11 Velocity curves before and after correcting by EEMD

圖11 表明EEMD 方法對于速度信號的零漂失真有明顯改善:隨著白噪聲系數(shù)的提高,修正后的速度曲線逐漸靠近零點(diǎn)。當(dāng)k 值大于0.60 時,EEMD 對于速度信號的零漂修正效果幾乎沒有變化,這也與圖7(a)中失信頻段的修正指數(shù)分布圖相契合。同時,這也表明通過修正指數(shù)分析確定信號A 的最優(yōu)白噪聲系數(shù)范圍[0.056,0.066]是合理的:將白噪聲系數(shù)范圍上調(diào)已經(jīng)不能明顯改善速度曲線零漂趨勢;而將白噪聲系數(shù)范圍下調(diào),速度曲線明顯隨之下滑,零漂趨勢加重。但是,圖11也表明了EEMD 修正后的速度曲線無法完全消除零漂趨勢,優(yōu)化的噪聲系數(shù)也很難讓速度曲線基線完全回歸零點(diǎn)。

3.3 進(jìn)一步討論

以上的分析表明,EEMD 對于積分后的速度信號的零漂趨勢去除得不夠完善,存在一定的局限性。僅是考慮讓積分后速度信號基線完全恢復(fù)到零點(diǎn)上,滑動平均法是一種簡單有效的方法。影響滑動平均法去除信號零漂趨勢的因素有兩個[8]:平滑點(diǎn)數(shù)和滑動階次。對加速度信號A 進(jìn)行平滑點(diǎn)數(shù)為20 和滑動階次為10 的滑動平均法處理,處理后的信號積分得到速度曲線如圖12 所示。圖11 和圖12 對比顯示,滑動平均法得到的速度曲線是優(yōu)于EEMD 的,更符合爆破速度曲線的振動特征。

圖 12 滑動平均法修正后速度曲線Fig. 12 Velocity curve modified by sliding average method

通過圖13 給出的原始信號以及兩種方法處理后信號的沖擊響應(yīng)譜,發(fā)現(xiàn)滑動平均法相對于EEMD 方法在各個頻段上頻域損傷均更大。而沖擊響應(yīng)譜的低頻部分,往往反映了巖石自身大的位移或者應(yīng)變,其主要是由于零漂誤差還是巖石自身的動力響應(yīng)造成,還難以界定。所以,盡管滑動平均法很好地去除加速度信號積分后的速度信號的零漂趨勢,但這也會使得修正信號的頻譜失真加重以及巖石自身的動力響應(yīng)削弱。

圖 13 兩種方法處理后信號與原始信號的沖擊響應(yīng)譜Fig. 13 Shock response spectra of original signals and processed signals by the two methods

4 結(jié) 論

本文結(jié)合花崗巖的室外小型爆破試驗(yàn),對EEMD 修正爆破零漂信號中最優(yōu)白噪聲系數(shù)范圍進(jìn)行了探析,得出主要結(jié)論如下:(1)利用EEMD 方法并結(jié)合高低頻處理能夠較好地去除爆破加速度信號的零漂趨勢,是一種自適應(yīng)且高效的處理方法;(2)基于沖擊響應(yīng)譜定義的修正指數(shù),能夠較好地反映信號修正前后在沖擊響應(yīng)譜不同頻段上的平均偏差幅度;隨著白噪聲系數(shù)增大,不同頻段上修正指數(shù)均不同程度上升,二者呈現(xiàn)相關(guān)性較高的冪指數(shù)函數(shù)關(guān)系;(3)通過修正指數(shù)并結(jié)合對現(xiàn)場測試傳感器的頻響范圍分析,確定出加速度零漂信號對應(yīng)的白噪聲系數(shù)的最優(yōu)取值范圍,優(yōu)化的噪聲系數(shù)使得EEMD 法在最大程度削弱信號零漂趨勢的同時盡量保留信號的頻譜信息;(4)優(yōu)化EEMD 法中的噪聲系數(shù)只能改善加速度積分后速度信號的零漂趨勢而無法將其徹底去除;相對應(yīng)的滑動平均法在速度信號的零漂修正效果上表現(xiàn)更好,但同時也帶來了修正前后頻譜差異加大以及動力響應(yīng)的過度削弱的問題。

猜你喜歡
噪聲系數(shù)頻段修正
Some new thoughts of definitions of terms of sedimentary facies: Based on Miall's paper(1985)
修正這一天
快樂語文(2021年35期)2022-01-18 06:05:30
gPhone重力儀的面波頻段響應(yīng)實(shí)測研究
地震研究(2021年1期)2021-04-13 01:04:56
脈沖多普勒火控雷達(dá)系統(tǒng)接收通道噪聲系數(shù)分析
合同解釋、合同補(bǔ)充與合同修正
法律方法(2019年4期)2019-11-16 01:07:28
功分器幅相不一致對多路合成網(wǎng)絡(luò)噪聲系數(shù)的影響分析
最佳噪聲系數(shù)的接收機(jī)系統(tǒng)設(shè)計?
軟件修正
推擠的5GHz頻段
CHIP新電腦(2016年3期)2016-03-10 14:07:52
TD—LTE在D頻段和F頻段的覆蓋能力差異
中國新通信(2015年1期)2015-05-30 10:30:46
五寨县| 长宁县| 杂多县| 阿勒泰市| 墨竹工卡县| 潞西市| 崇明县| 左贡县| 乌拉特前旗| 措美县| 泽州县| 黄山市| 新野县| 肇庆市| 武隆县| 黔南| 阿图什市| 石景山区| 南木林县| 新田县| 保靖县| 上高县| 宽城| 荆门市| 辽宁省| 潼关县| 宜春市| 阿勒泰市| 东台市| 石屏县| 哈尔滨市| 宕昌县| 大洼县| 海南省| 噶尔县| 调兵山市| 资溪县| 阜阳市| 呼和浩特市| 鄂尔多斯市| 武清区|