基于無人機航拍的近地風場類別確定方法

李正農(nóng) ，代向奎，蒲鷗，吳紅華

［建筑安全與節(jié)能教育部重點實驗室（湖南大學），湖南 長沙 410082］

目前，在進行結構抗風設計時，為了確定風荷載，首先要確定其風場類別，傳統(tǒng)的風場類別確定方法主要依據(jù)《建筑結構荷載規(guī)范》（GB 50009—2012）［1］規(guī)定的四類地貌.但李正農(nóng)等［2］研究發(fā)現(xiàn)，高層建筑受到的風荷載是受周邊環(huán)境影響的，該規(guī)范中不同地貌的風場類別判定無法完全涵蓋一個地區(qū)的所有位置，因此，對于城市地區(qū)，這種分類方法較為粗糙且界定標準不夠細化.空中風場受地貌影響會隨高度呈現(xiàn)一定的變化規(guī)律，我國通常采用指數(shù)律來擬合風速剖面并確定其風場類別，常用的測量方法包括測風塔［3］、風廓線雷達［4］和移動激光雷達［5］等.測風塔安裝難度大、易受場地限制，在城市地區(qū)進行多點測量時難度很大［6］；移動雷達空間分辨率較低且價格昂貴［7］.因此城市風場難以采用傳統(tǒng)方法進行測量.吳紅華等［8］和李正農(nóng)等［9-10］提出了采用兩臺搭載風速儀的無人機同步測量的方法進行近地風場的測量并識別其風場類別，該方法具有較高的空間靈活性且成本較低，但由于一些風向很少出現(xiàn)，無法長期進行現(xiàn)場實測，該方法無法判定所有方向的風場類別.

鑒于通過空中風場測量確定風場類別的局限性，地貌識別或?qū)⒊蔀楦泳_且應用廣泛的確定風場類別的新方法.識別地面地貌狀況的方法大致分為兩種，第一種是通過空中風場測量，根據(jù)其變化規(guī)律來進行地貌狀況的判斷，第二種是根據(jù)實際地面粗糙度進行地貌狀況的判斷.在地面粗糙度識別方面，張相庭［11-12］利用風耗散能原理對各地貌進行分析，并通過建立α值與風耗散能η的關系實現(xiàn)了按照實際地貌計算地面粗糙度指數(shù)α值，但隨著大量新建建筑的出現(xiàn)，該方法的精確度及適用性都大大降低.Alam 等［13］利用衛(wèi)星遙感分析技術對衛(wèi)星圖像進行波段處理，根據(jù)其反射率確定土地覆蓋類型，然后結合空間算法形成粗糙度等級圖，實現(xiàn)了無須現(xiàn)場測量的地面粗糙度分類方法，但該方法只進行了幾千平方千米區(qū)域的分類測試，在城市局部區(qū)域中的分類精確度并未進行研究.城市地區(qū)地面粗糙度主要由建筑物形成，研究者們根據(jù)建筑物高度和建筑物密度等因素提出了地面粗糙度的計算公式［14-16］，但城市地區(qū)建筑物眾多，粗糙度識別較為復雜.近年來通過無人機影像進行地面建筑物信息識別的技術得到廣泛應用，Vacca 等［17］和Koc 等［18］對基于無人機航拍得到的建筑物精度進行了研究，結果表明通過設置合適的航拍參數(shù)可以獲得較高精度的建筑物尺寸信息；Subrahmanya 等［19］根據(jù)無人機圖像結合人工神經(jīng)網(wǎng)絡算法提供了準確的建筑物識別信息；董彪等［20］和Gao 等［21］均基于Yolov3 算法進行改進，在遙感圖像的建筑物檢測中展現(xiàn)出了較高的檢測精度.基于無人機航拍可以獲取城市地區(qū)地面建筑物信息，但現(xiàn)階段尚未建立有效的依據(jù)實際地貌狀況進行城市地區(qū)風場類別識別的方法.

本文應用無人機航拍城市地區(qū)地面照片，結合圖像識別技術獲取建筑物信息并對地貌狀況進行分析，通過計算出地面粗糙度Z0并利用風剖面對比獲得與Z0對應的粗糙度指數(shù)α，進而建立一種應用無人機航拍地貌狀況獲取測區(qū)地面粗糙度指數(shù)α并判定其風場類別的方法.本文所提出的方法具有較強的普適性，實現(xiàn)了依據(jù)實際地貌判定城市地區(qū)風場類別，其結論可為城市地區(qū)建筑抗風設計提供參考.

