穩(wěn)健最小二乘法在風電塔筒垂直度分析中的應用

程 帥，米 珂，李曉博，舒 進，陳 倉，蘭 昊，汪俊波

（西安熱工研究院有限公司，陜西 西安 710054）

隨著風能開發(fā)利用技術(shù)的逐漸成熟，風電已成為目前最具規(guī)模的可再生能源之一，是我國能源轉(zhuǎn)型的重要方向[1-2]。風電塔筒作為整個風電機組的支撐結(jié)構(gòu)，具有重心高、承受水平力和傾覆彎矩較大等特點，其垂直度對機組安全及穩(wěn)定性具有較強的敏感性[3]，塔筒垂直度超標輕則導致風機塔筒驅(qū)動側(cè)加速度超限，達不到設(shè)計運行標準，重則造成塔筒傾覆、倒塔等嚴重后果[4-5]。因此，提高風電塔筒垂直度分析精度并定期監(jiān)測，對機組的安全穩(wěn)定運行具有重要意義。

風電塔筒垂直度監(jiān)測的主要方法有工程測量法、固定監(jiān)測儀器在線監(jiān)測法及三維激光掃描法等。郭倩倩[6]、曾群意[7]等利用工程測量法得到建筑物的特征點坐標，通過直接擬合法分析類似風電塔筒的高聳柱狀結(jié)構(gòu)垂直度，但該方法僅在人工剔除明顯偏差測點的情況下，用原始測量數(shù)據(jù)直接擬合得到最終結(jié)果，其計算精度完全取決于測量點的誤差大小。付曉敏[8]、李智峰[9]等分別發(fā)明了一種風電機組塔筒垂直度的在線監(jiān)測裝置，將靜力水準儀與激光垂準儀分別應用在塔筒內(nèi)部作為傾斜傳感器監(jiān)測風機塔筒的垂直度。丁克良[10]、王二民[11]等引入三維激光掃描技術(shù)測得建筑物表面的三維空間點云數(shù)據(jù)，經(jīng)過簡單去噪處理后，對高聳柱狀結(jié)構(gòu)進行垂直度分析，并與工程測量法結(jié)果進行了對比分析。

上述針對柱狀結(jié)構(gòu)垂直度監(jiān)測分析中，核心都是獲取表面特征點，擬合特征截面圓心后計算其偏移量，針對原始數(shù)據(jù)僅進行簡單的去噪處理或人工剔除存在明顯偏差的測點，計算精度對原始測量點的誤差依賴較大。然而，由于現(xiàn)場測量條件復雜，加之塔筒結(jié)構(gòu)受自然環(huán)境的侵蝕、安裝工藝誤差，以及人為損壞或其他異物撞擊的影響所帶來的表面損傷、結(jié)構(gòu)不均勻等問題，使得原始測值序列存在粗差等異常測點，且數(shù)據(jù)偏離正態(tài)分布[12]，諸多存在誤差的測點在去噪處理或人工剔除時不易發(fā)現(xiàn)，最終導致計算結(jié)果失真。上述處理方法中均未考慮此情況，在擬合圓心時不具備消除此類誤差的能力。故甄別粗差觀測值并消除其影響，是提高分析精度的關(guān)鍵所在。

穩(wěn)健估計作為一種具有較強抗干擾性、可識別觀測值中偶然誤差和粗差的數(shù)據(jù)處理方法，已廣泛應用于諸多行業(yè)。陶葉青等[13]將穩(wěn)健估計應用于測繪數(shù)據(jù)的處理中，解決了我國參心坐標系中控制點已知坐標含粗差、精度差等問題。陳開端[14]采用穩(wěn)健總體最小二乘擬合，建立了地鐵施工監(jiān)測數(shù)據(jù)預測模型，該模型精度較高。李剛等[15]將穩(wěn)健估計用于火電機組設(shè)備的狀態(tài)評估模型中，解決了模型參數(shù)的偏差污染問題，參數(shù)具有較強穩(wěn)健性。馬洪磊等[16]在鐵路平面線形擬合中以穩(wěn)健估計為基礎(chǔ)，采用選權(quán)迭代法識別并剔除非擬合線形內(nèi)的粗差測點，顯著提高了擬合參數(shù)的可靠性。此外，穩(wěn)健估計理論與最小二乘法相結(jié)合在信息檢索[17]、流程監(jiān)控[18]、沉降預測[19]等多個領(lǐng)域得到廣泛應用。

