球面最短距离的求法公式及其推导

2005高教社杯全国大学生数学建模竞赛（C题）雨量预报方法的评价模型学校名称：江西财经大学参赛队员：肖潇杨靖雅孺指导老师：数学建模小组2005年9月19日

雨量预报方法的评价模型肖潇杨靖雅孺摘要：本文使用方差分析模型、相关性检验模型和创建效用模型分析了气象部门两种雨量预测方法的准确性，利用!ATLAB软件强大的计算功能，对问题给出的大量数据进行分析，并得到了一致的结果。首先，建立了坐标转换公式，将问题中给出的网格点和观测站点的经纬度转化为平面直角坐标，解决了网格点不等距问题。并利用已知的网格点上雨量的预测数据通过二维Lagrange插值多项式，求出了91个观测站点处雨量的预测值；其次，使用方差分析的方法，2191(f(k)Z(k))2，通过比较两种预测方法下所有预nijiji1测数据与实际数据的总方差不、同观测站点处预测数据与实际数据的方差以及不同时段下预测数据与实际数据的方差，均得到预测方法一较预测方法二的准确性更高的一致结论并，对观测站点的位置进行分析发，现长江三角洲地区和沿海地区的观测站点更适用第二种预测方法；同时发现了五个数据异常的观测站点，并对其成因进行了分析；第三，使用二分类方法，对两种预测方法的雨量预测数据和实际数据之间分别作了相关系数检验(—检验、独立性检验2—检验和回归系数检验Q—检验，均同样得到了与方差分析方法一致的结果：预测方法一较预测方法二的准确性更高；U(k)

ij第四，对于问题2，U(k)

ij(k)f(k)Z(k)，U(k)eLi?，Uijijj'j '分别针对两种预测方法，求出了所有预测等级的总效用、各观测站点上的预测等级的效用以及不同的四个时段上预测等级的效用通，过对两种预测方法效用的比较，依然得出了预测方法一比预测方法二效用更高、准确性更好的结论；最后，作为文章的结束，阐述了本文模型的优点及缺陷同，时也指出了需要改进的地方。关键字：雨量预报、预测方法.Lagrange插值、误差分析、相关性分析、效用分析

一、问题的重述雨量预报对农业生产和城市工作和生活有重要作用，但准确、及时地对雨量作出预报是一个十分困难的问题，广受世界各国关注。我国某地气象台和气象研究所正在研究6小时雨量预报方法，即每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度、北纬32度附近的53x47的等距网格点上。同时设立91个观测站点实测这些时段的实际雨量，由于各种条件的限制，站点的设置是不均匀的。气象部门希望建立一种科学评价预报方法好坏的数学模型与方法。气象部门提供了41天的用两种不同方法的预报数据和相应的实测数据。小于0.1毫米视为无雨。（1） 请建立数学模型来评价两种6小时雨量预报方法的准确性；（2） 气象部门将6小时降雨量分为6等：0.1—2.5毫米为小雨，2.6—6毫米为中雨，6.1—12毫米为大雨，12.1—25毫米为暴雨，25.1—60毫米为大暴雨，大于60.1毫米为特大暴雨。若按此分级向公众预报，如何在评价方法中考虑公众的感受？二、模型的分析与假设根据预测点（网格点）和实测点的经纬度，画出预测点和实测点的位置图（图1），但由于经纬度坐标并不对应于直角坐标系，得到的并不是问题中给出的规则网格图。同时发35纬

度343332313035纬

度34333231302928:+**++47X53个预测点与生个实测点经纬度位置图+++++++++++♦+t**********二117 118 119 120 121 122 123 *124经度125图（一）现91个观测点不但分布不均匀,而且基本上没有落在预测点上，因此，我们考虑将预测点和实测点的经纬度坐标转化为平面直角坐标。坐标转换后91个实测点并没有在网格点上，因此我们将使用二维插值多项式计算出实测点处两种预测方法预报的降雨量的近似值，然后将其与实测点处降雨量的实际值进行比较，给出评价模型和问题的解答。鉴于本问题中给出的数据十分庞大，（不完全统计，数据量达80万个），所以在分析过程中，本文大量使用了MATLAB软件中提取文件中数据的dlmread函数。根据本文的需要，我们提出以下合理的假设：1．雨量用毫米做单位，小于0.1毫米视为无雨；2.地球为规则球体，地球半径取WGS84标准中的长半轴6378.137公里，短半轴6356.752公里的平均值6371公里；3.预测点位置位于东经120度、北纬32度附近的53x47的等距网格点上；不考虑人为因素所造成的测量值的误差；题目中所给经纬度坐标值和雨量数据为真实值，或者与真实情况相符合。三、符号说明如无特别说明，本文中的符号特指以下意义：（,）：经纬度坐标；ij（x,y）：修正后欧氏坐标；ij：球心角；：经度差；S：A,B两点的球面距离；ABz：降雨量；ijzk:表示第i个观测点第j天第k时段雨量的实际值；ijfk:表示第i个观测点第j天第k时段雨量的预测值;ij2（k）:表示用雨量的预测值与实际值的方差；U:表示效用值；四、模型建立与求解（一） 球面最短距离的求法公式及其推导设A点经纬度为（，）,B点经纬度为（，），R为地球半径，则A,B两点之AA BB间的球面距离S为ABSRarccos（sinsincoscoscos（））AB AB AB AB具体的球面距离公式的推导见附录［1］。（二） 预测点和观测点处经纬度转化为直角坐标我们以最左下的预测点为原点，通过球面距离公式，将预测点和观测点转化为直角坐标系下的点，具体位置如图（二）。从图中可以看出，整个预测区域已经化为一个近

似的矩形区域（关于将预测点和观测点的经纬度转换为直角坐标的具体过程见附录［2］,但是与问题中给出的预测区域是位于东经120度、北纬32度附近的53X47的等距网格点不相符,因此,我们再次调整预测点位置,将其位置微调到等距网格点上,具体位置如图（三）。在后面的分析中,我们将使用分为非等距网格点和等距网格点两种数据进行比较分析。onunuonunu876公里转换后的观测点与预测点的位置图调整在等距网格下观测点和预测点位置图50403020jrl里700600500400300200100:::::::::::::::::::::::::::::::::::♦+onunuonunu876公里转换后的观测点与预测点的位置图调整在等距网格下观测点和预测点位置图50403020jrl里700600500400300200100:::::::::::::::::::::::::::::::::::♦+♦+♦+♦+♦♦♦+♦+♦♦♦■♦■♦♦♦+++♦♦♦+♦♦♦+++♦::::::::::::::::::::::::::::::::::::：：：：：：：：：：：：：：：：：：：：：：：：：：：：：：：：:計………t:：：：：：：：：：：U：：：：：：：：：：：：：：：：：：：：：：：：：：：：：：：：：：rfetr::u::r::::::::::::::::::::::::::+J+*tttJt++*tttJ**t++*t+♦♦***++♦*：：：：：：：：：：：：：t：：u：m::♦+♦*******★+♦ *竝***★+++** ******************：：：：：:;::站::：：;;：：：：«：:；：：：：：：：：：：：：：：：：：：：：：：： :u::::::::::::::::::::::*******■：****************fc.++ *********************+r-*****♦<；■+***-♦**轉农*+<■_*****************莪注淞矇鋳密詩録"由舲由由宙辽************■******¥******+♦★★*************出技陰:也请昭iiWlH◎密也出由由::::::::::::::::::::::::::::::::::::::::::::::占"占"''测测观预10 400 500 600公里旳°::::::::::100200300::::::400500600公里700图（三）图（二）图（三）（三）观测点处降雨量预测值的计算由于各种条件的限制，91个观测站点的设置不是均匀的，并且也不在预测位置的网格点上，为了与观测点处的降雨量的实际观测值作比较，必须根据网格点上降雨量的预测值计算出观测点处的预测值。我们由于网格点的步长较短，大概在15公里以内，可以使用二维二次插值多项式来预测观测点处的降雨量。假设第k个观测点的位置在P（x（k）,y（k））处，其相邻的四个预测点（网格点）的位置为（x,i测点（网格点）的位置为（x,iy),(x,yji(x,y)i j1),(x,yj1 i1),(x,yj i1(x」y丿P(x(k),y(k：丄1 j1)(x,y)ij(x,y)i1jj1），如图（四）图（四）如果四个预测点的降雨量的预测值为f（x,y）,f（x,y）,f（x,y）,f（x,y），使用二ijij1i1ji1j1维Lagrange插值法［1,572~573页］,可以计算出观测点P（x（k）,y（k））处降雨量的预测值f(x(k),y(k))：

f(k) f(x(k),yf(k) f(x(k),y(k))(xx)(yy)厂十7七f^y.)(xx)(yy)iji1 . .1(xx)(yy)( 〒一 f(x]y.)(xx)(yy)i1,.i1i. .1(xi(xLi1(xx)(yy)x)(yy)ij1i1j1 jx.)yy.)f()—jf(x,y)x)(yij1k1,2,)(y.1)(yy) j-y)i1.1.,91由于问题给出的预测点位置并非等距网格数据（图（二）），为合问题本意，我们修正了预测点位置，并使其为等距网格数据（图（三）），因此根据两种预报方法和两种网格数据，可以得到91个观测点上的四种不同的雨量预测值，并分别与观测点上雨量实际值进行比较，进而综合评价预报方法的准确性。图（四）有雨时两种预测值与实际值的比较。图（四）图（四）（四）预测方法的评价根据上面的分析，可以根据两种雨量预报方法给出的预测点（网格点）上雨量的数据，求出91个观测点41天每天各4个时段的雨量的预测值，同时使用MATLAB软件的dlmread函数，从问题附件MEASURING文件夹中提取91个观测点41天每天各4个时段的实际降雨量，比较两种雨量预报方法的准确性。1．模型一：（方差分类分析法）假定根据某种方法预测的第i个观测点第j天k时段雨量为f（k），而实际雨量为i.Z（k）,使用方差来衡量方法的准确性：i.（1）预测方法的总方差： 2 1 44191（f（k）Z（k））244191 i.i.k1.1i1由非等距网格数据（图（二））和修正的等距网格数据（见图（三））以及两种不同的预报方法得到的雨量预测值与实际值的总误差见表（一）。