1 實測概況

1.1 實測設備

如圖1（a）所示，利用兩臺搭載風速儀的六旋翼無人機進行來流風場的實測，一臺無人機懸停在100 m 高度處作為參考，另一臺無人機在10～120 m高度范圍內(nèi)進行測量，本文采用10 個高度進行風剖面計算.風速儀安裝在無人機機身上部20.20 cm 處（0.53 倍旋翼直徑）以減小機翼轉(zhuǎn)動產(chǎn)生的不利影響.經(jīng)改裝后風速儀由無人機直接供電并將所測數(shù)據(jù)通過與其相連接的無線電臺實時傳輸?shù)接嬎銠C端.本次實測前已在風洞實驗室對風速儀進行標定并進行誤差修正，通過該方式測量風場的可行性已由Li等［22］進行了試驗驗證.

圖1 實測所用系統(tǒng)Fig.1 System used in the measurement

如圖1（b）所示，采用P4 multispectral無人機對測區(qū)內(nèi)各建筑物進行拍攝，本文分別在100 m 和120 m高度處對各建筑物依次進行正射航拍，航拍均選擇在光照條件較好時進行以避免因照片清晰度不足而影響建筑物信息提取，通過對照片的識別獲得各建筑物尺寸信息，進而獲取測區(qū)內(nèi)的地貌狀況.

1.2 場地概況

本試驗選定兩處場地分別進行風場實測和地面粗糙度計算，如圖2 所示，實測場地均為建筑物較為密集的城市地區(qū).測區(qū)一中以五角星標為目標區(qū)域判斷來流方向（箭頭所示）的地貌類別，該測區(qū)來流方向有較多形狀規(guī)則且相似的建筑物.測區(qū)二中以圓標為目標區(qū)域判斷來流方向（箭頭所示）的地貌類別，測區(qū)二來流方向有較多的均勻小區(qū)和工廠，部分建筑物由于地圖未及時更新，因此在圖2（b）中未顯示.按照《建筑結構荷載規(guī)范》（GB 50009—2012）［1］（以下簡稱規(guī)范）中關于地面粗糙度的分類標準可將兩處測區(qū)大致判斷為C 類地貌，地面粗糙度指數(shù)為0.22左右.

圖2 實測場地Fig.2 Test site

實測期間風場條件比較理想，選取風速較大且風向穩(wěn)定（變化幅度集中在7°以內(nèi)）的時段進行來流方向風場測量.在地面粗糙度計算過程中僅考慮建筑物為主要粗糙元，并且考慮到較遠處建筑物對目標建筑的影響較小，對測區(qū)一來流方向約800 m 范圍內(nèi)、測區(qū)二來流方向約1 000 m 范圍內(nèi)的建筑物進行了拍攝及計算，具體范圍如圖2所示.

2 實測數(shù)據(jù)前處理與風場類別的判定

2.1 實測數(shù)據(jù)前處理

在實際測量過程中，由于無人機自身的抖動以及機身傾斜等因素會對所采集的風場數(shù)據(jù)有一定的影響，結合風速儀在風洞實驗室中的標定結果，按照公式（1）對無人機所采集的數(shù)據(jù)進行修正，不同風速區(qū)間的修正系數(shù)見表1.

表1 無人機機身修正系數(shù)Tab.1 UVA fuselage correction coefficient

式中：u（t）為修正后的無人機風速時程數(shù)據(jù)；uw（t）為無人機原始風速時程數(shù)據(jù)；λ為修正系數(shù).

2.1.1 平均風速和平均風向角

現(xiàn)場實測風速儀直接采集的原始數(shù)據(jù)為風速、風向時程，按照公式（2）進行分解可得到x、y方向的風速分量（圖3），各工況下的平均風速U和平均風向角φ根據(jù)公式（3）～公式（6）計算.

圖3 風速時程分解示意圖Fig.3 Wind speed time-history decomposition