基于此，本文將穩(wěn)健估計應用至風電塔筒的垂直度分析中，結(jié)合最小二乘圓擬合，通過Huber 選權(quán)迭代法，逐步消除測值序列存在的粗差、偶然誤差等異常測點的影響，并結(jié)合實例對其分析過程和效果進行闡述，提高了擬合參數(shù)及垂直度分析的可靠性。

1 風電塔筒垂直度分析

1.1 垂直度監(jiān)測特點與原理

垂直度指結(jié)構(gòu)中心線在不同高度處相應底部點的偏移現(xiàn)象，又稱為傾斜度[20]。風電機組塔筒是由多節(jié)塔段組成的高聳圓臺狀結(jié)構(gòu)，在安裝階段法蘭平面度細微偏差、接觸面清理不到位、螺栓預緊力不均勻等，運行階段基礎(chǔ)不均勻沉降、長期風載或其他人為因素等情況都會造成其垂直度超標。

風電塔筒垂直度的工程測量法包括極坐標法和偏心測量法。其中，偏心測量法通過觀測偏心點距離和角度得到圓心坐標，因其多余觀測量少，測量結(jié)果易受偏心點選取、測距及仰角大小的影響，精度較低[5,10]。故本文以極坐標法測得塔筒特征截面原始坐標。

極坐標法是以測站點為中心，后視點定向，測定塔筒特征截面坐標點的方法（圖1），其計算原理見式(1)，坐標測量與計算過程此處不再贅述。

圖1 極坐標法測量原理示意Fig.1 Schematic diagram of polar coordinate measurement principle

式中，XP、YP、HP分別為待測點坐標，XS、YS、HS分別為測站點坐標，DP為斜距，VP、αP分別為測站點Sn至待測點P的天頂距和方位角，is、ip分別為儀器高和棱鏡高。

將原始測量數(shù)據(jù)處理，擬合得到各特征截面圓心坐標(XTi,YTi,HTi)及(XBi,YBi,HBi)，按照式(2)進行風電塔筒垂直度計算（圖2）。

圖2 垂直度計算示意Fig.2 Schematic diagram of perpendicularity calculation

式中：i為風電塔筒垂直度（傾斜率），φ為傾斜角度，Δh、δ分別為兩特征截面高度差及其圓心位移偏移量。

1.2 塔筒垂直度國家標準

針對風電機組塔架的垂直度限值，目前國內(nèi)標準尚未統(tǒng)一，《風力發(fā)電機組塔架》（GB/T 19072—2010）及《風力發(fā)電機組裝配和安裝規(guī)范》（GB/T 19568—2017）規(guī)定在安裝完成后其垂直度應不大于0.1%，而在《風力發(fā)電機組驗收規(guī)范》（GB/T 20319—2017）規(guī)定機組驗收時塔架總體傾斜度不大于8 mm/m。隨著風電行業(yè)的發(fā)展和相關(guān)技術(shù)的進步，標準規(guī)范不斷完善，最新《高聳結(jié)構(gòu)設(shè)計標準》（GB 50135—2019）中新增了風力發(fā)電塔塔架傾斜率（垂直度）最大允許值為0.4%，根據(jù)筆者對甘肅及內(nèi)蒙古地區(qū)風電場超過500 臺風機塔筒垂直度的統(tǒng)計，認為其最大允許值定為0.4%較為合適。