表（一）两种数据两种方法预测的总误差非等距网格数据等距网格数据方差标准差方差标准差方法16.22834.028416.00814.0010方法—17.27314.156117.08744.1337从表（一），不管利用非等距网格数据还是等距网格数据，预报方法一产生的方差都要小于方法二，即方法一准确性要高于方法二。2）每个观测点处的方差：2 (f(k)Z(k))2,i1,2, ,91i441 ij ijk1j1由非等距网格数据（图（二））和修正的等距网格数据（见图（三））以及两种不同的预报方法得到的91个观测点的4组雨量预测值与91个观测点雨量实际值的误差与标准差见表（二）。表（二）91个观测站点雨量预测值与雨量实际值的方差及标准差非等距网格等距网格预测方法一预测方法二预测方法一预测方法二序号观测站号误差标准差误差标准差比较观测站号误差标准差误差标准差比较1581384.4662.11334.69452.16671581384.02912.00734.78022.186412581391.18791.08991.08251.04040581390.85340.92380.77910.882703581410.62010.78750.64790.80491581411.08141.03991.16421.07914581431.00791.0040.81940.90520581430.68520.82780.51060.714605581462.14971.46622.31811.52251581462.29381.51452.46681.570616581473.81931.95433.32851.82440581473.7031.92433.77171.942117581480.4430.66550.24070.49060581480.75560.86920.38890.623608581500.90970.95380.99180.99591581502.29551.51512.15821.469109581542.67711.63621.97211.40430581542.50321.58222.03071.425010581588.66182.94318.26032.87410581587.97712.82448.23282.869311158230149.95812.2457161.87712.723115823081.60499.033570.43288.3924012582360.89370.94540.96870.98421582360.37870.61540.4430.66561135823812.79263.576712.96553.60081582382.52691.58962.59381.6105114582400.72750.85290.62130.78830582400.97230.9860.97140.9856015582410.59680.77250.96160.98061582410.56170.74940.88980.9433116582430.43290.6580.41510.64430582430.75260.86750.85810.9263117582450.96970.98480.45280.67290582451.13341.06460.82360.9075018582460.69670.83470.79420.89121582460.71240.84410.77570.8807119582490.90340.95050.81760.90420582490.87940.93780.75290.8677020582510.40670.63770.27510.52450582510.77140.87830.6120.7823021582524.40072.09784.02372.00590582523.87921.96963.50431.872022582540.86130.92810.87380.93481582541.42251.19271.38931.1787023582550.70960.84241.09291.04541582551.02491.01240.75920.8713024582640.17580.41930.48030.69311582640.78110.88380.77890.8826025582651.63591.2791.76461.32841582651.64381.28211.67221.2931126582692.73951.65512.55211.59750582692.61371.61672.27711.5090275833315.32263.914414.96263.868205833326.31655.1324.46684.9464028583344.86162.20495.58432.36311583340.9730.98640.78690.8871029583355.20782.28215.58222.36271583356.01252.4525.97052.4435030583363.99781.99954.20542.05071583364.39462.09634.22482.0554031583371.09321.04561.16121.07761583371.43991.21.5221.2337132583416.29992.517.65752.76721583416.18632.48727.46412.7321335834213.9423.733912.26773.502505834218.1244.257217.35184.1655034583431.30991.14451.383919980.9991.33611.15591355834410.32163.212711.06133.32591583449.73593.120211.18233.344136583451.02341.01160.92390.96120583450.4010.63330.5930.77011375834613.68833.699812.46293.530305834611.21253.348510.36953.2202038583491.41921.19131.45341.20561583493.98981.99743.92071.9801039583517.29442.70087.33572.70841583517.73392.7817.31152.704040583521.00411.0020.8630.9290583520.61080.78150.60520.7779041583545.02622.24197.16732.67721583541.52341.23422.60371.6136142583561.96351.40122.04891.43141583561.75291.3241.94491.3946143583581.19941.09521.69961.30371583583.60711.89925.22582.286144583593.29871.81624.5322.12891583592.62061.61882.9551.719145583600.5080.71270.47020.68570583600.58030.76180.46540.6822046583612.59441.61072.82081.67951583613.70431.92473.86421.9658147583623.11471.76483.24431.80121583622.90111.70332.98911.7289148583650578807608069808355158365059340770308551092471

