地理聚类与判别分析

第5章地理聚类与判别分析地理聚类与类型判别分析是地理学研究的重要内容这一，现代地理学中的定量聚类与判别分析更为重要。第1节地理聚类分析传统地理学分类法已不适应现代地理学对地理分类问题的需要，现介绍定性与定量相结合的现代地理分类法，即地理聚类分析法，主要介绍地理聚类分析法的基本概念、聚类统计量、系统聚类法的基本原理和方法步骤等。地理聚类分析的概述概念地理聚类分析是指根据地理变量的属性或特征值求算出能表明分类对象之间相似性或亲疏程度的聚类统计量，用数学方法把它们逐步分型划类，最后得到一个能反映分类对象之间亲疏关系的客观的分类系统，即，地理定量分类法结果，其含义下：1）地理聚类分析法的理论基础是数理统计原理；）地理聚类分析法的数据基础是原始地理数据，它是聚类结果形成的必要条件；）地理聚类分析法的主要方法是聚类统计量和聚类过程求算的n*!-、■夕*.亠数学方法；）地理聚类分析法的核心任务是聚类统计量求算和逐步聚类过程；）地理聚类分析法的最终结果是聚类谱系图，即，能表明分类对象之间亲疏关系的图示。地理聚类分析的一般过程是：明确问题，收集资料，求算聚类统计量，逐步系统聚类，地理分析等等。主要用途：主要用于解决地理分类问题，如区域规划、农业区划、地貌分类、气候分类、土壤分类等等。现已广泛应用，但是，要注意定性和定量分析相结合，才会得到更好的分类结果。特点地理聚类分析无需知道分类对象的分类结构，而只需要表示一批分类对象的属性或特征值的地理数据，然后，选好分类统计量和聚类方法，并按规定的方法步骤逐步计算聚类，最后便能自然而然地客观地得出一份完整的分类谱系图，而且还能根据求得的不同统计量尺度划分出多级分类结构，具体特点如下：1）系统性具有多要素、多属性、多层次结构；2）定量性以地理数据为基础，以求算得的分类统计量为标度，进行定量化聚类；3）自然性通过聚类过程自然而然形成聚类结果；4）直观性聚类过程和聚类系统谱系图均直观可视；5）实用性非常科学、快捷、使用。区别是指与传统分类法的区别1）传统分类法靠经验定性描述定性结果；

2）聚类分析法靠数学方法和计算机技术定量分析定量结果。数据处理（略）聚类统计量聚类统计量是指能度量分类对象间相似程度或相关程度的数量指标。常用的有相似系数和距离系数两种。它是聚类分析的依据，也是聚类分析的必要条件。相似系数1）概念相似系数是指用来度量分类对象间相似程度的数量指标，它又分为相关系数和夹角余弦两种。2）类型（1）夹角余弦COS9ij概念COS9COS9ij=ijx2・乙x2tikjkk二1k二1其中，i,j表示分类对象k表示i,j的第k个特征值性质取值范围：—1<COS9<1ij聚类准则：COS9越大,表明i与j两地点越相似,应先归为ij一类。当cos0=1时，则COS0=0。，表明两地点间完全相似，但向ijij量长度不一定相等。当COS0=0时，则Cos0=90。，表明两地点间完全不同。ijij当COS0<0时，表明两地点间为负相似。ij矩阵形式夹角余弦的矩阵形式是指将分类对象间两两夹角余弦求出后写成矩阵形式，即：/COS0COS0COS0•COS0)1112131nCOS0COS0COS0•COS02122232n0=COS0COS0COS0•COS0:31:32:33.•3nCOS0COS0COS0•COS0n1n2n3nn性质：为实对称矩阵，这样只需求出上三角或者下三角矩阵即可COS0=COS0=COS0=…=COS0=1112233nn注：数据处理后再求夹角余弦为宜。（2）相关系数略）说明：相关系数值大为聚类准则。距离系数1）概念距离系数是指在M维空间系统中，以两点间的抽象距离来表示分类对象间相似程度的统计量指标称为距离系数，即将分类对象的特征

值看作为M维空间的点，而求得M维空间两点间的距离，也可以说是抽象距离，不是地理距离。2）类型距离系数的种类很多，常用的有如下几种：（1）欧氏距离d=..'X（X-X）2ij[ikjk丫k=1其中，I,j二1,2,3,・„,n为空间任意两点;k=1,2,3,„,m为指标个数。（2）相对欧氏距离1yz、d=（x一x）2ijmikjkk=13）绝对值距离d=Xx一xijikjkk=11/p4）明科夫斯基距离1/p^mIdij乙x一dijikjkk=15）切比雪夫距离d=maxx一xijkikjk说明：聚类准则：距离系数值越小，则两点间相似程度越大；反之，则越小，故可根据值大小逐步聚类。3距离系数矩阵1）一般形式若将所有两点间距离系数求出，即可写成矩阵形式为：fdfddd111213ddd212223D=ddd313233dddn1n2n3d)1nd2nd3ndnn这样只需求出上三角或者下三角2）性质：这样只需求出上三角或者下三角为实对称矩阵，即，d=dijji矩阵即可。=d=0

nn对角线上的元素：d=d==d=0