1.3 誤差來源

在風電塔筒垂直度監(jiān)測中，原始特征點坐標測量存在的具有定量特性的系統(tǒng)誤差如儀器照準誤差、大氣折光差、溫濕度偏差等，可通過一定方法消除[20]。但由其他因素帶來的偶然誤差或粗差，無規(guī)律可循，其消除手段較為有限。如塔筒結(jié)構(gòu)在長期風力、雨水等環(huán)境的侵蝕下，出現(xiàn)表面保護層剝落導致表面凹凸不平；制造及安裝工藝問題導致法蘭交界面存在細微的臺階狀突變；塔筒受到人為損壞或其他異物撞擊，在特征截面上留下凹陷印記等現(xiàn)象，諸如此類由自然因素或人因素所帶來的塔筒表面損傷、結(jié)構(gòu)不均勻的部位會直接導致特征點坐標失真，從而帶來原始數(shù)據(jù)的粗差，最終導致特征截面擬合圓心的精度不穩(wěn)定。傳統(tǒng)最小二乘無法對原始坐標序列進行篩除，在擬合圓心時不具備消除此類誤差以及其他異常誤差的能力。

2 穩(wěn)健估計最小二乘圓擬合

2.1 最小二乘圓擬合原理

最小二乘圓擬合是以觀測數(shù)據(jù)為基礎(chǔ)采用最小二乘法擬合得到圓曲線方程，最小二乘法以觀測值殘差的平方和最小作為最優(yōu)估計準則，是迄今為止最主要、應用最廣泛的參數(shù)估計方法[20-21]。設(shè)B為自變量x的系數(shù)矩陣，β為參數(shù)矩陣，L為觀測值y的列向量，則誤差估值V的方程為

最小二乘法原理為

風電塔筒外壁測量在同高處符合圓形分布，為得到其特征截面圓心，本文采用最小二乘圓擬合法對一系列觀測數(shù)據(jù)進行預處理，根據(jù)圓參數(shù)方程構(gòu)建多變量的二次目標函數(shù)

式中：(a,b)為圓心坐標；r為圓半徑；(xi,yi)為所測特征截面圓弧采集點坐標；n為測點數(shù)量，測點數(shù)量即參與擬合的點數(shù)，為計算出參數(shù)a、b、r，測點數(shù)n應大于待估參數(shù)個數(shù)，即n≥3。

因目標函數(shù)u(a,b,r)中含有根號，計算極為復雜，令，則將此非線性的最小二乘問題轉(zhuǎn)換為式(6)求解(d2–r2)最小值的線性問題。

根據(jù)最小二乘法原理求相應參數(shù)，使得目標函數(shù)f(a,b,r)取得最小值，因平方和函數(shù)f(a,b,r)≥0，存在非負極小值，分別對目標函數(shù)式(6)的各參數(shù)求偏導，上述問題可轉(zhuǎn)化為

即可得到方程組

為求解式(8)，引入新參數(shù)c，令c=r2–a2–b2，代入簡化后，可將其優(yōu)化成自變量為(a,b,c)的一階線性方程組。

將各測量點(xi，yi)代入式(9)中，即可求得特征截面圓的擬合參數(shù)a、b、r。

2.2 穩(wěn)健估計理論

最小二乘法是基于原始觀測數(shù)據(jù)服從正態(tài)分布的前提，抵御粗差的能力有限。在實際觀測中，粗差及偶然誤差的出現(xiàn)無可避免，其概率最大約占觀測總數(shù)的10%[21]。原始觀測數(shù)據(jù)中存在此類誤差會使擬合參數(shù)失真，極大降低垂直度監(jiān)測精度。傳統(tǒng)最小二乘法用統(tǒng)計及殘差檢測方法來剔除此類誤差，具有很大局限性。針對最小二乘抗差性不足的缺陷，1953年GE.P.BOX 首先提出了Robustness概念，后經(jīng)過Huber 等人[22-23]的研究發(fā)展，奠定了穩(wěn)健估計理論的基礎(chǔ)。穩(wěn)健估計理論的原則是充分利用原始數(shù)據(jù)的有效信息，限制和排除有害信息，盡可能減弱實際觀測中存在的粗差和偶然誤差對參數(shù)估值的影響，進而得到最優(yōu)參數(shù)估值。

