欧美在线观看视频网站,亚洲熟妇色自偷自拍另类,啪啪伊人网,中文字幕第13亚洲另类,中文成人久久久久影院免费观看 ,精品人妻人人做人人爽,亚洲a视频

低復(fù)雜度的分布式麥克風陣列聲源定位系統(tǒng)的制作方法

文檔序號:11275192閱讀:357來源:國知局
低復(fù)雜度的分布式麥克風陣列聲源定位系統(tǒng)的制造方法與工藝

本發(fā)明涉及一種低復(fù)雜度的分布式麥克風陣列聲源定位系統(tǒng)。主要涉及專利分類號g01測量;測試g01s無線電定向;無線電導(dǎo)航;采用無線電波測距或測速;采用無線電波的反射或再輻射的定位或存在檢測;采用其他波的類似裝置g01s5/00通過確定兩個或更多個方向或位置線的配合來定位;通過確定兩個或更多個距離的配合進行定位g01s5/18應(yīng)用了超聲波、聲波或次聲波g01s5/20由多個分隔開的定向器確定的信號源位置。



背景技術(shù):

在服務(wù)機器人等應(yīng)用中,常需定位聲源。相比使用傳統(tǒng)規(guī)則結(jié)構(gòu)的麥克風陣列,分布式麥克風陣列不僅靈活性好,而且具有更大的空間孔徑,可以獲得更好的定位性能,同時,對麥克風節(jié)點失能也具有一定容錯性。

現(xiàn)有的基于可控波束形成、信號到達時間差(tdoa)和空間譜估計等定位方法總體而言精度受混響和噪聲影響較大,個別算法還受麥克風陣列幾何結(jié)構(gòu)影響,導(dǎo)致算法可移植性差,對陣列節(jié)點失能容錯性差。

特征定位算法因其很好地利用室內(nèi)信號復(fù)雜的反射折射等信息來生成特征庫,使得該類方法在室內(nèi)定位問題中被廣泛使用。但是,受硬件處理器運算能力限制,特征方法較高的復(fù)雜度成為其進一步發(fā)展的障礙。因而,在定位精度足夠情況下,降低室內(nèi)定位方法低復(fù)雜度的就成為一個發(fā)明課題。

現(xiàn)有主流技術(shù)之一,采用了最大似然估計器計算信號到達時延,進而完成聲源定位。最大似然估計器來估計tdoa漸進等價于非線性最小二乘法。首先根據(jù)每個麥克風陣列節(jié)點接收到的信號,估計出響應(yīng)信號到達時延;然后通過添加距離坐標分量建立一個四維坐標系,從而構(gòu)造一個自變量為x、y和z坐標以及距離坐標的目標函數(shù);接下來,借助錐面定位法并且通過與參考麥克風的位置相比較,最終得到聲源位置。

該方法在構(gòu)建關(guān)于聲源位置目標函數(shù)時,需要預(yù)先知道所有麥克風的準確坐標,并且需要一個參考麥克風,這在很多應(yīng)用場景中無法實現(xiàn)。與此同時,在每個麥克風節(jié)點內(nèi),都假設(shè)麥克風陣列為均勻線性陣列,以確定接收信號模型;其次,該方法最終通過在錐面搜索到聲源的最小距離來得到聲源位置,而這個過程的目標函數(shù)非線性,不僅所得結(jié)果存在非全局最優(yōu)解的可能,而且運算過程復(fù)雜。

現(xiàn)有技術(shù)中,也出現(xiàn)了先用分步格點劃分法將整個定位區(qū)域劃分成若干個格點,每個麥克風分別接收聲源發(fā)出的信號;其次,通過對接收信號加窗處理,并經(jīng)過兩次離散余弦變換(dct)提取語音信號特征;接下來,利用上一步得到的降維特征向量構(gòu)造稀疏定位模型,并增加模型失配修正環(huán)節(jié)以將除噪聲外的其他因素考慮在內(nèi);最后,借助字典學習和稀疏重構(gòu),實現(xiàn)定位目的。

該方法所采用的稀疏重構(gòu)在求解過程中,目標函數(shù)存在收斂困難問題,并且在信號信噪比低時會導(dǎo)致定位準確度下降。另外,該發(fā)明中字典學習與更新過程使得其復(fù)雜度較高,運算量較大。



技術(shù)實現(xiàn)要素:

本發(fā)明要解決的技術(shù)問題是一種低復(fù)雜度的分布式麥克風陣列聲源定位系統(tǒng),麥克風陣列中的每個麥克風節(jié)點中由至少2個麥克風組成,包括:

離線訓練單元,該單元包括:

分別計算p個待定位點到n個麥克風節(jié)點中所述兩個麥克風的時間差,即時延τ的時延估計模塊i,每個麥克風節(jié)點對應(yīng)p個時延τ1、τ2……τp;

根據(jù)所述每個麥克風節(jié)點對應(yīng)的p個時延,計算每個麥克風節(jié)點的時延分布概率密度函數(shù)的時延分布函數(shù)估計模塊;

所述每個麥克風節(jié)點的時延分布概率密度函數(shù)記錄在特征庫中;

定位單元,該單元包括:

時延估計模塊ii,分別計算待定位位置處的聲源到n個麥克風節(jié)點中所述兩個麥克風的時間差,即時延τ,每個待定位位置對應(yīng)n個時延τ1、τ2……τn;

信度值求取模塊,該模塊分析由所述時延估計模塊ii計算得出的待定位位置的n個時延τ對應(yīng)一個所述待定位點的信度值m;

針對一個待定位的音源點,系統(tǒng)中的n個麥克風節(jié)點對應(yīng)n個信度值,m1、m2……mn;

分布式數(shù)據(jù)融合模塊,該模塊根據(jù)系統(tǒng)內(nèi)每個待定位位置對應(yīng)的n個信度值,采用加權(quán)ds理論進行融合,得出系統(tǒng)總信度值bel;

定位模塊,該模塊比較選取的待定位點位置信度值,得到與p個指紋點位置信度值最接近的標記位置,即為該待定位置的所述標記。

作為優(yōu)選的實施方式,所述的時延估計模塊i采用廣義互相關(guān)法并配合相位變換加權(quán)函數(shù)(gcc-phat)來求得時延,過程如下:

定義一個麥克風節(jié)點中的兩個麥克風接收到的信號分別為x1(t)和x2(t),聲源到麥克風節(jié)點的時延τ12可由式(1)、式(2)求得,

其中,x1(f)和x2(f)分別是兩個麥克風接收到的信號x1(t)和x2(t)的傅里葉變換,借助fft算法近似求得兩時域信號的傅里葉變換結(jié)果,(·)*表示共軛操作,e表示自然對數(shù)的底,j為虛單位,j2=-1,f為頻率,τ12表示時延。

作為優(yōu)選的實施方式,所述的時延分布函數(shù)估計模塊計算過程如下:

計算該節(jié)點處的時延分布概率密度函數(shù)公式為

式中,p為每個麥克風節(jié)點處得到的時延個數(shù),也是位置標記數(shù)。h為帶寬(或稱窗口),決定估計函數(shù)的平滑度。k(·)為核函數(shù),滿足對稱性及歸一性(∫k(x)dx=1),本發(fā)明中核函數(shù)選擇高斯形式,xc為核函數(shù)中心,具體公式為:

最終得到所有n個麥克風節(jié)點處的時延分布的概率密度函數(shù)

作為優(yōu)選的實施方式,所述信度值求取模塊利用時延分布估計模塊得到的該麥克風節(jié)點處的時延分布概率密度函數(shù),當給定一個時延τ時,對應(yīng)得到其信度值m,如下公式:

當整個系統(tǒng)有n個麥克風節(jié)點時,對于一個待定位點,就可以得到n個信度值m1,m2,…,mn。

作為優(yōu)選的實施方式,所述的分布式數(shù)據(jù)融合模塊的加權(quán)ds理論的融合過程如下:

其中,wi和wj分別表示第i個麥克風節(jié)點和第j個麥克風節(jié)點的權(quán)重;⊕表示ds融合;權(quán)重為

式中max{x,y}表示取x,y中的最大值,n為麥克風節(jié)點總數(shù);針對n個麥克風節(jié)點的分布式數(shù)據(jù)融合,一致性算法

迭代過程如下:

其中∪表示對于麥克風節(jié)點i的所有通信節(jié)點處的信度值按上述加權(quán)ds融合規(guī)則進行融合,l=10為迭代次數(shù)。