nn112233注：原始数据处理后再求算距离系数为宜。3）举例以表7—2的数据为例，用相对欧氏距离系数公式来求距离系数矩阵。地点八、、原始数据自然对数变换标准差标准化数据海拔年平均最大冻工8级海拔年平均最大冻工8级海拔年平均最大冻工8级高度降水量土深度大风日数高度降水量土深度大风日数高度降水量土深度大风日数X1X2X3X4X1X2X3X4X1X2X3X4哈巴河532.6173.8150.061.86.2785.1585.0114.124-0.9101.0330.6580.856阿勒泰735.1191.5146.037.76.6005.2554.9843.630-0.2501.1670.5970.541克拉玛依427.0114.4197.075.46.0574.7405.2834.323-1.3630.4551.2710.996巴楚1116.541.661.07.67.0183.7284.1112.0230.607—0.946-1.372-0.511莎车1231.242.593.011.07.1163.7504.5332.3980.808-0.915-0.420-0.268于田1427.046.481.01.47.2633.8374.3940.3361.100—0.795—0.734-1.623六个站点：1，哈巴河；2，阿勒泰；3，克拉玛依；4，巴楚；5，沙车；6，于田。四个指标：X，海拔高度；X，年平均降水量；x，最大冻土深123度；x4，8级大风日数。d=121I.1q1v乙(X-X)2=乙(X-X)2'mikjk41k2kk—1k—1-,11一\4[X-X)2+(X-X)2+(X-X)2+(X-X)2]1121122213231424一山T4[—0.910—(—0.250))2+(1.033—1.167)2+(0.658—0.597)2+(0.865—0.541)2]=0.375其中，1=1,哈巴河；J=2，阿勒泰；K=1,2,3,4。同理可求得两两站点间的相对欧氏距离系数矩阵，于是可求得初始距离系数矩阵为：100.3750.4831.7491.5161.972]0.37500.7761.5961.3361.743D(0)—0.4830.77601.9261.6622.1541.7491.5961.92600.5010.6931.5161.3361.6620.50100.7141.9721.7432.1540.6930.7140J说明：相似系数和距离系数是系统聚类的基础,即在此基础上逐步聚类的。聚类分析方法1概述1）聚类对象类型根据分类目的不同,可分为两类（1）Q型Q型是指分类对象为地点或地区,即对地点或地区进行分类,例如：农业区划,农业自然指标,农业经济指标,农业人文指标。2）R型R型是指分类对象为地理要素，即对某地理要素进行分类，例如：土壤分类，气候分类，地貌分类，城市分类等等。2）聚类方法聚类方法有许多种，例如：系统聚类法，逐步聚类法，逐步分解聚类法，最优分割聚类法等等。现在主要介绍系统聚类法。2系统聚类法1）系统聚类法的基本原理系统聚类法的基本原理是：首先，将N个聚类对象看作各自成一类；其次，选择距离系数最小的或相似系数最大的一对合成为一新类，并根据预定的八种方法之一求算性类与未聚类者之间的距离系数或相似系数，以形成新的距离系数矩阵或相似系数矩阵；然后，根据新形成的矩阵再选择距离系数最小的或相似系数最大的一对合成为一新类，这样每合并一次就减少一类，直至N个分类对象都合并为一类为止。说明：（1）聚类准则以距离系数最小或相似系数最大为原则。注：每一步都是如此。（2）聚类过程是逐步自然而然形成的。（3）聚类结果是直观的聚类谱系图。（4）求D⑴，D（2）,D⑶，…，D（n）是表示新的距离系数矩阵，其组成是：一是未聚类者的距离系数；二是新类与未聚类者经求算而

得的距离系数。这里要注意的关键是：新类与未聚类者的距离系数如何求算呢？其求算方法有八种，即：类平均法，离差平方和法。最短距离法，最长距离法，类平均法，离差平方和法。可变法，可变平均法，重心法以上方法的通用公式为：D2=aD2+ad2+pd2+rD2一D2krpkpqkqpqkpkq其中，D:表示新类与旧类的距离系数；krr:表示新类；k:表示未聚类者，即旧类；p,q:表示合成新类的两类；a,a,p,r:表示为不同方法的系数；如表所示，(P):pq135n,n,n,n:表示为各类的样本个数。pqkr此公式为编程提供了方便。求R求R(1)，R(2)，R(3)，…,R(n)类同。系统聚类法的步骤例1根据地理聚类分析的原理和方法，现表7.2地理数据为例介绍用最短矩离法进行聚类分析的具体步骤：PP131原始数据原始数据自然对数变换标准差标准化数据原始数据原始数据自然对数变换标准差标准化数据地点海拔地点海拔年平均最大冻工8级海拔年平均高度降水量土深度大风日数高度降水量x1x2x3x4x1x2哈巴河532.6173.8150.061.86.2785.158阿勒泰735.1191.5146.037.76.6005.255克拉玛依427.0114.4197.075.46.0574.740巴楚1116.541.661.07.67.0183.728莎车1231.242.593.011.07.1163.750于田1427.046.481.01.47.2633.837最大冻28级海拔年平均最大冻28级土深度大风日数高度降水量土深度大风日数xxxxxx3412345.0114.124-0.9101.0330.6580.8564.9843.630-0.2501.1670.5970.5415.2834.323-1.3630.4551.2710.9964.1112.0230.607—0.946-1.372-0.5114.5332.3980.808-0.915-0.420-0.2684.3940.3361.100—0.795—0.734—1.623第一步数据处理将表7.2的原始地理数据先进行对数处理再进行标准差标准化处理如表7.2所示：第二步求D(0)，即初始距离系数矩阵根据相对欧氏距离公式求得D(0)G1G1G2G3/00.3750.4830.37500.776D(0)=0.4830.77601.7491.5961.9261.5161.3361.662、1.9721.7432.154GGG4561.7491.5161.972、(G)11.5961.3361.743G21.9261.6622.154G300.5010.693G40.50100.714G50.6930.7140GJ6丿第三步逐步聚类1)以d(0)聚类，并求d⑴由于d]2=0.3751为最小，故将G]与G2合并为新类G7,并求余类与新类g7的距离系数可得d(i)。这^时p=1,q=2,r=7,k=3,4,5,6d=mid=mi4dn,d}=m371323d=mi4dn,d}=m471424d=mi4dn,d}=m571525d=mi4dn,d}=m671626GGGGG34567r01.92601.66082.15220.4824、1.926000.50050.69191.5941D(1)=1.66080.500500.71311.33552.15220.69190.713101.7416、0.48241.59411.33551.74160丿urg3G4G5G6IG7TOC\o"1-5"\h\zi{仇.482Q.77兴=50.4824i417476.594=11.5941i41515B335=51.3355i41970B.741=61.7416于是可得d⑴:说明：1）最长距离法的求算方法是：d=max{d,d}=max{0.4824,0.7755}=0.7755TOC\o"1-5"\h\z371323d=max4/,d}=max6.7476,1.5941}=1.7476471424(2)中心距离法的计算方法是：d=mid4/,d}=mid{0.4824,0.7755}=—(0.4824+0.7755)=0.62903713232d=mid4/,d}=midh.7476,1.5941}=1(1.7476+1.5941)=1.6709471424余者方法类同。2）以D（1）聚类，并求D（2）由于d=0.4824为最小，故将G与G7合并为新类G8，并求余类37378与新类G8的距离系数可得D（2）。8这时=&k=4,5,6这时=&k=4,5,6d=mi4dn,Dd}=m483447d=mi4dn,Dd}=m583557d=mi4dn,Dd}=m683667d=mi4dn,Dd}=m6716263,q=7,ri41926D.594=H.5941i{ln.66,D.335=51.3355i42n.15，2D12.741=61.7416i41n.97，0D15.741=61.7416于是可得D（2）:G4G5/00.50050.50050D(2)=0.69190.7131、1.59411.3355GG680.69191.5941、rg)40.71311.3355GU501.7416G61.74160丿:G8丿3）以D（2）聚类，并求D⑶由于d45=0.5005为最小，故将G4与G5合并为新类G9,并求余类与新类G9的距离系数可得D⑶这时这时pd=mi4dn,Dd}=m694656d=mi4dn,Dd}=m894858=4,q=5,r=9,k=6,8i40n.69，1D09.713=10.6919i41n.59，4D11.335=51.3355于是可得D⑶:G6G8G9r01.74160.6919、rg)6D(3)=1.741601.3355UG8、0.69191.33550J:G9丿4）以D⑶聚类，并求D⑷由于d69=0.6919为最小，故将G6与G9合并为新类Gio,并求余696910类与新类Gi0的距离系数可得D⑷】p=6,q=9,r=10,k=8d=miid,d}=mi{]n.74116.335=51.33558106889于是可得D⑷：GG810/01.3355、rg)D(4)=U81.33550G‘丿'10丿5)以D⑷聚类由于d=1.3355为最小，故将G8与G合并为新类G,最后全81081011部聚为一类，即告结束。第四步求绘谱系图根据上述聚类过程，可得联结表。聚类过程联结表联结顺序联结法距离系数新类类另U1G7GG12d=0.3751d122G8GG37d=0.4824373G9GG45d=0.5005454G10GG69d=0.6919695G11GG810d=1.3355810由联结表可得聚类谱系图如图7—1所示。p137注意1聚类过程的联系实意理解2确定各类在横坐标位置具体方法是：从大类到小类逐层寻找即可，如：G11(G8G)；10G(GG)；1069G(8G3G)等7第五步地理分析由聚类谱系图可知！)当d=0.5时'可分为四类，即：G,G,G123G;G；G452)当d=0.7时'可分为三类，即：G,G,G；123G,G；G4563)当d=0.8时，可分为二类，即：G,G,G；G,G,G该研究对象为新疆的六个地区，从图和实际情况可知，该地区可分为两大类：一类为准噶尔盆地类，一类为塔里木盆地类。在准噶尔盆地类又分为两小类型：一为哈巴河—阿勒泰型，一为克拉玛依型。在塔里木盆地类型中也分为两小类型：一为巴楚---莎车类型，一为于田类型。聚类结果与当地的自然地理环境完全一致。例2根据地理聚类分析的原理和方法，现介绍用最短矩离法进行聚类分析的具体步骤。以表5.1中地理数据为例。表5.1原始地理数据\变xXxxX量\12345