穩(wěn)健估計理論主要分為M 估計、L 估計和R 估計3 種類型，本文采用最具實用價值、應用最為廣泛的M 估計法[23]，即極大似然估計準則。為估計參數(shù)矩陣β，進行多次觀測，得到觀測值列向量L，f為隨機向量L的密度函數(shù)，則有

對式(10)求導，可得

由式(10)—式(11)求出參數(shù)矩陣β的估值，即為M 估計。M 估計的方法有多種，本文選取易于編程實現(xiàn)、抗粗差能力強、應用最廣泛的選權(quán)迭代法[24]，其誤差方程與式(3)一致，權(quán)函數(shù)矩陣為式(12)，各權(quán)因子初始值為1，則

限制條件為

則法方程為

求解式(12)—式(14)，得到待求參數(shù)和殘差的第一次估值：

由初值V(1)確定各觀測值新的權(quán)因子，實現(xiàn)粗差定位，構(gòu)成新的等價權(quán)矩陣，再次解算算法方程式(14)，如此迭代直至滿足終止條件，得到最終的參數(shù)估值為：

整個迭代過程中，含有粗差的測值的權(quán)值會逐漸減小，通過各測點殘差值的大小可直觀了解粗差等異常誤差的大小及其測點位置。

2.3 Huber 迭代權(quán)函數(shù)

極大似然估計中的選權(quán)迭代法，其迭代權(quán)函數(shù)有多種選擇，本文選用相應的Huber 迭代權(quán)函數(shù)，選取不同ρ函數(shù)，其φ函數(shù)及相應權(quán)函數(shù)也隨之變化，且在迭代過程中，權(quán)函數(shù)也隨改正數(shù)的更新而變化，根據(jù)Huber 法，ρ函數(shù)和φ函數(shù)[19,21]為：

式中：v為誤差函數(shù)，c為調(diào)和系數(shù)，按照“2σ準則”一般取2σ，對應的權(quán)函數(shù)為

3 穩(wěn)健最小二乘法計算塔筒垂直度

測得塔筒各特征截面測點原始坐標后，用所有測點參與擬合，采用最小二乘法計算各截面的圓心坐標和半徑的擬合初值。由于各截面測點沿圓周均勻布置，高程坐標zi基本處于同一高程，為簡化計算，在擬合過程中僅采用平面坐標(xi,yi)得到圓心平面坐標(a,b)，高程坐標h由各擬合點的高程均值代替。取得圓心坐標及半徑的初值后，計算各測點至擬合圓距離即擬合誤差di，以此誤差初值確定各觀測值新的權(quán)因子，構(gòu)成新的權(quán)因子矩陣，如此迭代直至滿足迭代終止條件，即可實現(xiàn)逐步弱化和剔除粗差等異常誤差測點的影響，得到各截面最終擬合圓心(a,b,h)及半徑r的穩(wěn)健結(jié)果。圖3 為穩(wěn)健最小二乘法計算風電塔筒垂直度算法流程，該算法采用Excel 內(nèi)嵌VBA 語言編程實現(xiàn)。

圖3 穩(wěn)健最小二乘法計算風電塔筒垂直度算法流程Fig.3 Flow of least-square in perpendicularity calculation of wind turbine tower

穩(wěn)健估計選權(quán)迭代法的迭代終止條件一般以待估計參數(shù)平差值、觀測值殘差或單位權(quán)中誤差作為判斷基礎(chǔ)。本文在垂直度分析算法流程中以相鄰2 次待估參數(shù)平差值作為終止迭代條件，認為相鄰2 次參數(shù)估計值差異小于10–4便可停止迭代。