通過采用上述技術(shù)方案,本發(fā)明提出的一種低復(fù)雜度的分布式麥克風陣列聲源定位系統(tǒng),有益效果如下:在滿足精度要求的同時,本發(fā)明算法降低了復(fù)雜度和運算量,核心算法沒有復(fù)雜的非線性計算。對低信噪比高混響環(huán)境的適應(yīng)性較強,在該類型環(huán)境中仍然可以保證較高準確度。采用分布式數(shù)據(jù)融合算法后,針對分布式麥克風節(jié)點有損壞的場景,依然可以得到令人滿意的定位效果,增強了算法的可移植性和魯棒性。

附圖說明

為了更清楚的說明本發(fā)明的實施例或現(xiàn)有技術(shù)的技術(shù)方案,下面將對實施例或現(xiàn)有技術(shù)描述中所需要使用的附圖做一簡單地介紹,顯而易見地,下面描述中的附圖僅僅是本發(fā)明的一些實施例,對于本領(lǐng)域普通技術(shù)人員來講,在不付出創(chuàng)造性勞動的前提下,還可以根據(jù)這些附圖獲得其他的附圖。

圖1為本發(fā)明的系統(tǒng)模塊圖

圖2為本發(fā)明模擬實驗中會議室俯視平面圖

圖3為本發(fā)明標記判斷誤差定義示意圖

圖4為本發(fā)明中混響時間對定位準確度的影響示意圖

圖5為本發(fā)明中噪聲對定位準確度的影響

圖6為本發(fā)明中分布式麥克風節(jié)點邏輯等效圖

圖7為本發(fā)明中節(jié)點損壞與定位誤差的關(guān)系示意圖

圖8為本發(fā)明實施例中三種算法在不同混響條件下的對比示意圖

圖9為本發(fā)明實施例中三種算法在不同信噪比條件下的對比示意圖

具體實施方式

為使本發(fā)明的實施例的目的、技術(shù)方案和優(yōu)點更加清楚,下面結(jié)合本發(fā)明實施例中的附圖,對本發(fā)明實施例中的技術(shù)方案進行清楚完整的描述:

如圖1-圖9所示:

本發(fā)明的整體框圖如圖1所示,總體包括兩部分,即離線訓練模塊(10)和定位模塊(20)。

其中離線訓練模塊(10)中包括時延估計模塊(101)、時延分布函數(shù)估計模塊(102)和構(gòu)建指紋庫模塊(103);定位模塊(20)包括時延估計模塊(201)、信度值求取模塊(202)、分布式數(shù)據(jù)融合模塊(203)和定位模塊(204)。

其中時延估計模塊(101)與時延估計模塊(201)功能、原理和實現(xiàn)方式相同。

離線訓練模塊(10)

訓練階段,在室內(nèi)可允許人員走動的范圍內(nèi),每間隔1米設(shè)定一個標記位置。本實施例中每兩個麥克風構(gòu)成一個節(jié)點,多個節(jié)點構(gòu)成麥克風陣列,隨機分布在室內(nèi)。

本模塊主要目的是得到每個麥克風節(jié)點處關(guān)于時延分布的概率密度函數(shù),以此構(gòu)成指紋庫,具體流程如表1所示。

表1時延分布的概率密度函數(shù)估計流程

時延估計模塊(101)

時延估計是指求出聲源到一個麥克風節(jié)點中兩個麥克風的時間差。本發(fā)明采用廣義互相關(guān)法并配合相位變換加權(quán)函數(shù)(gcc-phat)來求得時延。

首先定義一個麥克風節(jié)點中的兩個麥克風接收到的信號分別為x1(t)和x2(t),聲源到麥克風節(jié)點的時延τ12可由式(1)、式(2)求得,

其中,x1(f)和x2(f)分別是兩個麥克風接收到的信號x1(t)和x2(t)的傅里葉變換,在本發(fā)明中借助fft算法近似求得兩時域信號的傅里葉變換結(jié)果,(·)*表示共軛操作,e表示自然對數(shù)的底,j為虛單位,j2=-1,f為頻率,τ12表示時延。

時延分布函數(shù)估計模塊(102)

該模塊估計每個麥克風節(jié)點處時延的分布概率密度,采用核函數(shù)估計法。

在訓練階段,得到所有位置標記點分別到每個麥克風節(jié)點的時延,即針對每個麥克風節(jié)點有一系列時延τ1,τ2,…,τp,計算該節(jié)點處的時延分布概率密度函數(shù)公式為