式中：ux（t）為E 方向分解風速時程；uy（t）為N 方向分解風速時程；u（t）為風速時程；θ（t）為風向時程.

式中：n為測量時距內(nèi)的采集數(shù)目為x方向分解風速平均值為y方向分解風速平均值；U、φ分別為各工況下測量時距內(nèi)的平均風速和平均風向角.

2.1.2 地面粗糙度指數(shù)α1實測原理

由于地表粗糙程度的影響，在大氣邊界層內(nèi)風速隨高度發(fā)生變化，根據(jù)規(guī)范有關內(nèi)容可按照指數(shù)律［公式（7）］進行風剖面擬合.本試驗采用兩臺搭載風速儀的六旋翼無人機同步測量，無人機B 分別測量10～120 m 范圍內(nèi)不同高度處的風速，無人機A 懸停在100 m 高度作為參考并獲得與無人機B 對應的10 個樣本數(shù)據(jù)，用參考文獻［10］的方法對無人機所測風場數(shù)據(jù)進行歸一化，并采用公式（8）擬合得到地面粗糙度指數(shù)α1.

式中：Z為離地高度；UZ為Z高度處的風速；U10為參考高度10 m 處的風速；UA-i和UB-i分別為機身修正后無人機A和無人機B在各工況下測量的平均風速，i=1，2，…，10；ZB-i為無人機B 在各工況下的對應高度；ZA為無人機A懸停參考高度，本文為100 m.

2.1.3 湍流度

湍流度是反映脈動風速相對強度的一種度量指標，通常采用平均時距內(nèi)脈動風速均方根與平均風速的比值來表示，即：

在我國規(guī)范中，湍流度沿高度的分布按照公式（11）進行計算：

式中：I10為10 m 參考高度處的名義湍流度，對應于A、B、C 和D 類地面粗糙度，可分別取0.12、0.14、0.23、0.39.

2.2 無人機航拍識別地面粗糙度

本研究以城市區(qū)域建筑物作為主要粗糙元進行地面粗糙度的計算，由于建筑物成像大小受航拍高度影響，本文采用兩個高度分別進行拍攝以實現(xiàn)建筑物尺寸的計算，其中航拍高度計算方法［23］如式（12）：

式中：H為航拍高度，m；f為鏡頭焦距，mm；GSD為地面采樣距離，cm/pixel，其實質(zhì)為用地面距離單位表示像素大小；a為像元尺寸，即圖像中單位像素點的尺寸（影像傳感器寬度/圖像寬度）.

2.2.1 相機標定

攝影成像的過程是將三維地物投影到二維平面上，即坐標轉(zhuǎn)換的過程，相機標定主要有兩個作用，一個是根據(jù)坐標轉(zhuǎn)換過程計算出相機的內(nèi)、外參數(shù)［如公式（13）］，另一個是獲取鏡頭畸變參數(shù)進行圖像校正.理想情況下由三維物理坐標系中一點（Xw，Yw，Zw）到二維像素坐標系中（m，n）的轉(zhuǎn)換公式如式（13）：

式中：Zc為物體與光心之間的距離；fx=f/dx，fy=f/dy，f為焦距；（m0，n0）為圖像的中心位置，該矩陣為相機的內(nèi)參矩陣；R為3×3的旋轉(zhuǎn)矩陣，T為3×1的平移向量，該矩陣為相機的外參矩陣.

在實際成像過程中機身姿態(tài)、光學透鏡固有的透視失真等因素會導致相機成像產(chǎn)生一定程度的變形，本文按照Zhang［24］標定方法采用Matlab對無人機相機進行標定.如圖4（a）所示，對標準7 行9 列形式的25 mm×25 mm正方形網(wǎng)格標定板進行不同位姿的拍攝（本文采用20 張照片進行標定），根據(jù)各網(wǎng)格角點的空間坐標和對應的像點坐標之間的變換關系進行相機內(nèi)參、外參及畸變參數(shù)的計算.圖4（b）給出了相機標定的重投影誤差（各網(wǎng)格角點的實際投影像點與根據(jù)相機標定結果計算出的理論投影像點之間的距離）直方圖，通常認為誤差在0.3 以內(nèi)即滿足精度要求，本次標定結果的平均誤差為0.14 Pixel 且最大誤差為0.16 Pixel，因此本次標定結果有效.