49583661.98831.41011.75221.32370583664.13442.03333.87351.9681050583670.9660.98280.69370.83290583673.37331.83662.96931.7232051583691.26951.12671.23381.11070583691.21611.10281.11481.0559052583703.41111.84693.66341.9141583702.2541.50132.17461.4747053583771.3411.1581.05371.02650583770.44070.66390.26440.5142054584266.21672.49335.70712.38890584268.68262.94669.5093.0837155584312.77141.66473.20631.79061584312.04411.42972.98291.7271156584323.72551.93023.96351.99091584322.85891.69082.19841.4827057584331.02011.011.2241.10631584331.46181.20912.03761.4274158584356.66962.58268.5272.92011584355.71942.39156.89852.6265159584360.36270.60220.66640.81631584360.12160.34870.33230.5765160584389.71913.117612.65223.5571584389.0483.0089.38483.0635161584414.49032.1196.03652.45691584414.37572.09186.23532.497116258442120.72710.9876147.16412.131115844282.92979.1066113.70510.6633163584432.96971.72332.73861.65490584433.16271.77842.98131.7266064584461.39261.18011.46891.2121584461.781.33421.95751.3991165584484.95352.22575.02712.24211584485.59492.36535.89712.42841665844913.88373.726116.98054.12071584497.36582.7149.07343.0122167584501.59891.26452.28661.51211584500.34610.58830.97670.98831685845154.02377.350152.27317.230584513.87811.96933.57021.889506958452653.34125.5605644.24125.3819058452750.127.3879778.02227.893170584532.35021.5332.43351.561584531.16221.0780.97780.98880715845410.23283.198913.40793.661715845410.7533.279215.20343.8992172584552.87111.69443.36481.83431584552.71491.64773.28781.813217358456102.80710.1394144.54812.0228158456100.53910.0269144.60312.0251174584571.57971.25691.59281.26211584573.06761.75153.04561.7452075584598.10352.84678.97312.99551584598.09712.84556.17452.4849076584600.51070.71470.45320.67320584600.37510.61240.47510.6893177584613.99931.99983.87831.96930584614.04932.01234.13172.0327178584620.74220.86150.62310.78940584620.93810.96860.75430.8685079584632.80381.67453.49071.86831584633.21431.79283.35091.8305180584641.58981.26091.78411.33571584642.46771.57092.39781.5485081584672.5451.59534.29552.07261584672.81991.67924.51582.1251182584681.42331.1931.3031381.17471.16471.0792083584720.91480.95650.90060.9490584720.41170.64170.3710.6091084584774.5372.134.50912.12350584770.60990.7810.4890.6993085584844.02982.00743.18681.78520584845.12362.26354.10522.02610865853084.48159.191483.20689.1218058530141.08611.878134.23811.58610875853126.8885.185432.09275.6651585313.18371.78436.50112.5497188585344.3232.07925.47932.340815853417.55454.189815.32623.9149089585420.37480.61220.80980.89991585420.52010.72120.75740.8703190585501.14461.06980.95360.97650585501.36661.1691.39881.182719158562176821329731309176945856211686108116989130341从表（二）可以得到，根据非等距网格数据，91个观测站点中有55个点方法一的方差更小，而根据等距网格数据，91个观测站点中有47个点方法一的方差更小，见表（三）综合数据可以，方法一的准确性高于方法二。表（三）91个、观测站点中方差较小的$站点数非等距网格数据等距网格数据方差小的观测站点数方差小的观测站点数预测方法5547预测方法二3644假设经纬度数据符合实际情况，根据经纬度作出对应的地图。地理位置是沿海地区，分析观测点的地理位置得到以下结论：1、对于91个观测点来说，方法一具有更普遍的适用性。2、如果考虑地理位置因素，两种方法的好坏又具有相对性，对于长江三角洲和沿海地区的观测点，预报方法二更优于预报方法一，而内陆地区预报方法一更优。3、在长江三角洲和沿海地区，预报方法一不如方法二，可能是由于预报方法一预测时受水环境的影响很大，从而导致测量值误差较大。如图（五）。另外，从表（二）还发现，有五个观测站点的方差异常，与其他站点预测数据的方差相比，两种方法的结果明显偏大，产生该结果的原因是多方面的，如这五个观测站点的地理位置，工作人员的责任心等因素都可能产生这种结果。具体数据见表（四）