式中,p為每個麥克風節(jié)點處得到的時延個數(shù),也是位置標記數(shù)。h為帶寬(或稱窗口),決定估計函數(shù)的平滑度。k(·)為核函數(shù),滿足對稱性及歸一性(∫k(x)dx=1),本發(fā)明中核函數(shù)選擇高斯形式,xc為核函數(shù)中心,具體公式為

以此類推,可以得到所有n個麥克風節(jié)點處的時延分布的概率密度函數(shù)這些概率密度函數(shù)即構(gòu)成本發(fā)明的指紋庫(103)。

定位模塊(20)

信度值求取模塊(202)

本模塊求聲源在某處時,到第i個麥克風節(jié)點時延τ對應(yīng)的信度值。本發(fā)明中使用時延分布概率密度函數(shù)作為信度函數(shù),利用時延分布估計模塊(102)得到的該麥克風節(jié)點處的時延分布概率密度函數(shù),當給定一個時延τ時,就可以對應(yīng)得到其信度值m,如下公式:

這樣,當整個系統(tǒng)有n個麥克風節(jié)點時,對于一個待定位點,就可以得到n個信度值m1,m2,…,mn。

分布式數(shù)據(jù)融合模塊(203)

本發(fā)明采用加權(quán)ds理論進行數(shù)據(jù)融合。質(zhì)函數(shù)(又稱信度函數(shù))可表示為

m:p(θ)→[0,1](6)

其中θ為一個元素相互獨立且有限的集合,在本發(fā)明中稱為定位問題的全域,由一系列可能的子命題構(gòu)成。式(6)的含義為給定全域θ的任一子集,都可以得到一個在[0,1]范圍內(nèi)的值,表征該子集對應(yīng)命題的信度值。

質(zhì)函數(shù)m需同時滿足條件于是,對于全域中一子命題a各麥克風節(jié)點所得信度值利用ds理論進行融合,得到總信度值。ds融合規(guī)則為:

式中,命題a、b、c集合滿足都為全域θ的子集;i、j表示麥克風節(jié)點序號。因此,對于某待定位位置,將在信度值求取模塊中得到的n個麥克風節(jié)點處觀測到的關(guān)于該位置的信度值m1,m2,…,mn按上述ds融合規(guī)則進行融合,最終得到所有麥克風節(jié)點對該位置的總信度m。

在本實施例中提出了一種更適應(yīng)實際需求的分布式加權(quán)ds數(shù)據(jù)融合方式,

其中,wi和wj分別表示第i個麥克風節(jié)點和第j個麥克風節(jié)點的權(quán)重;⊕表示ds融合。

在現(xiàn)有技術(shù)中給出了一種基于傳感器在相似條件下的歷史表現(xiàn)估計得到正確率來表示節(jié)點的權(quán)重。本實施例使用不基于歷史先驗信息的權(quán)重確定法(即按照無向圖思想),根據(jù)分布式麥克風陣列的通信半徑,得到與每個節(jié)點通信的節(jié)點個數(shù)d,在本實施例中稱為節(jié)點的度,則權(quán)重為

式中max{x,y}表示取x,y中的最大值,n為麥克風節(jié)點總數(shù)。

針對n個麥克風節(jié)點的分布式數(shù)據(jù)融合,本實施例根據(jù)一致性算法

提出基于所提加權(quán)ds融合規(guī)則的一致性迭代方法,

其中∪表示對于麥克風節(jié)點i的所有通信節(jié)點處的信度值按上述加權(quán)ds融合規(guī)則進行融合。l為迭代次數(shù),本實施例中l(wèi)取10。

定位模塊

該模塊利用數(shù)據(jù)融合模塊得到的所有麥克風節(jié)點對待定位位置的信度值來確定出該待定位位置屬于訓練模塊中的哪個位置標記。

首先,對于每一個待定位位置,構(gòu)造證據(jù)全域θ={a,r,n},其中a指的是命題“聲源在該標記處”,r表示命題“聲源不在該標記處”,n表示命題“不確定聲源是否在該標記處”。則命題a、r和n可信度的計算公式[4]分別為

a=(1-n)×mi(12)

r=(1-n)×(1-mi)(13)

表2ds數(shù)據(jù)融合全域子集集合關(guān)系

在定位階段,只考慮“聲源在該標記處”的總信度值,計算過程參考式(8),其中b、c的選取條件滿足b,c∈{a,r,n}且式(8)中的命題a、b和c的集合關(guān)系滿足如表2(以兩個麥克風節(jié)點為例),其中表示空集。

定位步驟如下:

(1)對于某一待定位位置,n個麥克風節(jié)點分別接收到音頻信號,計算出聲源到每個節(jié)點的時延τ1,τ2,…,τn。

(2)根據(jù)時延分布函數(shù)估計模塊中式(5)得到n個麥克風節(jié)點對于該位置的信度值m1,m2,…,mn。

(3)根據(jù)預(yù)設(shè)定的麥克風節(jié)點通信半徑,以及數(shù)據(jù)融合模塊中式(8)、式(9)和基于所提加權(quán)ds融合規(guī)則的一致性迭代算法式(11)進行分布式數(shù)據(jù)融合,得到n個麥克風節(jié)點對該位置數(shù)據(jù)融合后的信度值m'1,m'2,…,m'n。

(4)再將該n個數(shù)據(jù)融合后的信度值m'1,m'2,…,m'n按加權(quán)ds數(shù)據(jù)融合規(guī)則融合,得到最后該位置的總信度值bel。

(5)將bel與訓練階段得到的p個位置標記處的信度值做比較,得到與該待定位位置信度值最接近的標記位置信度值,即為該待定位位置的所屬標記。

實施例

本實施例中,采用matlab進行仿真,模擬會議室的平面圖如圖2所示,其中房間大小為8×6×3米;會議桌大小為5×2×0.8米;說話人的高度設(shè)定為1.6米;位置標記點的間隔為1米,共有24個標記位置;每個麥克風節(jié)點包含兩個麥克風,間距固定為30厘米,共15對麥克風分布式地散放在會議桌上。用image模型用來模擬房間脈沖響應(yīng)及各麥克風接收到的音頻信號,聲速為343米/秒。

本實施例中定義了“標記判斷誤差”來衡量算法準確度,其定義為待定位位置到所得標記位置的距離與理論所屬的位置距離之差,如下圖3所示。

在仿真實驗中,定義n表示montecarlo實驗次數(shù),采用均方根誤差(rmse)表征聲源位置所屬標記的估計與理論所屬標記的偏離程度。

本實施例方法的抗混響和噪聲性能

圖4給出了相同信噪比(snr=20db)、不同混響時間情況下,對本實施例算法的定位準確度的影響。其中,rt60表示混響時間;ms表示毫秒。圖5給出了相同混響時間rt60=0.3秒、不同信噪比(snr)的條件下,對本實施例算法的定位準確度影響。其中db表示分貝。

從圖4可見,定位誤差會隨著混響時間的加大而增加,但也可發(fā)現(xiàn)每增加100毫秒混響時間,各定位誤差區(qū)間對應(yīng)百分比減小不足5%,在混響為200ms和

300ms條件下,定位誤差有近90%小于0.5米,即使在500ms混響下,定位誤差也有85%小于1米;從圖5可見,噪聲嚴重,即snr=5db時,定位誤差雖有下降,但仍有近65%定位誤差小于1米,相較于位置標記而言,該誤差等價于將待定位的位置標記錯判給相鄰的一個標記;而對于信噪比為15db和25db情況下,分別有近85%和90%定位誤差小于0.5米,這等價于所屬特征標記錯判不到相鄰一個標記。從而說明本算法在給待定位位置進行所屬特征標記分類時,受混響和噪聲影響較小。

本實施例方案抵抗節(jié)點失能的性能

本實施例中,可將麥克風節(jié)點分布及通信情況等效成一種無向圖,如圖6所示。麥克風節(jié)點等效成圖結(jié)構(gòu)中的節(jié)點,圖結(jié)構(gòu)中的邊表示這兩個節(jié)點所代表的麥克風節(jié)點間可通信。并且定義與該節(jié)點通信的節(jié)點個數(shù)為該節(jié)點的度。黑色節(jié)點表示損壞節(jié)點,白色節(jié)點為正常工作節(jié)點,實邊表示正常通信,虛邊表示因節(jié)點損壞導(dǎo)致的通信失敗。

在節(jié)點失能對定位影響的仿真實驗中,選定室內(nèi)環(huán)境參數(shù)為信噪比20db,混響時間t60=0.3秒。首先仿真了無節(jié)點損壞情況下,算法的定位準確度;其次隨機選擇一個節(jié)點損壞,得到對應(yīng)的定位準確度,并做多次montecarlo仿真求均值,以得到損壞一個節(jié)點情況下,算法的定位準確度;最后以同樣的隨機方法,得到損壞兩個節(jié)點的情況下,算法的定位準確度。將三條曲線對比在圖7中顯示。從圖7可見,當有節(jié)點損壞時,本實施例所提方法定位準確度會有所下降,但仍可得到可接受的定位位置所屬標記準確度。

