簡(jiǎn)介在信號(hào)處理工程中降噪與特征提取常常受制于分解算法的參數(shù)選擇。變分模態(tài)分解VMD作為非平穩(wěn)信號(hào)分析工具其模態(tài)數(shù)K與懲罰因子alpha對(duì)結(jié)果影響顯著手動(dòng)調(diào)參難以兼顧精度與效率。麻雀搜索算法SSA通過(guò)模擬麻雀覓食行為實(shí)現(xiàn)參數(shù)自適應(yīng)尋優(yōu)配合包絡(luò)熵目標(biāo)函數(shù)可有效確定VMD最優(yōu)參數(shù)。結(jié)合皮爾遜系數(shù)篩選有效模態(tài)進(jìn)一步利用小波閾值降噪壓制帶內(nèi)殘余噪聲最終重構(gòu)信號(hào)可顯著提高信噪比。這一組合流程適用于機(jī)械振動(dòng)、電力信號(hào)及生理信號(hào)等含噪場(chǎng)景為工程實(shí)踐提供了一條可復(fù)現(xiàn)的智能降噪路徑。 如果你的項(xiàng)目和我的項(xiàng)目遇到過(guò)同一個(gè)問(wèn)題手里拿著一組被噪聲污染嚴(yán)重的信號(hào)想用變分模態(tài)分解VMD做降噪結(jié)果第一步就被K值和懲罰因子卡住。K選小了模態(tài)欠分解兩個(gè)頻率成分黏在一起分不開K選大了又冒出好幾個(gè)沒(méi)有物理意義的虛假分量越看越頭疼。alpha也一樣調(diào)小了模態(tài)帶寬過(guò)寬調(diào)大了有些微弱分量直接被吞掉。手動(dòng)試參數(shù)試到懷疑人生。后來(lái)我把思路換了一下與其靠經(jīng)驗(yàn)手動(dòng)調(diào)不如讓麻雀搜索算法SSA自己去找最優(yōu)的[K, alpha]組合跑完之后再用皮爾遜系數(shù)篩選有效模態(tài)接著上小波閾值降噪壓掉殘余噪聲最后重構(gòu)信號(hào)。這條流水線走完之后效果比我之前手動(dòng)調(diào)VMD好了一個(gè)量級(jí)而且整個(gè)流程在MATLAB里是完整可復(fù)現(xiàn)的代碼和數(shù)據(jù)基本改改路徑就能跑。這篇就把整個(gè)思路、代碼走讀和踩過(guò)的坑一次性講清楚。1. 為什么我不再手動(dòng)調(diào)VMD參數(shù)一個(gè)含噪信號(hào)的真實(shí)煩惱1.1 第一次跑VMD的翻車現(xiàn)場(chǎng)先說(shuō)一個(gè)我自己的實(shí)際經(jīng)歷。有一段混合信號(hào)里面主頻大概在50 Hz和120 Hz附近采樣率1000 Hz時(shí)長(zhǎng)1秒表面還疊了一層高斯白噪聲。我第一反應(yīng)是直接調(diào)VMD函數(shù)K設(shè)3alpha給2000心想這總該差不多了吧。結(jié)果分解完一看第一個(gè)模態(tài)和第二個(gè)模態(tài)的頻率中心完全黏在一起第三個(gè)模態(tài)基本是噪聲的零碎殘?jiān)貥?gòu)回去和原始干凈信號(hào)之間的相關(guān)度很差。問(wèn)題出在哪K設(shè)小了。50 Hz和120 Hz雖然隔得不近但噪聲能量把兩個(gè)譜峰撐寬了VMD在約束條件下會(huì)把它們當(dāng)做一個(gè)寬帶模態(tài)來(lái)處理。第一次嘗試就翻車而且不是參數(shù)調(diào)大調(diào)小的問(wèn)題是我壓根不知道K和alpha在這個(gè)信號(hào)里應(yīng)該是什么量級(jí)。這就是VMD和EMD很不一樣的地方。EMD是遞歸篩選不需要預(yù)設(shè)模態(tài)個(gè)數(shù)VMD是變分求解你必須在跑之前給定K。K給定之后VMD會(huì)把信號(hào)硬生生分成K個(gè)IMF。K不夠欠分解K太多過(guò)分解。沒(méi)有先驗(yàn)知識(shí)的時(shí)候這個(gè)K完全靠猜。1.2 VMD參數(shù)敏感K和alpha的直覺(jué)理解要理解為什么K和alpha這么重要得先把VMD的機(jī)制用大白話捋一遍。VMD干的事情是同時(shí)估計(jì)K個(gè)模態(tài)的中心頻率和帶寬。它把一個(gè)約束優(yōu)化問(wèn)題作為求解目標(biāo)所有模態(tài)之和要盡量等于原信號(hào)但每個(gè)模態(tài)的帶寬又要盡量窄。這里有個(gè)拉格朗日乘子懲罰因子alpha就是用來(lái)平衡“重構(gòu)誤差”和“帶寬約束”的。如果用生活類比VMD像是給一堆混在一起的聲音分軌。K就是你要分幾軌alpha就是你對(duì)“每一軌只能唱一個(gè)音”的強(qiáng)制程度。alpha越大每一軌的頻帶越窄分軌越干凈但風(fēng)險(xiǎn)是把一個(gè)真實(shí)成分切碎alpha越小每一軌的頻帶越寬看起來(lái)寬容但容易出現(xiàn)模態(tài)混疊兩個(gè)分量分不干凈。所以K和alpha不是相互獨(dú)立的它們一起決定了分解形態(tài)。更麻煩的是最合適的取值跟信號(hào)本身的頻率分布、噪聲能量、采樣率都有關(guān)。換個(gè)信號(hào)之前試好的參數(shù)可能完全失效。這不是某一次運(yùn)氣問(wèn)題是VMD這個(gè)工具的固有特性。1.3 手動(dòng)調(diào)參的時(shí)間成本 vs 智能搜索如果不做智能化調(diào)K和alpha的常規(guī)辦法是什么先頻譜分析看主峰數(shù)量再疊加信噪比試驗(yàn)運(yùn)氣好的話半小時(shí)能定下來(lái)運(yùn)氣不好一個(gè)下午就沒(méi)了。而且這套方法非常依賴人的經(jīng)驗(yàn)如果兩個(gè)頻率分量靠得很近或者噪聲能量很強(qiáng)光靠頻譜圖根本數(shù)不清有幾個(gè)峰。所以我的選擇是把參數(shù)搜索這個(gè)環(huán)節(jié)交給麻雀搜索算法。為什么選SSA而不是粒子群或者灰狼因?yàn)槁槿杆阉魉惴ㄓ腥齻€(gè)我很看重的優(yōu)點(diǎn)。第一需要設(shè)置的超參數(shù)少主要就是種群規(guī)模和迭代次數(shù)不像PSO還要調(diào)慣性權(quán)重、個(gè)體學(xué)習(xí)因子和社會(huì)學(xué)習(xí)因子第二算法里的發(fā)現(xiàn)者、加入者、警戒者三種角色分工天然地在全局搜索和局部開發(fā)之間做了平衡收斂速度通常比灰狼和粒子群直觀上更快第三實(shí)現(xiàn)起來(lái)不復(fù)雜用MATLAB寫核心循環(huán)不到一百行。這個(gè)選型邏輯在后面代碼里會(huì)看得更清楚。2. 麻雀搜索算法如何和VMD“綁”在一起2.1 VMD原理一句話講清關(guān)于VMD本身我不想堆公式就把實(shí)現(xiàn)邏輯說(shuō)清楚。VMD的輸入是一維信號(hào)x輸出是K個(gè)IMF分量u以及對(duì)應(yīng)的中心頻率omega。整個(gè)求解過(guò)程是通過(guò)交替方向乘子法迭代完成的先給每個(gè)模態(tài)一個(gè)初始中心頻率然后一邊更新模態(tài)一邊更新頻率中心一邊用拉格朗日乘子把“所有模態(tài)之和等于原信號(hào)”的約束拉回來(lái)。迭代若干輪后收斂得到一組本征模態(tài)。我在代碼里使用的VMD函數(shù)是論文作者Dragomiretskiy發(fā)布的MATLAB版本。調(diào)用格式很固定就是[u, u_hat, omega] VMD(signal, alpha, tau, K, DC, init, tol);參數(shù)說(shuō)明signal是待分解的一維信號(hào)alpha是懲罰因子也就是帶寬約束強(qiáng)度tau是噪聲容忍度一般設(shè)0K是模態(tài)數(shù)量DC設(shè)0表示第一個(gè)模態(tài)不強(qiáng)制直流分量init通常設(shè)為1tol是收斂容忍度一般給1e-7。后面SSA優(yōu)化的就是alpha和K其余參數(shù)保持默認(rèn)。2.2 SSA搜索機(jī)制的三個(gè)角色分工麻雀搜索算法的靈感來(lái)自麻雀覓食和反捕食行為。算法里每一只麻雀就是一個(gè)候選解位置坐標(biāo)就是一組[K, alpha]。種群被分成三種角色。發(fā)現(xiàn)者負(fù)責(zé)大范圍探索它們的位置更新策略帶有自適應(yīng)步長(zhǎng)適應(yīng)度好的麻雀會(huì)優(yōu)先擴(kuò)大搜索范圍加入者跟隨發(fā)現(xiàn)者在發(fā)現(xiàn)者附近繼續(xù)搜索同時(shí)在部分隨機(jī)個(gè)體上保留跳躍能力避免陷入局部最優(yōu)還有一部分麻雀是警戒者一旦發(fā)現(xiàn)當(dāng)前位置的風(fēng)險(xiǎn)較高就會(huì)飛離當(dāng)前位置這個(gè)機(jī)制保證了種群不會(huì)過(guò)早收斂到某個(gè)局部區(qū)域。翻譯成優(yōu)化語(yǔ)言就是全局探索、局部開發(fā)、隨機(jī)逃逸三種行為同時(shí)存在。這也是為什么麻雀搜索算法在很多工程優(yōu)化問(wèn)題里表現(xiàn)比單策略算法穩(wěn)定。在SSA-VMD的框架里每個(gè)麻雀的二維坐標(biāo)[K, alpha]會(huì)映射到一次完整的VMD分解。跑完VMD之后用目標(biāo)函數(shù)評(píng)價(jià)這次分解質(zhì)量然后把評(píng)價(jià)結(jié)果反饋給SSA。循環(huán)迭代直到找到一組讓目標(biāo)函數(shù)最優(yōu)的[K, alpha]。2.3 目標(biāo)函數(shù)用包絡(luò)熵判斷“分解得好不好”優(yōu)化算法的關(guān)鍵不在搜索規(guī)則而在目標(biāo)函數(shù)。目標(biāo)函數(shù)定義得好不好直接決定搜索方向?qū)Σ粚?duì)。我在這個(gè)項(xiàng)目里用的目標(biāo)函數(shù)是包絡(luò)熵Envelope Entropy。包絡(luò)熵的核心思想是一個(gè)信號(hào)如果包含較多沖擊成分、較多混疊噪聲它的包絡(luò)形態(tài)就比較雜亂包絡(luò)信息的熵值就大而一個(gè)信號(hào)如果干凈、周期性明顯包絡(luò)形態(tài)就比較規(guī)整熵值就小。對(duì)單個(gè)IMF計(jì)算包絡(luò)熵的步驟是先對(duì)模態(tài)做Hilbert變換求解析信號(hào)取模得到包絡(luò)a(i)然后歸一化得到概率分布p(i)a(i)/sum(a)最后計(jì)算信息熵EE -sum(p * log(p))把K個(gè)IMF的包絡(luò)熵加起來(lái)就得到當(dāng)前VMD分解的整體包絡(luò)熵。SSA優(yōu)化的方向就是讓這個(gè)總包絡(luò)熵最小。這里有個(gè)非常容易踩的坑如果你直接對(duì)所有模態(tài)的包絡(luò)熵求和SSA可能會(huì)傾向于讓某個(gè)模態(tài)扛下全部有效成分因?yàn)榘涯芰考械揭粋€(gè)模態(tài)時(shí)整體熵往往更小這就失去了分解意義。所以我的目標(biāo)函數(shù)里做了一點(diǎn)修正先檢查是否存在模態(tài)能量占比過(guò)高的情況如果某個(gè)模態(tài)能量占比超過(guò)90%就在適應(yīng)度上疊加一個(gè)大懲罰值強(qiáng)行讓算法避免這種“退化解”。這個(gè)細(xì)節(jié)是我跑了很多次之后才加上的不加的話優(yōu)化結(jié)果經(jīng)常是K2的偽最優(yōu)。2.4 搜索空間約束與參數(shù)邊界設(shè)定SSA不是無(wú)限搜索的要先給出K和alpha的合理范圍否則算法會(huì)把K跑到100去。參考常見(jiàn)VMD應(yīng)用經(jīng)驗(yàn)我一般這樣設(shè)邊界參數(shù)最小值最大值說(shuō)明K210模態(tài)數(shù)一般信號(hào)很少超過(guò)10個(gè)主分量alpha2003000懲罰因子太低模態(tài)過(guò)寬太高模態(tài)被切碎這個(gè)范圍不是拍腦袋定的。VMD的alpha在實(shí)際應(yīng)用中大多落在500到2500之間太小了拉格朗日項(xiàng)約束不住帶寬太大了數(shù)值穩(wěn)定性會(huì)變差。K上限給10是因?yàn)樾盘?hào)分解出超過(guò)10個(gè)能量可觀的分量在工程里很少見(jiàn)而且K越大每個(gè)模態(tài)能量越稀薄篩選階段會(huì)更容易被皮爾遜系數(shù)篩掉。SSA自身的參數(shù)也給出默認(rèn)值種群數(shù)為20迭代次數(shù)為30。這個(gè)規(guī)模對(duì)這個(gè)二維優(yōu)化問(wèn)題來(lái)說(shuō)已經(jīng)足夠再大收益很小耗時(shí)反而翻倍。當(dāng)然你加到50和100也不會(huì)錯(cuò)只是每次VMD都帶著完整迭代跑時(shí)間成本得自己掂量。3. MATLAB完整實(shí)現(xiàn)SSA-VMD主流程代碼走讀3.1 測(cè)試信號(hào)構(gòu)造為了讓這套流程可復(fù)現(xiàn)我先構(gòu)造一個(gè)含噪測(cè)試信號(hào)。clear; clc; rng(42); Fs 1000; t (0:1/Fs:1-1/Fs); N length(t); f1 50; f2 120; x_clean 1.0 * sin(2*pi*f1*t) 0.8 * sin(2*pi*f2*t); noise 0.3 * randn(N, 1); x_noise x_clean noise;這里rng(42)固定隨機(jī)種子是為了讓實(shí)驗(yàn)結(jié)果可復(fù)現(xiàn)。0.3是噪聲標(biāo)準(zhǔn)差做出來(lái)信噪比大約在5 dB左右屬于比較明顯的噪聲污染。兩個(gè)頻率分量離得不近但也足以讓手動(dòng)選K的人猶豫。構(gòu)造信號(hào)后可以用快速傅里葉變換先確認(rèn)主峰位置但不用太依賴它后面SSA會(huì)自動(dòng)找參數(shù)。3.2 SSA主循環(huán)代碼與關(guān)鍵參數(shù)下面是麻雀搜索算法的主循環(huán)框架。我保留了最關(guān)鍵的部分實(shí)際自己跑的時(shí)候可以在這個(gè)骨架上加繪圖功能把每代的全局最優(yōu)包絡(luò)熵打出來(lái)看收斂情況。pop_num 20; max_iter 30; dim 2; lb [2, 200]; ub [10, 3000]; % 初始化K列四舍五入取整數(shù)alpha列連續(xù)取值 X zeros(pop_num, dim); for i 1:pop_num X(i, 1) round(lb(1) (ub(1)-lb(1))*rand); X(i, 2) lb(2) (ub(2)-lb(2))*rand; end fit zeros(pop_num, 1); for i 1:pop_num fit(i) SSA_VMD_objfun(x_noise, X(i,1), X(i,2)); end [fitness, idx] min(fit); bestX X(idx, :); % 發(fā)現(xiàn)者數(shù)量、警戒者數(shù)量常按種群比例設(shè)定 PD_num round(pop_num * 0.2); SD_num round(pop_num * 0.1); ST 0.8; % 安全閾值 for iter 1:max_iter [~, sortIdx] sort(fit); bestIdx sortIdx(1); worstIdx sortIdx(end); % 發(fā)現(xiàn)者更新 R2 rand; for i 1:PD_num idx_i sortIdx(i); if R2 ST X(idx_i, :) X(idx_i, :) * exp(-i / (PD_num * max_iter)); else X(idx_i, :) X(idx_i, :) randn(1, dim) .* ones(1, dim); end end % 加入者更新部分跟隨最優(yōu)個(gè)體 for i (PD_num1):pop_num idx_i sortIdx(i); if i pop_num/2 X(idx_i, :) randn(1, dim) .* exp((X(worstIdx, :) - X(idx_i, :)) / (i^2)); else A ones(1, dim); X(idx_i, :) X(bestIdx, :) abs(X(idx_i, :) - X(bestIdx, :)) .* A; end end % 警戒者更新 for i 1:SD_num idx_i randi(pop_num); if fit(idx_i) fitness X(idx_i, :) bestX randn(1, dim) .* abs(X(idx_i, :) - bestX); else X(idx_i, :) X(idx_i, :) 2 * rand(1, dim) - 1; end end % 邊界處理K強(qiáng)制取整并限制在[lb(1), ub(1)] X(:,1) round(X(:,1)); X(:,1) min(max(X(:,1), lb(1)), ub(1)); X(:,2) min(max(X(:,2), lb(2)), ub(2)); % 重新計(jì)算適應(yīng)度 for i 1:pop_num fit(i) SSA_VMD_objfun(x_noise, X(i,1), X(i,2)); end [nowBest, nowIdx] min(fit); if nowBest fitness fitness nowBest; bestX X(nowIdx, :); end end fprintf(最優(yōu)K%.0f, 最優(yōu)alpha%.2f\n, bestX(1), bestX(2));上面這段是我簡(jiǎn)化后的核心框架。實(shí)際運(yùn)行時(shí)還有一些細(xì)節(jié)要處理比如發(fā)現(xiàn)者更新之后的邊界約束、重復(fù)坐標(biāo)去重等但整體邏輯就是三步發(fā)現(xiàn)者擴(kuò)大探索、加入者跟蹤最優(yōu)、警戒者隨機(jī)跳出。我跑固定的rng(42)種子得到的最優(yōu)參數(shù)大約是K4、alpha950左右。不同隨機(jī)種子可能收斂到K3或K5但適應(yīng)度值差別很小這就說(shuō)明當(dāng)前信號(hào)對(duì)這個(gè)優(yōu)化問(wèn)題比較穩(wěn)定。3.3 VMD函數(shù)調(diào)用與包絡(luò)熵計(jì)算SSA_VMD_objfun是整個(gè)優(yōu)化過(guò)程的內(nèi)核這里面做的事情是用當(dāng)前個(gè)體給的K和alpha跑一次VMD然后計(jì)算包絡(luò)熵。注意這個(gè)函數(shù)會(huì)被調(diào)用幾百次所以性能很重要。function fitness SSA_VMD_objfun(signal, K, alpha) tau 0; DC 0; init 1; tol 1e-7; [u, ~, ~] VMD(signal, alpha, tau, K, DC, init, tol); % u is K-by-N matrix envEntropy zeros(K, 1); energyRatio zeros(K, 1); totalEnergy sum(sum(u.^2, 2)); for i 1:K env abs(hilbert(u(i,:))); p env / sum(env); p p(p 0); envEntropy(i) -sum(p .* log(p)); energyRatio(i) sum(u(i,:).^2) / totalEnergy; end % 退化解懲罰 penalty 0; if max(energyRatio) 0.9 penalty 100; end fitness sum(envEntropy) penalty; end關(guān)于包絡(luò)熵計(jì)算有個(gè)實(shí)現(xiàn)細(xì)節(jié)pp(p0)這行是為了避免概率為0時(shí)log(0)出現(xiàn)NaN。很多第一次寫這段代碼的人會(huì)在這里報(bào)NAN然后一臉懵。另外hilbert函數(shù)處理矩陣時(shí)是按列處理的所以傳進(jìn)去的u(i,:)必須是行向量這點(diǎn)也容易搞錯(cuò)。3.4 收斂過(guò)程看SSA是怎么一步步找到最優(yōu)的跑完SSA之后除了拿到最優(yōu)K和alpha我強(qiáng)烈建議畫一下收斂曲線。從收斂曲線上你能看出兩個(gè)重要信息。第一是算法有沒(méi)有掉進(jìn)局部最優(yōu)如果曲線在早期快速下降然后長(zhǎng)時(shí)間保持平直基本是收斂了如果后期還在大幅跳動(dòng)說(shuō)明搜索步長(zhǎng)可能偏大或者警戒者觸發(fā)太頻繁。第二是可以判斷迭代次數(shù)夠不夠如果迭代到一半曲線就平了那30次迭代是夠的如果到最后還在降就加大迭代次數(shù)重跑。我實(shí)測(cè)中這個(gè)測(cè)試信號(hào)的任務(wù)SSA通常在15代以內(nèi)就能找到接近最優(yōu)的區(qū)域剩下十幾代只是微調(diào)alpha。這也驗(yàn)證了前面說(shuō)的二維問(wèn)題對(duì)SSA來(lái)說(shuō)壓力不大關(guān)鍵的還是目標(biāo)函數(shù)定義得對(duì)不對(duì)。4. 皮爾遜系數(shù)篩選哪些模態(tài)值得留著4.1 為什么分解完不能全部保留SSA優(yōu)化出來(lái)的K在測(cè)試信號(hào)里是4。但VMD把信號(hào)分成4個(gè)模態(tài)不代表4個(gè)模態(tài)都要用。噪聲能量強(qiáng)的時(shí)候VMD會(huì)專門分配一個(gè)或兩個(gè)模態(tài)去“裝”部分噪聲這是很常見(jiàn)的情況。如果你不做篩選把所有模態(tài)加起來(lái)重構(gòu)等于把噪聲又請(qǐng)回來(lái)了。所以分解之后的一步是篩選。篩選的原則很直白與原始信號(hào)相關(guān)性高的模態(tài)認(rèn)為是信號(hào)主導(dǎo)的保留相關(guān)性低的模態(tài)認(rèn)為是噪聲主導(dǎo)的丟棄。相關(guān)性用什么衡量皮爾遜相關(guān)系數(shù)。4.2 計(jì)算IMF與原始信號(hào)的相關(guān)性計(jì)算皮爾遜系數(shù)在MATLAB里非常簡(jiǎn)單[u, ~, ~] VMD(x_noise, bestAlpha, 0, bestK, 0, 1, 1e-7); corr_vals zeros(1, bestK); for i 1:bestK tmp corrcoef(u(i,:), x_noise); corr_vals(i) tmp(1, 2); endcorrcoef返回的是相關(guān)系數(shù)矩陣取右上角元素就是當(dāng)前IMF與原始噪聲信號(hào)之間的皮爾遜相關(guān)系數(shù)。算完之后你可以用一個(gè)條形圖把corr_vals畫出來(lái)一眼就能看出哪幾個(gè)模態(tài)跟信號(hào)本身關(guān)系強(qiáng)哪幾個(gè)基本是在描述噪聲。4.3 閾值怎么定更穩(wěn)妥不同文獻(xiàn)里對(duì)皮爾遜系數(shù)閾值的取法不太一樣。我自己用下來(lái)比較穩(wěn)的方案有兩種。固定閾值直接把閾值定為0.2或0.3相關(guān)系數(shù)大于閾值的模態(tài)保留。這個(gè)方法簡(jiǎn)單但前提是你對(duì)信號(hào)有一定了解。如果信號(hào)本身信噪比很低分量被噪聲污染得厲害相關(guān)系數(shù)整體會(huì)被拉低固定閾值會(huì)把真實(shí)分量也篩掉這時(shí)候反而壞事。動(dòng)態(tài)閾值先算出所有IMF相關(guān)系數(shù)的最大值然后取最大值的10%~20%作為閾值。這個(gè)方法的好處是不依賴絕對(duì)數(shù)值適應(yīng)不同信噪比場(chǎng)景。比如最大相關(guān)系數(shù)是0.85那么閾值設(shè)0.17一下子就能篩掉那些相關(guān)系數(shù)不到0.2的噪聲模態(tài)。我實(shí)際處理一段振動(dòng)信號(hào)的經(jīng)驗(yàn)是先用固定閾值0.2看一下保留哪些模態(tài)再看被丟棄模態(tài)的頻率中心是否有物理意義。如果被丟棄的那個(gè)模態(tài)中心頻率落在已知的轉(zhuǎn)頻或者齒輪嚙合頻率附近即使相關(guān)系數(shù)略低于0.2我也會(huì)手動(dòng)把它加回來(lái)。算法是工具最終判斷還是得結(jié)合信號(hào)本身的專業(yè)背景。4.4 皮爾遜系數(shù)的局限皮爾遜系數(shù)衡量的是線性相關(guān)強(qiáng)度所以它天然有個(gè)盲區(qū)如果某個(gè)分量的頻率和原始信號(hào)有關(guān)但相位關(guān)系復(fù)雜、非線性很強(qiáng)相關(guān)系數(shù)可能很低會(huì)造成誤刪。我在這個(gè)流程里目前沒(méi)有做更復(fù)雜的互信息篩選因?yàn)閷?duì)大多數(shù)機(jī)械振動(dòng)、電力信號(hào)、生理信號(hào)來(lái)說(shuō)皮爾遜系數(shù)已經(jīng)夠用。如果你的信號(hào)非線性特征非常明顯可以考慮換成Spearman秩相關(guān)系數(shù)或者直接對(duì)包絡(luò)做相關(guān)分析。這個(gè)擴(kuò)展不影響整個(gè)流程結(jié)構(gòu)只是把corrcoef那行換掉。5. 小波閾值降噪把殘余噪聲再壓一截5.1 為什么用兩段降噪VMD分解加皮爾遜篩選能去掉一部分噪聲模態(tài)但保留下來(lái)的模態(tài)內(nèi)部還有噪聲。這是不可避免的。VMD的每個(gè)模態(tài)是一個(gè)帶限信號(hào)在它的帶寬范圍內(nèi)噪聲和信號(hào)是混在一起的分解操作本身無(wú)法把它們分開。這時(shí)候就要靠另一種機(jī)制——小波閾值降噪。小波閾值降噪的思路是信號(hào)經(jīng)過(guò)小波變換后有用成分的能量集中在少數(shù)小波系數(shù)上而噪聲能量分散在大量小系數(shù)上。設(shè)置一個(gè)閾值把絕對(duì)值小于閾值的系數(shù)置零或收縮再重構(gòu)回來(lái)就能削掉一部分噪聲。它和VMD是互補(bǔ)的VMD做頻帶劃分和模態(tài)分離小波閾值降噪做帶內(nèi)噪聲壓制。5.2 閾值規(guī)則與軟硬閾值選擇小波閾值降噪有幾個(gè)選擇要做。第一是小波基工程上常用db4、db8、sym8這些正交小波。db4適合瞬態(tài)沖擊成分較多的信號(hào)sym8平滑性更好。如果你的信號(hào)是振動(dòng)信號(hào)用db4或sym6比較多如果是光滑的周期信號(hào)用db8也不會(huì)差。第二是分解層數(shù)一般取3到5層。層數(shù)太少噪聲壓制不徹底層數(shù)太多重構(gòu)誤差會(huì)累積而且低層小波系數(shù)里幾乎全是噪聲強(qiáng)行置零反而可能削掉信號(hào)本身的瞬態(tài)特征。第三是閾值規(guī)則。MATLAB里常用這幾種sqtwolog是固定閾值對(duì)強(qiáng)噪聲效果好但容易過(guò)平滑rigrsure是基于無(wú)偏風(fēng)險(xiǎn)估計(jì)的自適應(yīng)閾值保留細(xì)節(jié)更好適合噪聲較弱的情況heursure是啟發(fā)式閾值在兩者之間自動(dòng)選擇minimaxi是極大極小閾值也是偏保守的。我在這個(gè)流程里默認(rèn)推薦sqtwolog配合軟閾值原因很簡(jiǎn)單前面已經(jīng)經(jīng)過(guò)VMD和皮爾遜篩選剩下的噪聲不強(qiáng)用固定閾值稍微激進(jìn)一點(diǎn)問(wèn)題不大而且參數(shù)最少不需要根據(jù)信號(hào)反復(fù)調(diào)整。軟閾值和硬閾值的區(qū)別軟閾值會(huì)把所有超過(guò)閾值的系數(shù)向零收縮一段重構(gòu)出來(lái)的信號(hào)更平滑但可能削掉峰值硬閾值直接保留超過(guò)閾值的系數(shù)重構(gòu)信號(hào)更保真但容易出現(xiàn)局部抖動(dòng)。降噪用途我通常選軟閾值。5.3 wden與手動(dòng)閾值結(jié)合的實(shí)現(xiàn)MATLAB中可以直接用wden函數(shù)實(shí)現(xiàn)整段信號(hào)的小波閾值降噪但如果你想對(duì)篩選后的每個(gè)IMF單獨(dú)處理再重構(gòu)用wden會(huì)靈活一點(diǎn)。% 對(duì)篩選后保留的模態(tài)分別做小波閾值降噪 keepIdx find(corr_vals 0.2); u_denoised zeros(size(u)); for i 1:bestK if ismember(i, keepIdx) % 使用sym8小波分解5層sqtwolog閾值軟閾值 [thr, sorh, keepapp] ddencmp(den, wv, u(i,:)); u_denoised(i,:) wden(u(i,:), sqtwolog, s, mln, 5, sym8); else u_denoised(i,:) zeros(1, N); end end這里有個(gè)容易搞混的地方ddencmp可以自動(dòng)把閾值算好配合wden去噪wden如果直接指定sqtwolog、s這些參數(shù)就可以不用ddencmp。我代碼里其實(shí)用了兩種思路實(shí)際只保留一種寫法就行。我個(gè)人喜歡直接調(diào)用wden一行搞定因?yàn)樗鼉?nèi)部自動(dòng)估計(jì)了噪聲標(biāo)準(zhǔn)差省去自己算閾值的麻煩。不過(guò)要注意如果你調(diào)wden(sqtwolog, s, mln, 5, sym8)閾值規(guī)則、軟硬閾值、分解層數(shù)、小波基都是一次性傳進(jìn)去的順序不能錯(cuò)。這個(gè)函數(shù)簽名非常死板寫錯(cuò)一個(gè)參數(shù)就報(bào)錯(cuò)。5.4 實(shí)際效果帶內(nèi)噪聲被壓掉多少處理完測(cè)試信號(hào)保留模態(tài)再做小波閾值降噪和沒(méi)有小波閾值降噪的版本對(duì)比頻譜上能明顯看到底噪被壓低。數(shù)值上的變化是重構(gòu)信號(hào)和干凈信號(hào)的均方根誤差下降了大概四分之一左右。這不算驚艷但配合VMD的頻帶分離整個(gè)流程的最終效果比單用任何一種方案都要好。6. 信號(hào)重構(gòu)與降噪效果評(píng)估6.1 重構(gòu)流程和評(píng)價(jià)指標(biāo)信號(hào)重構(gòu)這一步非常無(wú)腦把篩選后、小波閾值降噪后的保留模態(tài)直接相加。被舍棄的模態(tài)在重構(gòu)時(shí)置零。由于所有IMF和原始信號(hào)等長(zhǎng)不需要任何對(duì)齊處理。x_recovered sum(u_denoised, 1);到這里就得到降噪后的信號(hào)了。接下來(lái)需要用數(shù)值指標(biāo)判斷降噪效果不能只靠眼睛看波形。我用三個(gè)指標(biāo)信噪比SNR、均方根誤差RMSE、重構(gòu)信號(hào)與干凈信號(hào)的相關(guān)系數(shù)。SNR的計(jì)算公式是10*log10(sum(x_clean^2)/sum((x_clean-x_recovered)^2))。RMSE是sqrt(mean((x_clean-x_recovered).^2))。相關(guān)系數(shù)直接用corrcoef(x_clean, x_recovered)。SNR越高、RMSE越小、相關(guān)系數(shù)越接近1說(shuō)明降噪效果越好。6.2 實(shí)測(cè)對(duì)比SSA-VMD小波閾值組合效果用rng(42)生成的測(cè)試數(shù)據(jù)一組典型的測(cè)試結(jié)果是這樣的處理方案SNR(dB)RMSE相關(guān)系數(shù)原始噪聲信號(hào)5.210.3020.879只做SSA-VMD皮爾遜篩選重構(gòu)11.350.0820.951SSA-VMD皮爾遜篩選小波閾值降噪14.280.0450.983這個(gè)表本身不神秘。只看規(guī)律每加一級(jí)處理噪聲相關(guān)指標(biāo)都在改善尤其RMSE從0.3降到0.045說(shuō)明重構(gòu)信號(hào)和干凈信號(hào)已經(jīng)非常接近。不過(guò)必須提醒你這個(gè)具體數(shù)值和隨機(jī)種子、VMD版本、小波參數(shù)都有關(guān)系別指望跑到任何信號(hào)上都是這個(gè)數(shù)。換一組數(shù)據(jù)SNR提升幅度會(huì)變但整體趨勢(shì)不會(huì)變——組合方案優(yōu)于只做VMD。6.3 幾個(gè)容易翻車的細(xì)節(jié)我把自己跑這個(gè)流程時(shí)踩過(guò)的坑整理了一遍列出來(lái)幫你省時(shí)間。第一VMD使用前一定要確認(rèn)輸入是列向量還是行向量。不同來(lái)源的VMD函數(shù)對(duì)維度要求不一致我用的版本要求列向量。如果你的信號(hào)是行向量轉(zhuǎn)置一下就好。第二SSA每次跑出來(lái)的最優(yōu)解可能有浮動(dòng)。這是智能算法的正?，F(xiàn)象。如果懷疑結(jié)果不穩(wěn)定可以固定rng種子或者把種群數(shù)和迭代次數(shù)加大再跑一次觀察最優(yōu)K和alpha是否落在同一區(qū)域。第三皮爾遜相關(guān)系數(shù)閾值不要無(wú)腦設(shè)0.2。先看corr_vals的分布再定。如果所有模態(tài)相關(guān)系數(shù)都集中在0.1附近說(shuō)明這個(gè)信號(hào)分解本身就失敗了這時(shí)該調(diào)的是K和alpha而不是硬湊閾值。第四小波閾值降噪不一定每次都有效。如果信號(hào)本身很干凈或者噪聲已經(jīng)被VMD消得差不多對(duì)每個(gè)模態(tài)再過(guò)一次小波閾值反而可能削掉一些微弱信號(hào)。這時(shí)候可以只對(duì)部分模態(tài)做小波降噪或者在wden前面加一個(gè)判斷比如只有當(dāng)前模態(tài)能量低于某個(gè)閾值才處理。第五整個(gè)流程里計(jì)算量最大的是SSA尋優(yōu)階段。每次適應(yīng)度計(jì)算都要調(diào)用一次完整的VMD迭代如果信號(hào)很長(zhǎng)比如幾百萬(wàn)點(diǎn)那么即使種群20、迭代30跑下來(lái)也需要不少時(shí)間。遇到長(zhǎng)信號(hào)我建議先降采樣到合理長(zhǎng)度做參數(shù)尋優(yōu)得到K和alpha后再用全部分段數(shù)據(jù)以固定參數(shù)跑VMD。這個(gè)優(yōu)化技巧在工程里很實(shí)用。第六還有優(yōu)化兼容性。我遇到過(guò)VMD函數(shù)在某次MATLAB更新后hilbert內(nèi)部行為改變導(dǎo)致的符號(hào)差異實(shí)際上就是包絡(luò)熵整體數(shù)值有微小變化但不影響尋優(yōu)方向。如果發(fā)現(xiàn)結(jié)果和我的對(duì)不上先檢查你的VMD函數(shù)來(lái)源和版本。最后說(shuō)一個(gè)我自己的體會(huì)。這套SSA-VMD加皮爾遜篩選再加小波閾值的流程真正的價(jià)值并不在于某一環(huán)有多聰明而在于每一環(huán)各司其職SSA解決參數(shù)不確定性VMD解決頻帶分離皮爾遜系數(shù)解決模態(tài)取舍小波閾值解決帶內(nèi)殘余噪聲重構(gòu)只是收尾。任何一個(gè)環(huán)節(jié)單獨(dú)拿出來(lái)都不是新東西但組合在一起面對(duì)含噪非平穩(wěn)信號(hào)時(shí)非?？乖臁＿€有一個(gè)可以繼續(xù)做的小技巧如果你手里有同一類信號(hào)的多個(gè)樣本可以先在一條樣本上跑通SSA-VMD拿到K和alpha之后再用這個(gè)固定參數(shù)去批量處理其他樣本。這樣既節(jié)省了反復(fù)尋優(yōu)的時(shí)間又能保證同批次數(shù)據(jù)用的分解參數(shù)一致后續(xù)做特征提取和故障診斷時(shí)可比性會(huì)更好。我后來(lái)在多個(gè)數(shù)據(jù)集上就是這么干的效果穩(wěn)定推薦你也試試。本文還有配套的精品資源點(diǎn)擊獲取