若原始觀測數(shù)值整體誤差較大，雖經(jīng)多次迭代后滿足終止條件。但此時觀測點權(quán)值接近0（小于1×10–3）的數(shù)目較多，有效測點數(shù)不足，最終估計值不可信。鑒于此，在分析算法流程中加入對整體測點中誤差及有效點數(shù)的限制條件。當滿足終止迭代條件后，考察整體測點中誤差是否滿足《工程測量規(guī)范》（GB 50026—2007）中關(guān)于高聳構(gòu)筑物結(jié)構(gòu)變形測量精度等級及相應測量點位中誤差要求，且有效測點數(shù)目是否大于20 個，滿足條件后方可計算最終風電塔筒垂直度。

4 應用實例

4.1 研究對象概況

以甘肅某風電場三期工程為研究對象，該工程安裝風電機組單機容量為2.5 MW，采用圓柱形鋼塔筒，由基礎(chǔ)環(huán)和5 節(jié)塔筒組成，輪轂高度120 m，葉片長度70 m，風輪直徑143 m，各節(jié)塔筒間的連接采用高強度螺栓連接。風電場地處典型的黃土高原地區(qū)，機組沿黃土梁頂布置，地質(zhì)較為疏松，部分區(qū)域水土流失嚴重，垂直度監(jiān)測尤為重要。

4.2 實例分析與驗證

采用全站儀均勻測得塔筒底部、中部及頂部特征截面圓周各點坐標，各截面測點數(shù)量不少于40 個。采用上述穩(wěn)健最小二乘法對原始數(shù)據(jù)進行處理，以其中某一臺風機為例。為直觀描述特征截面在穩(wěn)健最小二乘法下的迭代過程與收斂情況，選擇該風機的頂部截面，將各次迭代擬合得到的圓心坐標統(tǒng)一繪制在同一坐標系下（圖4），圖4 中序號為迭代次數(shù)，每次迭代得到的截面半徑及其變化過程如圖5所示。

圖4 截面圓心擬合過程Fig.4 The fitting process of section center

圖5 擬合半徑變化過程Fig.5 The change process of fitting radius

由圖4 可知，空心圓點為迭代終止前圓心擬合結(jié)果，經(jīng)過21 次迭代后滿足迭代終止條件，最終圓心位置如實心圓點所示，坐標為（907.835 7，664.857 5），可見穩(wěn)健估計過程中，圓心擬合結(jié)果在一定范圍內(nèi)不斷變化，隨著迭代次數(shù)增加，逐漸向最終圓心位置處靠攏。圖5 中截面半徑的擬合變化過程也同樣反映此規(guī)律，迭代初期受異常測點影響，擬合半徑在終值附近變動較大，隨著迭代次數(shù)增加，半徑的擬合結(jié)果逐漸趨于穩(wěn)定。底部和中部截面與此類似。

為描述穩(wěn)健估計最小二乘法在塔筒擬合過程中的精確性，繪制迭代過程中測點截面圓心擬合偏差均值與中誤差的變化（圖6），測點擬合最大偏差及平均偏差對比如圖7所示。

圖6 擬合偏差均值與中誤差變化Fig.6 Changes of the fitting deviation mean and its mean erro

由圖6—圖7 可知，在整個迭代過程中，測點擬合平均偏差及擬合中誤差呈均勻下降趨勢，平均偏差由2.1 mm 減小至0.5 mm，減小76%；中誤差由2.6 mm 減小至0.6 mm，減小77%。測點擬合最大偏差由6.9 mm 降低至1.2 mm，降低83%，最大偏差在迭代初期降低速率較高，后續(xù)降低速率逐漸減小。總之，隨著迭代次數(shù)的增加，測點擬合最大偏差、平均偏差及擬合中誤差都逐漸減小，最后收斂穩(wěn)定并達到迭代終止條件，得到誤差更小、更為精確的圓心坐標及半徑的最終擬合值，提高了截面圓心擬合精度。

