Patch近场声全息-声品质客观参量三维分布可视化方法

摘要

本发明公开了一种基于Patch近场声全息技术的声品质客观参量三维可视化方法：在大型结构声源辐射形成的存在来波和去波的声场中，利用小型麦克风阵列对大型复杂声源表面进行多次扫描测量，由本发明公开的Patch近场声全息-声品质耦合矩阵计算模型，得到声学量和声品质客观参量在整个空间声场中的分布信息，实现了声学物理量和声品质客观参量三维分布可视化。

说明书

技术领域：

本发明涉及大型复杂结构噪声源识别定位技术与声品质评价，更具体的说是基于Patch近场声全息理论的声品质客观参量三维分布可视化方法。 

背景技术：

近场声全息方法NAH(Near‐field Acoustic Holography)是利用靠近声源的全息测量面获得的声学近场信息，通过各种数学变换重构出三维空间内的声场物理信息(声压、振速等)的方法。NAH方法通过近场测量捕捉到声场中的倏逝波，得到物体的振动细节，理论上可以获得不受波长λ限制的高分辨率，极大地提高了声源识别定位的精度。目前，NAH方法已被广泛地应用于声源的识别和定位中，尤其适用于低频声源的特性判别、声散射测量和结构振动引起的声辐射测量等场合。但是，常规近场声全息方法要求麦克风阵列测量孔径至少是声源面积的两倍以上。因此，如果要测量与分析大尺寸以及复杂结构声源，就需要配置大量的麦克风，使得测量所需的设备价格昂贵。例如对0.4×0.6m²大小的板声源在40dB信噪比条件下采用傅立叶变换近场声全息算法重构辐射声场，若要获得0.05m分辨率，需要采用882个麦克风来测量，这在工程实际测量中基本无法实现，麦克风数量限制了近场声全息方法在大型复杂结构声源测量等场合的推广应用。 

而Patch近场声全息算法很好的克服这一问题。它通过将大尺寸及复杂结构声源分解为多个局部声源，采用同一麦克风阵列对各个局部声源依次进行扫描测量，通过Patch声全息孔径合成算法重构得到整个声源的声学量的空间分布，突破了常规近场声全息对于全息测量孔径的限制。 

另一方面，研究人员发现噪声对于人的影响不仅与声压级有关，还与声音自身频率的组成、人耳听觉系统的物理特性以及人的心理特性有关。相同声压级的两种声音由于各自频率组成的不同，会导致在人心理感觉上响度的巨大差异。因此需要引入能够反映人对声音主观感受的量‐声品质，作为对噪声评价的参考。Blauet给出的声品质定义为：“声品质”是指人通过人耳对声音事件的听觉感知，在主观上做出喜好判断的过程。而声品质好坏的描述可由一系列反应了噪声对人主观情绪上影响程度的指标(如尖锐度、粗糙度、响度等)来表示。 

综上所述，使用近场声全息方法来计算分析噪声能给出声压、声强等声学量的空间分布，但是近场声全息方法的分析过程没有把人的主观感受考虑在内，计算结果没有给出响度、尖锐度等声品质客观参量的空间分布结果，不能实现对人主观听觉感受影响最大声源位置的定位。而现有声品质分析方法都是针对若干个麦克风测得的声场信息进行计算、分析和评价，使用现有的声品质分析方法只能得到该空间局部测试点位置的声品质客观参量信息，不能作为整个空间声品质好坏的评价标准，反映不了人的主观听觉感受与噪声源空间分布之间的相互关系，无法获得对人主观听觉感受影响最大的声源位置信息。 

课题组先前研发的基于球面近场声全息方法的声品质三维分布可视化方法(申请号：201310261258.2)主要用于封闭空间声场(如汽车车厢内部)等声源的初步识别定位。在相同传声器数量和相同声场条件下，球面近场声全息方法的声源定位精度没有Patch声全息方法的精度高。同时Patch声全息也不能用于狭小封闭空间内部的噪声源定位，两种方法的声场重构的物理模型和数学公式均不同，是应用于不同类型声场测量分析的不同方法。 

