基于極限梯度提升和探地雷達時頻特征的水泥路面脫空識別

張 軍， 姜文濤， 張 云， 羅婷倚， 余秋琴， 楊 哲

（1.長安大學 公路養護裝備國家工程實驗室，陜西 西安 710064；2.長安大學 道路施工技術與裝備教育部重點實驗室，陜西 西安 710064；3.廣西北投公路建設投資集團有限公司，廣西 南寧530028；4.廣西交科集團有限公司，廣西 南寧 530007）

水泥路面具有強度高、耐久性好的特點，但存在水泥路面層與基層剛度不匹配問題。在溫度和車輛荷載的重復作用下，面層底部出現脫空區域并逐步擴展，最終導致路面出現斷板。脫空區域的橫向尺寸是影響路面結構承載力的重要參數［1］，因此準確識別水泥路面脫空的橫向尺寸對水泥路面的養護工作具有指導意義。

探地雷達（GPR）是利用高頻電磁波在介質中的反射和散射來實現淺層成像和定位的高分辨率深層無損探測技術，通過路面電磁特性變化定性或定量地識別地下目標體，因此被廣泛應用于地下空洞［2］、管線［3-4］和路面病害［5］的定位與識別。路面病害的識別集中于基于B-scan 的圖譜分析和基于A-scan 的機器學習。圖譜分析法是借助深度遷移學習對雷達B-scan圖譜中的病害區域進行識別［6-8］，已用于傳力桿［9-10］、管道［11］、裂縫［12］和其他目標的智能識別［13］。王輝等［14］設計了級聯結構的卷積神經網絡，用于雙曲線目標分類。針對脫空區域的圖譜特征識別，筆者采用淺層混合網絡Resnet18-YOLOv2 實現了機場跑道空洞區域的定位［6］。然而，深度遷移學習需要大量準確標記的樣本圖譜和大量取芯驗證，加上路面病害區域無固定的寬高比，難以在圖譜上準確標注病害區域，導致病害區域的邊界定位存在誤差，因此需要研究病害邊界的高精度識別方法。

相比圖譜識別，基于A-scan信號特征的病害識別在橫向尺寸上獲得更高的精度，目前常見的建模方法有支持向量機（SVM）［15-16］、人工神經網絡（ANN）［7］、隨 機 森 林（RF）［17］和 極 限 梯 度 提 升（XGBoost）等方法［18-19］。覃暉等［20］將GPR 信號分段并提取方差、標準絕對偏差和四階矩等3個特征，對隧道襯砌空洞病害的識別準確率為93.56%，但時域特征抗強干擾能力較差。為提高模型的抗干擾能力，周輝林等［21］采用2 個時域信號特征（幅值、平均絕對偏差）和4 個小波近似系數作為SVM 輸入建模，公路路基病害檢測識別的準確率為92.7%，但未明確病害類型。為獲得路面病害的敏感特征參數，筆者借鑒旋轉設備的故障識別方法，用28個時頻特征表征病害，通過主成分分析（PCA）降維和ANN構建了瀝青路面水損害識別模型，識別準確率達到92.4%［7］。杜豫川等［19］將每道A-scan 數據劃分為p段，分別提取能量、方差、偏度和對數功率譜等特征值，并將特征值輸入XGBoost用于模型訓練，對路面脫空病害的識別準確率為96%，但段數p需要根據經驗設置。現有研究為脫空病害識別提供了參考，但仍存在以下問題：①未考慮不同的GPR設備天線頻率和采樣頻率之間存在的差異；②模型的泛化能力不強，需要人工設置合理的模型參數；③模型中病害區域的數據樣本量小。因此，還需要研究高效、高精度的脫空病害識別方法。

以水泥路面脫空病害為研究對象，采用正演模擬、室內模型和現場試驗相結合的方法建立了標準脫空數據集，其中包含10 732 條脫空數據和10 251條正常數據。針對GPR 天線頻率和采樣參數不一致的問題，提出用重采樣和標準化方法統一數據采樣信息和幅值范圍。在時頻特征提取方法［22］的基礎上，增加了平均能量和偏度［19］2 個時域特征，將脫空數據集進行特征處理，構建了30 個時頻特征數據集。采用XGBoost 算法建立脫空識別模型，并通過取芯驗證了模型的準確性。