表（四）：五个方差异常的观测站点非等距网格等距网格预测方法一预测方法二预测方法一预测方法二序号观测站号误差标准差误差标准差比较观测站号误差标准差误差标准差比较1158230149.95812.2457161.87712.723115823081.60499.033570.4328&392406258442120.72710.9876147.16412.131115844282.92979.1066113.70510.663316958452653.34125.5605644.24125.3819058452750.127.3879778.02227.89317358456102.80710.1394144.54812.0228158456100.53910.0269144.60312.02511865853084.48159.191483.20689.1218058530141.08611.878134.23811.5861014191（3）每个时段的方差： （k）2 （f（k）z（k））2,k1,2,3,44191ijijj1i1由非等距网格数据（图（二））和修正的等距网格数据（图（三））以及两种不同的预报方法得到的四个时段观测点雨量实际值的误差与标准差见表（五）。表（五）：四个时段观测点雨量实际值的误差与标准差非等距网格等距网格预测方法一预测方法二预测方法一预测方法二时段误差标准差误差标准差时段误差标准差误差标准差时段129.46905.428530.54715.5269时段131.37095.601035.55805.9631时段222.36674.729323.95464.8943时段220.80374.561119.27494.3903时段39.51793.085110.83123.2911时段38.04052.83569.54523.0895时段43.55961.88673.75961.9390时段43.81731.95383.97141.9928从表（五）可以看出，四个时段预报方法一的方差均比方法二的方差要小，因此达到预报方法一的准确性高于预报方法二。根据前面的分析，利用方差分析的方法，针对各类数据进行方差分析，得到的结果都是预报方法一的准确性要高于预报方法二。（五）分级预报方法的评价1．模型二（利用相关系数与回归系数评价分类预报方法）如果将雨量预报划为二分类预报，把原降雨量划为两个等级赋值，将无雨、小雨、中雨、大雨雨量较小情况列为第一级别，级别赋值为0；把暴雨、大暴雨、特大暴雨雨量偏大情况列为第二类，级别赋值为1。设N为总预报次数，n.为N次预报中，实际ij值为第i类，但预报为第j类的预报次数。如表（六）表（六） 二分类列联表预报实况、第类/v/v——第—类合计第类n11n12n nn11112/v/v——第—类nnn n n212222122合计n nn11121n nn21222N根据[2424~425页]建立预测结果与实际情况的相关系数r检验：nnnnr ]]22 1221（nn）6n）6n）6n）%1121122211121222r的变化范围在-1到+1之间，r的越大越好。r1时，预报结果与实际情况完全正相关，即所有预报结果都是正确的；r1时，预报结果与实际情况完全负相关，即所有预报结果都是相反的；r0时，预报结果与实际情况完全无关。可以建立预测值与实际值独立性的2检验N(nnnn)22 11 22 12 21Nr2nnnn1212r0时，2 0,这种情况表示实际情况与预报结果之间完全独立。如果N次预报性质上全部报对或全部错误，2N。奥布霍夫提出用下列指标评定预报质量，即回归系数检验：Q1（）nn其中 -^2 , —21nn