圖4 相機標定Fig.4 Camera calibration

通常將鏡頭畸變分為徑向畸變和切向畸變，徑向畸變是成像過程中的主要畸變，它是以畸變中心為中心點，沿著鏡頭半徑方向產(chǎn)生的位置偏差；切向畸變是安裝過程中鏡頭本身與圖像平面不平行而產(chǎn)生的誤差.消除畸變誤差的計算公式如式（14）：

式中：m、n為原始圖像上的坐標位置；x、y為校正后的坐標位置；r為點到成像中心的距離；k1、k2、p1、p2和k3為5個畸變參數(shù)，通過相機標定得到.

2.2.2 建筑物識別定位

城市地區(qū)建筑物眾多且特征較為復雜，傳統(tǒng)的目標檢測方法無法滿足工程應用需求，本文采用基于端到端的Yolo 算法［25］進行建筑物檢測.Yolo 系列算法在不斷改進與發(fā)展［26-27］，Yolov5 是目前較為成熟的識別算法，其通過回歸能夠快速檢測出圖像中物體的類別和邊界框坐標，在檢測速度和精度方面均有很好的效果且在工業(yè)領域被廣泛使用.

Yolov5的網(wǎng)絡結構分為Input、Backbone、Neck和Head 四部分，相較于之前的版本進行了一定的優(yōu)化：在輸入端采用Mosaic 數(shù)據(jù)增強等方式來提高小目標的檢測效果；在Backbone 階段引入Focus 結構和CSP 結構提高網(wǎng)絡的學習能力；在Neck 階段引入FPN 和PAN 結構來增強網(wǎng)絡特征融合的能力；在輸出端采用CIOU loss 作為預測框的損失函數(shù).在檢測過程中Yolov5將輸入圖片劃分成三種不同尺度的網(wǎng)格，當物體中心點落在網(wǎng)格內(nèi)部時，由該網(wǎng)格負責預測該物體，每個網(wǎng)格都會預測出3 個邊界框和其所包含物體的類別概率，每個邊界框都包含其坐標信息和置信度.置信度是邊界框包含物體的概率與邊界框定位準確性的乘積［如公式（16）］，將邊界框的置信度與其所屬網(wǎng)格預測出的類別概率相乘可得到該邊界框包含某類物體的概率，通過閾值過濾和非極大值抑制篩選出置信度評分最高的邊界框作為輸出結果.

式中：c為邊界框的置信度；Pr（Object）代表邊框包含物體的概率；P表示預測邊界框；G 表示人工標注框；IOU越接近于1代表定位越準確.

預測框輸出結果為預測物體的置信度以及預測框的坐標，根據(jù)預測框?qū)φ鎸嵖虻念A測誤差構造損失函數(shù)，并通過梯度下降的方式不斷優(yōu)化獲得網(wǎng)絡結構在各階段的最優(yōu)權重參數(shù).Yolov5 算法的損失函數(shù)包含預測框定位損失、置信度損失和分類損失三部分，其中置信度損失和分類損失與之前版本一致都采用交叉熵損失［26］，預測框定位損失采用CIOU loss［公式（18）］，綜合考慮了預測框和標注框的重疊面積、中心點距離和長寬比等因素.

式中：Lobj為Yolov5 算法損失函數(shù)；Lbox為預測框定位損失；Lconf為置信度損失；Lclass為分類損失；λ1、λ2和λ3為權重系數(shù).

式中：b、h和w分別為預測框的中心點位置、高度和寬度；上標gt（ground truth）表示人工標注框；IOU 為預測框和標注框的交互比；ρ為預測框和標注框的中心點距離；d為預測框和標注框的最小外接框的對角線長度；s為預測框和標注框的寬高比的相似度.