1 材料與方法

1.1 水泥路面脫空病害

在實際服役中水泥路面受到溫度荷載和車輛荷載的耦合作用，出現由板底塑性變形導致的板底脫空，最終形成斷板。板底脫空包括充氣脫空和含水脫空，充氣脫空是板底脫空的早期形態。圖1 展示了板底脫空的部位。隨著時間的推移，脫空尺寸逐漸擴大，雨水滲透后形成含水脫空。

圖1 水泥路面板底脫空Fig.1 Void under cement pavement slab

1.2 脫空數據采集

以板底充氣脫空為例，簡化充氣脫空的形狀，使用16～300 mm不等的矩形脫空和圓形脫空進行室內試驗和正演模擬，并綜合現場試驗數據構建脫空數據集。正演模擬和室內試驗提供標準脫空特征，并為現場判斷脫空提供依據，現場試驗數據用于模型訓練和驗證模型的準確性。脫空數據集構成如圖2所示。

圖2 脫空數據集構建Fig.2 Void dataset construction

1.2.1 正演模擬

gprMax 正演模擬是以時域有限差分為基礎的雷達信號仿真算法，被廣泛應用于GPR 的病害模擬，以獲得病害區域的信號特征［23-24］。構建如圖3所示的脫空模型，C1～C4為圓形脫空，R1～R4為矩形脫空。根據實際路面的材料特性，設置脫空、水泥混凝土和半剛性基層的介電常數分別為1、7和9，水泥混凝土厚度為240 mm，仿真參數如表1所示。采用gprMax 3.1.5 進行正演模擬，獲得如圖4 所示的結果，脫空病害區域均表現為雙曲線特征，其寬度與病害橫向尺寸相關。

表1 gprMax仿真模型參數Tab.1 Parameters of gprMax simulation model

圖3 正演脫空模型（單位：mm）Fig.3 Simulated void model（unit：mm）

圖4 正演脫空模型的B-scan結果Fig.4 B-scan results of simulated void model

1.2.2 室內試驗

從仿真結果可知，脫空形狀并不影響脫空特征。為此，室內模型僅考慮不同尺寸的矩形脫空，構建了橫向尺寸分別為100、90、80、70、60、50、40、30、25、16 mm 的10 個矩形脫空區域A1～A10，脫空區域的中心深度和中心間隔分別為200 mm和150 mm。模型材料為型號C30 的水泥混凝土，尺寸為2 070 mm×400 mm，如圖5 所示。意大利IDS 公司的RIS 型900 MHz 和美國US 公司的1 000 MHz 天線沿著如圖5所示的雷達測線方向對模型中的脫空區域進行多次數據采集，具體雷達參數設置如表2所示。獲得的雷達圖譜如圖6 所示，分別為RIS 和US雷達實測剖面圖譜。

表2 探地雷達系統參數設置Tab.2 GPR system parameter settings

圖5 室內試驗水泥板模型Fig.5 Model of lab test

圖6 室內模型雷達實測剖面圖譜Fig.6 GPR results of lab model

探地雷達的縱向分辨率與雷達天線的頻率有關，具體關系如下所示：

式中：Δr為探地雷達的縱向分辨率；λe為探地雷達波在介質中傳播的波長。

以波速0.1 m·ns-1為例，900 MHz和1 000 MHz的電磁波在混凝土的波長分別為111 mm 和100 mm，根據式（1）可得最小可探測脫空尺寸分別約為28 mm 和25 mm，因此2 個天線難以分辨25 mm 和16 mm兩處脫空區域，導致這兩處區域的GPR圖譜不清晰，但可以使用更高天線頻率的GPR設備或者三維GPR檢測小尺寸的脫空區域。

1.2.3 現場試驗

