Accession number: 20101712882095
Title: The research of parallel lagrange interpolation method and its application in dew-point detection
Authors: Liu, Chunfeng1 ; Wang, Ai-Ling2 ; Yang, Aimin1
Author affiliation: 1 College of Science, Hebei Polytechnic University, Tangshan, Hebei Province, China
2 Department of Mathematics, Heze University, Heze, Shandong 274000, China
Corresponding author: Liu, C.
Source title: Proceedings of the International Symposium on Test and Measurement
Abbreviated source title: Proc Int Symp Test Meas
Volume: 2
Monograph title: ICTM 2009 - 2009 International Conference on Test and Measurement
Issue date: 2009
Publication year: 2009
Pages: 192-196
Article number: 5413076
Language: English
ISBN-13: 9781424447008
Document type: Conference article (CA)
Conference name: 2009 International Conference on Test and Measurement, ICTM 2009
Conference date: December 5, 2009 - December 6, 2009
Conference location: Hong Kong, Hong kong
Conference code: 79911
Sponsor: Institute ofElectrical and Electronics Engineers; IEEE Instrumentation and Measurement Society; Intelligent Information Technology; Application Research Association
Publisher: International Academic Publishers, 137 Chaonei Dajie, Beijing, 100010, China
Abstract: In this paper, a simple introduction about the interpolation method and the parallel algorithms is given, and Lagrange interpolation is carried out and applied in accordance with the idea of the parallel, analysing the developments and principles of the parallel, giving the application that the parallel algorithm is used in the Lagrange interpolation, then the parallel algorithms' advantages are showed in the numerical analysis according to the speedup ratio and efficiency. After careful scrutiny, parallel algorithm is applied to the interpolation reduces the time greatly, further, the fields of the parallel algorithm and Lagrange interpolation in practical applications are expanded. Finally it is applied to the low humidity dew point detection area with low humidity environment. The parallel Lagrange interpolation method is studied in the non-linear relationship which is shown by the measure of the capacitive sensors frequency and dew-point values. © 2009 IEEE.
Number of references: 23
Main heading: Interpolation
Controlled terms: Lagrange multipliers - Moisture - Parallel algorithms
Uncontrolled terms: Capacitive sensor - Dew points - Dewpoint - Interpolation method - Lagrange interpolation method - Lagrange interpolations - Low humidity - Non-linear relationships - Speedup ratio
Classification code: 443.1 Atmospheric Properties - 723 Computer Software, Data Handling and Applications - 801.4 Physical Chemistry - 921 Mathematics - 921.6 Numerical Methods
DOI: 10.1109/ICTM.2009.5413076
Database: Compendex
Compilation and indexing terms, © 2009 Elsevier Inc.
题目：平行拉格朗日插值方法及其在露点检测中的应用研究
摘要：本文对插值法和并行算法做了简要介绍，并给出发展和平行原则下拉格朗日插值的应用程序，分析了给并行算法的发展趋势和原则，是用在拉格朗日插值，那么'并行算法的优点是在数值分析表明根据加速比和效率。经过认真审议，并行算法应用到的时间大大减少插值，进一步的并行算法和拉格朗日插值在实际应用领域不断扩大。最后将它应用到低湿度，低湿度环境露点检测区域。平行拉格朗日插值方法研究了非线性的关系，这是由电容式传感器的频率和露点值测量显示。
1 中国古代数学的发展
在古代世界四大文明中，中国数学持续繁荣时期最为长久。从公元前后至公元14世纪，中国古典数学先后经历了三次发展高潮，即两汉时期、魏晋南北朝时期和宋元时期，并在宋元时期达到顶峰。
与以证明定理为中心的希腊古典数学不同，中国古代数学是以创造算法特别是各种解方程的算法为主线。从线性方程组到高次多项式方程，乃至不定方程，中国古代数学家创造了一系列先进的算法(中国数学家称之为“术”)，他们用这些算法去求解相应类型的代数方程，从而解决导致这些方程的各种各样的科学和实际问题。特别是，几何问题也归结为代数方程，然后用程式化的算法来求解。因此，中国古代数学具有明显的算法化、机械化的特征。以下择要举例说明中国古代数学发展的这种特征。
1.1 线性方程组与“方程术”
中国古代最重要的数学经典《九章算术》(约公元前2世纪)卷8的“方程术”，是解线性方程组的算法。以该卷第1题为例，用现代符号表述，该问题相当于解一个三元一次方程组：
3x+2y+z＝39
2x+3y+z＝34
x+2y+3z＝26
《九章》没有表示未知数的符号，而是用算筹将x�y�z的系数和常数项排列成一个(长)方阵：
1 2 3
2 3 2
3 1 1
26 34 39
“方程术”的关键算法叫“遍乘直除”，在本例中演算程序如下：用右行(x)的系数(3)“遍乘”中行和左行各数，然后从所得结果按行分别“直除”右行，即连续减去右行对应各数，就将中行与左行的系数化为0。反复执行这种“遍乘直除”算法，就可以解出方程。很清楚，《九章算术》方程术的“遍乘直除” 算法，实质上就是我们今天所使用的解线性方程组的消元法，以往西方文献中称之为“高斯消去法”，但近年开始改变称谓,如法国科学院院士、原苏黎世大学数学系主任P.Gabriel教授在他撰写的教科书[4]中就称解线性方程组的消元法为“张苍法”，张苍相传是《九章算术》的作者之一。
1.2 高次多项式方程与“正负开方术”
《九章算术》卷4中有“开方术”和“开立方术”。《九章算术》中的这些算法后来逐步推广到开更高次方的情形，并且在宋元时代发展为一般高次多项式方程的数值求解。秦九韶是这方面的集大成者，他在《数书九章》(1247年)一书中给出了高次多项式方程数值解的完整算法，即他所称的“正负开方术”。
用现代符号表达，秦九韶“正负开方术”的思路如下：对任意给定的方程
f(x)=a0xn+a1xn-1+……+an-2x2+an-1x+an=0 (1)
其中a0≠0，an<0,要求(1)式的一个正根。秦九韶先估计根的最高位数字，连同其位数一起称为“首商”，记作c，则根x＝c+h，代入(1)得
f(c+h)=a0(c+h)n+a1(c+h)n-1+……+an-1(c+h)+an=0
按h的幂次合并同类项即得到关于h的方程:
f(h)=a0hn+a1hn-1+……+an-1h+an=0 (2)
于是又可估计满足新方程(2)的根的最高位数字。如此进行下去，若得到某个新方程的常数项为0，则求得的根是有理数；否则上述过程可继续下去，按所需精度求得根的近似值。
如果从原方程(1)的系数a0,a1，…,an及估值c求出新方程(2)的系数a0,a1，…,an的算法是需要反复迭代使用的，秦九韶给出了一个规格化的程序，我们可称之为“秦九韶程序”， 他在《数书九章》中用这一算法去解决各种可以归结为代数方程的实际问题，其中涉及的方程最高次数达到10次，秦九韶解这些问题的算法整齐划一，步骤分明，堪称是中国古代数学算法化、机械化的典范。
1.3 多元高次方程组与“四元术”
绝不是所有的问题都可以归结为线性方程组或一个未知量的多项式方程来求解。实际上，可以说更大量的实际问题如果能化为代数方程求解的话，出现的将是含有多个未知量的高次方程组。
多元高次方程组的求解即使在今天也绝非易事。历史上最早对多元高次方程组作出系统处理的是中国元代数学家朱世杰。朱世杰的《四元玉鉴》(1303年)一书中涉及的高次方程达到了4个未知数。朱世杰用“四元术”来解这些方程。“四元术”首先是以“天”、“地”、“人”、“物”来表示不同的未知数，同时建立起方程式，然后用顺序消元的一般方法解出方程。朱世杰在《四元玉鉴》中创造了多种消元程序。
通过《四元玉鉴》中的具体例子可以清晰地了解朱世杰“四元术”的特征。值得注意的是，这些例子中相当一部分是由几何问题导出的。这种将几何问题转化为代数方程并用某种统一的算法求解的例子，在宋元数学著作中比比皆是，充分反映了中国古代几何代数化和机械化的倾向。
1.4 一次同余方程组与“中国剩余定理”
中国古代数学家出于历法计算的需要，很早就开始研究形如：
X≡Ri (mod ai) i=1,2,...,n (1)
(其中ai 是两两互素的整数)的一次同余方程组求解问题。公元4世纪的《孙子算经》中已有相当于求解下列一次同余组的著名的“孙子问题”：
X≡2(mod3) ≡3(mod5) ≡2(mod7)
《孙子算经》作者给出的解法，引导了宋代秦九韶求解一次同余组的一般算法——“大衍求一术”。现代文献中通常把这种一般算法称为“中国剩余定理”。
1.5 插值法与“招差术”
插值算法在微积分的酝酿过程中扮演了重要角色。在中国，早从东汉时期起，学者们就惯用插值法来推算日月五星的运动。起初是简单的一次内插法，隋唐时期出现二次插值法(如一行《大衍历》，727年)。由于天体运动的加速度也不均匀，二次插值仍不够精密。随着历法的进步，到了宋元时代，便产生了三次内插法(郭守敬《授时历》，1280年)。在此基础上，数学家朱世杰更创造出一般高次内插公式，即他所说的“招差术”。 朱世杰的公式相当于
f(n)=n△+ n(n�1)△2+ n(n�1)(n�2)△3
+ n(n�1)(n�2)(n�3)△4+……
这是一项很突出的成就。
这里不可能一一列举中国古代数学家的所有算法，但仅从以上介绍不难看到，古代与中世纪中国数学家创造的算法，有许多即使按现代标准衡量也达到了很高的水平。这些算法所表达的数学真理，有的在欧洲直到18世纪以后依赖近代数学工具才重新获得(如前面提到的高次代数方程数值求解的秦九韶程序，与1819年英国数学家W. 霍纳重新导出的“霍纳算法”基本一致；多元高次方程组的系统研究在欧洲也要到18世纪末才开始在E. 别朱等人的著作中出现；解一次同余组的剩余定理则由欧拉与高斯分别独立重新获得；至于朱世杰的高次内插公式，实质上已与现在通用的牛顿-格列高里公式相一致)。这些算法的结构，其复杂程度也是惊人的。如对秦九韶“大衍求一术”和“正负开方术”的分析表明，这些算法的计算程序，包含了现代计算机语言中构造非平易算法的基本要素与基本结构。这类复杂的算法，很难再仅仅被看作是简单的经验法则了，而是高度的概括思维能力的产物，这种能力与欧几里得几何的演绎思维风格截然不同，但却在数学的发展中起着完全可与之相媲美的作用。事实上，古代中国算法的繁荣，同时也孕育了一系列极其重要的概念，显示了算法化思维在数学进化中的创造意义和动力功能。以下亦举几例。
1.6 负数的引进
《九章算术》“方程术”的消元程序，在方程系数相减时会出现较小数减较大数的情况，正是在这里，《九章算术》的作者们引进了负数，并给出了正、负数的加减运算法则，即“正负术”。
对负数的认识是人类数系扩充的重大步骤。公元7世纪印度数学家也开始使用负数，但负数的认识在欧洲却进展缓慢，甚至到16世纪，韦达的著作还回避负数。
1.7 无理数的发现
中国古代数学家在开方运算中接触到了无理数。《九章算术》开方术中指出了存在有开不尽的情形：“若开方不尽者，为不可开”，《九章算术》的作者们给这种不尽根数起了一个专门名词——“面”。“面”，就是无理数。与古希腊毕达哥拉斯学派发现正方形的对角线不是有理数时惊慌失措的表现相比，中国古代数学家却是相对自然地接受了那些“开不尽”的无理数，这也许应归功于他们早就习惯使用的十进位制，这种十进位制使他们能够有效地计算“不尽根数”的近似值。为《九章算术》作注的三国时代数学家刘徽就在“开方术”注中明确提出了用十进制小数任意逼近不尽根数的方法，他称之为“求微数法”，并指出在开方过程中，“其一退以十为步，其再退以百为步，退之弥下，其分弥细，则……虽有所弃之数，不足言之也”。
十进位值记数制是对人类文明不可磨灭的贡献。法国大数学家拉普拉斯曾盛赞十进位值制的发明，认为它“使得我们的算术系统在所有有用的创造中成为第一流的”。中国古代数学家正是在严格遵循十进位制的筹算系统基础上，建立起了富有算法化特色的东方数学大厦。
1.8 贾宪三角或杨辉三角
从前面关于高次方程数值求解算法(秦九韶程序)的介绍我们可以看到，中国古代开方术是以�c+h n的二项展开为基础的，这就引导了二项系数表的发现。南宋数学家杨辉著《详解九章算法》(1261年)中，载有一张所谓“开方作法本源图”，实际就是一张二项系数表。这张图摘自公元1050年左右北宋数学家贾宪的一部著作。“开方作法本源图”现在就叫“贾宪三角”或“杨辉三角”。二项系数表在西方则叫“帕斯卡三角”�1654年 。
1.9 走向符号代数
解方程的数学活动，必然引起人们对方程表达形式的思考。在这方面，以解方程擅长的中国古代数学家们很自然也是走在了前列。在宋元时期的数学著作中，已出现了用特定的汉字作为未知数符号并进而建立方程的系统努力。这就是以李冶为代表的“天元术”和以朱世杰为代表的“四元术”。所谓“天元术”，首先是“立天元一为某某”，这相当于“设为某某”，“天元一”就表示未知数，然后在筹算盘上布列“天元式”，即一元方程式。该方法被推广到多个未知数情形，就是前面提到的朱世杰的“四元术”。因此，用天元术和四元术列方程的方法，与现代代数中的列方程法已相类似。
符号化是近世代数的标志之一。中国宋元数学家在这方面迈出了重要一步，“天元术”和“四元术”，是以创造算法特别是解方程的算法为主线的中国古代数学的一个高峰�。
2 中国古代数学对世界数学发展的贡献
数学的发展包括了两大主要活动：证明定理和创造算法。定理证明是希腊人首倡，后构成数学发展中演绎倾向的脊梁；算法创造昌盛于古代和中世纪的中国、印度，形成了数学发展中强烈的算法倾向。统观数学的历史将会发现，数学的发展并非总是演绎倾向独占鳌头。在数学史上，算法倾向与演绎倾向总是交替地取得主导地位。古代巴比伦和埃及式的原始算法时期，被希腊式的演绎几何所接替，而在中世纪，希腊数学衰落下去，算法倾向在中国、印度等东方国度繁荣起来；东方数学在文艺复兴前夕通过阿拉伯传播到欧洲，对近代数学兴起产生了深刻影响。事实上，作为近代数学诞生标志的解析几何与微积分，从思想方法的渊源看都不能说是演绎倾向而是算法倾向的产物。
从微积分的历史可以知道，微积分的产生是寻找解决一系列实际问题的普遍算法的结果�6�。这些问题包括：决定物体的瞬时速度、求极大值与极小值、求曲线的切线、求物体的重心及引力、面积与体积计算等。从16世纪中开始的100多年间，许多大数学家都致力于获得解决这些问题的特殊算法。牛顿与莱布尼兹的功绩是在于将这些特殊的算法统一成两类基本运算——微分与积分，并进一步指出了它们的互逆关系。无论是牛顿的先驱者还是牛顿本人，他们所使用的算法都是不严格的，都没有完整的演绎推导。牛顿的流数术在逻辑上的瑕疵更是众所周知。对当时的学者来说，首要的是找到行之有效的算法，而不是算法的证明。这种倾向一直延续到18世纪。18世纪的数学家也往往不管微积分基础的困难而大胆前进。如泰勒公式，欧拉、伯努利甚至19世纪初傅里叶所发现的三角展开等，都是在很长时期内缺乏严格的证明。正如冯·诺伊曼指出的那样：没有一个数学家会把这一时期的发展看作是异端邪道；这个时期产生的数学成果被公认为第一流的。并且反过来，如果当时的数学家一定要在有了严密的演绎证明之后才承认新算法的合理性，那就不会有今天的微积分和整个分析大厦了。
现在再来看一看更早的解析几何的诞生。通常认为，笛卡儿发明解析几何的基本思想，是用代数方法来解几何问题。这同欧氏演绎方法已经大相径庭了。而事实上如果我们去阅读笛卡儿的原著，就会发现贯穿于其中的彻底的算法精神。《几何学》开宗明义就宣称：“我将毫不犹豫地在几何学中引进算术的术语，以便使自己变得更加聪明”。众所周知，笛卡儿的《几何学》是他的哲学著作《方法论》的附录。笛卡儿在他另一部生前未正式发表的哲学著作《指导思维的法则》(简称《法则》)中曾强烈批判了传统的主要是希腊的研究方法，认为古希腊人的演绎推理只能用来证明已经知道的事物，“却不能帮助我们发现未知的事情”。因此他提出“需要一种发现真理的方法”，并称之为“通用数学”(mathesis universakis)。笛卡儿在《法则》中描述了这种通用数学的蓝图，他提出的大胆计划，概而言之就是要将一切科学问题转化为求解代数方程的数学问题：
任何问题→数学问题→代数问题→方程求解而笛卡儿的《几何学》，正是他上述方案的一个具体实施和示范，解析几何在整个方案中扮演着重要的工具作用，它将一切几何问题化为代数问题，这些代数问题则可以用一种简单的、几乎自动的或者毋宁说是机械的方法去解决。这与上面介绍的古代中国数学家解决问题的路线可以说是一脉相承。
因此我们完全有理由说，在从文艺复兴到17世纪近代数学兴起的大潮中，回响着东方数学特别是中国数学的韵律。整个17—18世纪应该看成是寻求无穷小算法的英雄年代，尽管这一时期的无穷小算法与中世纪算法相比有质的飞跃。而从19世纪特别是70年代直到20世纪中，演绎倾向又重新在比希腊几何高得多的水准上占据了优势。因此，数学的发展呈现出算法创造与演绎证明两大主流交替繁荣、螺旋式上升过程：
演绎传统——定理证明活动
算法传统——算法创造活动
中国古代数学家对算法传统的形成与发展做出了毋容置疑的巨大贡献。
我们强调中国古代数学的算法传统，并不意味中国古代数学中没有演绎倾向。事实上，在魏晋南北朝时期一些数学家的工作中，已出现具有相当深度的论证思想。如赵爽勾股定理证明、刘徽“阳马”�一种长方锥体 体积证明、祖冲之父子对球体积公式的推导等等，均可与古希腊数学家相应的工作媲美。赵爽勾股定理证明示意图“弦图”原型，已被采用作2002年国际数学家大会会标。令人迷惑的是，这种论证倾向随着南北朝的结束，可以说是戛然而止。囿于篇幅和本文重点，对这方面的内容这里不能详述，有兴趣的读者可参阅参考文献�3�。
3 古为今用，创新发展
到了20世纪，至少从中叶开始，电子计算机的出现对数学的发展带来了深远影响，并孕育出孤立子理论、混沌动力学、四色定理证明等一系列令人瞩目的成就。借助计算机及有效的算法猜测发现新事实、归纳证明新定理乃至进行更一般的自动推理……，这一切可以说已揭开了数学史上一个新的算法繁荣时代的伟大序幕。科学界敏锐的有识之士纷纷预见到数学发展的这一趋势。在我国，早在上世纪50年代，华罗庚教授就亲自领导建立了计算机研制组，为我国计算机科学和数学的发展奠定了基础。吴文俊教授更是从70年代中开始，毅然由原先从事的拓扑学领域转向定理机器证明的研究，并开创了现代数学的崭新领域——数学机械化。被国际上誉为“吴方法”的数学机械化方法已使中国在数学机械化领域处于国际领先地位，而正如吴文俊教授本人所说：“几何定理证明的机械化问题，从思维到方法，至少在宋元时代就有蛛丝马迹可寻，”他的工作“主要是受中国古代数学的启发”。“吴方法”，是中国古代数学算法化、机械化精髓的发扬光大。
计算机影响下算法倾向的增长，自然也引起一些外国学者对中国古代数学中算法传统的兴趣。早在上世纪70年代初，著名的计算机科学家D.E.Knuth就呼吁人们关注古代中国和印度的算法�5�。多年来这方面的研究取得了一定进展，但总的来说还亟待加强。众所周知，中国古代文化包括数学是通过著名的丝绸之路向西方传播的，而阿拉伯地区是这种文化传播的重要中转站。现存有些阿拉伯数学与天文著作中包含有一定的中国数学与天文学知识，如著名的阿尔·卡西《算术之钥》一书中有相当数量的数学问题显示出直接或间接的中国来源，而根据阿尔·卡西本人记述，他所工作的天文台中就有不少来自中国的学者。
然而长期以来由于“西方中心论”特别是“希腊中心论”的影响以及语言文字方面的障碍，有关资料还远远没有得到发掘。正是为了充分揭示东方数学与欧洲数学复兴的关系，吴文俊教授特意从他荣获的国家最高科学奖中拨出专款成立了“吴文俊数学与天文丝路基金”，鼓励支持年轻学者深入开展这方面的研究，这是具有深远意义之举。
研究科学的历史，其重要意义之一就是从历史的发展中获得借鉴和汲取教益，促进现实的科学研究，通俗地说就是“古为今用”。吴文俊对此有精辟的论述，他说：“假如你对数学的历史发展，对一个领域的发生和发展，对一个理论的兴旺和衰落，对一个概念的来龙去脉，对一种重要思想的产生和影响等这许多历史因素都弄清了，我想，对数学就会了解得更多，对数学的现状就会知道得更清楚、更深刻，还可以对数学的未来起一种指导作用，也就是说，可以知道数学究竟应该按怎样的方向发展可以收到最大的效益”。数学机械化理论的创立，正是这种古为今用原则的硕果。我国科学技术的伟大复兴，呼唤着更多这样既有浓郁的中国特色、又有鲜明时代气息的创新。
地层面的拟合是多源地质建模中最为重要的步骤。无论是通过Delaunay细分方法增加节点还是通过网格分级增加的节点，都需要进一步求取其高程值。因此，必须借助插值方法来对有限的数据点信息所形成的初始地层面进行精确、光滑处理。插值是指在根据已知的数据计算未知值的过程，结果是形成一个连续分布的数据场(Spragueetal.，2005)。目前存在的插值方法众多，但适合三维数据场表达的插值方法还是比较有限的(胡小红等，2007)，常用的有距离加权反比法(Inverse-Distance Weighting，IDW)和普通Kriging法，其中距离加权反比法属于一种确定性差值方法，Lu et al.(2008)对该方法进行扩展，根据样本的数量和分布密度特征，使其能够根据样本的特征来确定其参数的取值；而Kriging法属于一种不确定性差值方法，Jessell(2001)对其进行了深入研究，基于该方法提出一种势场(potential-field)的插值方法，能够处理存在断层的不连续数据场。这里，仅对这两种方法进行讨论。
5.3.3.1 距离加权反比法
距离反比加权法是最常用的地质数据插值方法之一。它首先由气象学家及地质学工作者提出，后由D.Shepard进行改进，故该方法被称为Shepard方法。
距离反比加权法的基本思想是将插值函数f(P)定义为各数据点函数fk的加权平均，它认为与待插值点距离最近的若干个已知采样点对待插值点的值贡献最大，其贡献与距离的某次幂成反比。
距离反比加权法的基本原理可用下式表示：
数字地下空间与工程三维地质建模及应用研究
式中：f(P)是待插值点P的估计值；fi是第i(i=1，…，n)个已知采样点Pi的样本值；di是第i个样本点Pi与待插值点P的广义距离；v是距离的幂，它显著影响内插的结果，它的选择标准是最小平均绝对误差。相关研究结果表明，幂越高，内插结果越具有平滑的效果(Lu et al.，2008)。在幂指数为2时，不仅能得出较满意的内插结果，而且具有容易计算的优点。在实际应用中，一般用距离平方反比法来求待估算值。
传统的距离反比加权法是非常简单和自然的，但应用于实际的地质特征插值时，却存在以下明显的缺陷：
(1)当数据点的数目非常庞大时，f(P)的计算量将变得十分巨大，计算量的庞大甚至可能导致该方法变得无法实现。
(2)该方法只考虑了从Pi到P的距离，而没有考虑其方向。事实上，只考虑距离的大小是不充分的，有的已知离散点虽然距离待插值点较近，但它对待插值点的影响可能会被其他点屏蔽掉。
(3)在已知采样点Pi的邻域内，由于di≈0计算的误差将变得非常敏感，尤其是当两个项形式占优而又符号相反时，计算的误差更是如此。
因此，在实际应用时，必须对传统的距离反比加权法加以改进。考虑到地质特征数据的空间相关性，对距离反比加权法可附加地质体结构和影响距离两个限制条件，以提高空间几何或属性数据插值的合理性和精度。具体改进如下：
(1)地质体几何结构限制条件。在空间特征插值过程中，在选择影响待估计值的原始样本数据时只选取同一地层岩性地质内的样本；不是同一地层岩性地质体内的样本，即便距离很近，也不采用。即按照地质构造划分空间单元，并只在同一地层岩性地质体内选取样本点。
(2)邻近样本点的选择条件。选择待插值点的邻近点时，可考虑三个原则：一是距离原则，即根据地质属性数据的特征给出一距离r，在该距离之外的样本点对待插值点的估算无影响；二是点数原则，即给定一数据m，以距离待插值点最近m个样本点进行估算；三是利用Voronoi图求取待插值点的邻近点。
(3)建立数据点索引表，提高待插值点周围样本点的搜索效率，从而大大减少大数据量队的计算量。
5.3.3.2 普通Kriging插值方法
设研究区域为A，区域化变量(即欲研究的物理属性变量)为 {Z(x)∈A}，x表示空间位置。Z(x)在采样点xi(i=1，2，…，n)处的属性值(或称为区域化变量的一次实现)为Z(xi)(i=1，2，…，n)，则根据普通Kriging插值原理，未采样点x0处的属性值Z(x0)估计值是n个已知采样点属性值的加权和，即：
数字地下空间与工程三维地质建模及应用研究
其中，λi(i=1，2，…，n)为待求权系数。
假设区域化变量Z(x)在整个研究区域内满足二阶平稳假设：
(1)Z(x)的数学期望存在且等于常数：E[Z(x)]=m(常数)。
(2)Z(x)的协方差Cov(xi，xj)存在且只与两点之间的相对位置有关。或满足本征假设：
(3)E[Z(xi)-Z(xj)]=0。
(4)增量的方差存在且平稳：Var[Z(xi)-Z(xj)]=E[Z(xi)-Z(xj)]2。
依据无偏性要求：E[Z*(x0)]=E[Z(x0)]。
推导可得： 。
在无偏条件下使估计方差达到最小，即：
min{ Var[Z*(x0)-2μ( (λi-1))}，其中μ为拉格朗日乘子。
可求得求解权系数λii(=1，2，…，n)的方程组：
数字地下空间与工程三维地质建模及应用研究
求出诸权系数λii(=1，2，…，n)后，就可以求出采样点x0处的属性值Z*(x0)。
上述求解权系数λii(=1，2，…，n)的方程组中协方差Cov(xi，xj)若用变异函数γ(xi，xj)表示时，形式为：
数字地下空间与工程三维地质建模及应用研究
变异函数的定义为：
数字地下空间与工程三维地质建模及应用研究
由Kringing插值所得到的方差为：
数字地下空间与工程三维地质建模及应用研究
或
数字地下空间与工程三维地质建模及应用研究
A.Kringing插值法中的变异函数
变异函数是Kringing插值法插值的基础。插值中需要首先确定所研究的区域化变量的变异函数。假设研究的区域为A，区域A中有一区域化变量Z(x)，它在位置xi(i=1，2，…，N)上的一次采样为Z(xi)(i=1，2，…，N)，则Z(x)的变异函数的定义为：
数字地下空间与工程三维地质建模及应用研究
一个空间变量的空间变异性是指这个变量在空间中如何随着位置的不同而变化的性质。变异函数通过其自身的结构及其各项参数从不同的角度反映空间变异性，确定变异函数的过程就是一个对空间变异性进行结构分析的过程。
设h是一个模为r=|h|，方向为a的向量，如果存在着被向量h所隔开的Nh对观测数据点，则在a方向上相应于向量h的实验变异函数γ*(h)可表示为如下形式：
数字地下空间与工程三维地质建模及应用研究
其中，z(xi+h)和z(xi)分别位于点xi+h和xi(i=1，2，…，Nh)上的观测数据。
B.变异函数理论模型
当获取实验变异函数值后，需要先选择变异函数理论模型，然后对所选择的变异函数理论模型进行参数拟合，这一过程被称为“结构分析”。
变异函数理论模型参数一般包括：变程(range，一般用a表示)、基台(sill，一般用C(0)表示)、拱高(一般用C表示)、块金常数(Nugget，一般用C0表示)，如图5.9所示。
图5.9 理论变差函数
变程a表示了从空间相关性状态(|h|＜a)向不存在相关性状态(|h|＞a)转变的分界线；变异函数在原点处的间断性称为块金效应，相应的常数C0= (h)称为块金常数；基台C(0)具有协方差函数：C(h)=C(0)-γ(h)的二阶平稳区域化变量Z(x)的先验方差：Var{Z(x)}=C(0)=γ(h)；拱高C为变异函数中基台C(0)与块金常数C0之差：C=C(0)-C0。
C.变异函数理论模型分类
变异函数理论模型一般分为有基台值和无基台值两大类。有基台值的变异函数理论模型包括球状模型、指数模型、高斯模型等(图5.10)。最常用的是球状模型。球状模型公式为：
数字地下空间与工程三维地质建模及应用研究
图5.10 变异函数中的球状模型、指数模型、高斯模型
D.变异函数对空间变异性结构的反映
变异函数作为定量描述空间变异性的一种统计学工具，通过其自身的结构及其各项参数，从不同角度反映了空间变异性结构。利用变异函数可以对空间变量的连续性、相关性、变量的影响范围、尺度效应、原点处的间断性、各向异性等要素进行描述。
E.变异函数理论模型参数拟合
变异函数理论模型参数拟合就是利用原始采样点数据或实验变异函数取值对所选定的理论模型参数以特定的方法进行估计。拟合方法一般采用手工拟合法。
手工拟合就是依据实验变异函数的取值，一方面通过观察实验变异函数图；另一方面对所研究的区域化变量进行必要的分析，采用肉眼观察来确定变异函数模型参数，并对参数反复进行交叉验证，最终确定模型参数。其拟合的大致过程如下：
(1)首先对所研究的区域化变量进行必要的结构、背景等方面的分析，结合专家经验，确定变异函数理论模型。
(2)利用实验变异函数散点图确定变异函数参数中的块金常数、基台值、变程、各向异性角度以及各向异性比值。
(3)交叉验证。
彩色图象的二维变形_电子通信论文
摘 要 该文讨论了彩色图像的变形扭曲技术，并针对二维变形给出了一个速度、精度均令人满意的算法。
一、引言
在图像处理的应用中，一般图像所覆盖区域边界是规则的矩形。为获得某种特殊效果，常常需要将图像变换到具有任意不规则边界的二维区域或映像到三维空间曲面，简单地说，这就是所谓的图像变形技术。本文重点讨论了其中的任意二维多边形区域的变形问题，并针对彩色图像给出一个切实可行的算法。而三维情况下，则属于计算机图形学中的纹理贴面范围，一般均会牵涉到立体图形消隐、明暗处理等技术，比较复杂，本文未作深入探讨。
二、变换原理
本文所要讨论的二维变形问题可以形式化说明如下:图像定义在矩形区域ABCD之上，源多边形区域P=p1p2…pnp1(Pi为顶点，i=1，2，…n)完全包含在ABCD内；变形就是通过变换f，将P上的图像变换到目的多边形区域Q=Q1Q2…QnQ1(Qi为顶点，i=1，2，…n)，其中，P与Q中的各顶点一一对应，即有:Qi=f(Pi)(i=1,2…n)。图1是变形的一个简单例子:图中的源多边形区域是矩形区域ABCD，目的多边形为任意四边形EFGH，阴影部分在变换前后的变化清楚地说明了变形的效果。
@@T5S13200.GIF;图1@@
那么，变换应该如何进行呢?
一种直接的思路是显式地求出变换f的表达式。而f的实施又分两种方法；其一为正向变换法，即用f将P内的任一像素点变换到Q内，取原像素值加以显示。由于P与Q所包含像素点的数目一般不相同，甚至相差很大，造成Q中的像素点或者未被赋值，形成令人讨厌的空洞，或者被多次赋值，浪费了时间，总的效果不理想；其二利用f的反变换f-1，将Q内的每一像素点反变换至P内的对应点，一般此点具有实数坐标，则可以通过插值，确定其像素值，这样，结果图像中的每一像素点均被赋值唯一的一次，既提高了精度，又可以避免不必要的赋值，使用效果较好。
上述显示求变换(或反变换)的表达式的思路，比较精确，但是这往往牵涉到复杂的多元方程求解问题，并非轻易可以完成。本文所给出的另外一条思路是:既然P与Q中各顶点一一对应，组成变换对，即源多边形P中的任一顶点Pi(i=1,2…n)经过变换f，得到目的多边形Q中的顶点Qi(i=1,2…n)，则Qi的反变换点也必为Pi。这样，对Q内(包括边界)的各像素点A，可以利用各顶点的反变换点的坐标值通过双线性插值技术近似求出其反变换点B；再用点B的坐标值在源图像中进行插值，最终求得结果像素值，用于显示A。
第二种方法在保留一定精度的前提下，避免了变换表达式的显式求解，实现简便。本文基于此思想，设计了一个快速变形算法；另外，算法中还借鉴了多边形区域扫描转换的扫描线算法的思路，以实现对Q内各像素点的高效扫描。以下，本文首先介绍了插值技术及增量计算技术，然后将给出二维变形算法的详细步骤。
三、插值技术
已知目的多边形Q各顶点Qi(i=1,2…n)的变换坐标值，如何求出Q内任一像素的反变换坐标呢?双线性插值法是一种简单快速的近似方法。具体做法是:先用多边形顶点Qi(i=1,2…n)的反变换坐标线性插值出当前扫描线与多边形各边的交点的反变换坐标，然后再用交点的反变换坐标线性插值出扫描线位于多边形内的区段上每一像素处的反变换坐标值用于以后的计算。逐条扫描线处理完毕后，Q内每一像素点的反变换坐标值也就均求出来了。如图2中所示，扫描线Y(纵坐标=Y)与多边形相交于点A和B两点，D则是位于扫描线上位于多边形内的区段AB上的任一点。已知多边形的3个顶点Qi(i=1,2，3)的反变换坐标为(RXi,RYi)；
又令A、B及D各点的反变换坐标分别是(RXa,RYa),(RXb,RYb)和(RXd,RYd)。则RXp可按以下公式求出:
RXa=uRX1+(1-u)RX2 式1
RXb=vRX1+(1-v)RX3式2
RXd=tRXa+(1-t)RXb 式3
其中，u=|AQ2|/|Q1Q2|，v=|BQ3|/|Q1Q3|，t=|DB|/|AB|，
称为插值参数。
RYd的值亦可完全类似地求出，甚至不必改变插值参数的计算。(Rxd,Ryd)即是D点在原图像中对应点的坐标近似值。
@@T5S13201.GIF;图2@@
上述的双线性插值过程可以通过增量计算方法提高速度。其中，在水平方向上，位于多边形内的各区段上的各像素的反变换坐标可以沿扫描线从左至右递增计算。仍以反变换的X坐标为例。如图2所示，在扫描线Y上，C与D是相邻两像素点，对C点，插值参数tc=|CB|/|AB|，对D点，td=|DB|/|AB|，则插值参数之差△t=|CD|/|AB|，由于C与D相邻，且在同一扫描线上，|CD|=1，即△t=1/|AB|，在AB区段上为常数。根据式1～式3，不难推得D点的反变换X坐标Rxd与C点的反变换X坐标Rxc之间的关系如下:
Rxd=Rxc+(Rxa-Rxb)·△t=Rxc+△Rxx
由于△Rxx在AB区段仍为常数，故AB区段上各像素点的反变换X坐标均可由A点的Rxa依次递增求得，而反变换Y坐标的递增求法亦是相同。这样，AB区段上各像素点的反变换坐标值的计算简化为两次加法，时间的节省是惊人的。事实上，在垂直方向，每条边也可在相邻扫描线间递增计算其与扫描线交点的反变换坐标。如图2中的Q1Q2边，其与相邻的两条扫描线Y与Y-1分别交于A点和E点。则两点的插值参数之差△u=|AE|/|Q1Q2|，而Q1Q2边与扫描线交角固定为θ，且A和E两点的Y坐标之差为1，则有:|AE|=1/Sinθ，对于Q1Q2边而言是常量，因此△u对此边也是常量，于是推得两点反变换X坐标关系如下:
Rxa=Rxe+(Rx1-Rx2)△u=Rxe+△Rxy
显然，△Rxy沿Q1Q2边亦是常数，故而可知，相邻扫描线与各边交点的反变换坐标也只要两次浮点加法的计算量。这样，区域内每一像素点的反变换均可通过增量计算高效地完成，这大大提高了整个变形算法的速度。
另外，前面提到，经过反变换后的点一般具有实数坐标，无法直接在原图像中获得颜色值。但我们知道，一幅所谓数字图像，其实质是对连续图像在整数坐标网格点上的离散采样，因而可以用插值的方法，得到区域内具有任意坐标的点的颜色值。插值即是对任意坐标点的颜色值，用其周围的若干像素(具有整值坐标值，颜色值确定)的颜色值按一定插值公式近似计算。一般有最近邻点法、双线性插值法及3次样条函数法等插值方法，出于精度与速度的折衷要求，选用双线性插值方 法对绝大多数的应用场合是适宜的。需特别指出的是，应该对颜色的3原色分量分别进行插值，而不要直接使用读像素点得到的颜色索引号。详细讨论见文献[1]。
四、算法细节
下面将要给出的彩色图像的二维变形算法以多边形区域扫描转化的扫描线算法为框架，且使用相仿的数据结构，对目的多边形区域高效地进行逐点扫描，同时实现前面讨论的各种技术。
首先给出的是用C语言描述的数据结构:
struct Edge {
float x; /*在边的分类表ET中表示边的下端点的x坐标；在边的活化链
表AEL中则表示边与扫描线的交点的x坐标；*/
float dx; /*边的斜率的倒数；即沿扫描线间方向X的增量值*/
int Ymax; /*边的上端点的y坐标*/
float Rx; /*在ET中表示边的下端点*/
float Ry; /*的反变换坐标；在AEL中则表示边与扫描线交点的反变换坐标*
/
表float dRx; /*沿扫描线间方向，反变*/
float dRy; /*换坐标(Rx,Ry)的增量值*/
struct Edge *next;/*指向下一条边的指针*/
}; /*多边形的边的信息*/
struct Edge *ET[YResolution];
/*边的分类表，按边的下端点的纵坐标Y对非水平边进行分类的指针数组。
下端点的Y值等于i的边归入第i类，同一类中，各边按X值及△X的值递增顺序排列；YResolution为扫描线数目*/
struct Edye *AEL;
表 /*边的活化链表，由与当前扫描线相交的所有多边形的边组成，记录了多边形边沿当前扫描线的交点序列。*/
struct Polygon {
int npts; /*多边形顶点数*/
struct Point *Pts;
/*多边形的顶点序列*/
}; /*多边形信息*/
struct Point {
int X;
int Y; /*顶点坐标*/
float Rx;
float Ry; /*顶点的反变换坐标*/
}; /*多边形各顶点的信息*/
注意以上注释中边的下端点指纵坐标值较小的一端，另一端即为上端点。
以下则为算法的详细步骤:
1.数据准备
对于每一条非水平边QiQi+1，设Qi与Qi+1的坐标分别为(Xi,Yi)
及(X
i+1,Yi+1)；其反变换坐标为(Rxi,Ryi)及(RXi+1,RYi+1)。
则按以下各式对此边的信息结构各域进行填写:
X=Xi,Yi＜Yi+1
Xi+1,Yi＞Yi+1
RX=RXi,Yi＜Yi+1
RXi+1,Yi＞Yi+1
RY=RYi,Yi＜Yi+1
RYi+1,Yi＞Yi+1
dx=(xi-xi+1)/(yi-yi+1)
Ymax=max(yi，yi+1)
dRx=(Rxi-Rxi+1)/(yi-yi+1)
dRy=(Ryi-Ryi+1)/(yi-yi+1)
然后将其插入链表ET[min(yi，yi+1)]中。活化边表AEL置空。
当前扫描线纵坐标y取为0，即最小序号。
2.扫描转换
反复作以下各步，直到y等于YResolution
(1)若ET[y]非空，则将其内所有边插入AEL。
(2)若AEL非空，则将其按X及dx的值从小到大排列各边，接(3)；否则转
(3)将AEL内各边按排列顺序两两依次配对。则沿当前扫描线Y组成若干水平区间[xLeft,xRight]，其左右端点的反变换坐标分别为:(lRx,lRy)，(rRx,rRy)。则对于每一个这样的区间作以下各步:
dRxx=(lRx-rRx)/(xleft-xRight)
dRyx=(lRy-rRy)/(xleft-xRight)
又设原图像已读入二维数组Image之中。令XX=xleft, Rxy=lRx, Ryx=lRy则对于每个满足xLeft≤xX≤xRight的坐标为(xx,y)的像素，其反变换坐标(Rxy,Ryx)可按下式增量计算:
Rxx=Rxx+dRxx
Ryx=Ryx+dRyy
用(Rxx,Ryx)在数组Image之中插值，(参见文献[1])，按所得颜色值显示该像素。然后边x=x+1，计算下一像素。
(4)将AEL中满足y=Ymax的边删去，然后按下式调整AEL中各边的信息。
X=X+dx
Rx=Ry+dRx
Ry=Ry+dRy
(5)y=y+1,重复下一点。
五、讨论
上述算法针对彩色图像的二维变形问题，给出了一个简单快速的实现方案。至于三维变形，由于一般会牵涉到隐藏面消除等问题，比较复杂。但在一些情况下，可以避开消隐问题，如目的曲面形状比较简单，投影到屏幕后，各部分均不发生重叠，也就没有必要使用消隐技术，直接投影就可以了。
摘要
本文采用k-NN法，从2491个网格点的预报雨量得到91个站点预报雨量的估计值，并与对应的实测雨量进行统计对比分析。我们建立了平均绝对误差、模糊评分、面雨量三个模型，对两种雨量预报方法进行比较评价，最后为了在评价方法中考虑公众的感受，我们还构造了一个新的评价模型。结果表明方法I和方法II具有较好的晴雨预报能力，但总体上方法I的准确性高于方法II。在上述两种雨量预报方法的基础上,我们还可采用动态权重系数法对其进行综合集成，形成一种新的预报方法，该方法的可靠性要好于每种单独的预报方法。
一． 问题重述
我国某地气象台和气象研究所正在研究6小时雨量预报方法，即每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度、北纬32度附近的53×47的等距网格点上。同时设立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,2,3,4； ）
: 表示用第一种雨量预报方法测得第 天，第 个时段，第 个观测站点雨量的
预测值（ =1,2,3,4； ）
: 表示用第二种雨量预报方法测得第 天，第 个时段，第 个观测站点雨量的
预测值（ =1,2,3,4； ）
：表示第 天，第 个时段，第 个观测站点的实测雨级（ =1,2,3,4； ）
：表示用第一种预报方法测得第 天，第 个时段，第 个观测站点的预报雨级（ =1,2,3,4； ）
：表示用第二种预报方法测得第 天，第 个时段，第 个观测站点的预报雨级（ =1,2,3,4； ）
: 表示对第 种雨量预报方法在第 天，第 个时段，第 个观测站点的预测值的模糊评分（ =1,2,3,4； ）
：表示公众对用第c种预报方法测得的第 天，第 个时段，第 个观测站点的预报雨级的满意度（ =1,2,3,4； ）
三． 问题分析与数据处理
要评价两种雨量预报方法的准确性，就要在相同地点，相同时段对雨量的实测值和预测值进行比较。由于已知数据给出的91个观测站点的地理位置并不在用来预报的2491个网格点位置上，为了使得数据采集的地理位置相同，有两种思路对数据进行处理。一是把91个站点的实测数据扩充到2491个网格点上；二是利用2491个网格点的预报值给出91个站点的预报值。显然，第一种数据处理方式损失的信息比较多，而且是把预报值和处理过的实测值进行比较，其结果难以令人信服。第二种数据处理方式在保持91个观测站点的实测数据不被处理（从而保持实测数据真实可靠）的前提下，利用2491个网格点处的预报值来估计91个站点处的预报值；再对两种不同预报方法给出的预报值和实测到的数据分别进行比较。第二种方式虽然只在91个地点进行比较，但要比第一种方式更加有说服力。所以我们采取第一种数据处理方式，首先要利用2491个网格点的预报值给出91个站点的预报值。
有以下三种方法可以用来给出91个站点处的预报值。
1．k-最近邻居法（k-NN法）[1][2]
给定某个观测站点的位置后，从2491个网格点中找出距离这个站点最近的k个网格点，把这k个网格点处的预报值的平均值作为这个站点处的预报值。k的大小可以根据实际需要进行调整。
2．邻域平均法
对每个观测站点取相同的球形邻域或者正方形邻域，把邻域内网格点处的预报值的平均值作为站点处的预报值（邻域内的网格点的数目不一定相同）。当然邻域的大小可以根据实际情况进行调整。
3．二维插值法
可参见计算数学方面的参考书。
需要说明的是：二维插值法计算程序相对复杂，计算量也较大。k-NN法与邻域平均法都是利用站点周围网格点处的预报值的均值作为这个站点处的预报值。这两种方法计算相对简洁，同时也具有很高的精度。如果假设预报值在整个预报区域内是连续的话，在适当的条件下，k-NN法给出的估计将收敛到真实预报值（参见文献[1]，[2]等）。
进一步，可以认为距离站点近的网格点对站点的影响要大些，距离站点远的网格点对站点的影响要小些，因此可以根据这些影响的不同而赋予不同的权重。
本文采用3-NN法计算，以距离第 个观测站点最近的3个网格点的预报值的加权均值作为该观测站点的预报值(程序见附录)。
四． 问题（1）模型建立及求解
1．模型I
为了比较两种预报方法的预报质量，我们对其进行绝对误差值检验分析。误差值为预测值与实测值之差，绝对误差即对误差取绝对值。
第一种雨量预测方法在第 时段的平均绝对误差 = ;
第二种雨量预测方法在第 时段的平均绝对误差 =
结果如表1所示：
表1 两种预报方法的平均绝对误差
平均绝对误差
预测方法I
预测方法II
时段1(21点至次日3点)
1.8780
1.8712
时段2(次日3点至9点)
2.4936
2.4880
时段3(9点至15点)
2.1391
2.1426
时段4(15点至21点)
2.0347
1.9124
由表中得到以下结论：
两种预报方法的平均绝对误差都在2.5mm以内，但方法I和II的预报质量差距并不太大，在第一、二、四时段，方法I的平均绝对误差略小于方法II的平均绝对误差。
进一步我们还给出了两种预报方法的预报-实测相关图(图1)。图中落在斜率为1的直线上的点为实测结果，预测点则落在该直线附近，其偏离直线越远表示预报误差越大。从图中可以看出，这两种雨量预测方法的共同点是：在各个时段对实测雨量较大的预报，大多数均变小。这说明实况出现大雨时预报水平较差。
图1 两种预报方法的预报-实测相关图
2．模型II
为了较客观地评定两种预报方法，我们利用模糊数学中的模糊综合评判方法。模糊数学的创立者 L. A.Zadeh 为了描述和处理事物的模糊关系，把“属于”关系进一步数量化，即集合A 中的某个元素ui 对A 不是要么“属于”要么“不属于”关系而是可以不同程度的“属于”和不同程度的“不属于”，这个程度叫做隶属度。隶属度的范围在0 与1 之间,即ui 的隶属度值域是[0 ,1 ]。“属于”关系用函数关系表示,将论域与值域相对应,故形成子集合A 唯一确定的一个映射,它们一一对应。其特点是在众多的“属于”关系的评价指标基础上进行加权平均,得出一个无量纲的综合评价值,然后比较综合评价值的大小,对受到多个因素制约的事物或对象作出一个总的评价,这就是所谓的综合评判问题[6 ,7 ] 。根据所给的条件,给每个对象赋予一个评判指标,称之为模糊评分。
第 天，第 个时段，第 个观测站点的预测值的模糊评分为
（1）
式(1) 中第一项是预报基础分，规定为60 分；第二项为强度(量级)预报的加权分。当预报雨量与实况一致时(即预报与实况误差为0)，该预报评分为100。当预报有误差时,按其误差大小给分，误差越大，分值越低，相反分值越高，预报值越接近于实测值。可以看出,根据误差大小计算的模糊评分,能够很好地表征预报贴近实况的程度,从而较好地检验两种预报方法的预报水平[3]。
为了便于比较，我们给出了在各个时段的模糊评分公式。方法I和II在第 时段的模糊评分为
=1,2,3,4
=1,2,3,4
结果见表2。由表2可以看出，两种方法的模糊评分都在80分以上，预报质量都比较稳定，但在第一、二、四时段，方法I的模糊评分都高于方法II的模糊评分。由此可见，在对雨量预报的准确程度上，方法I高于方法II。
表2 两种方法的预报模糊评分
方法
时段
时段1
时段2
时段3
时段4
预报方法I
84.2593
83.8684
82.3264
80.1272
预报方法II
84.1566
83.6940
82.2567
80.5580
3．模型III
由于91个观测站点的设置是不均匀的，它们较集中地分布在53×47的矩形网格的中央区域内，而面雨量能够更真实地反映平面区域降水的总状况，因此我们用其作为评价这两种方法预报好坏的另一个标准。面雨量是单位面积上的降水量, 实际上为某一特定区域或流域的平均降水状况，它有多种计算方法，如算术平均法、泰森多边形法、逐步订正格点法、三角法、等雨量线法等[4]。
这里我们采用算术平均法计算面雨量：其计算公式如下：
第 种预报方法在第 时段的面雨量 ，
第 时段的实测面雨量
结果见表3。由表3可以看出，无论哪个时段，方法I的面雨量都比方法II更接近实测
值，由此可见，方法I的预报效果好于方法II。
表3 实测及两种方法的面雨量
面雨量
方法I
方法II
实测情况
时段1
1.1167
1.1029
1.4827
时段2
1.5837
1.5498
2.1179
时段3
1.3144
1.3128
1.8007
时段4
1.2437
1.0887
1.6362
五． 对问题2建模及求解
1．模型I
雨量按无雨、小雨、中雨、大雨、暴雨、特大暴雨定义为0～ 6 级，见表4 ;
表4 雨量等级划分表
雨量（mm）
<0.1
0.1-2.5
2.6-6
6.1-12
12.1-25
25.1-60
>60.1
名称
无雨
小雨
中雨
大雨
暴雨
大暴雨
特大
暴雨
雨级
0
1
2
3
4
5
6
在此基础上，我们可将预测雨量与实测雨量转化为相应的雨级，并计算出第 天，第 个时段，第 个观测站点的预测雨级的模糊评分为
（2）
同样我们也可得到方法I和II在第 时段的雨级模糊评分，结果见表5。
表5 两种方法的预报模糊评分（雨级）
方法
时段
第1时段
第2时段
第3时段
第4时段
方法I
86.1192
85.8709
84.5196
82.5788
方法II
85.9831
85.6847
84.4693
82.88772
从表5中我们也发现总体上方法I的预报能力优于方法II。
2．模型II
在评定两种方法的预报质量同时，我们需要考虑公众的感受，由于公众对雨量预报的感受都具有一定的模糊性，可以用{很满意，比较满意，基本满意，不满意}来代替。对此我们有两种方法：
（1） 若预报雨级与实测雨级一致，观众很满意，评分为100；若预报雨级与实测雨级相差一级，观众比较满意，评分为80；若预报雨级与实测雨级相差两级，评分为60；若预报雨级与实测雨级相差两级以上，观众不满意，评分为0。
（2） 采用模糊数学中的隶属函数来处理公众的感受。在隶属函数的选取上，在此这里选用偏大型中升岭形分布函数为隶属函数来处理公众的感受。升岭形函数公式为：
其中 是预报雨级与实测雨级差的绝对值[5]。
用上述两种方法都可以求出公众对第 天，第 个时段，第 个观测站点的预测雨
级的满意度 ，在此我们选取了方法(2)。最后把模糊评分和公众满意度的权重分配分
别确定为0.6和0.4，得到一个新的评价指标 为
=0.6 +0.4
对第 时段关于41天、91个站点求和平均得到预报方法I和II在这个时段的评价指标 ，结果见表6。
表6 两种方法的预报模糊评分（雨级）
方法\时段
第1段
第2段
第3段
第4段
方法I
86.2715
85.2865
84.5318
83.1993
方法II
86.2219
85.2148
84.5096
83.5206
由表6可见，考虑公众感受以后，总体上仍然是方法I优于方法II。
你可以根据自己的水平加以改写