本文采用Yolov5 算法識別圖像中建筑物，從TorontoCity 數(shù)據(jù)集和中國典型城市建筑物實例數(shù)據(jù)集［28］以及本研究在多個城市區(qū)域所拍攝照片中挑選出2 000 張圖像，共包含建筑物近10 000 座，并使用Labelimg 標注工具對圖像中建筑物的矩形邊界進行標注（圖5），制作本文所需的數(shù)據(jù)集.為得到適合建筑物檢測的先驗框維度，在訓練之前使用k-means算法對標注框進行聚類分析，本文得到的先驗框尺寸為59×42、45×84、116×80、89×165、173×113、152×260、257×176、412×268 和269×440，并將數(shù)據(jù)集中的圖像按照7∶2∶1 的比例分別作為訓練集、測試集和驗證集.本文訓練過程中初始學習率設置為0.001，采用Adam 優(yōu)化器，動量為0.937，置信度閾值為0.5.圖6 給出了本文模型訓練階段的損失圖，可看出該模型表現(xiàn)出較好的收斂性，訓練迭代次數(shù)為300 時最終損失值穩(wěn)定在0.13 左右，此時平均準確率MAP達到80%.

圖5 標注界面Fig.5 Dimension interface

圖6 損失函數(shù)曲線Fig.6 Loss function curve

根據(jù)訓練得到的權重參數(shù)對校正后的圖像進行預測，預測框的輸出信息包括物體類別、置信度和預測框定位參數(shù)（中心點坐標、寬度和高度），其中一組輸出結果如圖7所示.

圖7 建筑物識別Fig.7 Building identification

2.2.3 建筑物尺寸計算

本文采用兩個高度分別對建筑物進行拍攝，圖8給出了無人機航拍獲得建筑物正射影像的示意圖，建筑頂面垂直于光軸且與像平面平行，以2.2.2 節(jié)中輸出的矩形預測框作為建筑物的平面邊界，由于建筑物大多為規(guī)則體且頂面為矩形，故以長方體作為建筑物的外邊界進行尺寸計算.圖8 中根據(jù)同一建筑物在兩個高度所拍攝照片中得到的建筑物定位框與建筑物實際邊界尺寸之間的幾何相似關系可得到各建筑物的尺寸信息，具體計算公式如式（19）：

圖8 建筑物測量Fig.8 Building measurement

式（19）中及圖8中，h、L和W分別為建筑物的實際高度、長度和寬度；P1、P2和P′1、P′2分別為像平面1、像平面2 中建筑物定位框的長度和寬度，單位為pixel；D1、D2分別為像平面1 和像平面2 與建筑物頂面之間的距離；O1、O2分別為兩個高度處的鏡頭中心；a為像元大小，本文中為 3 μm×3 μm；f為焦距，根據(jù)2.2.1節(jié)標定結果為5.83 mm；H1和H2為無人機離地高度，本文采用100 m和120 m.

經(jīng)驗證本文得到的建筑物尺寸誤差控制在5%以內(nèi)，可滿足本研究所需精度要求.本文中建筑物識別及尺寸計算具體流程如圖9所示.

圖9 建筑物識別流程Fig.9 Building recognition process

2.2.4 地面粗糙度Z0

對數(shù)律風剖面公式中，平均風速U可表示為地面以上高度Z的函數(shù)，若已知ZS高度處風速值US，其一般表達式［29］為：

式中：U為Z高度處的風速（Z為離地高度，本文僅對20 m 以上進行分析）；Z0為風速等于0 的高度，隨地面粗糙程度而變化，稱為地面粗糙度.

關于地面粗糙度Z0因?qū)Υ植谠目紤]因素不同而有不同的估算方法［30］，在本文中只考慮建筑物為粗糙元，主要影響因素為來流方向建筑物的高度、密度等，并且考慮到各個粗糙元的變化，本文采用Kondo等［15］給出的方法進行Z0的計算：

式中：Si代表粗糙元i占有的實際水平面積；hi代表該粗糙元的高度；S為計算范圍內(nèi)的總面積.

3 風場實測數(shù)據(jù)及地貌狀況的分析

3.1 測區(qū)一分析

3.1.1 風場實測α1

本次試驗采用兩臺搭載風速儀的無人機進行來流風場測量，無人機A 懸停在100 m 高度處作為參考，無人機B 在10～120 m 范圍內(nèi)依次測量10 個高度處的風速，測量過程中各工況下兩風速儀所采集風向角誤差均在3°以內(nèi)，故可直接進行比較.無人機所測風場數(shù)據(jù)經(jīng)機身校正、歸一化處理后按照公式（8）和公式（11）分別進行風剖面和湍流度剖面擬合，結果如圖10所示.