12一般，Q值在-1~1之间，Q值越大，预报结果与实际情况拟合越好，反之，越差。

以上三种检验可以单独使用，也可以结合使用，以加强检验效果，针对雨量预报我们得到的结果见表（七）。表（七）二分类检验表预报方法预报方法—n14241,n541112n 99,n 5302122N14924n14241,n541112n103,n 5262122N14924相关系数检验r0.8692r0.865542检验2 112752 11178回归系数检验0.0038, 0.1574Q0.83880.0038, 0.1638Q0.8325表（七）中的结果表明，预报方法一的各种检验结果都要优于预报方法二，这再一次验证了在模型一中方差检验的结果。2．模型三（在分级预报中考虑公众的感受的效用模型）气象部门常采用分级预报方式进行雨量的预报，表（八）中给出了分级预报的分级设置，一般，一个预报方法是否准确，只要预报的雨量等级与实际雨量等级相同即可认为预报准确，因此我们可以将预测的各观测站点的雨量和实际雨量转化为雨量等级再进行不同方法准确性的比较。表（八）：雨量分级预报等级设置（单位：mm）z<0.10.1〈z〈2.52.5〈z〈66<z<1212<z<2525<z<60z>60尢雨小雨中雨大雨暴雨大暴雨特大暴雨1234567仍然假设z（k）为第i站点第j天k时段的实际雨量等级，f（k）为第i站点第j天k时ij ij段的预测雨量级别，i1,2,,91;（k）f（k）Z（k），j1,2,,41;ij ij ijk1,2,3,4为第i站点第j天k时段的预测等级的级差。公众关注天气预报是要满足想知道未来准确天气状况的欲望，从而更好的安排生产生活。效用即对欲望的满足，具有很强的主观性。于是我们在模型中引入了效用的概念。当观测站的预测值与实际值同在一个阶段时，我们认为预报是准确的，公众欲望得到满足，效用为1。当观测站的预测值与实际值不在同一阶段时，公众从天气预报中获得的效用将减少，当实际值与预测值相差等级越大时，公众会感到预报不准确程度更高，所获得的效用就越小，这正符合边际效用递减规律。定义

