脈動風速隨機傅里葉譜模型的參數(shù)聚類分析

洪 旭， 彭勇波， 李 杰 ,

(1. 同濟大學 土木工程學院， 上海 200092； 2. 同濟大學 土木工程防災國家重點實驗室， 上海 200092)

作為數(shù)據(jù)分析方法的一種，聚類分析旨在將數(shù)據(jù)集合按相似程度進行分類[1].在上世紀40年代，Tryon在心理學研究中最早提出了聚類分析方法[2]；在當代，通過聚類分析，互聯(lián)網(wǎng)搜索引擎可以針對用戶的查詢需求將搜索結(jié)果分類分層顯示[3]；在地球物理科學中，聚類分析被用于影響氣候的主要因素的研究[4]；在分子生物學中，聚類分析在功能基因組的研究中起到重要作用[5]，等等.考察應用聚類分析的動機與結(jié)果不難發(fā)現(xiàn)，聚類分析的主要用途在于揭示數(shù)據(jù)之間的潛在聯(lián)系與內(nèi)部結(jié)構(gòu)，并根據(jù)這些特征結(jié)構(gòu)合理地預測其他屬性.由于數(shù)據(jù)之間特征組織方式的形成往往反映著特定物理背景的存在，因此聚類分析能夠為物理模型勾勒出大致輪廓.顯然，這種分析方法也可以應用于風工程中.

在風工程學科誕生之初，不少學者依據(jù)對自然風速觀測記錄的研究，建立了以功率譜密度為核心的統(tǒng)計分析方法[6-8].然而，盡管功率譜密度模型抓住了風速隨機過程的主要特征，卻并不能描述風中湍流的全部概率信息，也難以對非平穩(wěn)的風速時程做出合理的反映.

為了解決傳統(tǒng)結(jié)構(gòu)動力激勵建模中概率信息難以完備描述的問題，文獻[9-11]提出了物理隨機建模方法，明確指出物理規(guī)律是概率信息轉(zhuǎn)移的軌道.對于脈動風速時程，文獻[12-14]基于大氣邊界層能譜與零點演化時間概念給出了模擬脈動風速的隨機Fourier譜模型.研究證明，該模型對實測結(jié)果的概率統(tǒng)計特征具有良好的刻畫能力.

物理隨機系統(tǒng)的隨機性源自基本隨機變量，為了合理描述脈動風速乃至結(jié)構(gòu)風致響應的概率信息，需要仔細處理基本隨機變量及其概率分布.應該指出：在上述隨機Fourier譜模型中，對關(guān)鍵參數(shù)——分界波數(shù)的處理仍是粗糙的；已有工作的數(shù)據(jù)基礎(chǔ)局限于個別風場觀測臺陣，模型對不同地域的適用性仍是需要研究的問題.有鑒于此，本文利用廈門風場實測記錄，借助聚類分析方法，分析了隨機Fourier譜中基本隨機變量的分布，合理選擇數(shù)據(jù)并進行建模.

1 隨機Fourier譜模型

物理隨機系統(tǒng)的觀點認為，物理規(guī)律是隨機性傳播的途徑[10].如果記脈動風速沿時間t的變化為u(Θ,t)，對于給定的參數(shù)Θ，就能夠獲得一條風速時程u(Θ,t).將該樣本時程u(Θ,t)進行Fourier變換，可得到等價的描述方式

F(Θ,n)=|F(Θ,n)|eiφ(Θ,n)=

(1)

式中：n為自然頻率；T為脈動風速的總時長；|F(Θ,n)|稱為幅值譜，它說明了能量在不同頻率的分布；φ(Θ,n)稱為相位譜，它確定了不同頻率處能量在時間上的分布.

由于Θ為隨機變量，上述隨機Fourier譜(幅值譜與相位譜)也是隨機函數(shù).因此，只要給出幅值譜與相位譜的物理模型并確定隨機參數(shù)的概率信息，就能對脈動風速隨機過程進行合理的反映.