圖10 測區(qū)一風場擬合Fig.10 Fitting of wind field in test site 1

由圖10（a）可知，測區(qū)一中由風剖面擬合所的實測地面粗糙度指數(shù)α1為0.204 8，擬合優(yōu)度R2為0.968，擬合曲線介于B、C類地貌之間且與C類地貌吻合程度更好.圖10（b）給出了測區(qū)一的實測湍流度數(shù)據(jù)及規(guī)范中B類、C類地貌的湍流度剖面，可知實測湍流度與規(guī)范中C類地貌更貼近，同風剖面擬合結果較一致.

3.1.2 地面粗糙度Z0計算

測區(qū)一總面積為S1=520 104.152 6 m2，所選范圍內(nèi)共計38 處建筑物，根據(jù)2.2 節(jié)中方法識別圖像中的建筑物，對各建筑物進行信息提取如圖11 所示.根據(jù)公式（21）計算可得測區(qū)一的平均地面粗糙度為：

圖11 測區(qū)一建筑物信息Fig.11 Building information in test site 1

3.1.3 由地面粗糙度Z0計算所得α2

在各種地貌的梯度風高度以上區(qū)域，風的流動已不受地面粗糙度的影響，各處風速均屬相同，均為梯度風速［29］.在同一梯度風速下分別按照公式（7）和公式（20）繪制風剖面，并采用歐式距離公式［31］判斷出吻合程度最高的α2值.具體計算公式如下式（22）～（24）：

式中：v0為測區(qū)一所在地區(qū)的基本風速；vG為測區(qū)一所在地區(qū)的梯度風速值；U1、U2分別為根據(jù)指數(shù)律和對數(shù)律得到的不同高度處風速值；ZG為規(guī)范中定義的與α2相對應的梯度風高度，并根據(jù)不同的α2值按比例確定；d為歐式距離指標.

根據(jù)規(guī)范［1］中《全國基本風壓標準值表》查得測區(qū)一基本風壓為0.75 kN/m2，計算可知該地區(qū)基本風速為v0=34.64 m/s（考慮到城市地區(qū)下墊面粗糙度較大，一定高度以下風速隨高度的變化比較紊亂，不一定符合對數(shù)規(guī)律或指數(shù)規(guī)律［32］，故本文從20 m 高度處開始，并以200 m 高度范圍為例進行計算）.圖12為α2在0.150～0.250（精確度為0.001）范圍內(nèi)對應的d值變化曲線，對應于α2=0.198 時d取得最小值.圖13為d取得最小值時兩指標對應的風剖面對比圖.

圖12 測區(qū)一d值變化曲線Fig.12 Change curve of d value in test site 1

圖13 測區(qū)一風剖面對比圖Fig.13 Comparison of wind profiles in test site 1

3.1.4 誤差分析

根據(jù)地面粗糙度計算所得α2值與風場實測所得α1值比較接近，圖14 給出了兩種方法在同一梯度風速的風剖面誤差圖（近地面200 m 范圍內(nèi)）.由圖14可知，這兩種方法得到的風剖面曲線較為貼近，相對誤差為3.32%.分析可知測區(qū)一的風場類別介于B、C類之間，考慮到城市地區(qū)個別建筑物高度超過200 m，分別對近地面350 m（B類地貌梯度高度）和450 m（C 類地貌梯度高度）高度范圍進行計算及誤差分析，結果如表2所示.

表2 測區(qū)一不同高度誤差分析表Tab.2 Error analysis of different heights in test site 1

圖14 測區(qū)一風剖面誤差圖Fig.14 Wind profile error in test site 1

由表2 可知，在B 類地貌和C 類地貌梯度高度范圍內(nèi)根據(jù)地面粗糙度計算出的α2值均為0.196，與風場實測結果之間的相對誤差為4.30%.同200 m 高度范圍計算結果對比可知，地面粗糙程度對近地面一定高度范圍內(nèi)的風剖面變化規(guī)律有較顯著影響，采用的計算高度范圍不同計算結果略有差異，且計算高度達到一定值時其結果不再改變.

3.2 測區(qū)二分析

3.2.1 風場實測α1