對廣西桂林永福縣永鹿路部分路面進行數據采集，部分路面已經出現斷板（見圖7b）。由專業檢測人員和安全員協調下進行GPR 檢測，如圖7d 所示。試驗設備包括MALA ProEx 系統、800 MHz 天線和安裝了GroundVision 的個人電腦，采樣頻率為14.5 GHz，每道采樣點數為308，時窗觸發間隔為0.01 s。 典型路面圖譜如圖8所示。

圖7 現場試驗Fig.7 Field test

圖8 現場路面B-scan圖譜Fig.8 B-scan results of field test

采用以上3 種試驗方法，共獲取60 萬道A-scan數據，正演模擬5萬道，室內模型5萬道，現場檢測50萬道。使用ReflexW 軟件對上述數據進行分析，并對數據集中的脫空區域數據給出“1”的真值標簽，正常區域數據給出“0”的真值標簽。經過人工篩選后，選擇20 983 道A-scan 數據作為數據集，其中10 623道為室內和現場數據（脫空6 745道，正常3 878道），10 360道為正演數據（脫空3 987道，正常6 373道），合計脫空數據10 732道，正常數據10 251道。

1.3 GPR信號處理方法

對正反演GPR 數據進行分析，如圖9 所示。與現場數據相比，正演數據中沒有噪聲和干擾，而室內試驗與現場試驗數據均存在干擾和噪聲。因雷達設備的不同，室內模型數據的幅值明顯低于現場路面數據的幅值，但在脫空病害部位具有相似的回波形狀。因此，需要對GPR 數據進行預處理，才能將正反演數據結果都用于后期建模。

圖9 正演模擬、室內和現場試驗的A-scan數據對比Fig.9 Comparison of GPR signal among simulation,lab and field tests

通過對比，確定了如圖10所示的七步預處理方法以增強脫空特征，包括去直流漂移（去除直流漂移）、零點矯正（去除天線與地面距離的誤差）、能量增益（增強深層病害信號）、背景去除（減去平均道）、巴特沃斯帶通濾波（濾除高頻雜波干擾）、滑動平均和F-K 偏移（糾正脫空區域偏移成分）。同時，為了解決采樣頻率不一致，采用重采樣方法將原始數據轉為具有相同采樣頻率的數據。

圖10 雷達數據預處理流程Fig.10 Flow chart of radar data preprocessing

不同天線頻率以及數據源差異均導致電磁波幅值差異，因此采用標準化處理方法，將GPR 數據處理成均值為0 和標準差為1 的數據，具體方法如下所示：

式中：X為探地雷達采集的數據矩陣；為標準化后的數據矩陣；σ為整體數據方差；μ(X)為整體數據均值。

1.4 特征提取

參照水損害的特征提取方法［7］，對每一條Ascan數據提取了16個時域特征和12個頻域特征，另外增加平均能量和偏度（P18）2 個時域特征，累計30個特征（P1～P30）。P8、P17、P18的計算式如下所示：

式中：xi為該道A-scan第i個采樣點數據，i=1，2，…，n，n是每條A-scan的采樣點數。對提取的30個時頻域特征采用式（2）進行標準化處理。圖11為3種數據源中脫空區域和正常區域特征參數的對比結果。對比正常和脫空的信號特征，發現正演模型的特征中除P6、P9、P19、P27外，其余特征存在明顯的差異；在室內模型的特征數據中，除P19、P21和P22外，其余特征均存在明顯差異；在實際路面特征中，除P3和P5外，其余特征數據均存在明顯差異。結果表明，在這3種數據源中，正常與脫空數據的特征P2、P4、P7、P8、P10、P11、P13～P17、P20、P24～P26、P28和P30均存在明顯的差異，說明所選的時頻域特征可用于脫空病害表征。

圖11 30個時頻域特征值對比Fig.11 Comparison of 30 time-frequency domain feature values

