1.一种基于频谱分析的城市埋地管道早期泄漏诊断方法，其特征在于：包括以下步骤S01 在设定管道压力下采集管道频谱数据，得到管道非泄漏状态下的频谱图；

S02 在同S01中的管道压力状态下，开启泄漏球阀模拟管道泄漏状态，采集泄漏状态下管道频谱数据，得到管道泄漏状态下的频谱图；

S03 改变管道工况，重复步骤S01、S02，获取多种工况下管道非泄漏和泄漏状态下的频谱数据，建立管道非泄漏和泄漏状态下的频谱库；

S04 实际泄漏检测时，将管道实际泄漏采集获得的频谱与频谱库中的频谱数据进行特征比对，判断管道是否发生泄漏；

S05 对于有泄漏的管道，关闭下游末端阀门，使管道系统发出激励响应并采集管道系统的激励压力信号；然后以时间为横坐标，以瞬态压力与激励压力平均值的比值为纵坐标得出管道系统的时域压力信号图；

S06：将采集得到的时域压力信号进行CEEMD处理，去除信号噪声并识别信号奇异值拐点的位置，最后结合次声波波速计算泄漏点位置。

2.根据权利要求1所述的基于频谱分析的城市埋地管道早期泄漏诊断方法，其特征在于：所述的管道频谱数据包括通过次声波采集仪获取管道的次声波信号z(t)，通过智能压力变送器采集管道的瞬态压力p1(t)，通过涡轮流量计读取管道的瞬态流量q1(t)。

3.根据权利要求1所述的基于频谱分析的城市埋地管道早期泄漏诊断方法，其特征在于：通过所述管道频谱数据得到管道非泄漏状态下的频谱图的方法包括：a计算管道系统的传递函数ZD：已知，一维管道的连续及动量方程为：

式中，p为压力函数；q为流量函数；x为距管道上游段的距离；t为时间； 为单位长度感抗； 为单位长度容抗； 为单位长度阻抗；g为重力加速度；A为管道横截面积；a为压力波波速；D为管道内径；f为Darcy-Weibach摩擦系数；

频域中，压力函数和流量函数可表示为：

式中，为平均压力；p1(t)为瞬态压力；为平均流量；q1(t)为瞬态流量。

将式(2)代入式(1)，则有：

式中，p1(t)＝P(x)est；q1(t)＝Q(x)est； 为动态流阻。

管道进、出口处的压力和流量函数表达式为：

式中，PU为管道入口处压力；PD为管道出口处压力；QU为管道入口处流量；QD为管道出口处流量； 为传播常数； 为特性阻抗；

联立公式(3)(4)(5)，得到管道的输入、输出阻抗：

已知管道上游端为恒压源，即管道入口处压力的波动为零，于是将 代入式(6)，求得时域中管道系统的传递函数为：ZD＝-ZCtanhγl

式中，ZD为传递函数；Zc为特性阻抗；γ为传播常数；l为管道长度；b对传递函数ZD进行傅里叶变换，得到频域中传递函数的频谱图。

4.根据权利要求1所述的基于频谱分析的城市埋地管道早期泄漏诊断方法，其特征在于：所述管道工况包括管道压力、泄漏状态下的泄漏量和泄漏位置，非泄漏状态试验的工况是不同管道压力；泄漏状态涉及的工况有：不同管道压力、不同泄漏量、不同泄漏位置，记录每次试验的具体工况及时间，并将每次试验所对应的频谱图录入频谱库，建立起管道非泄漏和泄漏状态下的频谱库。

5.根据权利要求1所述的基于频谱分析的城市埋地管道早期泄漏诊断方法，其特征在于：所述步骤S04具体为：管道实际泄漏检测时，只需采集一次频谱数据，然后进行频谱分析并将所得的频谱图与频谱库中的频谱图进行特征比对，若实际泄漏检测频谱与频谱库中非泄漏状态下的频谱图相匹配，则判定管道未发生泄漏；若实际泄漏检测频谱图与频谱库中泄漏状态下的频谱图相匹配，则判定管道发生泄漏。

6.根据权利要求1所述的基于频谱分析的城市埋地管道早期泄漏诊断方法，其特征在于：步骤S5具体包括：对于有泄漏的管道，瞬间部分关闭管道下游末端阀门至1/4位置处，使管道末端产生激励响应；然后采集管道的激励压力信号p2(t)；

采集完毕后，将管道下游末端阀门恢复至完全开启状态，待管道平稳运行后重复上述操作，共采集n次激励压力信号并取平均值，得到激励压力的平均值最后，以时间t为横坐标，以瞬态压力p1(t)与激励压力平均值 的比值为纵坐标计算得出压力信号的时域分析图，得到管道系统的时域压力信号s(t)。

7.根据权利要求6所述的基于频谱分析的城市埋地管道早期泄漏诊断方法，其特征在于：所述n次为5次。

8.根据权利要求1所述的基于频谱分析的城市埋地管道早期泄漏诊断方法，其特征在于：步骤S6利用CEEMD对时域压力信号去噪并识别奇异值拐点的位置的方法具体包括：(1)在时域压力信号s(t)中加入高斯白噪声ωi(t)，对目标信号s(t)+ωi(t)进行第1次CEEMD分解，获得第一个固有模态函数IMF1式中，E1为定义的CEEMD运算符号，其角标“1”代表第一次迭代运算；ε1为每个阶段设定的信噪比系数，其角标“1”代表第1阶段；ωi为单位方差为零均值的高斯白噪声；

(2)计算得到一阶残差r1(t)

r1(t)＝s(t)-IMF1(t)            (8)(3)对于一阶残差r1(t)中添加的白噪声信号，构造信号r1(t)+ε1E1(ωi(t))并再次进行分解，获得第二个固有模态函数IMF2；

(4)以此类推计算第k个残量rk(t)rk，k＝2,3...Krk(t)＝rk-1(t)-IMFi(t)            (10)(5)再继续分解rk(t)+εkEk(ωi(t))，i＝1,2...n，将它的第1个固有模态函数分量作为CEEMD的IMFk+1(6)继续分解直至满足结束条件，否则返回至步骤(4)～(5)，最终结果表示为：式中，IMFk(t)为分解后的各阶IMF函数；

(7)选择CEEMD分解后的前三阶信号分量重构，得到去噪、重构后的有效信号x(t)，通过有效信号x(t)奇异点的位置，确定次声波反射到下游末端的时间，再结合次声波波速得到泄漏点位置。