測區(qū)二風場實測數(shù)據(jù)經(jīng)機身修正、歸一化處理后，按照指數(shù)律進行地面粗糙度指數(shù)α的擬合計算，結果如圖15所示.

圖15 測區(qū)二風場擬合Fig.15 Fitting of wind field in test site 2

由圖15（a）可知，測區(qū)二中實測風剖面擬合指數(shù)α1為0.210 4，擬合優(yōu)度為0.973 6，擬合曲線與C 類地貌吻合程度較高.根據(jù)圖15（b）可知測區(qū)二的實測湍流度數(shù)據(jù)與規(guī)范中C 類地貌非常貼近，同風剖面擬合結果較一致.

3.2.2 地面粗糙度Z0計算

測區(qū)二總面積為S2=1 039 631.443 2 m2，所選范圍內(nèi)共計86 處建筑物，對各建筑物進行識別和信息提取如圖16 所示.由于該測區(qū)存在較多的均勻小區(qū)和工廠，其外形相似性較高，故存在部分重疊區(qū)域.根據(jù)公式（21）計算可得測區(qū)二的平均地面粗糙度為：

圖16 測區(qū)二建筑物信息Fig.16 Building information in test site 2

3.2.3 由地面粗糙度Z0計算所得α2

由公式（23）～公式（25）可知，不同的梯度風速值不影響d值隨α2的變化規(guī)律，即不影響吻合程度最高的α2值判定結果，故本節(jié)以梯度風速vG=20 m/s 為例進行計算.

圖17 為同一梯度風速下α2在0.150～0.250（精確度為0.001）范圍內(nèi)對應的d值變化曲線，對應于α2=0.206 時d取得最小值.圖18 為d取得最小值時兩指標對應的風剖面對比圖（計算高度范圍同上）.

圖17 測區(qū)二d值變化曲線Fig.17 Change curve of d value in test site 2

圖18 測區(qū)二風剖面對比圖Fig.18 Comparison of wind profiles in test site 2

3.2.4 誤差分析

圖19 給出了分別根據(jù)地面粗糙度計算所得α2和風場實測得到的α1值在梯度風速均為vG=20 m/s時的風剖面誤差圖（近地面200 m 范圍內(nèi)）.由圖19 可知，這兩種方法得到的風剖面曲線較為貼近，相對誤差為2.09%.分析可知測區(qū)二的風場類別介于B、C類之間，分別以不同的梯度高度范圍進行計算及誤差分析，結果如表3所示.

表3 測區(qū)二不同高度誤差分析表Tab.3 Error analysis of different heights in test site 2

圖19 測區(qū)二風剖面誤差圖Fig.19 Wind profile error in test site 2

由表3 可知，在B 類地貌和C 類地貌梯度高度范圍內(nèi)根據(jù)地面粗糙度計算出的α2值均為0.203，與風場實測結果的相對誤差為3.52%.根據(jù)不同高度范圍計算結果得到的結論同測區(qū)一所得結論一致.

4 結論

本文基于無人機航拍和圖像識別方法對城市地區(qū)的建筑物尺寸信息進行提取，通過計算出測區(qū)內(nèi)地貌粗糙程度來判定其風場類別，并與空中風場實測結果進行對比，探討了基于無人機航拍識別城市地區(qū)近地風場類別的可行性，經(jīng)分析得到以下結論：

1）基于無人機航拍和圖像識別方法的應用，可快速、高精度地獲取建筑物的尺寸信息，進而準確地識別測區(qū)內(nèi)的地貌狀況.據(jù)此得到的粗糙度指數(shù)α值與空中風場實測得到的α值相差不大，兩者判定的風場類別一致.

2）按照無人機航拍識別風場類別的方法對近地面不同高度范圍進行了分析，采用不同的計算高度范圍得到的計算結果略有差異，且計算高度達到一定值時其結果不再改變.根據(jù)實際地貌計算得到的粗糙度指數(shù)α值與風場實測結果之間的相對誤差均小于5%，且計算高度范圍與風場實測高度范圍越接近兩者相對誤差越小.

3）本文基于無人機的輔助，提出了一套可靠度高且具有普適性的近地風場類別識別方法，該方法實現(xiàn)了基于實際地貌判定同目標建筑相關的風場類別，研究結果可為城市地區(qū)建筑抗風設計提供參考.