為了實現脫空識別，對數據集的20 983 道數據進行時域和頻域特征的提取，提取的時頻域特征組成一個30×N（N=20 983）的特征數據集矩陣。將1.2節中給出的對應每一條A-scan數據的真值標簽移植到時頻域特征數據集中，構成含有真值標簽的31×N特征數據集，這個數據集用于后續模型訓練。

1.5 特征重要性分析

對提取的30 個時頻域特征進行特征重要性分析，選擇最小的特征子集，可加快模型訓練速度和提升識別性能。特征重要性分析常用的方法有卡方檢驗、互信息、基于樹模型和費希爾信息的方法。由于各種特征分析方法存在評價指標的差異，因此采用基于多種方法融合的分析方法來選取重要特征。特征重要性分析流程如圖12所示。

圖12 特征重要性分析流程Fig.12 Flow chart of feature significance analysis

由于4種特征重要性分析方法的評價得分范圍存在差異，單純地將得分相加不能得到客觀的最終得分，因此對4種分析方法的得分進行歸一化處理，歸一化結果相加得到最終的重要性評價得分。4 種分析方法歸一化和相加的最終得分如圖13所示，時頻域特征重要性分析的最終得分排序如圖14所示。

圖13 特征重要性得分歸一化結果Fig.13 Normalization results of feature significance score

圖14 特征重要性得分排序Fig.14 Ranking of feature significance

選取綜合得分大于1的21個時頻域特征作為識別模型的輸入。進一步地去除31×N數據集中非重要的特征，得到含有真值標簽的22×N數據集，用于識別模型的訓練。

2 基于XGBoost的路面脫空識別算法

2.1 XGBoost算法

對于XGBoost算法，為避免模型過擬合，加快訓練速度，提升擬合及預測精度，損失函數應用了二階泰勒展開，在目標函數中加入正則項。將含有N個樣本m個特征的數據集分為2類，數據集

建立二分類模型，樣本Xi在第k輪預測值的表達式如下所示：

式中：K表示模型中總的決策樹數；F表示回歸樹空間；ω表示回歸樹葉子節點權重；q表示回歸樹結構，把每個樣本節點映射到對應葉子節點的索引；q(Xi)表示樣本Xi所在的葉子節點；T表示葉子節點數。

XGBoost算法的目標函數包括損失函數和正則項兩部分，定義如下：

2.2 橫向尺寸的權重設定

GPR 檢測過程中存在噪聲、干擾和天線跳動問題，連續的GPR數據中存在單道或者連續幾道的干擾信號，致使識別模型對此類數據給出錯誤的識別結果。小尺寸的脫空病害也會對識別結果產生干擾，降低識別結構的可讀性。采用如圖15所示的高斯權重系數，當連續的GPR數據識別結果中存在孤立點時，將雷達圖譜的識別結果中孤立點兩側的識別結果乘以高斯權重系數進行綜合考慮。應用時將圖15中的中心點（第i個識別結果）對應的數據作為該道識別結果的權重，其左右5 個權重系數分別賦給該道左右相鄰5 道，通過與設定的閾值進行比較來過濾孤立點。

圖15 高斯權重系數曲線Fig.15 Gaussian weight coefficient curve

2.3 識別模型建立流程

將GPR 信號數據集通過特征提取和歸一化后構建特征數據集，并分為訓練集和測試集，訓練集用于XGBoost 模型訓練，測試集用于模型的識別性能測試。預測模型建立流程如圖16所示。

圖16 XGBoost模型訓練流程Fig.16 Flow chart of XGBoost model training

（1） 將含有真值標簽的時頻特征數據集（22×N）隨機抽取70%作為XGBoost 模型訓練集樣本，30%作為模型測試集樣本。

（2） 將訓練集樣本輸入算法模型中進行訓練，并采用交叉驗證和網格搜索方法調整模型的max_depth（最大樹深）和n_estimators（子學習器個數）等超參數，獲得最優的預測模型。

（3） 將測試集樣本輸入預測模型，檢驗模型的預測準確率，并存儲預測模型。訓練后獲得的最優模型將用于后續路面脫空病害識別。

2.4 模型評價標準