因此，本发明提出一种基于Patch近场声全息和声品质联合(Patch Near‐field Acoustic Holography‐Sound Quality,Patch NAH‐SQ)分析方法。该方法突破了常规近场声全息对于全息测量孔径的限制，能够实现大尺寸以及复杂结构声源的声品质客观参量的空间分布可视化，揭示了主观听觉感受与大型复杂结构噪声源分布之间的相互关系，并高精度定位出了与人听觉感受最为密切的噪声源所在的位置。 

发明内容：

本发明要克服声全息阵列测量孔径的限制，实现大型复杂结构声源的测量分析与定位，以及单独使用声品质客观参量计算评价声场时不能给出整个空间中的声品质客观参量分布信息，以及单独采用近场声全息方法分析空间声场分布时没有考虑人的主观听觉因素的缺陷，实现根据人的主观听觉感受识别定位噪声源的目的。 

本发明提出一种基于Patch近场声全息方法的声品质客观参量的三维分布可视化方法，其计算结果实现了整个大型复杂结构声源表面 及其辐射声场内中低频噪声的可视化，提供了声品质客观参量的三维分布信息，并给出影响人主观听觉感受最大的噪声源位置信息，从而为大型结构的减振降噪提供最直接的指导。 

本发明所述的基于Patch近场声全息方法的声品质客观参量的三维可视化方法按如下步骤进行： 

1.采用同一个小尺寸高密度传声器阵列对大型复杂结构声源进行多次近场扫描测量，获得整个声源面的近场声信号信息。在每次测量的同时，记录参考点位置信号，作为每次测量参考点信号。 

2.根据1测量到的近场声信号信息，首先通过与参考麦克风信号进行互相干计算，获得所有测量面同步的麦克风复声压信号，然后采用Patch近场声全息方法获得声场声学量的三维空间分布信息(声压、粒子速度等)，计算结果以3D图像的形式给出。 

Patch近场声全息方法定义不同测量面的声压用如下公式计算： 

复声压相位：$\theta \left(\omega \right)=\arg \left[\frac{S_{r,h}\left(\omega \right)}{S_{r,r}\left(\omega \right)}\right]---\left(1\right)$

复声压幅值：$p\left(\omega \right)=\sqrt{S_{h,h}\left(\omega \right)}---\left(2\right)$

复声压大小：p(r)=p(ω)·e^jθ(ω)   (3) 

式中：S代表相干计算，下标r,h代表参考点和阵列测量面上麦克风信号。 

重构面的声场变换公式如下： 

p(r)=p^T(A⁺A+θ²I)^-1A⁺α   (4) 

式中的A⁺A和A⁺α由下式确定： 

${\left[A^{+}A\right]}_{{\mathrm{nn}}^{\prime }}=\frac{1}{{\mathrm{\pi k}}^2}{\int \int }_{-\infty }^{\infty }{e}^{j(k_z^{*}-k_z)d}{e}^{-j{K}_{\mathrm{xy}}(r_{\mathrm{xy},n}-r_{\mathrm{xy},n^{\prime }})}{\mathrm{dk}}_x{\mathrm{dk}}_y---\left(5\right)$

${\left[A^{+}\alpha \right]}_n=\frac{1}{{\mathrm{\pi k}}^2}{\int \int }_{-\infty }^{\infty }{e}^{j(k_z^{*}d-k_{z}z)}{e}^{-j{K}_{\mathrm{xy}}\left(r_{\mathrm{xy}}{-r}_{\mathrm{xy},n}\right)}{\mathrm{dk}}_x{\mathrm{dk}}_y---\left(6\right)$

式中：p为测量面上所有测量点的声压组成的列向量，p^T为p的转置；r为空间任意重构点的笛卡尔坐标，r=(x,y,z)，p(r)为重构点的声压；k为波数，k_z为重建面上z方向上的空间波数，*代表复数共轭；K_xy=(k_x,k_y)，k_x为重建面上x方向上的空间波数，k_y为重建面上y方向上的空间波数；r_xy=(x,y)为重构点的xy坐标，r_xy,n=(xn,yn)，r_xy,n′=(x_n′,y_n′)，r_xy,n和r_xy,n′分别表示麦克风测量点n和n′的xy坐标；d为测量面到声源面的距离。 