''地点11032.270030025221545.07002602523269.0130100102.2第一步数据处理将表5.1的原始地理数据进行标准差标准化处理得：表5.2标准差标准化数据\变量地点、x1X2x3x4X510.1299940.5773500.7559290.577350—0.57735220.9286470.5773500.3779640.577350—0.5773523—1.05864—1.15470—1.13389—1.154701.15470第二步求初始距离系数矩阵D(°)根据欧氏距离公式求得D(°)为:D（o）G1G2G3G100.8835743.73956G20.88357403.90321G33.739563.903210第三步逐步聚类以最短距离法进行逐步聚类，最长距离法、中心法等类同。(1)以d(°)求d(1)1111由于d=0.883574为最小，故将G与G合并为新类G，并求余21124类与新类G的距离系数，以最短距离法求之，即：4d=minid4331d}=min{3.73956,3.90321d=minid433132于是可得D（J为:D(1)G3G4G303.73956G43.739560（2）以D（1）聚类以D（1）可知，只可将G与G合并为新类G，则即可结束。D1345第四步求绘制聚类谱系图由上述聚类过程可得联结表联结顺序联结法距离系数新类类另1G4GG120.8835742G5GGG1233.739560图5－1聚类谱系图由联结表可绘制聚类谱系图图5－1聚类谱系图第五步地理分析由聚类谱系图可易知：与G；G。与G；G。23当d=3.739560时，可合并为一类，即G、G、G同属一类。123例3数据表根据上述聚类过程，可得联结表。聚类过程联结表联结顺序联结法距离系数新类类别1GGG0.21195343220142GGG0.224346633873GGG0.22695263418174GGG0.31404543523345GGG0.31714743635166GGG0.31988693732127GGG0.36433643825378GGG0.3688843393639GGG0.370276340383910GGG0.405868341302911GGG0.406433642274012GGG0.448102743151013GGG0.45408154433514GGG0.456624145444215GGG0.4644023461911

16GGG0.466426347224517GGG0.55038534844718GGG0.634460749434819GGG0.656208650282420GGG0.667750151134921GGG0.671997052465122GGG0.6762108539123GGG0.696616754505224GGG0.74276735565425GGG0.769757156315526GGG0.77226175725327GGG0.946027558415628GGG0.955551659585729GGG1.041591060215930GGG1.3054060612660由联结表可得聚类谱系图如图7—1所示。地理聚类分析的实习指导1.5.1实习目的1、巩固地理聚类分析的基本原理与方法步骤。2、掌握地理聚类分析程序的使用方法与技巧。3、求取地理分类结果及应用于与实际进行对比分析4、掌握地理聚类分析程序的变换应用方法。1.4.2实习内容1、标识符说明N对象个数M变量个数X(N，M)存放原始地理数据矩阵R(N,N)存放下三角距离矩阵L1存放数据处理方法代码当L1=1时，原始数据中心化当L1=2时，原始数据对数变换当L1=3时，原始数据极差标准化当L1=4时，原始数据标准差标准化Bl八种聚类方法参数，-丄<B1<04B2八种聚类方法参数，B2V1K存放聚类方法代码当K=1时，最短距离法当K=2时，最长距离法当K=3时，中间距离法当K=4时，可变法当K=5时，类平均法当K=6时，可变类平均法当K=7时，离差平方和法当K=8时，重心法2、程序5REM聚类分析程序10INPUT“对象数N,样本数M=”；N,M20DIMX(N,M),S(M)30FORI=1TON40FORJ=1TOM50READX(I,J)55PRINTX(I,J);60NEXTJ65PRINT70NEXTI75REM对原始数据进行标准化80INPUT“L1=”；L190IFL1=1THEN130100IFL1=2THEN270110IFL1=3THEN350120IFL1=4THEN500125GOTO650126REM中心化变换130FORJ=1TOM140A=0160FORI=1TON170A=A+X(I,J)180NEXTI190S(J)=A/N200NEXTJ210FORJ=1TOM220FORI=1TON230X(I,J)=X(I,J)－S(J)240NEXTI250NEXTJ260GOTO650265REM对数变换270RORJ=1TOM280FORI=1TON290IFX(I,J)<>0THEN310300X(I,J)=1E－10310X(I,J)=LOG(X(I,J))320NEXTI330NEXTJ340GOTO650345REM极差正规化变换350FORJ=1TOM360A=X(1,J)370B=X(1,J)380FORI=2TON

420430440450460470480490495500510520530540550560570580590600606610620625630635640645650700705710390IFX(I,J)<=ATHEN410400A=X(I,J)410IFX(I,J)>=BTHEN430B=X(I,J)NEXTIA=A－BFORI=1TONX(I,J)=(X(I,J)－B)/ANEXTINEXTJGOTO650REM标准差标准化FORJ=1TOMA=0B=0FORI=1TONA=A＋X(I,J)NEXTIA=A/NFORI=1TONB=B+(X(I,J)-A)f2NEXTIB=SQR(B/(N-1))PRINT“StandarXMatrix”FORI=1TONX(I,J)=(X(I,J)-A)/BPRINTX(I,J)；NEXTIPRINTNEXTJREM形成欧氏距离矩阵DIMR(N,N),W(N),Z(N)PRINT“DistanceDMatrix”PRINTFORI=2TON

72073074075076076577077578078580081082083084085086087088089090091092093094095096097010001010102010401050FORJ=1TOI－1R(I,J)=0FORK=1TOMR(I,J)=R(I,J)+(X(I,K)—X(J,K))f2NEXTKPRINTSQR(R(I,J))；NEXTJPRINT可变法的参数NEXTI可变法的参数REM计算最短距离法，最长距离法，中间距离法INPUT“B1,B2=”；B1，B2DIMP(8),Q(8),U(8),V(8)P(1)=.5Q(1)=.5U(1)=0V(1)=—.5P(2)=.5Q(2)=.5U(2)=0V(2)=.5P(3)=.5Q(3)=.5U(3)=B1V(3)=0P(4)=(1—B2)/2Q(4)=(1—B2)/2U(4)=B2V(4)=0INPUT“K=”；KFORI=1TONW(I)=1Z(I)=INEXTI1060FORT=1TON－11065REM找出距离矩阵中的最小元素

1070A=R(2,1)1080I0=21090J0=11100FORI=2TON－T+11110FORJ=1TOI－11120IFR(I,J)>=ATHEN11601130A=R(I,J)1140I0=I1150J0=J1160NEXTJ1170NEXTI1180PRINT“Step”；T1190PRINT“G”；N+T，“i=”；Z(I0)，“j=”；Z(J0)，“Distance”；SQR(A)1200IFK<>5THEN12301210GOSUB19001220GOTO14001230IFK<>6THEN12601240GOSUB19501250GOTO14001260IFK<>7THEN14001270GOSUB20001400FORI=J0+1TON－T+11405IFK<>8THEN14101406GOSUB20501410IFI>=I0THEN14401415REM每次并类后，形成新的距离矩阵1420R(I,J0)=P(K)*R(I,J0)+Q(K)*(I0,I)+U(K)*R(I0,J0)+V(K)*ABS(R(I,J0)－R(I0,I))1430GOTO14601440IFI<=I0THEN14601450R(I,J0)=P(K)*R(I,J0)+Q(K)*R(I,I0)+U(K)*R(I0,J0)+V(K)*ABS(R(I,J0)－R(I,I0))1460NEXTI1470FORI=1TOJ0－11475IFK<>8THEN1480

1476148014901495150015101520153015401545155015601570158015901600161016201630164016501660178017901800181018201830184018801890189519001910GOSUB2050R(J0,I)=P(K)*R(J0,I)+Q(K)*R(I0,I)+U(K)*R(I0,J0)+V(K)*ABS(R(J0,I)－R(I0,I))NEXTIIFI0>=N－TTHEN1545FORI=I0+2TON－T+1FORJ=I0TOI－2R(I,J)=R(I,J+1)NEXTJNEXTIIFI0>N－TTHEN1780FORI=I0TON－TFORJ=1TOI－1R(I,J)=R(I+1,J)NEXTJNEXTIPRINT“AfterChangingDMatrix”FORI=2TON－TFORJ=1TOI－1PRINTSQR(R(I,J))；NEXTJPRINTNEXTIW(Z(J0))=W(Z(J0))+W(Z(I0))W(Z(I0))=0IFI0>N－TTHEN1840FORI=I0TON－TZ(I)=Z(I+1)NEXTINEXTTPRINT“****************************”ENDrem计算类平均法，可变类平均法，重心法，离差平方和法参数的子程序P(5)=W(Z(J0))/(W(Z(J0))+W(Z(I0)))Q(5)=1－P(5)

1920U(5)=01930V(5)=01940RETURN1950P(6)=(1－B2)*W(Z(J0))/(W(Z(J0))+W(Z(I0)))1960Q(6)=(1－B2)*W(Z(I0))/(W(Z(J0))+W(Z(I0)))1970U(6)=B21980V(6)=01990RETURN2000P(7)=W(Z(J0))/(W(Z(J0))+W(Z(I0)))2010Q(7)=1－P(7)2020U(7)=－P(7)*Q(7)2030V(7)=02040RETURN2050A1=W(Z(I))+W(Z(I0))+W(Z(J0))2060P(8)=(A1－W(Z(I0)))/A12070Q(8)=(A1－W(Z(J0)))/A12080U(8)=－W(Z(I))/A12090V(8)=02100RETURN2210DATA1032.2,700,300,25,22220DATA1545,700,260,25,22230DATA269,130,100,10,2.21.5.3实习过程具体过程是：首先进入QBASIC状态，然后输入聚类分析程序；接着运行程序。运行程序，屏幕显示：对象数N,变量数M=?于“？”后输入3,5/屏幕显示：L1=？于“？”后输入4/屏幕显示！数据处理方法为标准差标准化B1,B2=？于“？”后输入－0.25，－0.5/屏幕显示：k=？于“？”后输入2/！聚类方法为最长距离法屏幕显示运行结果如下：RUN/L1=?4/Bl,B2=?-.25,—0.5/K=?2/StandardXMatrix.l29994.57735.755929.57735—.577352.928647.57735.377964.57735—.577352—l.05864—l.l547—l.l3389—l.l547l.l547DistanceDMatrix.8835743.739563.9032lSteplG4i=2j=1Distance0.883574Step2G5i=3j=1Distance3.903211.5.4程序的变换应用1、准备数据已知6组浅成岩均值和7组喷出岩均值数据如下表5.3所示，试用聚类分析程序将其分类并作出谱系图。表5.313块岩浆岩标本的8个氧化物含量序号岩石名称SiO2AgFe2O3FeOMgOCaONa?。K1玄武岩48.43514.7423.1997.3128.2819.2262.7640.9972辉长辉绿岩48.62316.0653.4886.1486.9778.1102.9971.9483安山玄武岩53.21116.5745.5254.0824.7427.6313.4311.4354闪长岩57.05617.3903.7473.6842.9545.8454.2502.0725安山岩59.51317.5212.6083.7353.4276.0412.8611.8796石英闪长岩61.62916.4691.5784.4262.8924.9173.4662.6817英安岩66.61915.2212.0221.9561.2103.1384.1472.2498花岗闪长玢LJU岩65.99915.1991.9172.4091.7882.7763.8543.4519流纹岩73.72713.2321.3160.9110.3401.0863.7423.97510花岗斑岩73.74213.5470.6841.6010.5790.9793.1634.50511粗面安山岩59.49117.6863.5982.3621.4413.8275.4814.30212粗面岩62.67017.1312.2961.9190.7832.4923.5726.92013正长斑岩63.91216.1892.0552.0371.2942.1224.2595.7982、编辑程序(1)对象数N,变量数M,数据处理方法L1聚类参数Bl、B2,聚类方法K等的值N=13M=8(2)DATA语句将DATA语句中数据均删除，重新输入所需表4.1数据，DATA语句的格式如：2210DATA48.435,14.742,3.199,7.312,8.281,9.226,2.764,0.9972220DATA48.623,16.065,3.488,6.148,6.977,8.110,2.997,1.9482230DATA2240DATA2250DATA3、运行程序运行程序，观察运行结果，并将两次运行结果作以比较。第2节地理判别分析地理判别分析对地理类型划分和区界判定等都具有十分重要的作用，现介绍其概念类型判别准则判别分析的基本原理和具体步骤2.1地理判别分析的概述1概念地理判别分析是指以某种判别准则为依据，根据某些地理对象的多种特征值求得判别函数和判别临界值，然后用来分辨和判别其类型归属问题的多变量数理统计方法，其含义如下：1）地理判别分析的理论基础是数理统计理论和方法；2）地理判别分析的关键任务是求得判别函数和判别临界值；3）地理判别分析的主要方法是最小二乘法；4）地理判别分析的充要条件是判别函数的显著性检验；5）地理判别分析的主要结果是定量结果，即类型归属表或类型判别表。2作用1）判别判别某地理对象的类型归属；2）定界确定区域间的界线；3）检验对原已确定类型的地理对象进行合理性检验；4）确定变量的贡献率求算判别地理对象的特征值的贡献率，即：各个变量在类型判别中的作用的大小。以示其作用大小来确定各特征变量在判别分析的作用，确定主要变量，优化判别函数。区别这是指与聚类分析的区别。1）从性质上聚类分析的属性是地理分类；判别分析的属性是判别与分类兼之，但重在判别其类型归属。2）从条件上聚类分析不必事先确定类型，而类型的形成是聚类分析的必然结果；判别分析必须事先知道分几种类型为前提，来判别未归类者应归属哪一类型。准则准则是指求算判别函数和判别临界值的原则，其种类有许多，常用的有费歇准则和贝叶斯准则。1）费歇准则费歇准则的要点是：类间或组间均值差尽量大，类内或组内离差平方和尽可能小，则其两者比值为最大，这样才能把地理类型区分得最清楚，故称为两组间最大分离原则。2）贝叶斯准则贝叶斯准则的要点是：把已知地理对象的地理数据分成几类或几组，然后计算未知地理类型或区域归属于各已知类型或组的概率值，它归属哪类的概率值最大，就把它划归为该类或该组。另外，还可以计算出划归各已知类的错分损失，其错分哪一类的平均损失为最小，就把它判归为该类，故称为条件概率原则。类型判别分析的类型有：两类判别分析，多类判别分析，逐步判别分析等等。现仅介绍两类判别分析iiii2.2两类判别分析1两类判别分析的基本原理1）一般表达式两类判别分析的判别函数的一般表达式为：y=ex+ex+ex+…+cx=乙ex112233mmkk其中，c,c,c,...,c123m其中，c,c,c,...,c123m为判别函数系数；m为变量个数。怎样从已知的地理数据求出具体的判别函数呢？实质上就是求c,c,c,...,c123mc,c,c,...,c123m判别函数的系数。2）求判别函数根据费歇准则可知，要使判别函数Y能充分地反映两类地理类型的差异，就得使：两类间均值差：［y（a）-y（b）］2尽可能大;两类内的离差平方和：送［y（A）-y（a）］+艺［y（B）-y（B）］尽可能小I这样其比值I才能最大，从而将两类清楚地分开，即表达式为:[y([y(a)-y(b)]2nmaxnmaxest》[y(A)-y(A)1+为[y(B)-y(B)]2nminii=1ii=1根据极值原理，可得判别函数的正规方程组为：1111lc+lc+lc+•…+lc=d=二x-x1111221331mm11A1Blc+lc+lc+••-+lc=d=x—x2112222332mm22A2Blc+lc+lc+••-+lc=d=x--x3113223333mm33A3Blc+lc+lc'+…+lc=d=x—xm11m22m33mmmmmAmB其中，lii=l(A)+liiii(B)l(A)=E(x其中，lii=l(A)+liiii(B)l(A)=E(xiiiA-xiA)2=Ex2-(ExiAiA)2/nA(B)=E(xiiiB-x)2iB二工x2-（工x2iBx)2/niBBl=l=l(A)+l(B)ijjiijijl(A)=E(x-x)(xijiAiAjA-l(A)=E(x-x)(xijiAiAjA-xjA)=ExxiAjA-ExExiAjA/nAl(B)=E(x-x)(xijiBiBjB-xjB)=xxiBjB-ExExiBjB/nB解此正规方程组求得c,c,c,...,c123m，即得判别函数。3）求判别指标求判别指标是指求判别临界指，即：判别类型归属的数量指标。当c,c,c,...,c123m

后，即可求得：乙cx(A)kkk=1乙cx(B)kkk=1于是用加全平均法可求出判别指标（y）cy=y=cn+nAB其中，叮为A类样本个数;n：为B类样本个数。By求出后，若判别函数又通过显著性检验，就可进行类型归属c判别，具体方法是：若一-若y(A)>y(B)当y>y时，则应归属于A类；ic当y<y时，则应归属于B类。ic若_若y(A)<y(B):当y<y时，则应归属于A类；ic当y>y时，则应归属于B类。ic4）几何意义几何意义是指判别函数与判别临界值的直观几何坐标图示。141（1）作图已知A,B两类，两类特征值：x,x及原始地理数据。12以为横轴，以为纵轴作直角坐标系，并将A，B两类原始地理数据点绘于坐标图上得判别分析的几何解析图。（2）问题如何最大限度的将A，B两类清楚分类，即使重合部分最少。（3）方法从图中可知，无论将其投影到x轴x轴上，都很难将其清楚分12类。若作直线时，A，B两类分别投影到其上，则可以清楚看到,重合部分显著减少，这样就能最大限度将其清楚分类。直线就相当于所求的判别函数直线，这就是上述介绍的判别分析原理的几何意义，即最大分离原则和最小二乘法原理。若在从A，B两类最外椭圆交点处作一直线垂直于，则交点C

就是判别指标值y，即A,B两类的分界点。c5）显著性检验n+n一mn+n一m一1、AB—m：综合距离函数；nnA_B(n+n)(n+n一2)ABAB其中，D2=(n+n-2)£cdABkkk=1M：变量个数;求F(f"),a若F>F时，则显著;a若F<F时，则不显著。a6）求贡献率（乂%）kcdx%=—k（n+n一2）x100%kD2AB由其式求得的结果可知各变量在判别分析中贡献作用的大小2地理判别分析的步骤例1根据地理判别分析的原理和方法，现介绍两类判别分析的具体步骤。以表7・11和表7・12的地理数据为例。P142一143已知：A类为塔城，特征变量x为大风日数，x为沙暴日数；1A2AB类为莎车，特征变量x为大风日数，x为莎暴日数；1B2B待判类：D类为库车，特征变量x为大风日数，x为沙暴日数。1D2D第一步求判别函数1、求算判别函数的正规方程组根据其基本原理和表7.11中的地理数据可先求算有关项，即可求得其正规方程组。(1)求工项A类：工x=583.01A工x2=37185.01A工x-x=109021A2AB类：工x=160.01B工x2=2814.01B工x-x=3051.0工x=17.02A工x2=73.02A工x=179.02B工x2=3593.02B1B2B⑵求x、x项(i=l，iAiBA类：云=58.31AB类：x=16.01Bx=1.72Ax=17.92B⑶求dk项d=x—x=58.3—16.0=42.311A1Bd=x—x=1.7—17.9=—16.222A2B⑷求L、L=L项G、j=1，2）iiijjiA类：l（A）=工x2—乞x1「n=37185111A1AA.0—583.02..10=3196.1同理求得：l（A）=44.122工x.-n2AAl(A)=l(A)=工x-x-工x.-n2AA12211A2A1A=1092.0-583.0x17.0/10=100.9B类：25.04l(B)=Ex2-E2.10=25.041B1BB同理求得：l(B)=388.922l(B)=l(B)=工x-x-工x-x■■n211B2B1B2BB=3051.0—160.0x179.0/10=187.0于是可得：=l(A)+l(B)=3196.1+254.0=345101111l二l(A)+l(B)二44.1+388.9二433.0222222l=l=l(A)+l(B)=10.90+18.07=28.9712211212⑷求d项(i=1，2)id=x—x=58.3—16.0=42.31A1Bd=x—x=1.7—17.9=—16.22A2B于是可得判别函数的正规方程组为：f3450.1c+287.9c=42.3<12I287.9c+433.0c=—16.2122、求判别函数解上述正规方程组可得：fc=0.0176<1Ic=—0.05212则判别函数为：—0.0521x2y=cx+cx—0.0521x211221第二步求判别指标yJc1、求yy(B)y(A)Cx第二步求判别指标yJc1、求yy(B)y(A)CxiiA0.0176x58.3—0.052117=0.9369y(B)CxiiB0.0176x16.0—0.0521x17.9=—0.6512i=122、求yCyC(ALB10yC(ALB1020.9369+102(—0.6512)=0.1429n+nAB10+10第三步显著性检验对判别函数进行显著性检验。1、求D2,第三步显著性检验对判别函数进行显著性检验。1、求D2,即综合距离函数。D2=(Cd+Cd)(n+n-2)

1122AB=[0.0176x42.3+(—0.0521x(—16.2)]x(10+10—2)=28.58542、求F-2)((n+n一m一1、B-mxD210x10(10x10(10+10—2—1AGo+10)(10+10—2)x28.585467.4933查表可得：3、查表可得：+n一m一1=17Ba=0.01,+n一m一1=17B12A

F(fl,f2)=F(2,17)=6.11a0.01因为F=67.4933〉F「7)=6>11，所以判别函数显著。0.01第四步判别分析1、判别准则由于y()=0.9369〉y()=-0.6512，若y〉y=0.1429，则应归属A类;(A)(B)C若yVy=0.1429，则应归属B类。C2、检验即检验已知A类、B类的正确率。将原始数据分别代入判别函数，可求得y将原始数据分别代入判别函数，可求得y、yA汴的序列值，并与yC比较,即可得检验结果,如表7.12所示。由该表可知yAyB均大于y,贝U皆应判归A类，正确率达100%；yAyBC均小于y,贝皓应判归B类，正确率达100%。C例如：y,y分别将原始地理数据代入判别函数即可求得。ABy二工Cx=0.0176X88-0.0521x0=1.547921AiiA因为y=,4792〉yC=0.1429,所以属于A类。1ACTOC\o"1-5"\h\zy二工Cx=0.0176X77—0.05X0L=1.35442AiiA因为y]35443>yC=0.1429,所以也属于A类。2ACy二工Cx=0.0176X13—0.0521x16=—0.60491BiiB因为y一0.60493VyC=0.1429,所以属于B类。1BCy二工Cx=0.0176X24—0.0521x27=—0.98456BiiB因为y一0.98454<yC=0.1429,所以也属于B类。6BC

同理可得检验结果如表5.5所示。表5.5检验结果序号A类B类1122123124125126127128129121012止确率100%100%注：1表示A类；2表示B类3、判别判别是指对待判类来判别其类型归属只要将待判类yD的数据序列分别代入判别函数求得订数据序列,并与yC比较，即可判别其类型归属，如表7.12。p43y二工Cx=0.0176X23-0.0521x10=-0.116431DiiD因为y一0.H643VyC=0.1429,所以属于B类。1DCy二工Cx=0.0176X26-0.0521x12=-0.16786DiiD因为y一0.16786<yC=0.1429,所以也属于B类。6DC判断结果均归为B类，这表明B类与D类同属一类，其结果符合该地区的自然地理特点。第五步求贡献率是指求各x的贡献率。iX］的贝献率=(n+n-2)x100%D2AB0.0176x42.3()=Qo+10-2丿X100%28.5854=46.85%X2的贝献^=C2d2(n+n-2)X100%D2AB-0.0521X(-16.2)=(10+10—2)x100%28.5854=53.15%由此可知，x2的判别贡献率大于*的判别贡献率。这表明在该类判别分析中沙暴日数X。的作用大于大风日数x的作用。2x1例2根据地理判别分析的原理和方法，现介绍两类判别分析的具体步骤。以表5.4的地理数据为例。

表5.4地理数据类另I」样本号x1X2X3A类1952116.8181956139.8141958-113.8181962111.519196310.0191965216.121197019.3191972128.916均值0.917.018B类1953-538.4151954-295.4171955-129.7121957-254.9211959-333.0171960133.3161961-129.3161964-240.0181966154.4151967012.315

196804.2131969-149.9171971-115.215均值-1.237.716已知：A类为塔城，特征变量x为大风日数，x为沙暴日数；TOC\o"1-5"\h\z1A2AB类为莎车，特征变量x为大风日数，x为莎暴日数；1B2B待判类：D类为库车，特征变量％为大风日数，x为沙暴日数。1D2D第一步求判别函数1、求算判别函数的正规方程组根据其基本原理和表6.3中的地理数据可先求算有关项，即可求得其正规方程组。(1)求工项A类：工x=71A工x2=111A工x=13622A工x2=336.92A工x-x=124.7工x-x=1291A2A1A3AB类：工x=-161B工x2=521B工x•x1B2B工x=490.02B工x2=24977.52B-708工x-x=-2671B3B工x=1443A工x2=26243A工x-x=2303.72A3A工x=2073B工x2=33573B工x-x=812.32B3B⑵求x、x项(i=1，2，3)iAiBA类：x=0.9x=17x=181A2A3AB类：x=-1.2x=37.4x=161B2B3B⑶求L、iiL=Lij项(i、j=1,ji2，3)A类：L(A)=L(A)=Ex2111A-(Ex1Al/n1A=11-72.；8=4.8750同理求得：L(A)=1051.095022

L33(A)=32・0000L(A)=L(A)=Ex12211A・x2Ax・Ex-n=12.74-7x13.2681A2A•A'=5.5250同理求得：L13(A)=L13(A)=L31(A)=3・°000L23(A)=L32(A)=-147.9000B类：L(B)L(B)=Ex2-111Bx1B=52-(-16)2.1332.3077同理求得：L(B)=6508.269222L33L(B)=6508.269222L33(B)=60.9231L(B)=L(B)=E1221x1B・x2B-Ex-xn1B2B•B—708-(—16)x490.13=－104.9231同理求得：L(B)=L(B)=318.99232332L(B)=L(B)=318.992323321331于是可求得：L=L(A)+L(B)=4.8750+32.3077=37.1827111111ii=1ii=1同理可得：L22=7.53964L33L22=7.53964L33=92-9231L=L1221L(A)+12L12(B)=5.525—104.9231=—99.3981同理可得：L13=L31—9.2708L13=L31—9.2708L23=L32=171-0923(4)求d项(i=1.di2，3)0.9—(—1.2)=2.11B同理可得：d=—20.72于是可得判别函数的正规方程组为：|37.1827C—99.3981C—9.2708C=2.1TOC\o"1-5"\h\z123<—99.3981C+7559.3642C+171.0923C=—20.7123—9.2708C+171.0923C+92.9231C=21232、求判别函数解上述正规方程组可得：「C=0.05761<C=—0.00272C=0.03313则判别函数为：y=0.0576x—0.0027x+0.0331x123第二步求判别指标yC♦求入、♦求入、(B)_3_y(A)=ZCx=0.0576x0.9—0.0027x17+0.0331x18iiA=0.6017y(B)=fCx=0.0576x(—1.2)—0.0027x37.7+0.0331x16y(B)=i=1=0.358722、求ycy()+ny()8X.6017+13X0.35878+13★3B3==8+13n+nAB第三步显著性检验对判别函数进行显著性检验。1、求D2,即综合距离函数。D2=(Cd+Cd+Cd)(n+n—2)112233AB=(0.0576x2.1+(—0.0027)x(—20.7)+0.0331x2)6+13—2)=4.61802、-2)((n+n一m一1、AB-mXD28X13(8+13)(8+13—2)、x4.6180丿3、=6.8222查表可得：F(f1，f2)=F(3,17)=5.18a0.01因为F=6.8222〉f(3,17)=5.18，0.01所以判别函数显著。第四步判别分析1、判别准则由于y=0.6017>y=0.3587，右y>y=0.4513，则应归属A类；（A）（B）C若yVy=0.4513，则应归属B类。C2、检验即检验已知A类、B类的正确率。将原始数据分别代入判别函数，可求得y、y的序列值，并与yayAByC比较，即可得检验结果。例如：y=ECx=0.0576X1-0.0027x16.8+0.0331x181AiiA=0.6080因为y=0.6080>yc=0.4513,所以属于A类。1Ay=ECx=0.0576X1-0.0027x39.8+0.0331x142AiiA=0.4135因为y=0.4135VyC=0.4513,所以属于B类。2Ay=ECx=0.0576X（-5）—0.0027x38.4+0.0331x15TOC\o"1-5"\h\z1BiiB=0.1048因为y=0.1048VyC=0.4513,所以属于B类。1BCy=ECx=0.0576x1-0.0027x33.3+0.0331x166BiiB=0.4973因为y=0.4973>yC=0.4513,所以属于A类。6BC同理可得检验结果如表5.5所示。

表5.5检验结果序号A类B类11222231241251261171281292101112122132正确率87%84.6%注：1表示A类；2表示B类3、判别判别是指对待判类来判别其类型归属只要将待判类yD的数据序列分别代入判别函数求得y数据序DD列，并与yC比较，即可判别其类型归属。第五步求贡献率是指求各x的贡献率。ix的贡献率=CA(n1D2A+n—2)x100%B0-0576x2-1(8+13—2)x1°°%4.6180=49.77x2的贡献率"d+n—2)x100%B=—°.°027X(—2°・7)(8+13—2)x100%4.6180=23.00的贡献率=。dx^^(n+n—2)x100%3D2AB=0.0331x2(8+13—2)x100%4.6180=27.23由此可知，X］的判别贡献率大于x3的判别贡献率大于x2的判别贡献率。2.3地理判别分析的实习指导2.3.1实习目的1、巩固地理判别分析的基本原理及方法步骤。2、掌握地理判别分析程序的使用方法与技巧。3、求取判别分析结果并应用于与实际的对比分析4、掌握地理判别分析程序的变换应用方法。2.3.2实习内容1、标识符说明N1第一类样本数第二类样本数N2N3待判类样本数P变量数A1(P，N1)存放第一类原始数据A2(P，N2)存放第二类原始数据A3(P，N3)存放待判类原始数据S(P，P+1)存放正规方程组系数B(P)存放判别函数系数Y0判别临界值FF检验统计量D・F1F统计量的第一自由度D・F2F统计量的第二自由度2、程序10REM两判别分析程序20INPUT“变量数P,—类，二类，待判类样本数Nl,N2,N3=”；P,Nl,N2,N3100PRINTP,N1,N2,N3110'QOIsTheMaximumofNl,N2,N3l20Q0=Nl130IFQ0>N2GOTO170140Q0=N2150IFQ0>N3GOTO170160Q0=N3170DIMA1(P,N1),A2(P,N2),A3(P,N3),C(P),B1(N1+1),B2(N2+1),B3(N3),S1(P,P),S2(P,P),X(P),M1(P),M2(P),A(P,Q0),S(P,P+1),B(P)180REM输入原始数据200FORI=1TON1210FORJ=1TOP220READA1(J,I)230PRINTA1(J,I)；240NEXTJ250PRINT260NEXTI280FORI=1TON2290FORJ=1TOP300READA2(J,I)310PRINTA2(J,I)；320NEXTJ330PRINT

340350360380390400410420430440450460470480490500510520540550560570590600610620630640650660680700710730740750760NEXTIIFN3=0GOTO450REMFORI=1TON3FORJ=1TOPREADA3(J,I)PRINTA3(J,I)；NEXTJPRINTNEXTIN=N1FORI=1TON1FORJ=1TOPA(J,I)=A1(J,I)NEXTJNEXTIREM计算正规方程组系数GOSUB1840FORI=1TOPM1(I)=X(I)FORJ=1TOPS1(I,J)=S(I,J)NEXTJNEXTIN=N2FORI=1TON2FORJ=1TOPA(J,I)=A2(J,I)NEXTJNEXTIGOSUB1840FORI=1TOPM2(I)=X(I)FORJ=1TOPS2(I,J)=S(I,J)NEXTJNEXTI

780REM790FORI=1TOP800B(I)=(M1(I)－M2(I))*(N1+N2－2):S(I,P+1)=B(I)820FORJ=1TOP830S(I,J)=S1(I,J)+S2(I,J)835PRINTS(I,J);850NEXTJ855PRINTS(I,P+1)，870NEXTI890REM求解正规方程组900GOSUB2080910REM判别系数920FORI=1TOP930C(I)=S(I,P+1)940NEXTI950Y1=0960Y2=0970FORI=1TOP980Y1=Y1+C(I)*M1(I)990Y2=Y2+C(I)*M2(I)995NEXTI1000Y0=(N1*Y1+N2*Y2)/(N1+N2)1010REMF统计量1020D=01030FORI=1TOP1040D=D+ABS(C(I)*(M1(I)－M2(I)))1050NEXTI1060D1=P1070D2=B1+N2－P－11080REM1100F=N1*N2*D2*D/(P*(N1+N2))1110FORI=1TON1+11120B1(I)=01130NEXTI1140FORI=1TON2+11150B2(I)=01160NEXTI

1170118011901200121012201230124012501260127012801290130013101320133013401350136013701380139014001410142014301440145014601470148014901500151015201530IFN3=0GOTO1220FORI=1TON3B3(I)=0NEXTIREM原样本判别FORJ=1TON1FORI=1TOPB1(J)=B1(J)+C(I)*A1(I,J)NEXTIIF(B1(J)－Y0)*(Y1－Y2)>0GOTO1290B1(J)=2GOTO1310B1(J)=1B1(N1+1)=B1(N1+1)+1NEXTJFORJ=1TON2FORI=1TOPB2(J)=B2(J)+C(I)*A2(I,J)NEXTIIF(B2(J)－Y0)*(Y1－Y2)>0GOTO1400B2(J)=2B2(N2+1)=B2(N2+1)+1GOTO14