模型識別結果的評價指標可以從混淆矩陣中導出。具體的分類評價標準及其公式和含義如表3所示。表3中，N表示測試集樣本總數，Mii表示模型對類別ci準確識別的樣本個數，M*i表示類別ci樣本個數，Mi*表示模型識別為類別ci的樣本個數。上述指標越接近1，代表模型的分類效果越好。

表3 準確性評價指標Tab.3 Accuracy evaluation index

3 試驗結果分析

3.1 模型優化

為提高模型的準確率，采用五折交叉驗證和網格搜索方法調優模型的關鍵參數［25］。XGBoost模型中max_depth和n_estimators對預測模型的過擬合性敏感，并影響模型收斂的穩定性。根據XGBoost 模型超參數調參時的經驗，將max_depth 調優范圍設置為［2，20］，n_estimators 調優范圍設置為［10，2 000］，比較在不同的最大樹深和不同的子學習器個數的情況下模型的收斂速度與識別精度。根據圖17a 中模型max_depth 參數調優曲線，將max_depth設置為11，模型的準確率為94.12%。根據圖17b中n_estimators調優曲線，將n_estimators設置為1 000，模型的準確率為94.14%。對預測模型其他參數采取同樣的調優操作，最終確定模型的最優參數組合，如表6所示。

圖17 參數調優曲線Fig.17 Parameter tuning curve

3.2 模型性能和對比

使用表4中調優后的模型超參數數值和70%的數據集進行XGBoost 模型訓練。為比較模型性能，相同的數據集還用于RF和ANN模型訓練。模型建模平臺均使用Python 語言，在CPU-AMD R7 5800H（3.20 GHz）、16 GB RAM 和GPU-NVIDIA GeForce RTX 3060 Laptop（6 GB）計算機上進行。將驗證集（30%的數據集）輸入建立的3種識別模型中，計算3 種識別模型在驗證集上的評價指標。由表5可見，3種模型在測試集上的識別準確率排序為XGBoost （99.63%） ＞ANN （99.30%） ＞RF（99.28%）。從識別準確率看，XGBoost模型比另外2 種模型更有效地區分雷達數據中的脫空病害和正常數據，對測試集中的脫空數據產生的誤判率更低。

表4 模型參數調優值Tab.4 Parameter values of optimal model

表5 模型識別性能評價Tab.5 Performance evaluation of recognition models

為了確定所提方法的有效性，在測試數據集（950 道數據）上進行對比，結果如圖18 所示。可以看出，3 種識別模型結果中都存在如圖18 中E 標識所指的單道（孤立點）A-scan誤判結果。XGBoost模型的識別結果中存在的單道A-scan 誤判情況明顯優于RF和ANN模型的識別結果。統計3種模型的識別準確率和消耗時間，如表6 所示。 XGBoost 模型對水泥路面雷達數據的識別準確率可達98.10%，識別過程共耗時0.141 s；RF 算法的識別準確率為93.17%，識別過程共耗時0.196 s； ANN 模型的識別準確率為95.10%，識別過程共耗時0.203 s。與RF 和ANN 模型相比，XGBoost 模型具有更高的預測準確率和更快的識別速度 。

表6 模型識別性能Tab.6 Performance of recognition models

圖18 XGBoost、RF和ANN模型識別結果Fig.18 Performance of XGBoost, RF and ANN models

4 模型適用性

4.1 正演模型數據

將圖4中的正演模型數據采取與訓練集相同的預處理流程，提取雷達數據的時頻域特征，輸入XGBoost 識別模型中，獲得如圖19 所示的預測結果。可以看出，XGBoost 模型對正演圖譜中大多數脫空區域都給出了“1”的識別結果，準確給出圓形脫空的雙曲線特征范圍，對矩形脫空的橫向尺寸給出了較高準確率的識別。圖20 為XGBoost 正演模型識別脫空尺寸和實際脫空尺寸的對比。從圖20 可見，XGBoost 正演模型對矩形脫空和圓形脫空的橫向尺寸的誤差在8 mm以內。造成此現象的原因是：F-K偏移對仿真數據中的圓形脫空雙曲線特征產生了誤差。