實際上，充分發(fā)展的湍流由一系列不同尺度的渦旋組成.大尺度渦旋在其與主流剪切的相互作用中吸取能量；而各尺度的渦旋又不斷破碎形成更小尺度的渦旋直到黏性起主導作用的Kolmogorov尺度，在此過程中能量逐級傳遞；最終在黏性的作用下湍流能量耗散為大氣的內(nèi)能.研究指出，在與主流相互作用強烈的含能子區(qū)和能量傳遞達到平衡的慣性子區(qū)中，能譜密度隨波數(shù)的變化呈不同冪律關(guān)系[15].據(jù)此，文獻[12]提出了雙線性Fourier幅值譜模型

(2)

圖1 能量級串與Fourier幅值譜

在文獻[13]中，分界波數(shù)可以通過渦量比值確定，但這一途徑應用了比較粗糙的經(jīng)驗擬合方法，本質(zhì)上是統(tǒng)計均值.本文的統(tǒng)計結(jié)果表明，分界波數(shù)具有較大的變異性，將其視為獨立的隨機變量是更合理的處理方式.即對于隨機Fourier幅值譜，剪切波速u*與分界波數(shù)kc應分別作為獨立隨機變量.

在相位的研究方面，文獻[14]建立了演化相位譜，它認為不同頻率的渦旋振動速度不同，因而相位變化速度不同，通過合理猜想，這些渦旋具有相同的相位演化起點，并將共同起點時刻到實際觀測時刻經(jīng)歷的演化時間定義為零點演化時間Te，通過幅值譜可以得到渦旋的振動速度為

(3)

式中：Δn是頻率寬度.于是，對應頻率n的演化相位為

Φ(Θ,n,Te)=v(Θ,n)k(n)Te

(4)

注意到零點演化時間在物理上表征了大氣中旋渦相位的演化時間，它在時間尺度上一般較大，且受局部場地風環(huán)境影響較小，因此文獻[14]的結(jié)果是合理的.

值得指出的是，作為風場平均風的描述，風剖面對刻畫風場的宏觀特性十分重要.在諸多風剖面模型中，對數(shù)律的風剖面因為有一定的理論基礎(chǔ)且與實際結(jié)果較吻合而得到廣泛的應用[16].因此，本文采用對數(shù)律作為風場平均特性的度量，即

(5)

式中：U(z)為高度z處的10 min平均風速；z0為地面粗糙度.

根據(jù)式(5)不難理解，風剖面的確定僅需要2個獨立參數(shù)，本文采用剪切波速u*與地面粗糙度z0這一組合作為風剖面基本參數(shù).在物理意義上加以考察，這種選擇是合理的：剪切波速u*本質(zhì)上反映了渦動黏性力在大氣邊界層中對平均風速的滯變效應，它建立起風的平均大尺度特性與高頻細部結(jié)構(gòu)之間聯(lián)系的橋梁[6]；地面粗糙度z0則將邊界條件引起的底部渦旋以特征長度的方式計入風剖面，它實質(zhì)上是邊界條件產(chǎn)生的阻力對平均風速的影響[17].這樣，加上分界波數(shù)kc，本文所研究的風速模型有3個基本隨機變量.

2 參數(shù)識別

2.1 廈門風場觀測臺陣

本研究小組于2014年在廈門大嶝島(北緯24°33′9″, 東經(jīng)118°19′47″)建立了一座強風觀測臺陣.這一臺陣共有4座觀測塔P1、P2、P3、P4，它們的空間布置如圖 2所示.在P1塔與P4塔的10 m與20 m處安裝有三維超聲風速儀，在P2塔與P3塔的10 m、20 m、30 m及40 m處安裝有三維超聲風速儀，采樣頻率均為10 Hz.觀測臺陣周圍以農(nóng)業(yè)用地與灌木為主，在其四周100～200 m遠處分布著多個建筑密集區(qū)，多為10 m以下民房，總體上該地區(qū)的地貌與《建筑結(jié)構(gòu)荷載規(guī)范》[18]的B類地貌或C類地貌相當.在本文的分析中，采用了2015年6月28日至2015年8月19日之間的數(shù)據(jù)，經(jīng)過檢查，所采用數(shù)據(jù)都處于良態(tài)風環(huán)境.經(jīng)過對比發(fā)現(xiàn)：由其他各塔的數(shù)據(jù)得出的結(jié)論與P3塔相一致，因此這里只討論P3塔的結(jié)果.