3.根据2中得到的空间中的声场信息(声压、粒子速度等)以及本发明公开的空间中声场信息与声品质客观参量算法的耦合计算模型，得到空间各位置处的声品质信息(如响度、尖锐度等)三维分布信息，实现声品质客观参量的三维分布可视化，并据此识别出与人主观听觉感受密切相关声源的所在位置。 

声品质响度的计算模型：根据Moore的响度定义，心理声学客观参量的响度计算以ERB(Equivalent Rectangular Width—等效矩形宽度)尺度为基础，在人耳0.05kHz～15kHz的听觉频率范围内建立了372个听觉滤波器，其中临界频带宽度ERB和频率f(kHz)的对应关系可近似为： 

ERB=24.673(0.004368f+1)   (7) 

式中：f为带宽的中心频率。 

0.05kHz～15kHz频率范围内人耳听觉滤波器的中心频率可由下式得到： 

ERB-number=21.366lg(0.004368f+1)   (8) 

式中：f为滤波器的中心频率。 

而滤波器的输出激励可由下式求得： 

$E_i=\frac{\mathrm{\Sigma W}\left(g_{\mathrm{ij}}\right)P_j^2}{P_0^2}{E}_0---\left(9\right)$

式中：E_i为第i滤波器的输出激励，W(g_ij)为第i个滤波器对频率j处输入的响应，为输入信号的有效值功率，P₀为参考声压2×10^-5帕。在获得滤波器输出信号基础上，可求得特征响度N'。 

最终响度是对372个滤波器求得的特征响度之和，如下式： 

$N=\underset{i=1}{\overset{372}{\Sigma }}{N}_i^{\prime }---\left(10\right)$

尖锐度模型：在尖锐度计算过程中，尖锐度用S表示，其计算公式为： 

$S=0.1043\times \frac{{\int }_0^{24}{N}^{\prime }g\left(z\right)\mathrm{dz}}{{\int }_0^{24}{N}^{\prime }\mathrm{dz}}---\left(11\right)$

式中：N'为Zwicker响度模型中的特征响度，而 

$g\left(z\right)=\left(\begin{array}{cc}1& (z<16)\\ 0.066e^{0.171z}& (z\ge 16)\end{array}\right)---\left(12\right)$

根据2中计算得到的声场中声压的空间分布结果，并结合空间中单点声品质客观参量的计算模型，建立空间中声压场与响度场的耦合三维矩阵映射模型，即： 

${{[p(x,y,z,{\omega }_1)··p(x,y,z,{\omega }_m)]}_{1\times m}\left(\begin{array}{ccc}{w}_1& ··& w_{372}\\ w_1& ··& w_{372}\\ ·& ··& ·\\ w_1& ··& w_{372}\end{array}\right)}_{m\times 372}={[N_1^{\prime }(x,y,z)··N_{372}^{\prime }(x,y,z)]}_{1\times 372}---\left(13\right)$

或简写为： 

P_iW=N′   (14)式中：p_i为声场中某位置点处任一频率下的声压，P_i为声场中某位置处各频率下的声压重构值所组成的矢量，W为由372个滤波器w_i组成的听觉滤波器矩阵，表示人耳对可听频带内所有频率的响应，N′为特征响度矢量。由公式(10)对特征响度矢量N′中的各项求和便可获得声场指定点的响度，而对声场三维空间节点重复这一计算过程便可获得声场响度三维分布结果。其它声品质客观参量三维分布结果也可采用与响度计算类似的分析流程获得。 

通过上述方法能够计算出声场各个位置的声品质客观参量值，并以三维图像形式给出其空间分布，实现了声品质客观参量的空间分布可视化，识别出对人主观听觉感受影响最大的噪声源。 

与以往技术相比，本发明的有益效果： 

1.本发明通过小型高密度阵列面对大型复杂结构声源进行多次近场扫描测量，突破原先近场声全息阵列孔径的限制，实现大型复杂结构声场的重构。 

2.本发明能获得不受声波波长限制的高精度声场重构结果，并能获得含粒子振速、声强三维分布结果等的所有声场信息。 