本實施例方案與兩種常見的聲源特征定位方法的對比

參與對比的兩種方法如下:對比方法1:采用聲源到各麥克風節(jié)點的時延作為參考量,并借助支持向量機(svm)進行分類訓練來構(gòu)成特征庫。在定位測試階段,用測得的待定位聲源位置到各麥克風節(jié)點的時延通入到已訓練好的svm分類器中得到所屬位置標記。對比方法2:采用聲源到各麥克風節(jié)點的時延作為參考量直接構(gòu)成特征庫,在定位測試階段,用測得的待定位聲源位置到各麥克風節(jié)點的時延與特征庫中時延做比較,得到與其最接近的時延序號所對應(yīng)的標記,即為待定位聲源位置所屬標記。為了對比需要,本實施例與上述兩文獻中方法對比時,均只取到兩方法特征分類階段的結(jié)果。

圖8和圖9分別顯示了三種算法針對不同混響和信號噪聲比的情況下的性能。圖8比較三種算法的抗混響能力,固定信噪比為20db;圖9比較三種算法的抗噪能力,混響時間t60=0.3秒。橫坐標分別為混響時間(rt)和信噪比(snr),縱坐標表示在定位誤差小于1米(即本實施例中相鄰標記點的間隔)所對應(yīng)的比例。由圖8、圖9可見,三種方法在低混響、高信噪比環(huán)境下性能相近,但對高混響時間或低信噪比環(huán)境,本實施例方法的性能顯著優(yōu)于其它兩種對比方法。三種方法主要運算量的對比。

此處為方便起見,定義了一些特定符號來表示除加乘基本運算外的復(fù)雜操作,具體如表3所示。同時設(shè)一致性算法融合迭代次數(shù)為c;任意一麥克風節(jié)點的度為ns;麥克風節(jié)點總數(shù)為m。則表4給出三種方法分別的主體計算量。

表3運算量計算的相關(guān)符號

表4計算量對比

分析表4可知,對比方法2(即直接以時延作為參考量構(gòu)建特征庫)的運算量最低,與此同時對比方法1(即以時延作為參考量并借用支持向量機構(gòu)建特征庫)的運算量最高。算法一運算量主要高在svm訓練階段權(quán)重迭代更新時的非線性運算,而且與訓練時迭代次數(shù)正相關(guān),事實上支持向量機往往需要多次迭代,因此其運算量最高。而本算法在基本的加乘運算量上比其他兩種方法高出m2,這主要是因為ds融合時需要先算出全域的“不確定”子集對應(yīng)的信度,雖與m2成正比,但往往麥克風節(jié)點數(shù)在10數(shù)量級上,所以相比而言運算量不會劇烈增加。其次本算法在其他運算量上遠小于機器學習相關(guān)方法,較算法二只增加一個概率密度函數(shù)估計的計算量。而由式(3)可知,該估計計算量與位置標記數(shù)正相關(guān),即φ1~o(p),因此實際增加的運算量較小。

以上所述,僅為本實施例較佳的具體實施方式,但本實施例的保護范圍并不局限于此,任何熟悉本技術(shù)領(lǐng)域的技術(shù)人員在本實施例揭露的技術(shù)范圍內(nèi),根據(jù)本實施例的技術(shù)方案及其發(fā)明構(gòu)思加以等同替換或改變,都應(yīng)涵蓋在本實施例的保護范圍之內(nèi)。

當前第1頁1 2 
網(wǎng)友詢問留言 已有0條留言
  • 還沒有人留言評論。精彩留言會獲得點贊!
1
元江| 屏山县| 阳西县| 清镇市| 栾城县| 华安县| 陕西省| 莱州市| 五指山市| 延津县| 特克斯县| 湘阴县| 普兰县| 微博| 高邮市| 钟祥市| 武川县| 临汾市| 资兴市| 滨州市| 伊宁市| 阜阳市| 宿迁市| 光泽县| 大足县| 榕江县| 高阳县| 揭西县| 汉寿县| 成安县| 特克斯县| 庆安县| 广德县| 商南县| 开江县| 龙里县| 黑河市| 苏尼特左旗| 泽州县| 安徽省| 巴中市|