圖2 風場觀測臺陣的布置(單位：m)

2.2 參數(shù)識別方法與結(jié)果

本文采用最小二乘法進行參數(shù)識別.這一方法實質(zhì)上是在合適的距離定義下尋找實測樣本的最佳逼近.考慮實測樣本y與模型y(Θ)的二范數(shù)誤差

(6)

式中：y是觀測值；y(Θ)是物理模型.參數(shù)識別結(jié)果ΘId應使式(6)取得最小值.如果將樣本y與模型y(Θ)分別取為實測脈動風速的Fourier幅值譜與雙線性Fourier幅值譜(式(2))，就能夠識別出其基本參數(shù)剪切波速u*與分界波數(shù)kc；同理，可以進行對數(shù)律風剖面的識別.

應當指出，根據(jù)脈動風速Fourier幅值譜(式(2))與風剖面(式(5))都可以識別得到剪切波速.如前所述，在Fourier譜中，剪切波速反映了考察位置處渦動黏性引起的主流剪切；而對數(shù)律中剪切波速表征了豎向剪切的平均.因此，本文先識別了各高度處Fourier幅值譜中的剪切波速；然后取其平均值作為風剖面的剪切波速進一步識別地面粗糙度.圖 3和圖 4是同一10 min間隔內(nèi)的典型脈動風速Fourier幅值譜及風剖面的識別結(jié)果.可以看到，各高度處的剪切波速十分接近，說明這一識別方法是合理的.

將所有識別結(jié)果作統(tǒng)計分析，u*、z0及kc的結(jié)果分別見圖 5、 圖 6與圖 7.與以前研究結(jié)果相比[19]，關(guān)于u*識別結(jié)果的統(tǒng)計分布沒有本質(zhì)差別.但大量z0集中在10 m左右，顯然不符合實際地面粗糙度分布范圍；此外3個高度kc的識別結(jié)果中出現(xiàn)了明顯的雙峰特征，這表明分界波數(shù)識別結(jié)果中可能存在特有的數(shù)據(jù)結(jié)構(gòu).

為初步檢驗由上述識別結(jié)果得到的概率分布對風速時程基本概率信息的影響，可以考察它們所導出的功率譜密度函數(shù)并與經(jīng)典譜對比.這里取Kaimal譜作為參照

(7)

式中：f=nz/U；S(z,n)為脈動風速的功率譜.注意到功率譜是樣本Fourier譜的平均

(8)

采用MATLAB中概率分布擬合的核密度估計方法，給出各高度分界波數(shù)的概率密度函數(shù)，如圖 7所示.在此基礎(chǔ)上，結(jié)合式(8)計算了由隨機Fourie譜模型導出的量綱一功率譜，如圖 8所示.顯然，在高頻段模型結(jié)果與Kaimal譜基本一致，但在中、低頻段很大范圍內(nèi)，由分界波數(shù)直接識別結(jié)果導出的功率譜密度函數(shù)顯著低于Kaimal譜.

a 10 m

b 20 m

c 30 m

d 40 m

圖4 典型風剖面與識別結(jié)果

圖5 剪切波速統(tǒng)計結(jié)果

圖6 地面粗糙度統(tǒng)計結(jié)果

為深入研究上述問題， 在圖 9的各子圖中作出不同高度分界波數(shù)的關(guān)系，其中各點的橫、縱坐標為相應時刻不同高度分界波數(shù).顯然，任意兩個高度處分界波數(shù)組成的點集都涇渭分明地分為兩簇.這一特征及其原因可以通過聚類方法進行分析.