3.与现行的声品质单点测量技术相比，本发明通过小型高密度 阵列面多次近场扫描测量，给出整个空间的声品质客观参量的三维空间分布图，使工程师能够“看到”声品质客观参量的空间分布。 

4.本发明可以识别出大型复杂结构声源中对人主观听觉感受影响最大的噪声源位置信息，为大型结构减振降噪工作提供直接指导。 

附图说明：

图1小孔径平面麦克风阵列示意图 

图2麦克风阵列测量布置及双声源声场示意图 

图3(a)声源1(800Hz,79dB)和声源2(1000Hz,84dB)两个点声源声场的响度和声压计算结果对比 

图3(b)声源1(800Hz,79dB)和声源2(1000Hz,84dB)两个点声源声场的尖锐度和声压计算结果对比 

图4(a)声源1(800Hz,82dB)和声源2(1000Hz,84dB)两点声源声场的响度和声压计算结果对比 

图4(b)声源1(800Hz,82dB)和声源2(1000Hz,84dB)两点声源声场的尖锐度和声压计算结果对比 

具体实施方式：

参照附图： 

本发明按如下步骤进行： 

1.采用同一个小尺寸高密度传声器阵列对大型复杂结构声源进行多次近场扫描测量，获得整个声源面的近场声信号信息。在每次测量的同时，记录参考点位置信号，作为每次测量参考点信号。 

2.根据1测量到的近场声信号信息，首先通过与参考麦克风信号进行互相干计算，获得所有测量面同步的麦克风复声压信号，然后采用Patch近场声全息方法获得声场声学量的三维空间分布信息(声压、粒子速度等)，计算结果以3D图像的形式给出。 

Patch近场声全息方法定义不同测量面的声压用如下公式计算： 

复声压相位：$\theta \left(\omega \right)=\arg \left[\frac{S_{r,h}\left(\omega \right)}{S_{r,r}\left(\omega \right)}\right]---\left(15\right)$

复声压幅值：$p\left(\omega \right)=\sqrt{S_{h,h}\left(\omega \right)}---\left(16\right)$

复声压大小：p(r)=p(ω)·e^jθ(ω)   (17)式中：S代表相干计算，下标r,h代表参考点和阵列测量面上麦克风信号。 

重构面的声场变换公式如下： 

p(r)=p^T(A⁺A+θ²I)^-1A⁺α   (18) 

式中的A⁺A和A⁺α由下式确定： 

式中：S代表相干计算，下标r,h代表参考点和阵列测量面上麦克风信号。 

${\left[A^{+}A\right]}_{{\mathrm{nn}}^{\prime }}=\frac{1}{{\mathrm{\pi k}}^2}{\int \int }_{-\infty }^{\infty }{e}^{j(k_z^{*}-k_z)d}{e}^{-j{K}_{\mathrm{xy}}(r_{\mathrm{xy},n}-r_{\mathrm{xy},n^{\prime }})}{\mathrm{dk}}_x{\mathrm{dk}}_y---\left(19\right)$

${\left[A^{+}\alpha \right]}_n=\frac{1}{{\mathrm{\pi k}}^2}{\int \int }_{-\infty }^{\infty }{e}^{j(k_z^{*}d-k_{z}z)}{e}^{-j{K}_{\mathrm{xy}}\left(r_{\mathrm{xy}}{-r}_{\mathrm{xy},n}\right)}{\mathrm{dk}}_x{\mathrm{dk}}_y---\left(20\right)$

式中：p为测量面上所有测量点的声压组成的列向量，p^T为p的转置；r为空间任意重构点的笛卡尔坐标，r=(x,y,z)，p(r)为重构点的声压；k为波数，k_z为重建面上z方向上的空间波数，*代表复数共轭；K_xy=(k_x,k_y)，k_x为重建面上x方向上的空间波数，k_y为重建面上y方向上的空间波数；r_xy=(x,y)为重构点的xy坐标，r_xy,n=(x_n,y_n)，r_xy,n′=(x_n′,y_n′)， r_xy,n和r_xy,n′分别表示麦克风测量点n和n′的xy坐标；d为测量面到声源面的距离。 