圖7 最大偏差及平均偏差變化Fig.7 Variations of the maximum and average deviation

利用傳統(tǒng)最小二乘法與穩(wěn)健最小二乘法分別計算該風電塔筒垂直度，結(jié)果見表1。

表1 塔筒垂直度計算結(jié)果對比Tab.1 Result contrast of the perpendicularity calculation for the tower

由表1 可知，采用穩(wěn)健最小二乘法進行數(shù)據(jù)處理后塔筒擬合圓心坐標及半徑都有不同程度的差異，由于底部存在粗差點較多，2 種方法坐標擬合結(jié)果相對較大。對比塔筒整體垂直度（頂部相對于底部）計算結(jié)果，穩(wěn)健最小二乘法所得偏移量及垂直度與傳統(tǒng)最小二乘法的相差約19%，可見，穩(wěn)健最小二乘法對截面圓心的擬合精度有較大提高。

4.3 權(quán)函數(shù)簡化與粗差點識別

根據(jù)迭代權(quán)函數(shù)式(19)，經(jīng)過多次迭代后，含有粗差或偶然誤差的奇異測值的權(quán)值逐漸減小，直至趨近于0，則可將權(quán)函數(shù)簡化為

式中，di為各觀測測點的殘差值，即測點距擬合圓的距離。

以權(quán)函數(shù)式(20)進行迭代，在迭代過程中權(quán)值為0 的點即為粗差點或奇異測點，即可識別并剔除誤差點，剔除時采用最大殘差逐步剔除，即每次迭代后在di＞2σ的所有測點中，剔除di最大的測點后繼續(xù)迭代直至滿足終止條件。

對比簡化權(quán)函數(shù)計算結(jié)果，進一步說明穩(wěn)健最小二乘在識別和剔除粗差測值的優(yōu)越性。以上述風機頂部截面為例，統(tǒng)計有效測點（權(quán)值小于1×10–3的測點）數(shù)目為24 個，在此24 個有效測點中隨機選取3 個點加入1～2 dm 的粗差。采用簡化權(quán)函數(shù)及殘差逐步剔除法，經(jīng)過4 次迭代，逐步將粗差測點賦權(quán)值為0，各次迭代結(jié)果見表2。由表2 可知，權(quán)函數(shù)簡化后可快捷有效識別并剔除粗差測點，對最終計算結(jié)果無影響，也進一步說明穩(wěn)健最小二乘圓擬合法基本不受粗差、偶然誤差等異常測值的影響。

表2 簡化權(quán)函數(shù)各次迭代結(jié)果Tab.2 The results of each iteration of the simplified weight function

5 結(jié)語

穩(wěn)健最小二乘擬合法在截面圓心坐標和半徑的擬合過程中能有效避免粗差、偶然誤差等異常測點的影響，隨著算法迭代次數(shù)的增加擬合結(jié)果逐漸趨于穩(wěn)定，測點擬合最大偏差、平均偏差及擬合中誤差都逐漸減小，有效提高了擬合參數(shù)及塔筒垂直度的計算精度。

簡化權(quán)函數(shù)可快速識別并剔除原始測值序列中存在粗差、偶然誤差等異常誤差的測點，由自然或人為因素所帶來的塔筒表面損傷、不均勻部位所導致的個別測點坐標失真的情況對最終計算結(jié)果無影響，基于穩(wěn)健最小二乘法的垂直度分析具有較強的抗差能力。

本文方法易于編程實現(xiàn)，可推廣應用于火電廠冷卻塔、水塔、煙囪等其他類似高聳塔筒結(jié)構(gòu)的垂直度分析中，對處理相關(guān)監(jiān)測數(shù)據(jù)，提高監(jiān)測精度具有一定參考價值。