3 分界波數(shù)的聚類分析與建模

3.1 聚類分析的基本原理

已經(jīng)指出，聚類分析能夠利用樣本屬性的相似性對樣本集合進行分類.本文采用歐幾里得距離作為樣本相似性的描述，即：將n維點集X= {xi,i= 1,2,…,n}中的點分成N類{Ci,i= 1, 2, …,N}，使得總偏差最小，總偏差的定義為

a 10 mb 20 m

c 30 md 40 m

圖7 分界波數(shù)統(tǒng)計結(jié)果

c 30 md 40 m

圖8由原始識別結(jié)果擬合分布導出的功率譜

Fig.8Powerspectrumdensityderivedfromstatisticaldistributionoforiginalidentification

(9)

一般可以使用經(jīng)典的遞歸算法實現(xiàn)上述目標：

(1) 首先適當?shù)剡x取N個中心點，作為各類的特征點，將集合X中的點按照距所選擇的中心點的遠近將其歸類，接下來重復(2)與(3)的操作，直至分類最終穩(wěn)定；

(2) 按上一次分類的結(jié)果重新計算N個中心點，即每一類所有點的平均值；

(3) 按新的中心點集進行分類.

根據(jù)對圖 9的分析，本文分類數(shù)N取為2即可.

a 橫軸-10 m, 縱軸-20 mb 橫軸-10 m, 縱軸-30 m

c 橫軸-10 m, 縱軸-40 md 橫軸-20 m, 縱軸-30 m

e 橫軸-20 m, 縱軸-40 mf 橫軸-30 m, 縱軸-40 m

圖9分界波數(shù)識別結(jié)果

Fig.9Clusterresultsofcut-offwavenumber

3.2 分界波數(shù)的聚類結(jié)果

將同一時間段內(nèi)、4個高度識別出的分界波數(shù)看作R4空間(每一維度代表某一高度的分界波數(shù)值)的一個點，則所有記錄得到的4維分界波數(shù)數(shù)組全體構(gòu)成了待聚類的點集.本文利用MATLAB中k-means函數(shù)進行聚類，聚類結(jié)果如圖 9所示.顯然，通過聚類分析，上節(jié)所述的分類特征得以細致地展現(xiàn).需要說明的是，圖 9中“+”代表實測風剖面與對數(shù)律風剖面相對誤差大于5.5%的樣本，而聚類分析只考慮相對誤差小于5.5%的樣本.

3.3 聚類現(xiàn)象原因分析

在極坐標圖 10中作出各高度10 min平均風速U10與風向的關(guān)系，并用相應的記號表現(xiàn)上述聚類結(jié)果，對圖 10的考察可以得出以下規(guī)律：

(1) 風剖面受顯著干擾的樣本幾乎完全集中在0°(正北)方向以及附近區(qū)域；

(2) 與其他方向比較，平均風速在0°(正北)方向有明顯的凹進，說明它在此處被顯著削弱；

綜合上述結(jié)果，可以推測：當風從正南方吹來時，受到某種因素的擾動，風速被削弱，且風速受到的擾動隨著風向與正北方的夾角增大而減小.考察風場觀測臺陣設(shè)備設(shè)施，注意到各高度處的三維超聲風速儀都安裝在由觀測塔朝正北方向伸出約1 m左右的構(gòu)件上，而且觀測塔直徑也在1 m左右，則有理由相信，觀測塔是擾動實測記錄的主要因素.

a 10 mb 20 m

c 30 md 40 m

圖10分界波數(shù)聚類結(jié)果與風向的關(guān)系(圖中半徑大小代表10m高、10min平均風速，單位:m·s-1)

Fig.10Relationofclusterresultstowinddirection(Themagnitudeofradiusdesignate10minmeanwindspeedataheightof10m,unit:m·s-1)

在已有的風場數(shù)據(jù)分析文獻中，鮮有考慮觀測設(shè)備設(shè)施對實測記錄影響的工作.本文工作說明，為了獲得正確合理的結(jié)論，應當重視觀測系統(tǒng)自身帶來的誤差.