3根据2中得到的空间中的声场信息(声压、粒子速度等)以及本发明公开的空间中声场信息与声品质客观参量算法的耦合计算模型，得到空间各位置处的声品质信息(如响度、尖锐度等)三维分布信息，实现声品质客观参量的三维分布可视化，并据此识别出与人主观听觉感受密切相关声源的所在位置。 

声品质响度的计算模型：根据Moore的响度定义，心理声学客观参量的响度计算以ERB(Equivalent Rectangular Width—等效矩形宽度)尺度为基础，在人耳0.05kHz～15kHz的听觉频率范围内建立了372个听觉滤波器，其中临界频带宽度ERB和频率f(kHz)的对应关系可近似为： 

ERB=24.673(0.004368f+1)   (21) 

式中：f为带宽的中心频率。 

0.05kHz～15kHz频率范围内人耳听觉滤波器的中心频率可由下式得到： 

ERB-number=21.366lg(0.004368f+1)   (22) 

式中：f为滤波器的中心频率。 

而滤波器的输出激励可由下式求得： 

$E_i=\frac{\mathrm{\Sigma W}\left(g_{\mathrm{ij}}\right)P_j^2}{P_0^2}{E}_0---\left(23\right)$

式中：E_i为第i滤波器的输出激励，W(g_ij)为第i个滤波器对频率 j处输入的响应，为输入信号的有效值功率，P₀为参考声压2×10^-5帕。在获得滤波器输出信号基础上，可求得特征响度N'。 

最终响度是对372个滤波器求得的特征响度之和，如下式： 

$N=\underset{i=1}{\overset{372}{\Sigma }}{N}_i^{\prime }---\left(24\right)$

尖锐度模型：在尖锐度计算过程中，尖锐度用S表示，其计算公式为： 

$S=0.1043\times \frac{{\int }_0^{24}{N}^{\prime }g\left(z\right)\mathrm{dz}}{{\int }_0^{24}{N}^{\prime }\mathrm{dz}}---\left(25\right)$

式中：N'为Zwicker响度模型中的特征响度，而 

$g\left(z\right)=\left(\begin{array}{cc}1& (z<16)\\ 0.066e^{0.171z}& (z\ge 16)\end{array}\right)---\left(26\right)$

根据b中计算得到的声场中声压的空间分布结果，并结合空间中单点声品质客观参量的计算模型，建立空间中声压场与响度场的耦合三维矩阵映射模型，即： 

${{[p(x,y,z,{\omega }_1)··p(x,y,z,{\omega }_m)]}_{1\times m}\left(\begin{array}{ccc}{w}_1& ··& w_{372}\\ w_1& ··& w_{372}\\ ·& ··& ·\\ w_1& ··& w_{372}\end{array}\right)}_{m\times 372}={[N_1^{\prime }(x,y,z)··N_{372}^{\prime }(x,y,z)]}_{1\times 372}$

                                          (27) 

或简写为： 

P_iW=N′   (28)式中：p_i为声场中某位置点处任一频率下的声压，P_i为声场中某位置处各频率下的声压重构值所组成的矢量，W为由372个滤波器w_i组成 的听觉滤波器矩阵，表示人耳对可听频带内所有频率的响应，N′为特征响度矢量。由公式(24)对特征响度矢量N′中的各项求和便可获得声场指定点的响度，而对声场三维空间节点重复这一计算过程便可获得声场响度三维分布结果。其它声品质客观参量三维分布结果也可采用与响度计算类似的分析流程获得。 

通过上述方法能够计算出声场各个位置的声品质客观参量值，并以三维图像形式给出其空间分布，实现了声品质客观参量的空间分布可视化，识别出对人主观听觉感受影响最大的噪声源。 

下面通过具体实施例子对本发明作进一步的描述。 

本实施例子中，均以小孔径平面麦克风阵列作为测量阵，如图1所示，平面麦克风阵列上均匀分布64个麦克风，麦克风之间的间距相等为0.03m，阵列孔径0.21m×0.21m。 