i1,2,,91;U(k)eL(ikj)，j1,2,,41;ijk1,2,3,4为该次预报给公众的感受，成为公众的效用，其中L表示公众对未来天气的关注程度，0L1表示公众不太关注未来天气，L1表示公众关注未来天气。也就是说，当L取不同的值时针对的是对天气变化敏感度不同的人群，L越大，效用越小。例如城镇居民与农民相比受天气变化影响较小。U（＜）是（k）的单调递减下凸函数，效用函数如图（六）。ijij0 0 1 F———_I 1 ' o1 2 3 4 5级差6同样可以定义：图（六）总效用：U44191U44191(k) eLij(k)ijk1j1i11Jk1j1i1每个观测站点的效用：Ui441U(k)ij441eL(k)，iji1,2,,91Jk1j1k1j1每个时段的效用：U4191U(k)ij4191eL(ikj)，ijk1,2,3,4Jj1i1j1i1分别取L=0.5和L=2,针对不同的预报方法，我们计算出各种效用，如表（九）。表（九）效用表L=0.5L=2预测方法一预测方法二预测方法一预测方法二总效用13679136691225312231时段一3440.33443.33102.33108.5时段二3444.33444.93122.73124.3时段三3422.53417.83070.53061.4时段四3371.93362.92957.52937

