管道中多个泄漏的识别：线性化模型、最大似然和超分辨率定位研究（Matlab代码实现）

 💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

文献来源：

摘要：

本文考虑了基于逆瞬态波理论的充水管道中多次泄漏的识别问题。该问题的解析解决方案涉及各种泄漏之间的非线互项。本文从分析和数值上表明，这些非线性项的泄漏大小与幂二和的量级相同;因此，可以忽略不计。由于这种简化，制定并测试了分别识别泄漏位置和泄漏大小的最大似然 （ML） 方案。结果表明，ML估计方案在噪声方面具有高效和鲁棒性。此外，ML方法是一种超分辨率泄漏定位方案，因为它的可解析泄漏距离（大约0.15λ，其中λ是最小波长）低于奈奎斯特-香农采样定理极限（0.5λ）。此外，还推导了 Cramér-Rao 下界 （CRLB） 并用于显示 ML 方案估计的效率。ML估计器的方差近似于CRLB，证明ML方案属于泄漏定位方法的最佳无偏估计器类。

关键词：

泄漏检测 缺陷检测 瞬态波 最大似然性 超分辨率 Cramér-Rao下限 

📚2 运行结果

部分代码：

% this figure shows the accuracy of linearization
figure; 
plot(omega/omega_th,abs(h_x2)); hold on;
plot(omega/omega_th,abs(h_x2_approximation),'r--'); 
xlabel('\omega/\omega_{th}','Interpreter','Tex'); ylabel('Head magnitude [m]');
legend('Transfer matrix','Linearized model')
title('Frequency Response Diagram');

h_NL = zeros(length(omega),M);   % Eq. (14)
for mm = 1:M
    h_NL(:,mm) = -Z.'.*sinh(xM(mm)*mu.').*q_x1 + cos(xM(mm)*mu.').*h_x1;
end

noise_method = 2;
SNR = 10;
sigma = 1/(10^(SNR/20));

h = h_hydro - h_NL;
h_MFP = real(h)+sigma*35*randn(size(h)) + 1j*imag(h)+sigma*13*randn(size(h));
h_MFP = h_MFP(:);
h_xM_up_esti = real(h_xM_up)+sigma*mean(mean(abs(real(h))))*randn(size(h_xM_up))+1j*imag(h_xM_up)+sigma*mean(mean(abs(imag(h))))*randn(size(h_xM_up));

q_x1_esti = (-h_xM_up_esti.'./(Z.*sinh(xM_up*mu))).';  % Eq. (17)

%% leak localization using MFP (1D search)
% For MFP, please refer to: X. Wang and M. S. Ghidaoui, Pipeline leakage detection using the matched-field processing method, Journal of Hydraulic Engineering, 144(6), 04018030, 2018

xL_test = 1:1:L-1; % discrete points for plotting cost function along the pipe
cost_MFP = zeros(size(xL_test));
S_MFP = zeros(size(xL_test));
G_all = zeros(length(xL_test),length(omega)*M);
for xx = 1:length(xL_test)
    G = zeros(length(omega),M);
    for mm = 1:M
        if xL_test(xx)>xM(mm)
            G(:,mm) = zeros(J,1);
        else
            G(:,mm) = -sqrt(9.8/2/(H1-(H1-H2)/L*xL_test(xx)))*Z.'.*sinh((xM(mm)-xL_test(xx))*mu.').*(Z.'.*sinh(xL_test(xx)*mu.').*q_x1_esti);
        end
    end
    G = G(:);
    G_all(xx,:) = G;
%     if noise_method == 3   % whitening blue noise
%         G = 1./sigma_blue'.*G;
%     end
    cost_MFP(xx) = h_MFP'*G*G'*h_MFP/(G'*G);   % Eq. (19) in the paper
    S_MFP(xx) = G'*h_MFP/(G'*G);    % Eq. (20) in the paper
end
CdAl_MFP = S_MFP;
[~,x_max] = max(cost_MFP);
if no_L==1
    disp(['Conventional MFP: the estimated leak position is ',num2str(xL_test(x_max)),' m and the estimated leak size is CdAl=',num2str(abs(CdAl_MFP(x_max)))]);
    disp(['Conventional MFP: the localization error is ',num2str(abs(xL_test(x_max)-xL(1))),' m ']);
end

figure;
plot(xL_test,abs(cost_MFP)); hold on; % plot(xL_test,abs(S_MFP),'r--');
for ii = 1:length(xL)
    plot([xL(ii) xL(ii)],[min(abs(cost_MFP)) 1.1*max(abs(cost_MFP))],'r--');
end
plot(xM,ones(size(xM))*min(abs(cost_MFP)),'rx','MarkerSize',8,'Linewidth',1.5);
plot(xM_up,min(abs(cost_MFP)),'rx','MarkerSize',8,'Linewidth',1.5);
xlabel('x [m]'); ylabel('|B|^2','Interpreter','Tex')
ylim([min(abs(cost_MFP)) 1.1*max(abs(cost_MFP))]);

🎉3 参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。

🌈4 Matlab代码实现