FCM的MATLAB实现
基于FCM聚类和漏失模拟的给水管网压力监测点布设
第４０卷 第６期Ｖｏｌ．４０ Ｎｏ．６徐志强基于ＦＣＭ聚类和漏失模拟的给水管网压力监测点布设徐志强，刘崎峰 ，刘佳峰，单黎明，裴文思（内蒙古大学生态与环境学院，内蒙古呼和浩特０１００２１）收稿日期：２０２１－０１－２５；修回日期：２０２１－０３－２１；网络出版时间：２０２２－０６－０７网络出版地址：ｈｔｔｐｓ：／／ｋｎｓ．ｃｎｋｉ．ｎｅｔ／ｋｃｍｓ／ｄｅｔａｉｌ／３２．１８１４．ＴＨ．２０２２０６０６．１６１１．０２０．ｈｔｍｌ基金项目：内蒙古自治区成果转化引导项目（ＣＧＺＨ２０１８０９５）第一作者简介：徐志强（１９９３—），男，河北承德人，硕士研究生（１３９４６４３２８７＠ｑｑ．ｃｏｍ．），主要从事智慧水务及管网漏失控制研究．通信作者简介：刘崎峰（１９７８—），男，内蒙古呼和浩特人，教授，博士生导师（ｎｄｌｑｆ＠ｉｍｕ．ｅｄｕ．ｃｎ），主要从事智慧水务及管网漏失控制、水污染控制技术研究．摘要：为了提高对供水管网压力监测的可靠性，以北方某城镇供水独立计量分区为模拟研究对象，提出一种结合模糊Ｃ均值（ＦｕｚｚｙＣ－Ｍｅａｎｓ，ＦＣＭ）聚类和对模拟漏失事件感知数量最大化的给水管网压力监测点优化布设方法．首先，基于水量变化的节点压力敏感系数矩阵，利用ＦＣＭ聚类算法对管网节点进行聚类，在各子类中预选出备选压力监测点；然后，随机模拟可能发生的漏失事件测试符合条件的不同优化组合的备选监测点对漏失事件引起压力变化的感知情况，确定优化的监测点布设方案．案例研究结果表明，该方法确定的压力监测点布设方案能够在保证至少２个监测点感知到同一个漏失事件的前提下，以最少３个监测点感知到９６．３６％以上的漏失事件．该方法在提高监测可靠性的同时可以节约监测点布设成本．关键词：压力监测点；优化布设；模糊Ｃ均值聚类；漏失模拟；给水管网中图分类号：ＴＵ９９１ 文献标志码：Ａ 文章编号：１６７４－８５３０（２０２２）０６－０６４２－０７Ｄｏｉ：１０．３９６９／ｊ．ｉｓｓｎ．１６７４－８５３０．２１．００１７ 徐志强，刘崎峰，刘佳峰，等．基于ＦＣＭ聚类和漏失模拟的给水管网压力监测点布设［Ｊ］．灌溉机械工程学报，２０２２，４０（６）：６４２－６４８． ＸＵＺｈｉｑｉａｎｇ，ＬＩＵＱｉｆｅｎｇ，ＬＩＵＪｉａｆｅｎｇ，ｅｔａｌ．ＯｐｔｉｍｉｚｅｄｌａｙｏｕｔｏｆｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｉｎｗａｔｅｒｓｕｐｐｌｙｎｅｔｗｏｒｋｂａｓｅｄｏｎＦＣＭｃｌｕｓｔｅｒｉｎｇａｎｄｌｅａｋａｇｅｅｖｅｎｔｓｓｉｍｕｌａｔｉｏｎ［Ｊ］．Ｊｏｕｒｎａｌｏｆｄｒａｉｎａｇｅａｎｄｉｒｒｉｇａｔｉｏｎｍａｃｈｉｎｅｒｙｅｎｇｉｎｅｅｒｉｎｇ，２０２２，４０（６）：６４２－６４８．（ｉｎＣｈｉｎｅｓｅ）ＯｐｔｉｍｉｚｅｄｌａｙｏｕｔｏｆｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｉｎｗａｔｅｒｓｕｐｐｌｙｎｅｔｗｏｒｋｂａｓｅｄｏｎＦＣＭｃｌｕｓｔｅｒｉｎｇａｎｄｌｅａｋａｇｅｅｖｅｎｔｓｓｉｍｕｌａｔｉｏｎＸＵＺｈｉｑｉａｎｇ，ＬＩＵＱｉｆｅｎｇ，ＬＩＵＪｉａｆｅｎｇ，ＳＨＡＮＬｉｍｉｎｇ，ＰＥＩＷｅｎｓｉ（ＣｏｌｌｅｇｅｏｆＥｃｏｌｏｇｙａｎｄＥｎｖｉｒｏｎｍｅｎｔ，ＩｎｎｅｒＭｏｎｇｏｌｉａＵｎｉｖｅｒｓｉｔｙ，Ｈｏｈｈｏｔ，ＩｎｎｅｒＭｏｎｇｏｌｉａ０１００２１，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：Ｉｎｏｒｄｅｒｔｏｉｍｐｒｏｖｅｔｈｅｒｅｌｉａｂｉｌｉｔｙｏｆｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｉｎｔｈｅｗａｔｅｒｓｕｐｐｌｙｎｅｔｗｏｒｋ，ａｌａｙｏｕｔｍｅｔｈｏｄｏｆｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｃｏｍｂｉｎｉｎｇｆｕｚｚｙＣ－ｍｅａｎｓ（ＦＣＭ）ｃｌｕｓｔｅｒｉｎｇｗｉｔｈｔｈｅｍａｘｉｍｕｍｐｅｒｃｅｉｖｅｄｎｕｍｂｅｒｏｆｌｅａｋａｇｅｅｖｅｎｔｓｗａｓｐｒｏｐｏｓｅｄｂａｓｅｄｏｎａＤｉｓｔｒｉｃｔＭｅｔｅｒｅｄＡｒｅａ（ＤＭＡ）ｉｎａｔｏｗｎｉｎｎｏｒｔｈｅｒｎＣｈｉｎａ．Ｆｉｒｓｔｌｙ，ｔｈｅＦＣＭａｌｇｏｒｉｔｈｍｗａｓｕｓｅｄｔｏｃｌｕｓｔｅｒａｌｌｎｏｄｅｓｂａｓｅｄｏｎｔｈｅｎｏｄｅｓｐｒｅｓｓｕｒｅｓｅｎｓｉｔｉｖｉｔｙｃｏｅｆｆｉｃｉｅｎｔｍａｔｒｉｘｏｆｗａｔｅｒｑｕａｎｔｉｔｙｃｈａｎｇｅｉｎｔｈｅｐｉｐｅｎｅｔｗｏｒｋ，ａｎｄｔｈｅａｌｔｅｒ ｎａｔｉｖｅｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｉｎｅａｃｈｓｕｂｃｌａｓｓｗｅｒｅｐｒｅｓｅｌｅｃｔｅｄ．Ｔｈｅｎ，ｔｈｅｐｏｓｓｉｂｌｅｌｅａｋａｇｅｅｖｅｎｔｓｗｅｒｅｓｉｍｕｌａｔｅｄｒａｎｄｏｍｌｙｔｏｔｅｓｔｔｈｅｐｅｒｃｅｐｔｉｏｎｏｆｔｈｅａｌｔｅｒｎａｔｉｖｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｗｉｔｈｄｉｆｆｅｒｅｎｔｏｐｔｉ ｍｉｚｅｄｃｏｍｂｉｎａｔｉｏｎｓｆｏｒｔｈｅｐｒｅｓｓｕｒｅｃｈａｎｇｅｓｃａｕｓｅｄｂｙｔｈｅｌｅａｋａｇｅｅｖｅｎｔｓ，ａｎｄｔｈｅｏｐｔｉｍａｌｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓ′ｌａｙｏｕｔｓｃｈｅｍｅｗａｓｄｅｔｅｒｍｉｎｅｄ．Ｔｈｅｓｉｍｕｌａｔｉｏｎｒｅｓｕｌｔｓｏｆｔｈｉｓｃａｓｅｓｈｏｗｔｈａｔｍｏｒｅｔｈａｎ９６．３６％ｏｆｔｈｅｌｅａｋａｇｅｓｗｅｒｅｄｅｔｅｃｔｅｄｗｉｔｈａｔｌｅａｓｔｔｈｒｅｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｏｎｔｈｅｐｒｅｍｉｓｅｅｎｓｕｒｉｎｇｔｈａｔａｔｌｅａｓｔｔｗｏｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｃａｎｐｅｒｃｅｉｖｅｔｈｅｓａｍｅｌｅａｋａｇｅｅｖｅｎｔ．Ｔｈｉｓｍｅｔｈｏｄｃａｎｉｍｐｒｏｖｅｔｈｅｍｏｎｉｔｏｒｉｎｇｒｅｌｉａｂｉｌｉｔｙａｎｄａｔｔｈｅｓａｍｅｔｉｍｅｓａｖｅｔｈｅｃｏｓｔ．Ｋｅｙｗｏｒｄｓ：ｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓ；ｏｐｔｉｍｉｚｅｄｌａｙｏｕｔ；ｆｕｚｚｙＣ－ｍｅａｎｓｃｌｕｓｔｅｒｉｎｇ；ｌｅａｋａｇｅｓｉｍｕｌａｔｉｏｎ；ｗａｔｅｒｓｕｐｐｌｙｎｅｔｗｏｒｋ 随着“智慧水务”的发展，给水管网中压力和流量数据能被实时采集和传送，但在复杂管网中对所有节点和管段进行监控是不现实的．对压力监测点优化布设可以在帮助供水机构节约投资的同时实时掌控管网运行状态、提升服务质量，对应对水资源短缺问题也具有重要意义［１－４］．在工程中对压力监测点布设主要通过经验选择和理论计算分析实现．经验法一般依据地形和用户性质在管网最不利末端、给水边界处、地形高差变化急剧以及需水大用户等重要位置设置监测点，用于监测指导管网的运行调度．但不足在于其位置由工程师人为指定、合理性和准确性差，这种过度依赖设计人员经验的方法也会使工程投资较大［５－７］．理论计算分析法是基于管网拓扑学，依靠数学分析及水力模型在建立监测点间的数学关系和水力联系的前提下对监测点进行优化．方潜生等［４］以节点坐标、影响度和压力标准差构建节点特征矩阵，对节点进行聚类分析，优化布设压力监测点．梁建文等［８］以局部漏损对压力的影响度为目标来实现压力监测点的优化布置，选出灵敏度高的监测点，实现了对给水管网健康状态的监测．张宏伟等［９］通过建立管网节点压力的灵敏度矩阵和方程，利用灵敏度分析法解决压力监测点的布设选址问题，该方法建立了监测点间的水力联系和数学联系．刘书明等［１０］通过建立反映管网节点水量和压力变化的节点覆盖水量和相关水压的多目标优化算法，对管网压力监测点进行选择，结果表明，监测点数量占总节点数的０．３％时就可以覆盖４５％以上的管网节点．金溪等［１１］采用遗传算法求解测压点优化布置模型，实现压力监测点的自动优化布置，并利用５个监测点对管网运行压力进行实时监测分析．此外，刘倍良［１２］以漏失定位为目的，利用水力模型结合压力驱动节点流量的漏失定位方法对某给水管网分区的压力监测点进行了优化布设，选出对爆管漏失感知灵敏度高的节点进行监测．ＳＯＲＯＵＳＨ等［１３］将地理位置统计和遗传优化结合，通过求解地理位置目标函数对压力监测点进行优化，保证监测点的分布更加均匀．ＲＡＹＡＲＯＴＨ等［１４］运用一种基于随机决策树分类器的混合蛙跳优化技术对监测点的压力数据优化，选取出对异常压力检测精度高的节点作为监测点．目前，对于监测点优化的研究都是在寻找区域内受水量摄动压力变化敏感性强、能够代表临近区域内水压情况的节点作为监测点，但在布设过程中对不同位置的监测点对漏失事故的感知情况考虑较少，用１个监测点对相应区域监测时，压力波动以及监测设备的测量误差和稳定性问题会影响监测的可靠性．为此文中以某给水ＤＭＡ管网为研究对象，结合ＦＣＭ聚类和水力模型产生的漏失事件对监测点进行优化，在保证１个漏失事件至少同时被２个监测点感知到的前提下确定可靠性和经济性较好的监测点布设方案．１ 研究案例概况案例研究对象为内蒙古某旗县供水管网中的一个给水ＤＭＡ，该区域是在２００８年开发建设并逐渐入住的综合商业住宅区，有１７６３户常住居民以及幼儿园、宾馆、餐饮饭店、银行、电玩城、物流快递、五金加工、汽修汽配和汽贸租赁公司等营业性网点７２３家．服务区域地势西高东低，地面标高在１３４１．０～１３７２．４ｍ．管网结构呈东西走向的矩形分布，平均埋深约为２ｍ，有节点３７９个，管道３８５段为２０１６—２０１７年铺设的管径ＤＮ６３到ＤＮ３００不等的ＰＥ管，服务面积约０．３ｋｍ２，共有需水用户２６４５户．经过精简优化后保留节点１３２个，给水干管长度约６．５ｋｍ，平均供水流量为６．０２Ｌ／ｓ，其中最大供水流量为９．０３Ｌ／ｓ，最小供水流量为３．６２Ｌ／ｓ，利用Ｅｐａｎｅｔ软件建立该区域的管网水力模型并结合ＭＡＴＬＡＢ计算平台进行研究分析．２ 监测点优化布设分析２．１ 对管网节点聚类划分通过聚类将水力变化性质相近的节点划分为一类，然后再逐类进行分析选取代表监测点，具体聚类步骤如下：①在正常情况下对管网模型进行水力模拟，计算各节点的水压值Ｈ０＝｛Ｈ０１，Ｈ０２，Ｈ０３，…，Ｈ０１３２｝，Ｈ０ｉ为正常情况下节点ｉ的压力值；②在第ｊ６４３个节点中增加部分水量作为漏失流量，进行水力模拟，计算假设ｊ节点漏失后各节点的压力值Ｈ１ｊ＝｛Ｈ１ｊ，１，Ｈ１ｊ，２，Ｈ１ｊ，３，…，Ｈ１ｊ，１３２｝，Ｈ１ｊ，ｉ为添加漏失流量后节点ｉ的压力值；③计算各节点的压力变化值△Ｈ＝｛Ｈ０１－Ｈ１ｊ，１，Ｈ０２－Ｈ１ｊ，２，Ｈ０３－Ｈ１ｊ，３，…，Ｈ０１３２－Ｈ１ｊ，１３２｝，Ｈ０ｉ－Ｈ１ｊ，ｉ为ｊ点发生漏失后ｉ点的压力变化值；④计算节点ｊ漏失对节点ｉ压力变化的影响系数Ｘ（ｉ，ｊ）＝△Ｈ／（Ｈ０ｊ－Ｈ１ｊ）；⑤清除添加给ｊ节点的漏失水量，选择下一节点，重复上述②－④步骤，直到遍历全部的１３２个节点，得到管网节点的压力影响系数特征矩阵Ｘ为Ｘ＝Ｈ０１－Ｈ１１，１Ｈ０１－Ｈ１１Ｈ０２－Ｈ１１，２Ｈ０１－Ｈ１１…Ｈ０１３２－Ｈ１１，１３２Ｈ０１－Ｈ１１Ｈ０１－Ｈ１２，１Ｈ０２－Ｈ１２Ｈ０２－Ｈ１２，２Ｈ０２－Ｈ１２…Ｈ０１３２－Ｈ１２，１３２Ｈ０２－Ｈ１２Ｈ０１－Ｈ１１３２，１Ｈ０１３２－Ｈ１１３２Ｈ０２－Ｈ１１３２，２Ｈ０１３２－Ｈ１１３２…Ｈ０１３２－Ｈ１１３２，１３２Ｈ０１３２－Ｈ１１３２．（１）由相关研究可知，当监测点数量占总节点数０．３％时，即可监测到４５％以上的节点，且一般要求监测范围覆盖８０％以上的管网节点和水量［１０，１５］．本研究区域管网节点总数为１３２，结合实际供水服务面积情况，理论上设置６个备选监测点即可达到监控整个管网要求．因此，初步设置６个监测点，同时对不同组合的监测位置进行优化，即对全部节点初步聚类为６类．为了保证同一个聚类分组中的需水节点之间有良好的相关性和相似性，有时候在对节点进行聚类时不能根据计算简单地将其划分到哪一类中，需要在分类计算的基础上结合管网的实际拓扑结构对聚类结果进行局部调整，保证同一聚类分组中的需水节点之间有直接通水路径．根据ＦＣＭ聚类划分中每个样本不是严格划分为某一类，而是以一定的隶属度属于某一类的聚类原则［１６］，采用ＦＣＭ算法对节点聚类，可利用ＭＡＴＬＡＢ中的“ｆｃｍ”函数进行分析计算，得到表１的聚类计算结果，表中ＩＤ为节点编号．第２子类中的ＳＴＪ２和第３子类中的Ｐ９７，Ｐ９９，ＰＲＶ１节点与各自子类中的其他节点都没有直接的水力连通性，而与第４子类的节点相邻，为保证每类中的节点都有直接相连的通水路径和相对该类较大的隶属度，在聚类计算的基础上对节点ＳＴＪ２，Ｐ９７，Ｐ９９，ＰＲＶ１重新划分到第４子类中．表１ 各节点聚类分布表Ｔａｂ．１ Ｃｌｕｓｔｅｒｄｉｓｔｒｉｂｕｔｉｏｎｏｆｐｏｉｎｔｓ类别ＩＤ１Ｐ１，Ｐ２，Ｐ３，ＳＴＪ６，Ｐ４，Ｐ４１，Ｐ８，Ｐ４３，Ｐ１２，Ｐ２０，Ｐ２６，Ｐ２８，Ｐ３０，Ｐ３２，Ｐ５４，Ｐ３４，Ｐ５，Ｐ４９，Ｐ４８，Ｐ４７，Ｐ４６，Ｐ５９，Ｐ５６，Ｐ４４－１，ＳＴＪ４，ＳＴＪ５，Ｐ４１－１，Ｐ６８，ＳＴＪ３，Ｐ６５－１，Ｐ６３－１，Ｐ６１－１，Ｐ５９－１，Ｐ５５－１，ＦＭ３－１，ＳＴＪ，ＳＴＪ７，Ｐ８３，Ｐ８２２Ｐ７９，Ｐ７７，Ｐ７５，Ｐ７３，Ｐ７０，ＳＴＪ２，Ｇ－０１，Ｇ－０２，Ｇ－０３，Ｇ－０４，Ｇ－０５，Ｇ－０６３Ｐ９７，Ｐ９９，ＰＲＶ１，Ｐ１２１，Ｐ１２３，Ｐ１２５，Ｐ１２７，Ｐ１２９，ＳＴＪ８，Ｐ１５０，Ｐ１３１，Ｐ１３３，Ｐ１３５，Ｐ１３７，Ｐ１３９，Ａ２，ＤＳ３４ＰＲＶ２，Ｐ１０２，Ｐ１０４，ＳＴＪ１，Ｇ－０７，Ｐ１１１，Ｐ１１２，Ｐ１１５，Ｐ１１７，Ｐ１１９５Ｐ１５２，Ｐ１５４，Ｐ１５６，Ｐ１５８，Ｐ１６０，Ｐ１６２，Ｐ１４１，Ｐ１４３，Ｐ１４５，Ｐ１４７，ＳＴＪ９，ＳＴＪ１０，Ｐ１９０，Ｐ１７１，Ｐ１７５，Ｐ１８３，Ｐ１８５，Ｐ１８７６ＤＳ５，ＤＳ７，ＤＳ１０，ＤＳ１２，ＤＳ１４，ＤＳ１６，ＤＳ２０，Ａ４，ＤＥ３，ＤＥ５，ＤＥ７，ＤＥ９，ＤＥ１２，ＤＥ１４，ＤＥ１７，ＤＥ１９，ＤＥ２１，ＤＥ２２，ＤＥ２４，ＤＥ２６，ＤＥ２９，Ｄ３０，ＤＥ３２，ＤＥ３７，ＤＥ３９，ＤＥ４４，Ａ５，Ａ９，Ａ１１，Ａ１４，ＤＹ２，ＤＹ７，ＤＹ１２，ＤＹ１４，ＤＹ１５，ＤＹ１６２．２ 选取备选监测点通过聚类把由水量漏失引起的压力变化相近的节点归为一类，在ＦＣＭ聚类中每类中隶属度ｕ越大的节点越具有所属子类的代表性［１６］，选取子聚类中隶属度最大的节点Ｐ４７，Ｇ－０１，ＳＴＪ１，ＳＴＪ８，ＳＴＪ１０和ＤＥ２１作为备选监测点，如表２所示，表中ｕ为隶属度．表２ 子类中隶属度最大节点Ｔａｂ．２ Ｐｏｉｎｔｓｗｉｔｈｔｈｅｌａｒｇｅｓｔｍｅｍｂｅｒｓｈｉｐｉｎｅａｃｈｓｕｂｃｌａｓｓ ＩＤｕＩＤｕＰ４７０．９９４３２０ＳＴＪ８０．９３４３２０Ｇ－０１０．９７７３０２ＳＴＪ１００．９６８４６６ＳＴＪ１０．９８２９７５ＤＥ２１０．９７４１９６图１ 管网节点聚类结果Ｆｉｇ．１ Ｃｌｕｓｔｅｒｉｎｇｒｅｓｕｌｔｓｏｆｐｉｐｅｎｅｔｗｏｒｋｐｏｉｎｔｓ２．３ 模拟漏失事件建立备选监测点压降集合由ＷＵ等［１７］的研究可知，当漏失发生时，距离漏失位置越近的监测点其压降越大，可以利用水力模型模拟相应的漏失事故，对监测点进行测试分析．研究中对漏失事件的模拟和建立相应的监测点压降集合主要有以下几个过程：①对水力模型进行模拟分析，得到正常情况下６个备选压力监测点压力６４４值Ｈ０＝｛Ｈ０１，Ｈ０２，Ｈ０３，Ｈ０４，Ｈ０５，Ｈ０６｝．②模拟漏失事件，计算备选监测点的漏失压力值．本案例中假设每个节点可能平均发生５次漏失事件，则利用蒙特卡罗算法随机生成６６０个漏失，模拟漏失水量如果太大，用户就会感到用水异常而发生热线报警，而模拟漏水量太小则因压力变化过小而无法被感知．在本研究中，拟设置漏失水量为平均供水总流量的２０％～６０％，考虑到管网中用户水量时刻变化的特点，分别取需水量最小、最大和中位数３个时段进行漏失事件模拟，总计生成１９８０个不同的随机漏失事件．计算每个漏失事件发生后备选监测点的压力值Ｈ１ｊ＝｛Ｈ１ｊ，１，Ｈ１ｊ，２，Ｈ１ｊ，３，Ｈ１ｊ，４，Ｈ１ｊ，５，Ｈ１ｊ，６｝，Ｈ１ｊ，ｉ为漏失事件ｊ发生后第ｉ个备选监测点的压力值．③计算备选压力监测点在每个漏失事件发生前后压力下降值△Ｈ＝｛△Ｈｊ，１，△Ｈｊ，２，△Ｈｊ，３，△Ｈｊ，４，△Ｈｊ，５，△Ｈｊ，６｝，△Ｈｊ，ｉ为备选压力监测点ｉ在漏失事件ｊ发生前后的压降值，即△Ｈ＝Ｈ０１，１－Ｈ１１，１Ｈ０１，１－Ｈ１１，１…Ｈ０１，１－Ｈ１１，１Ｈ０２，１－Ｈ１２，１Ｈ０２，１－Ｈ１２，１…Ｈ０２，１－Ｈ１２，１Ｈ０１９８０，１－Ｈ１１９８０，１Ｈ０１９８０，１－Ｈ１１９８０，１…Ｈ０１９８０，１－Ｈ１１９８０，１．（２）④对备选监测点压降数据进行处理．由于管网日常运行中存在一定的压力波动以及监测设备的测量误差和稳定性问题，需要对监测到的压降数据进行进一步处理．添加１个压力波动阈值ζ，当压降值小于ζ时为０，表示该监测点没有感知到此次漏失事件；当压降值大于等于ζ时为１，表示该监测点成功监测到此次漏失事件．阈值ζ需根据管网的实际运行情况选取，如果太大会导致不能及时监测到漏失，而太小又会增加监测的误报率．依据运行经验，本案例中ζ设置为０．２ｍ．ｓｊ，ｉ＝０，△Ｈｊ，ｉ＜０．２，１，△Ｈｊ，ｉ≥０．２，{（３）式中：ｓｊ，ｉ为处理过的备选压力监测点ｉ在漏失事件ｊ发生前后的压降值．根据式（３），将备选监测点的全部压降数据换算成０－１分布的形式，即未感知到和感知到漏失．表３以案例中前１０次模拟漏失事件的压降数据进行简要说明．按式（３）的换算方式对表３中的压力下降值ΔＨ进行处理，得到表４记录的各备选监测点的压力下降情况，表中ΔＨ′为换算后备选监测点压力下降值．表３ 备选监测点压力下降值Ｔａｂ．３ Ｐｒｅｓｓｕｒｅｄｒｏｐｖａｌｕｅｓｏｆａｌｔｅｒｎａｔｉｖｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｍ漏失事件ΔＨＰ４７Ｇ－０１ＳＴＪ１ＳＴＪ８ＳＴＪ１０ＤＥ２１１０．２８０．２９０．３００．３００．３００．３１２０．２５０．２３０．２１０．２２０．１８０．２１３０．１７０．２２０．１７０．１７０．１７０．１７４０．３００．３００．３１０．３２０．３２０．３３５０．２６０．２５０．２３０．２３０．１９０．１８６０．２００．２００．２１０．２１０．２１０．２１７０．１９０．１９０．２３０．２３０．２３０．２４８０．１９０．１９０．１９０．１９０．１９０．１８９０．３００．３００．３２０．３３０．３５０．３２１００．２６０．３５０．２６０．２６０．２６０．２６ 表４中漏失事件８引起６个备选监测点的压力下降值全部为０，即由于漏失事件８引起的压降值较小，没有被感知到，将诸如漏失事件８引起的所有备选监测点压力下降值为０的漏失事件定义为无效漏失，并从整个漏失事件集合中去除，最后得到１５６６个用于测试压力监测点优化布设的模拟漏失事件集合矩阵△Ｓ△Ｓ＝ｓ１，１ｓ１，２…ｓ１，６ｓ２，１ｓ２，２…ｓ２，６ ｓ１５６６，１ｓ１５６６，２…ｓ１５６６，６ ．（４）矩阵△Ｓ中，行元素表示某个漏失事件引起的压降值被全部备选监测点的感知情况，列元素则表示某个备选监测点对所有模拟漏失事件的感知情况．表４ 换算后备选监测点压力下降值Ｔａｂ．４ Ｐｒｅｓｓｕｒｅｄｒｏｐｖａｌｕｅｓｏｆａｌｔｅｒｎａｔｉｖｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｃｏｎｖｅｒｔｅｄｍ漏失事件ΔＨ′Ｐ４７Ｇ－０１ＳＴＪ１ＳＴＪ８ＳＴＪ１０ＤＥ２１１１１１１１１２１１１１０１３０１００００４１１１１１１５１１１１００６１１１１１１７００１１１１８００００００９１１１１１１１０１１１１１１２．４ 建立监测点优化布设目标函数将矩阵△Ｓ中列元素写成向量的形式Ｓｉ＝［ｓ１，ｉ，ｓ２，ｉ，ｓ３，ｉ，…，ｓ１５６６，ｉ］Ｔ，Ｓｉ为第ｉ个监测点感知到的６４５漏失事件向量．案例要求对同一个漏失事件至少同时被２个备选监测点感知来提高监测的可靠性，可以将问题理解为求将向量Ｓｉ（ｉ＝１，２，３，４，５，６）中相同位置元素相加，和向量中不小于２的元素个数的最大化，故最少可选出２个、最多可选择６个备选监测点．如以前１０个漏失事件为例，假设选取Ｐ４７和Ｇ－０１作为监测点，去除无效事件８后还剩余９个模拟漏失事件，即向量Ｓ１＝［１，１，０，１，１，１，０，１，１］Ｔ和Ｓ２＝［１，１，１，１，１，１，０，１，１］Ｔ，Ｓ１＋Ｓ２＝［２，２，１，２，２，２，２，０，２］Ｔ，将２个向量对应位置元素相加后得到的新向量中不小于２的元素有７个，说明该布设方案能够感知到７个漏失事件．据此，监测点优化选择的目标函数计算方法可表示为ＭＡＸ：ｓｕｍ｛［∑Ｃｎ６（Ｓｉ）］≥２｝，（５）式中：Ｃｎ６为从６个备选监测点中任选ｎ个的组合；ｓｕｍ表示对目标值求和．２．５ 对优化目标函数求解确定监测点布设方案案例中可将优化目标函数的求解看作是排列组合问题，通过对集合矩阵△Ｓ中列元素的求和，得到６个备选监测点能够感知到漏失事件数分别为１３８１，１３３５，１３９５，１４３１，１３５６和１４０５．每个备选压力监测点都能感知大多数的漏失事件，但对漏失事件的感知情况又各不相同，可从备选监测点中选取２，３，４，５，６个作为最终结果．由排列组合公式ＣｎＮ可知，上述方案分别有１５，２０，１５，６，１种不同的组合方式，对全部５７种组合方式按目标函数计算，得到表５中不同ｎ时感知漏失事件最少和最多２种布设方案及其漏失事件感知情况，表中Ｍ为感知到漏失的数量．表５ 不同监测点布设方式及感知漏失事件数量Ｔａｂ．５ Ｎｕｍｂｅｒｏｆｐｅｒｃｅｉｖｅｄｌｅａｋａｇｅｅｖｅｎｔｓｏｆｌａｙｏｕｔｗｉｔｈｄｉｆｆｅｒｅｎｔｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｎ选取备选监测点的组合方式Ｍ备注２Ｇ－０１，ＳＴＪ１０１１５５最少ＳＴＪ１，ＳＴＪ８１３０７最多３Ｐ４７，Ｇ－０１，ＳＴＪ１０１４３３最少Ｐ４７，ＳＴＪ８，ＤＥ２１１５０９最多４Ｐ４７，Ｇ－０１，ＳＴＪ１，ＳＴＪ１０１５００最少Ｐ４７，ＳＴＪ８，ＳＴＪ１０，ＤＥ２１１５３０最多５Ｐ４７，Ｇ－０１，ＳＴＪ１，ＳＴＪ８，ＳＴＪ１０１５２８最少Ｐ４７，Ｇ－０１，ＳＴＪ１，ＳＴＪ８，ＤＥ２１１５３０最多６Ｐ４７，Ｇ－０１，ＳＴＪ１，ＳＴＪ８，ＳＴＪ１０，ＤＥ２１１５３２— 在至少有２个监测点能同时作用的前提下，当拟设置２个监测点时，选取节点Ｇ－０１，ＳＴＪ１０感知到的漏失事件最少为１１５５个，当选取节点ＳＴＪ１，ＳＴＪ８感知到的漏失事件最多为１３０７个；当拟设置３个监测点时，选取节点Ｐ４７，Ｇ－０１，ＳＴＪ１０感知到的漏失事件最少为１４３３个，选取节点Ｐ４７，ＳＴＪ８，ＤＥ２１感知到的漏失事件最多为１５０９个；当拟设置４个监测点时，选取Ｐ４７，ＳＴＪ８，ＳＴＪ１０，ＤＥ２１感知到的漏失事件最多为１５３０个；当拟设置５个监测点时，选取Ｐ４７，Ｇ－０１，ＳＴＪ１，ＳＴＪ８，ＤＥ２１感知到的漏失事件最多为１５３０个；当选取全部备选节点时能够感知到的漏失事件为１５３２个．可见，当拟设置的监测点在４个以内时，随着监测点数量的增多，满足目标函数的被感知到的漏失数量明显增多；当拟设置的监测点为５或６时，虽然感知到的漏失事件有所增加，但增幅较小．经对比，当选取Ｐ４７，ＳＴＪ８，ＤＥ２１为监测点时最大可以保证２个以上监测点同时感知到漏失事件占全部漏失事件集合的９６．３６％以上；当设置４个监测点时却为９７．７０％．一般在工程中增加１个监测点，设备成本在２０００～４０００千元不等，安装费约５００元，后期更换电源的维护费以及数据收集处理的工作则相对漏失事件覆盖率仅增加１％．在考虑成本最小化和对可能发生的漏失事件的感知率最大化的基础上，设置３个监测点对管网运行压力进行实时监测最为科学．图２ 监测点优化布设位置Ｆｉｇ．２ Ｏｐｔｉｍａｌｌｏｃａｔｉｏｎｏｆｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓ３ 讨 论相对于Ｋ均值聚类等对样本非此即彼划分的硬聚类，ＦＣＭ聚类具有以一定隶属度将性质相近样本归为一类的特性，使得同一类内的样本具有较高同质性，类间样本具有较高异质性．因此，ＦＣＭ算法的软聚类特性对解决给水管网这种复杂的非线性系统中节点聚类问题十分适用．以管网中节点压力变化的敏感系数为样本节点的特征，利用ＦＣＭ聚类对节点聚类分析，得到节点分别属于各子类的隶属度．在计算的基础上，结合管网拓扑结构和隶属度大小对子类中离散的样本节点进行局部归类调整，最后确定的同一类中的节点既保证了以往单纯聚类６４６方法［５］中节点间紧密的数学联系，也保证了同一类节点间都有直接的通水路径，同时根据隶属度大小可从各子类中预选出备选监测点．利用Ｅｐａｎｅｔ软件基于ＭＡＴＬＡＢ平台可以直观模拟管网中不同位置可能发生的各种突发爆管漏失事件，随机产生的漏失事件可有效覆盖管网不同位置节点在不同供水流量下的漏失状态．利用目标函数测试不同的监测方案，经计算，当选取节点Ｐ４７，ＳＴＪ８，ＤＥ２１为监测点时的漏失事件监测覆盖率相对于只设置２个监测点时提高１０％以上，达到９６．３６％；设定的３个压力监测点对漏失事件感知率和相关研究中要求８０％以上的较高保证率［１５］结果一致；但当设置４个监测点时，监测覆盖率增大效果却并不明显，也验证了研究中假设６个备选监测点可有效覆盖到管网各部位的合理性．相比以往大多偏向于考虑监测点自身的特性和分布均匀性的布设方法，该方法更能突出设置的监测点与漏失事件感知与否间的数量关系，目标函数保证至少２个监测点同时作用，提高了对管网异常压力波动监测的可靠性，减少了高误报率给管网检修人员带来的额外工作量．４ 结 论１）通过ＭＡＴＬＡＢ和Ｅｐａｎｅｔ水力模型模拟可能发生的漏失事件，对ＦＣＭ聚类预选的备选监测点进行漏失感知测试和评价，该方法能直观地体现不同布设方案对可能发生的漏失事件的感知情况，最终确定的布设方案可以在保证监测点间良好数学关系和水力联系的前提下，利用Ｐ４７，ＳＴＪ８，ＤＥ２１这３个监测点监测到９６．３６％以上的漏失，同时方案要求至少２个监测点能感知到同一个漏失事件，提高了监测可靠性．２）案例管网为同一时期铺设且服务区域较小，在监测点的布设中没有考虑节点重要性以及地形变化的影响，对于大规模管网可在此基础上对部分节点进行单独监测．３）此外，由于给水管网是一个非线性的复杂系统，在未来的研究中需根据服务区域用户特征，对供水流量及压力波动变化等特征进行长期监测总结，建立更准确的监测点布设水力模型．参考文献（Ｒｅｆｅｒｅｎｃｅｓ）［１］ 尹兆龙，信昆仑，项宁银．给水管网压力监测点布置的实用方法［Ｊ］．中国给水排水，２０１４，３０（２）：１９－２３．ＹＩＮＺｈａｏｌｏｎｇ，ＸＩＮＫｕｎｌｕｎ，ＸＩＡＮＧＮｉｎｇｙｉｎ．Ａｐｒａｃｔｉｃａｌｍｅｔｈｏｄｆｏｒｌａｙｏｕｔｏｆｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｉｎｗａｔｅｒｄｉｓｔｒｉｂｕｔｉｏｎｎｅｔｗｏｒｋ［Ｊ］．Ｃｈｉｎａｗａｔｅｒａｎｄｗａｓｔｅｗａｔｅｒ，２０１４，３０（２）：１９－２３．（ｉｎＣｈｉｎｅｓｅ）［２］ 何忠华，袁一星．基于分区模型的城市给水管网压力监测点布置［Ｊ］．哈尔滨工业大学学报，２０１４，４６（１０）：３７－４１．ＨＥＺｈｏｎｇｈｕａ，ＹＵＡＮＹｉｘｉｎｇ．Ｌａｙｏｕｔｏｆｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｉｎｕｒｂａｎｗａｔｅｒｄｉｓｔｒｉｂｕｔｉｏｎｓｙｓｔｅｍｂａｓｅｄｏｎｄｉｓｔｒｉｃｔｍｏｄｅｌ［Ｊ］．ＪｏｕｒｎａｌｏｆＨａｒｂｉｎＩｉｎｓｔｉｔｕｔｅｏｆＴｅｃｈｎｏｌｏｇｙ，２０１４，４６（１０）：３７－４１．（ｉｎＣｈｉｎｅｓｅ）［３］ 张清周．基于水力模型的供水管网压力管理与漏损区域识别研究［Ｄ］．哈尔滨：哈尔滨工业大学，２０１７．［４］ 黄雅芳，卫海，管慧玲，等．聚类分析法选择给水管网压力监测点［Ｊ］．净水技术，２０１４，３３（Ｓ２）：７４－７６．ＨＵＡＮＧＹａｆａｎｇ，ＷＥＩＨａｉ，ＧＵＡＮＨｕｉｌｉｎｇ，ｅｔａｌ．Ｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｓｅｌｅｃｔｉｏｎｂｙｃｌｕｓｔｅｒｉｎｇａｎａｌｙｓｉｓｍｅｔｈｏｄｉｎｗａｔｅｒｓｕｐｐｌｙｎｅｔ［Ｊ］．Ｗａｔｅｒｐｕｒｉｆｉｃａｔｉｏｎｔｅｃｈｎｏｌｏｇｙ，２０１４，３３（Ｓ２）：７４－７６．（ｉｎＣｈｉｎｅｓｅ）［５］ 方潜生，张兆祥，谢陈磊，等．基于特征聚类的给水管网压力监测点优化布置［Ｊ］．安徽大学学报（自然科学版），２０１７，４１（４）：５５－６２．ＦＡＮＧＱｉａｎｓｈｅｎｇ，ＺＨＡＮＧＺｈａｏｘｉａｎｇ，ＸＩＥＣｈｅｎｌｅｉ，ｅｔａｌ．Ｏｐｔｉｍａｌｌｏｃａｔｉｏｎｏｆｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｉｎｗａｔｅｒｓｕｐｐｌｙｎｅｔｗｏｒｋｂａｓｅｄｏｎｆｅａｔｕｒｅｃｌｕｓｔｅｒｉｎｇ［Ｊ］．ＪｏｕｒｎａｌｏｆＡｎｈｕｉＵｎｉｖｅｒｓｉｔｙ（ｎａｔｕｒａｌｓｃｉｅｎｃｅｅｄｉｔｉｏｎ），２０１７，４１（４）：５５－６２．（ｉｎＣｈｉｎｅｓｅ）［６］ 孙征，王玉琼，涂杰，等．基于ＦＣＭ聚类的供水管网压力监测点的布置方法［Ｊ］．灌溉排水学报，２０１８，３７（Ｓ１）：４４－４６．ＳＵＮＺｈｅｎｇ，ＷＡＮＧＹｕｑｉｏｎｇ，ＴＵＪｉｅ，ｅｔａｌ．ＡｒｒａｎｇｅｍｅｎｔｍｅｔｈｏｄｏｆｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｏｆｗａｔｅｒｓｕｐｐｌｙｎｅｔｗｏｒｋｂａｓｅｄｏｎＦＣＭｃｌｕｓｔｅｒｉｎｇ［Ｊ］．Ｊｏｕｒｎａｌｏｆｉｒｒｉｇａｔｉｏｎａｎｄｄｒａｉｎａｇｅ，２０１８，３７（Ｓ１）：４４－４６．（ｉｎＣｈｉｎｅｓｅ）［７］ 陶涛，颜合想，信昆仑，等．基于ＳＣＡＤＡ压力监测的爆管定位分析［Ｊ］．供水技术，２０１６，１０（４）：１１－１４，２８．ＴＡＯＴａｏ，ＹＡＮＨｅｘｉａｎｇ，ＸＩＮＫｕｎｌｕｎ，ｅｔａｌ．ＰｉｐｅｂｕｒｓｔｌｏｃａｔｉｏｎａｎａｌｙｓｉｓｂａｓｅｄｏｎＳＣＡＤＡｓｙｓｔｅｍｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇ［Ｊ］．Ｗａｔｅｒｔｅｃｈｎｏｌｏｇｙ，２０１６，１０（４）：１１－１４，２８．（ｉｎＣｈｉｎｅｓｅ）［８］ 梁建文，肖笛，张宏伟，等．供水管网健康监测的压力监测点优化布置［Ｊ］．防灾减灾工程学报，２０１３，３３（Ｓ１）：５１－５７．ＬＩＡＮＧＪｉａｎｗｅｎ，ＸＩＡＯＤｉ，ＺＨＡＮＧＨｏｎｇｗｅｉ，ｅｔａｌ．Ｏｐｔｉｍａｌｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｉｎｗａｔｅｒｄｉｓｔｒｉｂｕｔｉｏｎｓｙｓｔｅｍｆｏｒｈｅａｌｔｈｍｏｎｉｔｏｒｉｎｇ［Ｊ］．Ｊｏｕｒｎａｌｏｆｄｉｓａｓｔｅｒｐｒｅｖｅｎｔｉｏｎａｎｄｍｉｔｉｇａｔｉｏｎｅｎｇｉｎｅｅｒｉｎｇ，２０１３，３３（Ｓ１）：５１－５７．（ｉｎＣｈｉ６４７ｎｅｓｅ）［９］ 张宏伟，张丽，梁建文．给水管网压力监测点的选址方法［Ｊ］．中国给水排水，２００３，１９（３）：５２－５５．ＺＨＡＮＧＨｏｎｇｗｅｉ，ＺＨＡＮＧＬｉ，ＬＩＡＮＧＪｉａｎｗｅｎ．Ｓｔｕｄｙｏｎｔｈｅｌａｙｉｎｇｏｆｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｎｏｄｅｓｏｆｗａｔｅｒｓｕｐｐｌｙｎｅｔｗｏｒｋ［Ｊ］．Ｃｈｉｎａｗａｔｅｒａｎｄｗａｓｔｅｗａｔｅｒ，２００３，１９（３）：５２－５５．（ｉｎＣｈｉｎｅｓｅ）［１０］ 刘书明，王欢欢，徐鹏，等．多目标大规模给水管网监测点的优化选址［Ｊ］．清华大学学报（自然科学版），２０１３，５３（１）：７８－８３．ＬＩＵＳｈｕｍｉｎｇ，ＷＡＮＧＨｕａｎｈｕａｎ，ＸＵＰｅｎｇ，ｅｔａｌ．Ｍｕｌｔｉｏｂｊｅｃｔｉｖｅｇｅｎｅｔｉｃａｌｇｏｒｉｔｈｍｆｏｒｏｐｔｉｍａｌｍｏｎｉｔｏｒｉｎｇｓｔａｔｉｏｎｉｎｌａｒｇｅｗａｔｅｒｄｉｓｔｒｉｂｕｔｉｏｎｓｙｓｔｅｍ［Ｊ］．ＪｏｕｒｎａｌｏｆＴｓｉｎｇｈｕａＵｎｉｖｅｒｓｉｔｙ（ｎａｔｕｒａｌｓｃｉｅｎｃｅｅｄｉｔｉｏｎ），２０１３，５３（１）：７８－８３．（ｉｎＣｈｉｎｅｓｅ）［１１］ 金溪，曾小兵，高金良，等．利用遗传算法进行供水管网压力监测点优化布置［Ｊ］．给水排水，２００７，３３（Ｚ１）：３４６－３４９．ＪＩＮＸｉ，ＺＥＮＧＸｉａｏｂｉｎｇ，ＧＡＯＪｉｎｌｉａｎｇ，ｅｔａｌ．Ｔｈｅｏｐｔｉｍａｌｐｌａｃｅｍｅｎｔｓｏｆｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｉｎｗａｔｅｒｓｕｐｐｌｙｎｅｔｗｏｒｋｂｙｕｓｉｎｇｔｈｅｇｅｎｅｔｉｃａｌｇｏｒｉｔｈｍ［Ｊ］．Ｗａｔｅｒａｎｄｗａｓｔｅｗａｔｅｒｅｎｇｉｎｅｅｒｉｎｇ，２００７，３３（Ｚ１）：３４６－３４９．（ｉｎＣｈｉｎｅｓｅ）［１２］ 刘倍良．城市供水管网爆管预警定位模型研究［Ｄ］．长沙：湖南大学，２０１９．［１３］ ＳＯＲＯＵＳＨＦ，ＡＢＥＤＩＮＩＭＪ．Ｏｐｔｉｍａｌｓｅｌｅｃｔｉｏｎｏｆｎｕｍ ｂｅｒａｎｄｌｏｃａｔｉｏｎｏｆｐｒｅｓｓｕｒｅｓｅｎｓｏｒｓｉｎｗａｔｅｒｄｉｓｔｒｉｂｕｔｉｏｎｓｙｓｔｅｍｓｕｓｉｎｇｇｅｏｓｔａｔｉｓｔｉｃａｌｔｏｏｌｓｃｏｕｐｌｅｄｗｉｔｈｇｅｎｅｔｉｃａｌｇｏｒｉｔｈｍ［Ｊ］．Ｊｏｕｒｎａｌｏｆｈｙｄｒｏｉｎｆｏｒｍａｔｉｃｓ，２０１９，２１（６）：１０３０－１０４７．［１４］ ＲＡＹＡＲＯＴＨＲ，ＳＩＶＡＲＡＤＪＥＧ．Ｒａｎｄｏｍｂａｇｇｉｎｇｃｌａｓ ｓｉｆｉｅｒａｎｄｓｈｕｆｆｌｅｄｆｒｏｇｌｅａｐｉｎｇｂａｓｅｄｏｐｔｉｍａｌｓｅｎｓｏｒｐｌａｃｅｍｅｎｔｆｏｒｌｅａｋａｇｅｄｅｔｅｃｔｉｏｎｉｎＷＤＳ［Ｊ］．Ｗａｔｅｒｒｅｓｏｕｒｃｅｓｍａｎａｇｅｍｅｎｔ，２０１９，３３：３１１１－３１２５．［１５］ 彭畅，彭森，吴卿，等．基于ＮＳＧＡ－Ⅱ的给水管网压力监测点多目标优化布置［Ｊ］．中国给水排水，２０１９，３５（１）：５８－６２．ＰＥＮＧＣｈａｎｇ，ＰＥＮＧＳｅｎ，ＷＵＱｉｎｇ，ｅｔａｌ．ＭｕｌｔｉｏｂｊｅｃｔｉｖｅｏｐｔｉｍｉｚａｔｉｏｎｏｆａｒｒａｎｇｅｍｅｎｔｏｆｐｒｅｓｓｕｒｅｍｏｎｉｔｏｒｉｎｇｐｏｉｎｔｓｉｎｗａｔｅｒｄｉｓｔｒｉｂｕｔｉｏｎｎｅｔｗｏｒｋｂａｓｅｄｏｎＮＳＧＡ－Ⅱ［Ｊ］．Ｃｈｉｎａｗａｔｅｒａｎｄｗａｓｔｅｗａｔｅｒ，２０１９，３５（１）：５８－６２．（ｉｎＣｈｉｎｅｓｅ）［１６］ 李柏年，吴礼滨．ＭＡＴＬＡＢ数据分析方法［Ｍ］．北京：机械工业出版社，２０１２．［１７］ ＷＵＬｉｐｅｎｇ，ＬＩＵＳｈｕｍｉｎｇ，ＷＡＮＧＸｉａｏｔｉｎｇ．Ｄｉｓｔａｎｃｅ ｂａｓｅｄｂｕｒｓｔｄｅｔｅｃｔｉｏｎｕｓｉｎｇｍｕｌｔｉｐｌｅｐｒｅｓｓｕｒｅｓｅｎｓｏｒｓｉｎｄｉｓｔｒｉｃｔｍｅｔｅｒｉｎｇａｒｅａｓ［Ｊ］．Ｊｏｕｒｎａｌｏｆｗａｔｅｒｒｅｓｏｕｒｃｅｓｐｌａｎｎｉｎｇａｎｄｍａｎａｇｅｍｅｎｔ，２０１８，１４４（１１）：０６０１８００９１．（责任编辑 朱漪云）６４８。
30个智能算法matlab代码
30个智能算法matlab代码以下是30个使用MATLAB编写的智能算法的示例代码： 1. 线性回归算法：matlab.x = [1, 2, 3, 4, 5];y = [2, 4, 6, 8, 10];coefficients = polyfit(x, y, 1);predicted_y = polyval(coefficients, x);2. 逻辑回归算法：matlab.x = [1, 2, 3, 4, 5];y = [0, 0, 1, 1, 1];model = fitglm(x, y, 'Distribution', 'binomial'); predicted_y = predict(model, x);3. 支持向量机算法：matlab.x = [1, 2, 3, 4, 5; 1, 2, 2, 3, 3];y = [1, 1, -1, -1, -1];model = fitcsvm(x', y');predicted_y = predict(model, x');4. 决策树算法：matlab.x = [1, 2, 3, 4, 5; 1, 2, 2, 3, 3]; y = [0, 0, 1, 1, 1];model = fitctree(x', y');predicted_y = predict(model, x');5. 随机森林算法：matlab.x = [1, 2, 3, 4, 5; 1, 2, 2, 3, 3]; y = [0, 0, 1, 1, 1];model = TreeBagger(50, x', y');predicted_y = predict(model, x');6. K均值聚类算法：matlab.x = [1, 2, 3, 10, 11, 12]; y = [1, 2, 3, 10, 11, 12]; data = [x', y'];idx = kmeans(data, 2);7. DBSCAN聚类算法：matlab.x = [1, 2, 3, 10, 11, 12]; y = [1, 2, 3, 10, 11, 12]; data = [x', y'];epsilon = 2;minPts = 2;[idx, corePoints] = dbscan(data, epsilon, minPts);8. 神经网络算法：matlab.x = [1, 2, 3, 4, 5];y = [0, 0, 1, 1, 1];net = feedforwardnet(10);net = train(net, x', y');predicted_y = net(x');9. 遗传算法：matlab.fitnessFunction = @(x) x^2 4x + 4;nvars = 1;lb = 0;ub = 5;options = gaoptimset('PlotFcns', @gaplotbestf);[x, fval] = ga(fitnessFunction, nvars, [], [], [], [], lb, ub, [], options);10. 粒子群优化算法：matlab.fitnessFunction = @(x) x^2 4x + 4;nvars = 1;lb = 0;ub = 5;options = optimoptions('particleswarm', 'PlotFcn',@pswplotbestf);[x, fval] = particleswarm(fitnessFunction, nvars, lb, ub, options);11. 蚁群算法：matlab.distanceMatrix = [0, 2, 3; 2, 0, 4; 3, 4, 0];pheromoneMatrix = ones(3, 3);alpha = 1;beta = 1;iterations = 10;bestPath = antColonyOptimization(distanceMatrix, pheromoneMatrix, alpha, beta, iterations);12. 粒子群-蚁群混合算法：matlab.distanceMatrix = [0, 2, 3; 2, 0, 4; 3, 4, 0];pheromoneMatrix = ones(3, 3);alpha = 1;beta = 1;iterations = 10;bestPath = particleAntHybrid(distanceMatrix, pheromoneMatrix, alpha, beta, iterations);13. 遗传算法-粒子群混合算法：matlab.fitnessFunction = @(x) x^2 4x + 4;nvars = 1;lb = 0;ub = 5;gaOptions = gaoptimset('PlotFcns', @gaplotbestf);psOptions = optimoptions('particleswarm', 'PlotFcn',@pswplotbestf);[x, fval] = gaParticleHybrid(fitnessFunction, nvars, lb, ub, gaOptions, psOptions);14. K近邻算法：matlab.x = [1, 2, 3, 4, 5; 1, 2, 2, 3, 3]; y = [0, 0, 1, 1, 1];model = fitcknn(x', y');predicted_y = predict(model, x');15. 朴素贝叶斯算法：matlab.x = [1, 2, 3, 4, 5; 1, 2, 2, 3, 3]; y = [0, 0, 1, 1, 1];model = fitcnb(x', y');predicted_y = predict(model, x');16. AdaBoost算法：matlab.x = [1, 2, 3, 4, 5; 1, 2, 2, 3, 3];y = [0, 0, 1, 1, 1];model = fitensemble(x', y', 'AdaBoostM1', 100, 'Tree'); predicted_y = predict(model, x');17. 高斯混合模型算法：matlab.x = [1, 2, 3, 4, 5]';y = [0, 0, 1, 1, 1]';data = [x, y];model = fitgmdist(data, 2);idx = cluster(model, data);18. 主成分分析算法：matlab.x = [1, 2, 3, 4, 5; 1, 2, 2, 3, 3]; coefficients = pca(x');transformed_x = x' coefficients;19. 独立成分分析算法：matlab.x = [1, 2, 3, 4, 5; 1, 2, 2, 3, 3]; coefficients = fastica(x');transformed_x = x' coefficients;20. 模糊C均值聚类算法：matlab.x = [1, 2, 3, 4, 5; 1, 2, 2, 3, 3]; options = [2, 100, 1e-5, 0];[centers, U] = fcm(x', 2, options);21. 遗传规划算法：matlab.fitnessFunction = @(x) x^2 4x + 4; nvars = 1;lb = 0;ub = 5;options = optimoptions('ga', 'PlotFcn', @gaplotbestf);[x, fval] = ga(fitnessFunction, nvars, [], [], [], [], lb, ub, [], options);22. 线性规划算法：matlab.f = [-5; -4];A = [1, 2; 3, 1];b = [8; 6];lb = [0; 0];ub = [];[x, fval] = linprog(f, A, b, [], [], lb, ub);23. 整数规划算法：matlab.f = [-5; -4];A = [1, 2; 3, 1];b = [8; 6];intcon = [1, 2];[x, fval] = intlinprog(f, intcon, A, b);24. 图像分割算法：matlab.image = imread('image.jpg');grayImage = rgb2gray(image);binaryImage = imbinarize(grayImage);segmented = medfilt2(binaryImage);25. 文本分类算法：matlab.documents = ["This is a document.", "Another document.", "Yet another document."];labels = categorical(["Class 1", "Class 2", "Class 1"]);model = trainTextClassifier(documents, labels);newDocuments = ["A new document.", "Another new document."];predictedLabels = classifyText(model, newDocuments);26. 图像识别算法：matlab.image = imread('image.jpg');features = extractFeatures(image);model = trainImageClassifier(features, labels);newImage = imread('new_image.jpg');newFeatures = extractFeatures(newImage);predictedLabel = classifyImage(model, newFeatures);27. 时间序列预测算法：matlab.data = [1, 2, 3, 4, 5];model = arima(2, 1, 1);model = estimate(model, data);forecastedData = forecast(model, 5);28. 关联规则挖掘算法：matlab.data = readtable('data.csv');rules = associationRules(data, 'Support', 0.1);29. 增强学习算法：matlab.environment = rlPredefinedEnv('Pendulum');agent = rlDDPGAgent(environment);train(agent);30. 马尔可夫决策过程算法：matlab.states = [1, 2, 3];actions = [1, 2];transitionMatrix = [0.8, 0.1, 0.1; 0.2, 0.6, 0.2; 0.3, 0.3, 0.4];rewardMatrix = [1, 0, -1; -1, 1, 0; 0, -1, 1];policy = mdpPolicyIteration(transitionMatrix, rewardMatrix);以上是30个使用MATLAB编写的智能算法的示例代码，每个算法都可以根据具体的问题和数据进行相应的调整和优化。
共空间模式CSP和滤波器组共空间模式FBCSP的MATLAB实现
共空间模式CSP和滤波器组共空间模式FBCSP的MATLAB实现首先，我们来介绍CSP的MATLAB实现。
CSP的目标是找到一组空间滤波器，能够使目标类别的变异性最大化，同时使非目标类别的变异性最小化。
下面是CSP的MATLAB实现的步骤：1. 输入数据准备：将原始EEG数据划分为两个类别的epochs，比如注意力集中和注意力放松时的脑电数据。
每个epoch都是一个时间段内的多个时刻信号。
2. 计算协方差矩阵：对于每个类别的epoch，计算该类别的协方差矩阵。
协方差矩阵是一个n x n的矩阵，其中n是脑电通道的数量。
3.特征提取：将协方差矩阵进行一个特征分解，得到特征向量和特征值。
根据特征值的大小，将特征向量按照升序排列。
4.空间模式提取：选择特定数量的最大和最小特征向量，相应地得到一组最佳空间滤波器。
空间滤波器可以通过对协方差矩阵进行特征投影来实现。
5. 转换特征：对于每个类别的epoch，将其信号与最佳空间滤波器进行投影，得到新的特征向量。
这些特征向量将用于比较不同类别的分离度。
6.分类器训练：使用提取的特征向量，可以使用各种分类算法，如线性判别分析(LDA)、支持向量机(SVM)等训练一个分类器。
接下来，我们讨论FBCSP的MATLAB实现。
FBCSP是CSP的扩展版本，它通过将原始EEG信号分解为一组频率子带，将CSP应用于每个子带，并合并各个子带的CSP特征向量来实现更好的分类性能。
FBCSP的MATLAB实现包括以下步骤：1. 输入数据准备：将原始EEG数据划分为两个类别的epochs，同样是根据任务要求分为不同的条件。
2. 将每个epoch的信号分解为一组频率子带：使用滤波器组将每个epoch的信号分解为一组频率子带。
可以使用一组带通滤波器，如Butterworth滤波器，将原始信号分解为不同频带。
3.对每个频率子带应用CSP：对于每个频率子带，应用CSP算法来提取空间模式。
可以使用CSP的MATLAB实现来实现这一步骤。
K-means算法matlab的实现
算法流程图三、实验源代码1、主程序clear allclc[FH FW]=textread('C:\Users\lenvo\Desktop\н¨Îļþ¼Ð\','%f %f'); [MH MW]=textread('C:\Users\lenvo\Desktop\н¨Îļþ¼Ð\','%f %f'); Data(1:50,1)=FH;Data(51:100,1)=MH;Data(1:50,2)=FW;Data(51:100,2)=MW;C=input('ÊäÈëC£º')；[U,P,Dist,Cluster_Res,Obj_Fcn,iter]=fuzzycm(Data,C)plot(Data(:,1), Data(:,2),'o');hold on;maxU = max(U);index1 = find(U(1,:) == maxU);index2 = find(U(2,:) == maxU);line(Data(index1,1),Data(index1,2),'marker','*','color','g');line(Data(index2,1),Data(index2,2),'marker','*','color','r');plot([P([1 2],1)],[P([1 2],2)],'*','color','k')hold off;2、子程序'function [U,P,Dist,Cluster_Res,Obj_Fcn,iter]=fuzzycm(Data,C,plotflag,M,epsm) if nargin<5epsm=;endif nargin<4M=2;endif nargin<3plotflag=0;end[N,S]=size(Data);m=2/(M-1);iter=0;、Dist(C,N)=0; U(C,N)=0; P(C,S)=0;U0 = rand(C,N);U0=U0./(ones(C,1)*sum(U0));¨while trueiter=iter+1;Um=U0.^M;P=Um*Data./(ones(S,1)*sum(Um'))';'for i=1:Cfor j=1:NDist(i,j)=fuzzydist(P(i,:),Data(j,:));endendU=1./(Dist.^m.*(ones(C,1)*sum(Dist.^(-m))));if nargout>4 | plotflagObj_Fcn(iter)=sum(sum(Um.*Dist.^2));endif norm(U-U0,Inf)<epsmbreak）endU0=U;endif nargout > 3res = maxrowf(U);for c = 1:Cv = find(res==c);Cluster_Res(c,1:length(v))=v;endendif plotflag^fcmplot(Data,U,P,Obj_Fcn);endfunction [U,P,Dist,Cluster_Res,Obj_Fcn,iter]=fuzzycm2(Data,P0,plotflag,M,epsm) if nargin<5epsm=;endif nargin<4M=2;endif nargin<3*plotflag=0;end[N,S] = size(Data); m = 2/(M-1); iter = 0;C=size(P0,1);Dist(C,N)=0;U(C,N)=0;P(C,S)=0;while trueiter=iter+1;for i=1:Cfor j=1:NDist(i,j)=fuzzydist(P0(i,:),Data(j,:));endend！U=1./(Dist.^m.*(ones(C,1)*sum(Dist.^(-m)))); Um=U.^M;P=Um*Data./(ones(S,1)*sum(Um'))';if nargout>4 | plotflagObj_Fcn(iter)=sum(sum(Um.*Dist.^2));endif norm(P-P0,Inf)<epsmbreakendP0=P;end-if nargout > 3res = maxrowf(U);for c = 1:Cv = find(res==c);Cluster_Res(c,1:length(v))=v;endendif plotflagfcmplot(Data,U,P,Obj_Fcn);end.function f = addr(a,strsort)if nargin==1strsort='ascend';endsa=sort(a); ca=a;la=length(a);f(la)=0;for i=1:laf(i)=find(ca==sa(i),1);ca(f(i))=NaN;endif strcmp(strsort,'descend'))f=fliplr(f);endfunction ellipse(a,b,center,style,c_3d)if nargin<4style='b';endif nargin<3 | isempty(center)center=[0,0];end:t=1:360;x=a/2*cosd(t)+center(1);y=b/2*sind(t)+center(2);if nargin>4plot3(x,y,ones(1,360)*c_3d,style)elseplot(x,y,style)endfunction fcmplot(Data,U,P,Obj_Fcn)~[C,S] = size(P); res = maxrowf(U);str = 'po*x+d^v><.h';figure(1),plot(Obj_Fcn)title('Ä¿±êº¯ÊýÖµ±ä»¯ÇúÏß','fontsize',8)if S==2figure(2),plot(P(:,1),P(:,2),'rs'),hold onfor i=1:Cv=Data(find(res==i),:);plot(v(:,1),v(:,2),str(rem(i,12)+1))ellipse(max(v(:,1))-min(v(:,1)), ...max(v(:,2))-min(v(:,2)), ...^[max(v(:,1))+min(v(:,1)), ...max(v(:,2))+min(v(:,2))]/2,'r:')endgrid on,title('2D ¾ÛÀà½á¹ûͼ','fontsize',8),hold off endif S>2figure(2),plot3(P(:,1),P(:,2),P(:,3),'rs'),hold on for i=1:Cv=Data(find(res==i),:);plot3(v(:,1),v(:,2),v(:,3),str(rem(i,12)+1))ellipse(max(v(:,1))-min(v(:,1)), ...，max(v(:,2))-min(v(:,2)), ...[max(v(:,1))+min(v(:,1)), ...max(v(:,2))+min(v(:,2))]/2, ...'r:',(max(v(:,3))+min(v(:,3)))/2)endgrid on,title('3D ¾ÛÀà½á¹ûͼ','fontsize',8),hold off endfunction D=fuzzydist(A,B)D=norm(A-B);：function mr=maxrowf(U,c)if nargin<2c=1;endN=size(U,2);mr(1,N)=0;for j=1:Naj=addr(U(:,j),'descend');mr(j)=aj(c);end四、实验结果1、FEMALE 和 MALE\(1)C=2, Z1(1)=(173,53)T, Z2(1)=(168,57)T。
模糊C均值聚类
主 单
讲：周润景 教授 位：电子信息工程学院
目 录
模糊C均值聚类应用背景 模糊C均值算法 模糊C均值聚类的MATLAB实现 模糊C均值聚类结果分析
一.模糊C均值聚类应用背景
传统的聚类分析是一种硬划分（Crisp Partition），它把每个待辨识的对 象严格地划分到某类中，具有“非此即彼”的性质，因此这种类别划分的界限 是分明的。然而实际上大多数对象并没有严格的属性，它们在性质和类属方面 存在着中介性，具有“亦此亦彼”的性质，因此适合进行软划分。Zadeh提出 的模糊集理论为这种软划分提供了有力的分析工具，人们开始用模糊方法来处 理聚类问题，并称之为模糊聚类分析。模糊聚类得到了样本属于各个类别的不 确定性程度，表达了样本类属的中介性，建立起了样本对于类别的不确定性的
三.模糊C均值聚类的MATLAB实现
426.31 3105.29 2057.8 1507.13 1556.89 1954.51 343.07 3271.72 2036.94 2201.94 3196.22 935.53 2232.43 3077.87 1298.87 1580.1 1752.07 2463.04 1962.4 1594.97 1835.95 1495.18 1957.44 3498.02 1125.17 1594.39 2937.73 24.22 3447.31 2145.01 1269.07 1910.72 2701.97 1802.07 1725.81 1966.35 1817.36 1927.4 2328.79 1860.45 1782.88 1875.13]; [center,U,obj_fcn] = fcm(data,4); plot3(data(:,1),data(:,2),data(:,3),'o');
fimplicit在matlab中的运用
fimplicit在matlab中的运用函数语法如下：其中，fun是用户定义的函数句柄，[xmin,xmax,ymin,ymax]定义了绘图的范围，n设置了网格大小，ax是指定图形绘制的轴。
在这个函数中，最重要的参数是fun，它定义了用户想要绘制的二元函数。
函数应该返回浮点数结果。
例如，如果想要绘制一个简单的二元函数y=x^2，则可以定义如下函数：function z = myfun(x, y)z=x^2-y;end然后，可以使用fimplicit函数来生成图形：这将绘制一个曲面，其中y轴表示y的值，x轴表示x的值，z轴表示myfun(x,y)的值。
另外，可以使用[f,g]=fun(x,y)的形式来返回辅助函数f和g的值。
这在一些情况下非常有用，因为fimplicit会计算并显示每个点上的交叉线。
通过检查辅助函数的值，可以确定函数在每个点上是否交叉。
为了更好地理解fimplicit的运用，下面是一个较为复杂的例子：1. 创建一个文件名为"myfun.m"的函数文件，内容如下：function z = myfun(x, y)z = sin(x) - cos(y);end2.在MATLAB中运行以下命令：title('fimplicit示例')xlabel('x轴')ylabel('y轴')zlabel('z轴')这将绘制出一个参数函数x^2+y^2-2=0的图像。
总而言之，fimplicit函数在MATLAB中的运用非常灵活，可以绘制各种二元函数的图形，并且可以通过调整参数来控制图像的范围和分辨率。
这对于理解和可视化数学函数在二维平面上的行为非常有用。
matlabfeedback函数用法
matlabfeedback函数用法Matlab中的feedback函数可以被用于控制系统分析和设计的简化。
feedback函数可以将系统的控制信号和系统输出信号进行操作合成，因此我们可以看到系统对于输入信号的控制效果。
feedback函数的函数原型如下：```sys_cl = feedback(sys, H, sign)```sys是一个包含系统状态空间模型的对象；H是包含反馈传递函数的对象；sign是反馈信号的方向。
在此函数中，反馈信号的方向是一个选择性参数，其值可以为正号或负号，可以控制反馈信号是系数为1还是-1。
如果sign参数定义为'+'，则表示反馈信号的系数为1；如果sign参数定义为'-'，则表示反馈信号的系数为-1。
反馈传递函数H是一个包含反馈传递函数的对象。
反馈传递函数可以指定控制系统的反馈路径，以便更好地控制系统。
反馈函数的传递函数通常是一个比例或积分传递函数。
有时也可以加入带通或低通滤波器等传递函数。
feedback函数的返回值是一个包含新系统状态空间模型的对象。
这就是我们可以用来控制系统的新模型。
在系统设计和分析中，feedback函数非常有用。
回归到feedback函数的工作原理，反馈函数的目的是将系统的输出信号反馈输入到系统中。
当反馈追踪路径较短时，即反馈信号可以追踪到控制信号之前，系统输出可以被控制。
控制系统在系统输出中引入反馈信号后，可以使系统变得更加稳定。
在电子电路中，我们常常用反馈来控制流经系统的电流、电压或功率等。
现在，我们将这个函数在一些实际的例子中实现。
让我们看看第一个实例：假设我们要设计一个标准的控制系统，用于控制各种机械设备的运动。
我们可以将系统的状态模型表示为：$$\begin{aligned}\dot{\mathbf{x}}(t)&=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t)\\\mathbf{y}(t)&=\mathbf{C}\mathbf{x}(t)+\mathbf{D}\mathbf{u}(t)\end{aligned}$$A、B、C、D是系数矩阵，x(t)是状态向量，u(t)是控制输入，y(t)是输出向量。
用matlab模拟FEC和交织两种方式
用matlab模拟FEC和交织两种方式FEC（前向纠错方式）：交织>>s1=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,3 4,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56];>>x1=(reshape(s1,8,7))';>>x1(2,2)=0;x1(2,3)=0;>>x1(2,4)=0;x1(2,5)=0;>> s2=reshape(x1,1,56);>>x2=reshape(s2,7,8);>>x2(2,2)=10;x2(2,3)=11;>> s3=reshape(x2',1,56);>> a=[s1,s2,s3];>>plot(s1,s2);x1 =1 2 3 4 5 6 7 89 0 0 0 0 14 15 1617 18 19 20 21 22 23 2425 26 27 28 29 30 31 3233 34 35 36 37 38 39 4041 42 43 44 45 46 47 4849 50 51 52 53 54 55 56L=1000;M=4;%每个符号的比特数N=2^M-1;%编码后码字长度K=N-4;%信息长度MSG=randint(L,1);%随机产生L比特信号TP=gftuple([-1:N-1]',M);%产生加罗华域元素PG=rspoly(N,K);%产生生成式[CODE,ADDED]=rsenco(MSG,TP,K);%编码NOI=rand(length(CODE)/M,1)<03;%加入3%的噪声NOI=(NOI*ones(1,M))';%产生突发错误NOI=NOI(:);CODE_NOI=rem(CODE+NOI,2);%噪声加入信号[DEC,ERR,CCODE,ERR_C]=rsdeco(CODE_NOI,TP,K);%译码MSG=[MSG;zeros(ADDED,1)];%调整长度max(abs(DEC-MSG));%比较X=[1:length(NOI)];Z=[1:M*N:length(NOI)];Y=zeros(1,length(Z));Z=[Z;Z];Y=[Y+min(ERR_C);Y+max(ERR_C)];subplot(211);plot(X,NOI,'yo',X,ERR_C,'rx',Z,Y,'g-');title('Error Detection Record');xlabel('o--placed error;x--detected error;vertical bar: RS-DECO section.'); axis([1,length(NOI),min(ERR_C),max(ERR_C)]);X=[1:length(MSG)];Z=[1:M*K:length(MSG)];Y=zeros(1,length(Z));Z=[Z;Z];Y=[Y;Y+max(MSG)];subplot(212);plot(X,MSG,'yo',X,DEC,'rx',Z,Y,'g-');title('Message and Decoded Signal Comparison');xlabel('o--original message;x--decoded result.');axis([1,length(MSG),min(min(MSG)),max(max(MSG))]);st1 = 27221; st2 = 4831; % States for random number generator n = 7; k = 4; % Parameters for Hamming codemsg = randint(k*500,1,2,st1); % Data to encodecode = encode(msg,n,k,'hamming/binary'); % Encoded data% Create a burst error that will corrupt two adjacent codewords.errors = zeros(size(code)); errors(n-2:n+3) = [1 1 1 1 1 1];% With Interleaving%------------------inter = randintrlv(code,st2); % Interleave.inter_err = bitxor(inter,errors); % Include burst error.deinter = randdeintrlv(inter_err,st2); % Deinterleave.decoded = decode(deinter,n,k,'hamming/binary'); % Decode.disp('Number of errors and error rate, with interleaving:');[number_with,rate_with] = biterr(msg,decoded) % Error statistics% Without Interleaving%---------------------code_err = bitxor(code,errors); % Include burst error.decoded = decode(code_err,n,k,'hamming/binary'); % Decode.disp('Number of errors and error rate, without interleaving:'); [number_without,rate_without] = biterr(msg,decoded) % Error statisticsNumber of errors and error rate, with interleaving:number_with =rate_with =Number of errors and error rate, without interleaving:number_without =4rate_without =0.0020>>msg=randint(k*500,1,2,st1);>>code = encode(msg,n,k,'hamming/binary');>>errors = zeros(size(code)); errors(n-2:n+3) = [1 1 1 1 1 1];>> inter = randintrlv(code,st2);>>inter_err = bitxor(inter,errors);>>deinter = randdeintrlv(inter_err,st2);>> decoded = decode(deinter,n,k,'hamming/binary');>>disp('Number of errors and error rate, with interleaving:');Number of errors and error rate, with interleaving:>>code_err = bitxor(code,errors);>> decoded = decode(code_err,n,k,'hamming/binary');>>disp('Number of errors and error rate, without interleaving:');Number of errors and error rate, without interleaving:>> [number_without,rate_without] = biterr(msg,decoded);L=1000;M=4;%每个符号的比特数N=2^M-1;%编码后码字长度K=N-4;%信息长度MSG=randint(L,1);%随机产生L比特信号TP=gftuple([-1:N-1]',M);%产生加罗华域元素PG=rsgenpoly(N,K);%产生生成式[CODE,ADDED]=rsenco(MSG,TP,K);%编码R=[0:0.01:1];%错误率for j=1:length(R) %循环修改错误率NOI=rand(length(CODE)/M,1)<R(j);%加入R%的噪声NOI=(NOI*ones(1,M))';%产生突发错误NOI=NOI(:);CODE_NOI=rem(CODE+NOI,2);%噪声加入信号[DEC,ERR,CCODE,ERR_C]=rsdeco(CODE_NOI,TP,K);%译码MSG=[MSG;zeros(ADDED,1)];%调整长度error=0;%错误码计数初值Q=length(DEC);%码长度for i=1:Q %循环寻找错误码if(MSG(i)~=DEC(i))error=error+1;endendp(j)=error/Q;%误码率。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
m文件1/7: function [U,P,Dist,Cluster_Res,Obj_Fcn,iter]=fuzzycm(Data,C,plotflag,M,epsm) % 模糊 C 均值聚类 FCM: 从随机初始化划分矩阵开始迭代 % [U,P,Dist,Cluster_Res,Obj_Fcn,iter] = fuzzycm(Data,C,plotflag,M,epsm) % 输入: % Data: N×S 型矩阵,聚类的原始数据,即一组有限的观测样本集, % Data 的每一行为一个观测样本的特征矢量,S 为特征矢量 % 的维数,N 为样本点的个数 % C: 聚类数,1% plotflag: 聚类结果 2D/3D 绘图标记,0 表示不绘图,为缺省值 % M: 加权指数,缺省值为 2 % epsm: FCM 算法的迭代停止阈值,缺省值为 1.0e-6 % 输出: % U: C×N 型矩阵,FCM 的划分矩阵 % P: C×S 型矩阵,FCM 的聚类中心,每一行对应一个聚类原型 % Dist: C×N 型矩阵,FCM 各聚类中心到各样本点的距离,聚类中 % 心 i 到样本点 j 的距离为 Dist(i,j) % Cluster_Res: 聚类结果,共 C 行,每一行对应一类 % Obj_Fcn: 目标函数值 % iter: FCM 算法迭代次数 % See also: fuzzydist maxrowf fcmplot if nargin<5 epsm=1.0e-6; end if nargin<4 M=2; end if nargin<3 plotflag=0; end [N,S]=size(Data);m=2/(M-1);iter=0; Dist(C,N)=0; U(C,N)=0; P(C,S)=0; % 随机初始化划分矩阵 U0 = rand(C,N); U0=U0./(ones(C,1)*sum(U0)); % FCM 的迭代算法 while true % 迭代计数器 iter=iter+1; % 计算或更新聚类中心 P Um=U0.^M; P=Um*Data./(ones(S,1)*sum(Um'))'; % 更新划分矩阵 U for i=1:C for j=1:N Dist(i,j)=fuzzydist(P(i,:),Data(j,:)); end end U=1./(Dist.^m.*(ones(C,1)*sum(Dist.^(-m)))); % 目标函数值: 类内加权平方误差和 if nargout>4 | plotflag Obj_Fcn(iter)=sum(sum(Um.*Dist.^2)); end % FCM 算法迭代停止条件 if norm(U-U0,Inf) break end U0=U; end % 聚类结果 if nargout > 3 res = maxrowf(U); for c = 1:C v = find(res==c); Cluster_Res(c,1:length(v))=v; end end % 绘图 if plotflag fcmplot(Data,U,P,Obj_Fcn); end
m文件2/7: function [U,P,Dist,Cluster_Res,Obj_Fcn,iter]=fuzzycm2(Data,P0,plotflag,M,epsm) % 模糊 C 均值聚类 FCM: 从指定初始聚类中心开始迭代 % [U,P,Dist,Cluster_Res,Obj_Fcn,iter] = fuzzycm2(Data,P0,plotflag,M,epsm) % 输入: Data,plotflag,M,epsm: 见 fuzzycm.m % P0: 初始聚类中心 % 输出: U,P,Dist,Cluster_Res,Obj_Fcn,iter: 见 fuzzycm.m % See also: fuzzycm if nargin<5 epsm=1.0e-6; end if nargin<4 M=2; end if nargin<3 plotflag=0; end [N,S] = size(Data); m = 2/(M-1); iter = 0; C=size(P0,1);Dist(C,N)=0;U(C,N)=0;P(C,S)=0; % FCM 的迭代算法 while true % 迭代计数器 iter=iter+1; % 计算或更新划分矩阵 U for i=1:C for j=1:N Dist(i,j)=fuzzydist(P0(i,:),Data(j,:)); end end U=1./(Dist.^m.*(ones(C,1)*sum(Dist.^(-m)))); % 更新聚类中心 P Um=U.^M; P=Um*Data./(ones(S,1)*sum(Um'))'; % 目标函数值: 类内加权平方误差和 if nargout>4 | plotflag Obj_Fcn(iter)=sum(sum(Um.*Dist.^2)); end % FCM 算法迭代停止条件 if norm(P-P0,Inf) break end P0=P; end % 聚类结果 if nargout > 3 res = maxrowf(U); for c = 1:C v = find(res==c); Cluster_Res(c,1:length(v))=v; end end % 绘图 if plotflag fcmplot(Data,U,P,Obj_Fcn); end m文件3/7: function fcmplot(Data,U,P,Obj_Fcn) % FCM 结果绘图函数 % See also: fuzzycm maxrowf ellipse [C,S] = size(P); res = maxrowf(U); str = 'po*x+d^v><.h'; % 目标函数绘图 figure(1),plot(Obj_Fcn) title('目标函数值变化曲线','fontsize',8) % 2D 绘图 if S==2 figure(2),plot(P(:,1),P(:,2),'rs'),hold on for i=1:C v=Data(find(res==i),:); plot(v(:,1),v(:,2),str(rem(i,12)+1)) ellipse(max(v(:,1))-min(v(:,1)), ... max(v(:,2))-min(v(:,2)), ... [max(v(:,1))+min(v(:,1)), ... max(v(:,2))+min(v(:,2))]/2,'r:') end grid on,title('2D 聚类结果图','fontsize',8),hold off end % 3D 绘图 if S>2 figure(2),plot3(P(:,1),P(:,2),P(:,3),'rs'),hold on for i=1:C v=Data(find(res==i),:); plot3(v(:,1),v(:,2),v(:,3),str(rem(i,12)+1)) ellipse(max(v(:,1))-min(v(:,1)), ... max(v(:,2))-min(v(:,2)), ... [max(v(:,1))+min(v(:,1)), ... max(v(:,2))+min(v(:,2))]/2, ... 'r:',(max(v(:,3))+min(v(:,3)))/2) end grid on,title('3D 聚类结果图','fontsize',8),hold off end
m文件4/7: function D=fuzzydist(A,B) % 模糊聚类分析: 样本间的距离 % D = fuzzydist(A,B) D=norm(A-B); m文件5/7: function mr=maxrowf(U,c) % 求矩阵 U 每列第 c 大元素所在行,c 的缺省值为 1 % 调用格式: mr = maxrowf(U,c) % See also: addr if nargin<2 c=1; end N=size(U,2);mr(1,N)=0; for j=1:N aj=addr(U(:,j),'descend'); mr(j)=aj(c); end
m文件6/7: function ellipse(a,b,center,style,c_3d) % 绘制一个椭圆 % 调用: ellipse(a,b,center,style,c_3d) % 输入: % a: 椭圆的轴长(平行于 x 轴) % b: 椭圆的轴长(平行于 y 轴) % center: 椭圆的中心 [x0,y0],缺省值为 [0,0] % style: 绘制的线型和颜色,缺省值为实线蓝色 % c_3d: 椭圆的中心在 3D 空间中的 z 轴坐标,可缺省 if nargin<4 style='b'; end if nargin<3 | isempty(center) center=[0,0]; end t=1:360; x=a/2*cosd(t)+center(1); y=b/2*sind(t)+center(2); if nargin>4 plot3(x,y,ones(1,360)*c_3d,style) else plot(x,y,style) end
m文件7/7: function f = addr(a,strsort) % 返回向量升序或降序排列后各分量在原始向量中的索引 % 函数调用:f = addr(a,strsort) % strsort: 'ascend' or 'descend'