3.4 基于聚類分析的分界波數(shù)建模

前文分析指出，將分界波數(shù)處理為隨機變量是一種更合理的建模方式.為了正確地統(tǒng)計其概率分布，應當將上文討論的外界擾動排除在外.據(jù)此，本文將風向在正北方向附近一定方位角內(nèi)及風剖面與對數(shù)律相對誤差大于5.5%的樣本記錄排除.由此得到的地面粗糙度的直方圖見圖 11，以及各高度處分界波數(shù)的對數(shù)直方圖見圖 12.可見，圖 6的地面粗糙度異常分布及圖 7分界波數(shù)的雙峰分布不再存在.

圖11 數(shù)據(jù)篩選后的地面粗糙度統(tǒng)計直方圖

a 10 mb 20 m

c 30 md 40 m

圖12分界波數(shù)直方圖與建模

Fig.12Histogramandmodelingofcut-offwavenumber

進一步考察樣本篩選是否能夠改善功率譜密度的模擬.采用對數(shù)正態(tài)分布對分界波數(shù)進行建模，建模結(jié)果如表 1所示.分析可見，kc對數(shù)均值μln kc與對數(shù)標準差σln kc與高度有明顯的線性關(guān)系(見圖 13與圖 14).

μln kc=-0.016 3z-0.809 8

σln kc=0.003 6z+0.186 2

式中：高度z單位為m；kc單位為m-1.

表1 分界波數(shù)建模結(jié)果

圖13 分界波數(shù)對數(shù)均值與高度的關(guān)系

圖14 分界波數(shù)對數(shù)標準差與高度的關(guān)系

按2.2節(jié)方法計算由隨機Fourier譜模型導出的功率譜，并與Kaimal譜對比，如圖 15所示.顯然，排除異常數(shù)據(jù)后，模型功率譜與Kaimal譜在一般工程所關(guān)心的頻率范圍內(nèi)幾乎一致，證明聚類分析與樣本篩選的有效性.同時，上述分析進一步驗證了隨機Fourier譜模型對脈動風速細部特征進行合理刻畫的能力.

4 結(jié)論

本文研究了脈動風速隨機Fourier譜模型基本隨機變量及其分布，得到以下結(jié)論：

a 10 mb 20 m

c 30 md 40 m

圖15數(shù)據(jù)篩選后的參數(shù)統(tǒng)計分布導出的功率譜

Fig.15Powerspectrumdensityderivedfromstatisticaldistributionafterdatasieving

(1) 脈動風速隨機Fourier譜模型的基本隨機變量包括剪切波速、地面粗糙度、分界波數(shù)以及零點演化時間.利用廈門地區(qū)實測風速記錄，識別了剪切波速、地面粗糙度、分界波數(shù).

(2) 初始識別結(jié)果中地面粗糙度分布異常；分界波數(shù)呈現(xiàn)雙峰特征.初始識別結(jié)果的統(tǒng)計分布導出的功率譜密度函數(shù)明顯偏離Kaimal譜.

(3) 研究發(fā)現(xiàn)，上述參數(shù)識別結(jié)果的異常源自風速樣本的聚類特征，且與風向緊密聯(lián)系.這一現(xiàn)象與觀測裝置對風場的擾動有關(guān)，因此在分析實測風場記錄時，應重視風場設(shè)備設(shè)施對觀測效果的影響.

(4) 依據(jù)聚類分析結(jié)果，合理地選擇樣本集合，進行了分界波數(shù)統(tǒng)計建模.由此導出的模型功率譜與Kaimal譜對比良好，這也證明了聚類分析的有效性.

參考文獻:

[1] JAIN A K. Data clustering: 50 years beyond K-means[J]. Pattern Recognition Letters, 2010, 31(8): 651.

[2] 方開泰, 潘恩沛. 聚類分析[M]. 北京：地質(zhì)出版社, 1982.