圖20 脫空橫向尺寸對比Fig.20 Comparison of lateral dimension

4.2 現場驗證

任取2 段存在脫空區域的數據，考慮單個脫空和多處脫空的狀態，現場檢測GPR 結果如圖21a、b所示。按照圖10預處理流程進行了特征提取后，將特征數據輸入XGBoost 模型，識別結果如圖21 所示。圖21c、e為模型直接輸出的識別結果。可知，脫空區域被準確識別，但識別結果中存在多道孤立點，具體原因如下：

圖21 實際路面脫空病害識別Fig.21 Identification of actual pavement void disease

（1） XGBoost 模型給出連續的“1”的識別結果比孤立的“1”的識別結果更能說明此處路面脫空病害的存在，在對路面脫空病害尺度的評估中占有更大權重。

（2） XGBoost 模型識別結果中存在稀疏的孤立“1”的結果，可能不是由路面脫空病害造成的，如數據集不全面；在進行雷達數據采集時存在干擾造成的雷達信號波動；存在早期的脫空病害（橫向尺寸小于2 cm），但尺寸過小難以驗證。

實際項目中主要關注大范圍脫空區域，為減少孤立點對識別結果的干擾，引入識別結果的后處理方法。當一條A-scan的XGBoost模型識別結果為“1”，左右相鄰的連續10條雷達數據的識別結果與高斯判斷系數相乘并相加的和等于或大于5時，此條A-scan識別結果才為“1”。此外，如果一個測試段（10條Ascan）雷達數據中有90%的A-scan都被XGBoost模型識別為“1”，整個測試段就被視為一個連續的脫空區域，XGBoost模型的結果都修改為“1”。

對XGBoost 模型識別結果進行后處理，得到如圖21d、f 所示的后處理結果，后處理可以有效地去除孤立點的影響。在圖21d 中，已經形成了長尺寸的脫空區域，連續脫空區域會在車輛交變荷載等作用下發展為更大的脫空，需要及時采取補救措施。GPR 專家可以識別圖21f 中的脫空區域。然而，在對比度較低的雷達圖譜中脫空區域很難評價。在0～100 條和300～400 條A-scan 之間存在兩處小尺寸的脫空區域，說明此脫空區域處于早期階段；在200～300 條A-scan 存在相鄰很近的兩處脫空區域，說明這兩處脫空會發展為連續的脫空區域。

上述驗證結果表明，所構建的模型可有效地確定路面脫空區域。為進一步驗證模型的準確性，在封閉交通的條件下，對模型識別出脫空的位置（見圖22a）采用破碎鎬進行破碎，并通過針孔攝像頭對破碎孔進行觀察，結果如圖22b所示。針孔成像結果驗證了脫空病害的存在，進一步驗證了算法的準確性。

圖22 取芯驗證Fig.22 Verification by coring

5 結論

（1）分析了GPR在脫空和正常區域的時域和頻域響應差異，提出了18個時域和12個頻域的時頻域特征，用以表征脫空和正常區域的GPR數據。通過正演模擬、室內試驗和現場試驗構建脫空數據集來訓練XGBoost模型，實現了脫空病害的自動識別。

（2） 提出了重采樣和標準化的GPR數據處理方法，統一不同GPR 數據的采樣頻率，解決了不同采樣頻率和不同天線頻率對識別模型的影響，可以解決不同源的數據使用和識別問題。

（3） 對 比 了XGBoost、RF 和ANN 的 模 型，XGBoost 模型最優，脫空識別準確率達到了98.10%，高于ANN（95.10%）和RF（93.17%）。

作者貢獻聲明：

張 軍：論文構思，論文撰寫和修改。

姜文濤：模型構建，程序設計，論文撰寫及修改。

張 云：開展室內和現場試驗。

羅婷倚：開展現場試驗和數據分析。

余秋琴：試驗數據分析。

楊 哲：孤立點的權重系數設計。