1.在空间中布置两个单极子点声源：目标声源1的参数设置为800Hz、79dB，放置在空间直角坐标系的x负半轴上0.12m处(‐0.12m,0,0)；目标声源2的参数设置为1000Hz、84dB，放置在空间直角坐标系的x正半轴上0.12m处(0.12m,0,0)，如图2所示。两声源间的距离为0.24m，而平面麦克风阵列孔径为0.21m，显然声源间距大于阵列孔径。Patch测量面距离声源面0.06m，且两个Patch测量面的边界刚好重合。采用本方法重构与声源面的距离为0.03m的平面处的声压、响度和尖锐度的空间分布图。 

图3(a)是800Hz、79dB的声源1和1000Hz、84dB的声源2同时 存在时的声压和响度的空间分布计算结果的比较图。图3(b)是800Hz、79dB的声源1和1000Hz、84dB的声源2同时存在时的声压和尖锐度的空间分布计算结果的比较图。 

2.图4(a)和图4(b)同样采用1所述的双声源声场模型，但声源1的参数设置为800Hz、82dB，放置在空间直角坐标系的x负半轴上0.12m处(‐0.12m,0,0)；目标声源2的参数设置为1000Hz、84dB，放置在空间直角坐标系的x正半轴上0.12m处(0.12m,0,0)，平面麦克风阵列(如图1所示)距离声源面0.06m，如图2(b)所示。采用本方法重构与声源面的距离为0.03m的平面处的声压、响度和尖锐度的空间分布图。 

图4(a)是800Hz、82dB的声源1和1000Hz、84dB的声源2同时存在时的声压和响度的空间分布计算结果的比较图。图4(b)是800Hz、82dB的声源1和1000Hz、84dB的声源2同时存在时的声压和尖锐度的空间分布计算结果的比较图。 

图3(a)是800Hz(79dB)和1000Hz(84dB)两声源同时存在时声压级和响度的空间分布图，图中声源2的声压级比声源1的声压级大，其响度值也较大；图4(a)中是800Hz(82dB)和1000Hz(84dB)两声源同时存在时声压级和响度的空间分布图，图中声压级仍然是声源1的声压级较大，但响度值却是声源2的值较大。对比图3(a)和图4(a)，可以看出当采用不同的度量识别噪声源位置时，得到的主要噪声源位置信息可能相反。因此不能根据声压级等声学物理量信息来识别内声场人听觉感受差的噪声源位置，应采用根据评价人听觉感觉的度量指标来 确定引起人恼怒的噪声源位置。 

本专利采取的方法能够获得空间中超出单个阵列测量孔径大范围区域的不同频率声场的声品质客观参量的空间分布信息。同时与传统的根据声压来识别声源的声全息方法不同，本方法能够识别出与人主观听觉感受密切相关的噪声源的位置信息。当当前时间在规定的高峰时间内，且满足以下两个条件时，则系统接受该租赁点的服务请求，并将其放入到集合A⁺中，将接受调度车辆的优先服务； 

2.3.1车锁比值α≤0.1或α≥0.9； 

2.3.2发出请求租赁点周围500m（居民可接受的最大步行距离）以内的租赁点车锁比值满足α≤0.1或α≥0.9条件的个数与该区域内租赁点的总个数的比值r≥0.5； 

2.4二级调度服务窗口A^-

当时间在预先规定的高峰时间内，且满足以下两个条件时，则系统接受该租赁点的服务请求，并将其放入到集合A^-中，准备接受调度车辆的服务； 

2.4.1车锁比值α≤0.1或α≥0.9； 

2.4.2发出请求租赁点周围500m（居民可接受的行走最大距离）以内的租赁点车锁比值α≤0.1或α≥0.9条件的个数与该区域内租赁点的总个数的比值0.3≤r<0.5； 

或者当时间在预先规定的平峰时间内时，且满足以下两个条件时，则系统接受该租赁点的服务请求，并将其放入到集合A^-中，准备接受调度车辆的服务； 

2.4.3车锁比值α≤0.2或α≥0.8； 

2.4.4发出请求租赁点周围500m以内的租赁点车锁比值满足条件α≤0.2或α≥0.8的个数与该区域内租赁点的总个数的比值r>0.5； 

2.5非调度服务窗口B 