FANG Kaitai, PAN Enpei. Cluster analysis[M]. Beijing: Geological Publishing House, 1982.

[3] 蒼宏宇, 譚宗穎. 聚類搜索引擎發(fā)展現(xiàn)狀研究[J]. 圖書情報工作, 2009, 53(2):125.

CANG Hongyu, TAN Zongying. Research on the development situation of clustering search engine[J]. Library & Information Service, 2009, 53(2):125.

[4] TAN P N, STEINBACH M, KUMAR V. Introduction to data mining [M]. Boston: Pearson Education Inc, 2006.

[5] BALDI P. Bioinformatics: the machine learning approach[M]. Boston: MIT Press, 2001.

[6] SIMIU E, SCANLAN R H. Wind effects on structures[M]. New York: John Wiley & Sons, 1996.

[7] COUNIHAN J. Adiabatic atmospheric boundary layers: a review and analysis of data from the period 1880-1972[J]. Atmospheric Environment, 1975, 9(10): 871.

[8] 閻啟, 謝強, 李杰. 風場長期觀測與數(shù)據(jù)分析[J]. 建筑科學與工程學報, 2009, 26(1): 37.

YAN Qi, XIE Qiang, LI Jie. Long-term observation and data analysis of wind field[J]. Journal of Architecture & Civil Engineering, 2009, 26(1): 37.

[9] 李杰. 隨機動力系統(tǒng)的物理逼近[J]. 中國科技論文在線, 2006, 1(2): 95.

LI Jie. A physical approach to stochastic dynamical systems[J]. Sciencepaper Online, 2006, 1(2): 95.

[10] 李杰. 物理隨機系統(tǒng)研究的若干基本觀點 [R]. 上海: 同濟大學, 2006.

LI Jie. Basic views of the stochastic physical system research[R]. Shanghai: Tongji University, 2006.

[11] LI J, CHEN J. Stochastic dynamics of structures[M]. Singapore: John Wiley & Sons, 2009.

[12] 李杰, 閻啟. 結(jié)構(gòu)隨機動力激勵的物理模型: 以脈動風速為例[J]. 工程力學, 2009 (A02): 175.

LI Jie, YAN Qi. Physical models for the stochastic dynamic excitations of structures: in the case of fluctuating wind speed[J]. Engineering Mechanics, 2009 (A02): 175.

[13] 李杰, 閻啟. 脈動風速隨機 Fourier 波數(shù)譜研究[J]. 同濟大學學報 (自然科學版), 2011, 39(12): 1725.

LI Jie, YAN Qi. Research on stochastic Fourier wave-number spectrum of fluctuating wind speed[J]. Journal of Tongji University(Natural Science), 2011, 39(12):1725.

[14] LI J, PENG Y, YAN Q. Modeling and simulation of fluctuating wind speeds using evolutionary phase spectrum[J]. Probabilistic Engineering Mechanics, 2013, 32: 48.

[15] KATUL G, CHU C R. A theoretical and experimental investigation of energy-containing scales in the dynamic sublayer of boundary-layer flows[J]. Boundary-Layer Meteorology, 1998, 86(2): 279.

[16] BLACKADAR A K, TENNEKES H. Asymptotic similarity in neutral barotropic planetary boundary layers[J]. Journal of the Atmospheric Sciences, 1968, 25(6): 1015.

[17] DYRBYE C, HANSWN S O. Wind loads on structures[M]. Chichester: John Wiley & Sons, 1996.

[18] 中華人民共和國住房和城鄉(xiāng)建設(shè)部. 建筑結(jié)構(gòu)荷載規(guī)范：GB 50009—2012[S]. 北京：中國建筑工業(yè)出版社,2012.

Ministry of Housing and Urban-Rural Development of the People’s Republic of China. Load code for the design of building structures: GB 50009—2012[S]. Beijing: China Architecture & Building Press,2012.

[19] LI Q S, ZHI L, HU F. Boundary layer wind structure from observations on a 325m tower[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2010, 98(12): 818.