针对91个观147.6291147.2356128.4317127.5671测站的效用150.2286149.8351133.7367132.8721140.9502140.5567114.4642113.5996142.6688143.0623117.9387118.8034146.2939146.0552125.0901124.9731146.997147.3905127.45128.3147155.3437154.9502144.9774144.1127152.1959152.5894138.0601138.9247150149.7513134.3515133.5027143.1562142.7627118.9046118.0399158.4076158.0141152.5253151.6607153.3763153.3763140.6541140.6541137.4227137.5775115.2149115.9625149.203148.416131.8904130.1611149.8351148.6547132.8721130.2781149.4416149.8351132.0074132.8721154.9502154.9502144.1127144.1127146.997147.7839127.45129.1794147.6291147.2356128.4317127.5671152.1959151.8025138.0601137.1954147.8549147.8549130.5576130.5576141.9529141.3207117.5877116.606152.5894152.9829138.9247139.7894157.7045158.098150.1654151.03136.2968135.2712108.2072106.3609152.7442152.3507139.6724138.8077146.8422146.0552126.7024124.9731149.7613149.7613134.2345134.2345148.8805148.8805132.404132.404149.6166148.8297134.2187132.4893156.2854155.105147.4544144.8604148.0226146.8422129.2964126.7024153.2925152.5056141.2847139.5553139.6759139.9145111.7691111.8861146.21146.21125.7207125.7207153.1377153.9246140.537142.2664141.0111141.4046115.1107115.9754145.9004145.5069124.2254123.3608140.7853140.7853112.9848112.9848143.9331144.3265119.9021120.7668155.3437154.9502144.9774144.1127148.7257148.7257131.6563131.6563142.9075142.1205118.0558116.3264154.3282155.902143.8628147.3215152.9829152.5894139.7894138.9247155.7472155.3537146.5738145.7092151.409150.2286136.3307133.7367154.1633153.7698142.3834141.5187154.9502154.5567144.1127143.248159.6718159.6718154.4887154.4887152.5894151.8025138.9247137.1954154.8725153.6921144.8424142.2484152.5894152.9829138.9247139.7894146.0552148.0226124.9731129.2964149.4416149.0482132.0074131.1427149.203148.5709131.8904130.9087152.9829152.5894139.7894138.9247149.1192149.7513132.521133.5027155.3437153.7698144.9774141.5187149.5965149.9899132.7551133.6197148.6547147.8678130.2781128.5487147.4743148.2612127.6841129.4134143.4558143.8492119.6681120.5327153.9956153.8408143.6446142.897148.1774148.8095130.044131.0257146.997148.8095127.45131.0257142.9075144.0879118.0558120.6497156.9377157.7246149.8997151.629153.7698154.1633141.5187142.3834147.9387149.5126129.927133.3857151.1703151.1703136.2137136.2137158.2628156.5342152.5095148.3032150.2319150.2319136.5765136.5765153.9246153.5311142.2664141.4017148.1774146.997130.044127.45157.4658157.4658150.0483150.0483154.3181154.4729143.131143.8787157.7045157.7045150.1654150.1654155.7371154.9502145.842144.1127156595115620161489496148085149.9061149.9061134.2503134.2503154.7116154.5567143.9957143.248152.1959152.9829138.0601139.7894152.1959151.1703138.0601136.2137147.4743147.0808127.6841126.8194149.4416149.0482132.0074131.1427155.4208155.5756146.4547147.2023155.5186155.5186147.1886147.1886157.311156.9176149.3007148.436150.1447151.3251134.3674136.9614151.409150.7769136.3307135.349根据表（九），在效用模型中，因为天气预报的三、四时间段（9:00—21:00）正处于人们正常的出行，工作期间，所以公众对三、四时间段的关心程度远高于一、二时间段（21:00—9:00）。在一二时间段两种预测方法效用接近的前提下，应以三、四时间段两种方法的效用比较来评价预测方法的优越性更合理。对L在区间（0,1],1,）上分别取L 0.5和L2来计算总效用，各时段效用，各观测站效用。在L两种取值下预测方法一的总效用大于预测方法二大，且两种预测方法在一、二时间段效用基本相同的前题下，方法一在三、四时间段的效用值大于方法二。这说明预测方法一能提供给公众更高的满意度，既预测方法一在公众眼中更准确。五、模型的评价与改进（一）模型优点：1、 本文通过球面距离公式将经纬度坐标转换为直角坐标（欧式坐标），得到了相对规则的矩形区域，在此基础上作的插值结果更为精确。2、 模型细致的分析了91个观测点，比较每个观测点上两种方法的准确度，在矩形区域经纬度所在范围内的地图上把每个观测点的不同方法区别标记，分析地理位置和地理环境，得到很有实践性的结论即考虑不同方法各自的优点特性，在不同区域选用不同的方法。3、 模型中采用了多个指标，建立多个模型来评价两种预测方法的准确性，通过比较不同模型的结果，验证结论的合理性。4、 模型中采用了气象学中标准的预报质量评价指标，定性相关系数与回归系数，简易而不失一般性。（二） 模型的不足：1、 本文采用的经纬度坐标转换直角坐标的公式简单易行，但是不具有普遍适用性。2、 本文建立的模型适合于比较各种方法的准确度，不适合对一个方法的独立评价。3、 模型需要处理大量的数据，比较繁杂，但这受到了问题本身的限制。（三） 模型的改进：计算预测值和观测值的相关系数和回归系数时采用的是二分类预报，可以推广成为多分类例如题中所给的七分类情况时的相关系数与回归系数，适用于不同分类数要求的情况，而且比二分类方法更加精准。参考文献：（苏）米凯拉德杰，《数学方法的数值方法》，北京市：科学技术出版社，1957。章基嘉、葛玲，《中长期天气预报基础》，北京市：气象出版社，1983美]保罗•萨缪尔森、威廉•诺德豪斯，萧琛等译，《经济学》北京市：华夏出版社、麦格劳•希尔出版公司，1999。《最新实用中国地图册》，北京市：中国地图出版社，1993。张志勇、徐彦琴等，《MATLAB教程》北京市：北京航空航天出版社，2001。周义仓、赫孝良，《数学建模实验》，西安市：西安交通大学出版社，1999。甘筱青，《数学建模教育及竞赛》，南昌市：江西高教出版社，2004。吴江，《经济数学模型及其应用》，南昌市：江西高教出版社，1995。谭冠日，《气象站数理统计预报方法》，北京市：科学出版社，1980。屠其璞、王俊德、丁裕国、史慧敏，《气象应用概率统计学》，北京市：气象出版社，1982。附录［1］：球面距离公式的推导球面距AB是指过球面A、B和球心O的大圆劣弧长AB。经度差算法：I I，若 |180;1212360 ，若 1801212其中若、 分别是东经和西经的度数时，视其中一个的度数为负值。此题情况，A、ABB两点分居同半球不同纬度、，经度差为 | |。12121在这种情况下，所求AB两点的球面距离，便是要求过AB两点的球大圆上劣弧AB的长根据平面几何知识，要求AB劣弧，只需求出弦AB所对的圆心角ZAOB（它不超过180°）,在这里，就是只要求出AB弦关于球心的张角=ZAOB即可，其范围仍限定在0°—180°，即0 180，并称此角为球心角（单位是弧度）。为了计算崩长，我们只需求出球心角=ZAOB即可。证明过程：Rt00证明过程：Rt00A中，0011Rt00B中，00220A0BR00000012210Asin，A00AcosA1A0Bsin，B00BcosB2BRsinRsinBA0101A与SB是异面直线，它们的公垂线是OU精度差 aB，利用异面直线上两点间距离公式：AB0两点间距离公式：AB0A2 0B21200220A0Bcos1212A0B中,cos0A20B2A0B中,cos0A20B2AB220A0BAB的球面距离为R综合整理得到球面距离公式：SabRarccos(sinsinABcoscoscosAB A附录［2］：经纬度坐标转换为大地坐标公式由于经纬度坐标不对应大地坐标，得到不规则的图形，通过把经纬度坐标值转换为大地坐标来修正图形。具体方法：把经纬度坐标平面最左下角的点(,)设为新坐标系下的原点(0,0)，经纬度坐标ij平面上最下端的一系列点('')作为新坐标系的横坐标轴上(x,l)点，经纬度坐标平i1i面上最左端的一系列点(",〃)坐落在新坐标系的纵坐标上1,y)。任意一点的坐标1jj(xy)的坐标值，横坐标为该点到对应同行边界点(",")的球面距离，纵坐标值为i,j1j该点到同列边界点(',')的球面距离。i1xRarccos(sinsin"coscos"cos("))ijjjji1yRarccos(sinsin'coscos'cos('))ij1j1ii附录［3］：91个观测点的在修正后的新坐标直角坐标系下的坐标求法：图1al(i),bl(i)表示观测点在新坐标系下对应的横纵坐标值，L1,L2,L3表示观测点到A,B两点间的距离和A,B的球面距离。(i从1到91)修正后的欧氏坐标下的图见图3L1arccos{siny1(i)/180siny(1,1)/180cosy(1,1)/180cosy(1,1)/180cosx1(i)x(1,1)/1806370L2arccos{siny1(i)/180siny(1,47)/180cosy1(i)/180cosy(1,47)/180cosx1(i)x(1,47)/1806370L3arccos{siny(1,1)/180siny(1,47)/180cosy(1,1)/180cosy(1,47)/180cosx(1,1)x(1,47)/1806370a1(i)(k2L12L22)/(2k)b1(i)sin{arcco[s(k2L12L22)/(2L3L1)]}L1附录[4]：程序functionshicedian()f=dlmread('020618.six');x1=f(:,3);y1=f(:,2);x=dlmread('lon.dat');y=dlmread('lat.dat');plot(x1,y1,'ro',x(:,1:47),y(:,1:47),'.k')legend观测点'预测点')axis([117125 27.535])pausea(1,1)=0;b(1,1)=0;fori=1:1:53forj=1:1:47a(i,j)=acos(sin(y(i,1)/180*pi)*sin(y(i,j)/180*pi)+cos(y(i,1)/180*pi)*cos(y(i,j)pi)*cos((x(i,j)-x(i,1))/180*pi))*6370;endendfori=1:1:53forj=1:1:47b(i,j)=acos(sin(y(1,j)/180*pi)*sin(y(i,j)/180*pi)+cos(y(1,j)/180*pi)*cos(y(i,j)pi)*cos((x(i,j)-x(1,j))/180*pi))*6370;endend%a=linspace(a(1,1),a(1,47),47);%b=linspace(b(1,1),b(53,1),53);%[a,b]=meshgrid(a,b);k=acos(sin(y(1,1)/180*pi)*sin(y(1,47)/180*pi)+cos(y(1,1)/180*pi)*cos(y(1,47)/180*pi)*co(x(1,1)-x(1,47))/180*pi))*6370;fori=1:1:91l1=acos(sin(y1(i)/180*pi)*sin(y(1,1)/180*pi)+cos(y1(i)/180*pi)*cos(y(1,1)/180*pi)*c((x1(i)-x(1,1))/180*pi))*6370;l2=acos(sin(y1(i)/180*pi)*sin(y(1,47)/180*pi)+cos(y1(i)/180*pi)*cos(y(1,47)/180*p)*cos((x1(i)-x(1,47))/180*pi))*6370;al(i)=(k“2+l「2—l2"2)/(2*k);endfori=1:1:911=acos(sin(y1(i)/180*pi)*sin(y(1,1)/180*pi)+cos(y1(i)/180*pi)*cos(y(1,1)/180*pi)*cos(1(i)-x(1,1))/180*pi))*6370;l2=acos(sin(y1(i)/180*pi)*sin(y(1,47)/180*pi)+cos(y1(i)/180*pi)*cos(y(1,47)/180*pi)*c((x1(i)-x(1,47))