当租赁点不满足以上条件时时，则将其放入到集合B中，该集合中的租赁点不接受服务请求。 

2.6二级调度服务窗口的租赁点插入到调度路径中需要满足这样的条件：尽量不增加车辆的调度路径，该租赁点的插入可以调整调度车辆上的自行车数量，更好的满足下一个租赁点需求条件，并减少或不增加当前使用的调度车辆数目。 

本发明的技术构思：调度指令生成之前对调度区域进行划分，避免了调度车辆行驶过程中的不合理路径；在自然调度的基础上，根据各个租赁点及区域在不同时间段的车锁比，将租赁点调度服务分为不同优先级调度，减少了调度车辆的盲目运行，同时提高了调度效率和群众满意度。 

本发明优点：可有效提高了公共自行车调度管理的科学性，并减少调度路径、节约降低调度成本，提高公共自行车利用率，缓解“租车难、还车难”的矛盾。 

图1基于自流动性的调度区域划分流程图 

参照附图： 

本发明所述的公共自行车系统的自然混合调度方法，具体步骤如下： 

步骤1、划分调度区域 

结合城市公共自行车系统的调度需求，因此可以利用公共自行车在租赁点间的这种自流动性，结合OD（起点与终点）运行数据，对城市公共自行车系统租赁点进行区域划分；进而结合各区域中租赁点潮汐方向、租赁点类型等特征属性，对区域进行二次划分，最终实现城市公共自行车系统动态调度的区域划分，是对公共自行车系统进行合理调度的前提。其流程图如图1所示： 

步骤2、生成调度指令 

2.1租赁点调度服务等级划分，针对公共自行车自流动性而导致租赁点服务需求量时候发生变化的特点，为减少调度车辆盲目的运行，同时又确保调度的效果。将公共自行车系统租赁点分为三类： 

2.3.1一级调度服务窗口A⁺，即优先享有服务的租赁点集合； 

2.3.2二级调度服务窗口A^-，即准备接受服务的租赁点集合 

2.3.3非调度服务窗口B，即尚未有服务需求的租赁点和已经接受过服务的租赁点集合。 

2.2高峰时间与平峰时间的定义，根据对杭州市公共自行车系统运行数据时空分析得出的系统运行规律，规定： 

高峰时间为：

平峰时间为：其余时间。 

2.3一级调度服务窗口A⁺

当当前时间在规定的高峰时间内，且满足以下两个条件时，则系统接受该租赁点的服务请求，并将其放入到集合A⁺中，将接受调度车辆的优先服务； 

2.3.1车锁比值α≤0.1或α≥0.9； 

2.3.2发出请求租赁点周围500m（居民可接受的最大步行距离）以内的租赁点车锁比值满足α≤0.1或α≥0.9条件的个数与该区域内租赁点的总个数的比值r≥0.5； 

2.4二级调度服务窗口A^-

当时间在预先规定的高峰时间内，且满足以下两个条件时，则系统接受该租赁点的服务请求，并将其放入到集合A^-中，准备接受调度车辆的服务； 

2.4.1车锁比值α≤0.1或α≥0.9； 

2.4.2发出请求租赁点周围500m（居民可接受的行走最大距离）以内的租赁点车锁比值α≤0.1或α≥0.9条件的个数与该区域内租赁点的总个数的比值0.3≤r<0.5； 

或者当时间在预先规定的平峰时间内时，且满足以下两个条件时，则系统接受该租赁点的服务请求，并将其放入到集合A^-中，准备接受调度车辆的服务； 

2.4.3车锁比值α≤0.2或α≥0.8； 

2.4.4发出请求租赁点周围500m以内的租赁点车锁比值满足条件α≤0.2或α≥0.8的个数与该区域内租赁点的总个数的比值r>0.5； 

2.5非调度服务窗口B 

当租赁点不满足以上条件时时，则将其放入到集合B中，该集合中的租赁点不接受服务请求。 

2.6二级调度服务窗口的租赁点插入到调度路径中需要满足这样的条件：尽量不增加车辆的调度路径，该租赁点的插入可以调整调度车辆上的自行车数量，更好的满足下一个租赁点需求条件，并减少或不增加当前使用的调度车辆数目。