CN1245638C - 作物性能分析的方法和系统 - Google Patents

Abstract

使用一个线性混合模型评价一个农作物品种广阔区域性能的方法和系统，该模型内含地质统计部分并包括固定效应、随机效应和协方差的参数(J20)。构建一个广阔区域数据库，包括共变量数据、一个或多个农作物品种的试验点空间坐标、试验点所在的地理区域和一个或多个农作物品种的性能特征值(H25，H30)。利用广阔区域数据库中的数据拟合线性混合模型，估计固定效应、随机效应和协方差的参数。使用参数估计结果，可以估计农作物品种的长期期望性能。使用参数估计结果，对于一个指定的时间阶段，可以预测农作物品种的平均性能。

Description

作物性能分析的方法和系统

对于1999年10月15日提交的、标题为“SYSTEM FOR PLANTPERFORMANCE ANALYSIS”、序列号为No.60/159,802的临时专利申请，本申请书要求其优先权。这里专门引用上述公开文件的全文作为参考，没有放弃权利要求。

技术领域

本发明涉及一个根据广阔区域中试验数据评价农作物品种性能的过程。本发明的一个过程根据一个统计混合效应模型，采用空间估计和空间预测来对比品种性能。

背景技术

在一个农作物新品种的培育中，要对该品种和与之竞争的其它品种采集性能数据。这些性能数据包括与指定农作物有关的多种农艺学特性的测量值；例如，对于Zea mays，要测量其粮食产量、粮食湿度和作物倒伏。

评价一个新的农作物品种(下文中称为“品种”)潜在的商业价值时，要将其农艺学性能与其它品种的农艺学性能进行对比。其它对照品种包括，培育该品种的公司中已商业化和未商业化的品种，以及竞争对手公司的商业化品种。注意，对现有的商业化品种也会进行这种同样类型的评价，以确定它们是应当保留在市场上还是应当被培育中的、更新的品种替代。

新品种和对照品种的农艺学性能数据来自多个试验点。这些试验点通常广泛分布于试验所包括品种的适应区域。覆盖这些试验点的适应区域往往很大，其范围达数百平方英里。例如，一个新的Zea mays培育品种的试验区可能从衣阿华州西部到密执安州东部，从威斯康星州中部到伊利诺伊州南部。

由于试验项目中的差异，对于一个指定的品种和它的竞争者，数据往往相当“不平衡”，指的是并非指定组中的所有品种都出现在所有的试验点上。考虑只有一对品种的试验数据，也就是新品种和单一的竞争者，某些试验点两个品种都有，而其它点只有这两个品种之一。

对这些性能数据进行分析，以便确定与对照品种相比，新品种具有足够大性能优势的地理范围，从而证明它应当引入这些地区的市场。在理想情况下，在整个试验区域中，与所有对照品种相比，正在培育的品种都具有显著的性能优势。然而在某些情况下，一个品种可能仅在一个地区具有性能优势，但是它仍然可以在该地区中满足相当大的市场需求。所以，不仅要在整个试验区域将指定品种与其它对照品种相比确定性能特征，而且也要在各个地区内这样做，这一点很重要。

不同地理位置或地区两个品种相对性能的差异来自所谓的“环境影响的基因型”(Sprague&Eberhart，1976)。环境影响的基因型是由不同品种对环境条件的反应差异造成的。这些环境条件可能包括如昼长、温度、土壤湿度和病虫害压力。注意，术语“环境”可以指指定一年的不同位置，或者指定位置的不同年度，或者位置与年度的某种组合。

Bradley等人(1988)介绍了品种性能评价的传统统计分析方法。这些传统方法通常是基于“位置匹配的”数据，也就是，对于一个指定品种相对于一个对照品种，在分析中使用的数据仅仅来自两个品种共同种植的试验点。采用一次配对的t试验，对两个品种没有性能差异的假设进行检验。此外，为了获得指定地区相对性能的有关结论，分析中只使用该地区内试验点的数据。

对于位置匹配的数据采用t试验的传统分析效率很低，至少有五个理由。第一，它没有使用全部数据；它使用的数据仅仅来自两个品种共同种植的试验点。第二，为了进行地区性能对比，它不使用所关注地区以外的邻近区域的数据。第三，它不利用与所关注之性能特征有关的共变量(如灌溉或土系)来帮助解释和预测这些差异。第四，在分析模型中它只使用年内的数据；如果一个模型使用多年度的数据，可以获得更稳健的结论。第五，它是基于一个不正确的假设：一个试验点的观测结果独立于其它点的结果。

有两个概括的理由，说明以上列举的低效为什么限制了广阔区域中的农作物评价。第一是数据的使用效率。要作出结论时，实验人员自然希望使用全部数据。在当前的情况下，这包括一个品种在没有种植其它品种之处的数据，包括与所关注地区类似区域的数据，包括有关共变量的数据，还包括跨年度的数据。统计方法应当争取充分利用所有可用的信息。传统分析受到限制的第二个理由是它基于数据独立的经典假设。品种性能数据几乎毫无例外地不满足这种独立性假设，使统计的结论无效，当数据并不真正表明差异时，往往使人推断出存在品种差异。

现有文献综述：

现有文献的简要综述也强调现有技术中潜在的不足，本发明争取解决这些问题。

对农场地块的实验应用空间统计的领域中，关系最密切的论文之一、可能是最值得关注的工作是由Zimmerman和Harville(1991)完成的。作者认为观测结果是一个随机田的表现，引入了所谓的随机田线性模型(RFLM)。在这个模型中，趋势是由一个平均结构所模拟的，而小尺度的相关性通过空间自相关结构模拟。参数估计是通过一种似然方法进行。通过实际数据分析，作者试图表明他们的模型优于最近邻域分析(NNA)方法。注意，他们的研究专门针对小区域估计的情形，空间相关性的范围限于一个试验点。

在空间预测中使用共变量的情况下，另一篇值得关注的论文是由Gotway和Hartford(1996)完成的。对粮食产量作为一个共变量的数据，通过应用协克里金法预测土壤硝酸盐含量，作者提出在空间预测中使用辅助或次要变量。将他们的方法应用于一个试验点的数据后，他们表明了其方法的效益超过更传统的外表漂移方法。他们研究的范围也是仅限于点上的预测。

涉及多地点产量试验的近期论文之一是由Cullis等人(1998)完成的。在这篇论文中，作者提出了一个多环境初代品种试验中空间分析的方法。这个方法使用最优线性无偏预测(BLUP)求取基因型效应和环境影响的基因型效应，使用REML求取空间参数和方差分量。不过，该方法对每个试验单独模拟协方差结构，也就是不考虑试验之间的相关性。

在广阔的陆地区域将地质统计学应用于土壤化学研究的领域中，Yost，Uehara和Fox(1982a，1982b)是农业科学研究的先驱。在两篇相继发表的论文中，他们论述了在夏威夷岛上进行土壤化学空间预测的研究结果。不过，他们的研究限于只使用克里金方法，而不考虑任何协克里金的方式。此外，在他们的研究中，除了使用选定的土壤化学变量以外，没有使用其它类型的变量。在土壤科学领域值得一提的另一项类似却相对更近期的工作是由Ovalles和Collins(1988)完成的。在整个佛罗里达西北部大约覆盖380公里×100公里的区域中，作者使用通用克里金方法研究选定土壤性质的空间变化，记载的自相关范围大约为40km。在他们的研究中没有试图进行任何空间估计或协克里金。

涉及农作物产量空间预测的论文之一是由Bhatti，Mulla和Frazier(1991)完成的。作者使用一块面积大约为655米×366米的实验田，研究小麦产量与土壤有机质和土壤磷含量的关系。预测产量时使用了克里金和协克里金。然而，研究限于单一的实验田，所以不包含任何模拟大尺度趋势的方式，而该趋势在广阔区域的试验中往往存在。农作物产量空间分析领域的另一篇论文是由Brownie，Bowman和Burton(1993)完成的。在这篇论文中，作者对比了三种可选的空间方法：趋势分析、最近邻域Papadakis分析和相关误差分析，以研究玉米(Zea mays)和大豆(Glycine max)产量数据的空间变化。与农作物产量的其它现有研究的情况相同，他们的研究也是仅在点上的，也就是不考虑各实验点之间的空间相关性。事实上，作者在其论文的结论中也注意到，虽然分析中结合了来自多点的数据，文献中还没有农作物产量数据的跨点空间分析。

有关农作物产量空间分析的其它论文在文献中也是存在的。Bhatti，Mulla，Kooehler和Gurmani(1991)使用了半变量图来识别农作物产量的空间自相关性。他们通过研究应用NNA前后的半变量图，显示了NNA消除空间可变性的效果。他们研究的范围也是仅限于点上分析，空间自相关的最大范围大约是20米。此外，他们的研究中不包括空间估计或预测的方式。Wu等人(1997)就谷类培育试验中消除产量数据的空间变化，对比了变量误差一阶差分(FD-EV)法(Besag和Kempton，1986)与更传统的Papadakis最近邻域法和经典的随机化完整区块(RCB)分析。不过，他们的方法不需要预先指定的趋势和空间自相关结构模型。另外，他们的研究仅限于点上空间变化。

Stroup，Baenziger和Mulitze(1994)使用了若干育种区的数据，对于有效地消除空间可变性产生的干扰，就处理效果比较，对比了传统的RCB分析、NNA和随机场线性模型分析(Zimmerman andHarville，1991)。自然，他们的空间相关性研究限于每个育种区之内的数据。

Gotway和Stroup(1997)的一篇论文的独特之处在于，作者把广义线性模型理论拓展到包括离散和分类的空间变量的空间估计和预测。他们将其拓展应用于两个数据集，一个有关作物虫害，另一个有关杂草计量。不过，与所有其它研究的情况相同，他们研究的范围及其应用限于每个实验点之内的数据。

在以上几段中，回顾了与本发明有关领域中现有的文献。这些领域中研究工作的总体主要可以划分为两个范畴：(1)地质统计学和空间统计学应用于土壤科学领域的大区域和小区域试验，(2)地质统计学和空间统计学应用于农田内小区域试验的农作物反应分析。现有文献缺乏在大区域试验情况下农作物反应分折(如Zea mays的粮食产量)领域的研究，其中空间相关性的范围要扩展到超出单个的实验点。

在农作物品种的性能评价中，根本问题是在多种环境下，也就是覆盖多个试验点的广阔地理范围内得出结论。在文献综述中也应当注意，现有文献的任何研究中都没有涉及空间估计和空间预测这二者的使用。相反，正如下面将要详细讨论的，本发明涉及解决品种评价问题的一种新方法，它不仅是基于多环境、大区域试验，它还具有两个部分来回答两个明显的问题：(1)一个估计部分，回答的问题是一个品种的长期平均性能或两个品种的长期平均性能差异，(2)一个预测部分，回答的问题是在若干指定的时间(年度)点上和/或若干指定的空间(地理)点上，一个品种的性能或两个品种的性能差异。通过通用克里金方法和通用区块克里金方法，本发明背后的方法考虑了大尺度的趋势，而且通过协克里金方法容易使用共变量。现有文献中的工作没有把以上的全部功能结合在一个统一的方法中，以便研究农作物品种的性能。

以上列举的问题并非面面俱到，对于目前已知的作物分析技术中有损于有效性的许多问题，只是涉及了一部分。其它值得关注的问题也可能存在；然而，上述问题应当足以表明现有技术中采用的方法总而言之还是不能令人满意。

发明内容

本发明的若干实施例采用一种称为线性混合模型的统计模型和地质统计方法，由广阔区域的试验数据评价某个农作物品种的性能。评价农作物性能可以通过测量商业上的重要特性，比如(但不限于)产量、粮食湿度或作物倒伏。混合模型能够使用共变量，比如年度、土壤类型、灌溉等等，作为固定效应，而且能够使用随机效应，比如试验位置，这样做有助于解释农作物的性能。固定效应和随机效应不能解释的剩余变化采用地质统计方法模拟。地质统计模型考虑了数据中的空间自相关性，能够获得有效的置信区间，以评价估计和预测中的不确定性。

本发明的实施例具有两个明显的部分：空间估计和空间预测。估计部分可用于：

(1)估计一个品种在一个指定空间位置的平均性能(点估计)，然后使用这些点估计结果构造

(a)一个品种在广阔区域(如多国或多州的集合)的平均性能曲面(跨越若干年度)，以及

(b)两个品种的性能在广阔区域(如多国或多州采集)的平均差异曲面(跨越若干年度)，

(2)估计一个品种在一个规定地理区域比如一个市场区(性能的区块估计)跨越若干年度的平均性能，以及

(3)估计两个品种之间在一个规定地理区域比如一个市场区(相对性能的区块估计)跨越若干年度的平均性能差异。

预测部分可用于：

(1)预测一个品种在一个指定年度一个指定空间位置的平均性能(点预测)，然后使用这些点预测结果构造

(a)一个品种在广阔区域(如多国或多州的集合)每一年度的性能预测曲面，

(b)两个品种在广阔区域(如多国或多州的集合)每一年度的性能差异预测曲面，

(2)预测一个品种在一个规定地理区域比如一个市场区(区块性能预测)每一年度的平均性能，以及

(3)预测两个品种之间在一个规定地理区域比如一个市场区(区块相对性能预测)每一年度的平均性能差异。

决定选用新培育的、与其它品种相比具有最优性能的品种时，需要如上所述的品种性能估计和预测。对于一个候选品种是否进入商业化状态需要作出决定，如果进入，将它置于适当的地区，使之与竞争对手相比表现良好。

从营销的角度看，需要进行一个品种相对性能的定量评价有两个明显的理由，首先，为了向市场引入一个新品种而将该品种推进到商业化状态，其次，为了决定由一个相对性能表现更好的新品种取代一个现有的商业化品种。

以上的决定所依据的最重要的准则之一，是一个品种的期望长期性能优势。评价品种在市场中的商业价值时需要这一点，其中度量一个品种在市场中的价值要跨越它在市场中的整个寿命。通过评价一个品种在一个指定位置上或一个指定地理区域中的长期相对性能，本发明的估计部分提供了这个问题的答案。

在一个指定地点上或一个指定地理区域中，本发明的预测部分可以评价在一个指定年度，一个品种的性能以及品种之间的性能差异。确定一个品种跨越空间和/或时间情况下性能或相对性能的一致性时，需要这种特定空间和/或特定时间的性能评价。对于一个商业化品种来说，性能的一致性——在作物培育术语中称为“品种稳定性”——是一个极为必要的属性。

在广阔区域的试验中，本发明优于传统方法的长处包括以下几点。第一，它并不要求数据是位置匹配的(数据集内的每个位置都包含两个品种的观测结果)，而是使用全部试验点的数据，只要该点至少具有所研究的品种之一。第二，为了进行地区性能评价，不限制它只使用本地区内的数据，而是也使用周围地区的试验点上的数据。第三，它易于加入共变量的信息，这些信息与所关注的主要特性有关(例如当产量是所关注的主要特性时，有关土系的数据)。第四，它使用的模型加入年度作为一个模型因素，因此可以容纳多年度的数据。第五，它并不假设不同试验点的观测结果具有独立性，而是利用不同试验点之间的空间相关性以提供改善的统计结论。

本发明的一个实施例中的步骤可归纳如下：

1.构建农作物在广阔区域中性能的一个数据库，包括试验点的空间坐标、每个试验点上品种性能的测量值和每个试验点上有关共变量的数据。

2.选择两个品种进行对比，例如，一个正准备推进到商业化状态的品种“头领”和一个已经在市场中的品种“竞争者”。

3.通过目测以及采用基于帽子矩阵的统计试验和Cook距离方法，从两个品种的数据中排除异常值。

4.为线性混合模型选择一个空间协方差模型。

5.估计空间协方差参数以及线性混合模型中其它固定效应和随机效应的参数。

6.使用估计的参数

(a)估计跨年度平均值或按年预测，

(b)计算在地理区域中连续曲面或区块平均值，

(c)计算每个品种性能特征或性能特征的差异，以及

(d)求取估计结果或预测结果的标准差。

7.使用估计结果和预测结果及其标准差评价品种性能，例如评价头领品种的相对性能，以决定将其推进到商业化状态。

一方面，本发明是一种在广阔区域中评价农作物品种性能的方法，它使用的线性混合模型包含地质统计部分，还包括固定效应、随机效应和共变量的参数。“广阔区域”仅仅表明一个试验品种假定的区域。例如，“广阔区域”可能是指一个试验品种的适应区域，它可能包括许多试验点。“农作物品种”是指培育一个指定的植物种类以及业内周知的任何其它用途。构建一个广阔区域的数据库，其中包括一个或多个农作物品种的试验点空间坐标以及一个或多个农作物品种的性能特征值。“数据库”是指数据的任何集合。例如，一个“数据库”可能指的是可搜索数据的一个电子形式的集合。“性能特征值”是指与特定品种相关联的任何受关注的农艺学特征。比如“性能特征值”可能指的是业内周知特征的任何数字。例如，它可能指的是比如Zea mays的粮食产量、粮食湿度和作物倒伏。根据广阔区域数据库中的数据拟合线性混合模型，从而估计固定效应、随机效应和共变量的参数。在一个或多个指定空间位置中的每一个，使用参数估计结果，估计农作物品种的长期期望性能。“估计长期期望性能”意为在多个试验期间(如年度)估计一个品种农艺学特征的一个数值，以及业内周知的该词组的正常用途。例如，它可能指的是对于该品种从产生直到它退出市场的整个时间周期，估计一个特征的一个数值，不过它也可以理解为该时间周期可随估计过程的需要而改变，而且它可能比从产生到退出的时间短得多。

在其它方面，数据库进一步可能也包括共变量数据。估计长期期望性能可能包括估计该农作物品种与另一个农作物品种之间的长期期望性能差异。估计这些参数的过程可能包括约束最大似然方法。估计长期期望性能的过程可能包括广义最小平方方法。本方法可能也包括在参数估计之前，使用一种倍率或Cook距离的方法排除数据。本方法可能也包括计算与长期期望性能相关联的一个标准差。本方法可能也包括形成长期期望性能的一个输出。该输出可能包括文本输出。该输出可能包括图形输出。该图形输出可能包括一个等值线图，表示长期期望性能的连续曲面。

另一方面，本发明是一种在广阔区域中评价农作物品种性能的方法，它使用的线性混合模型包含地质统计部分，还包括固定效应、随机效应和共变量的参数。构建一个广阔区域的数据库，其中包括一个或多个农作物品种的试验点空间坐标、试验点所处的地理区域以及一个或多个农作物品种的性能特征值。根据广阔区域数据库中的数据拟合线性混合模型，从而估计固定效应、随机效应和共变量的参数。在一个或多个指定地理区域中的每一个，使用参数估计结果，估计农作物品种的长期期望性能。

另一方面，本发明是一种在广阔区域中评价农作物品种性能的方法，它使用的线性混合模型包含地质统计部分，还包括固定效应、随机效应和共变量的参数。构建一个广阔区域的数据库，其中包括一个或多个农作物品种的试验点空间坐标以及一个或多个农作物品种的性能特征值。根据广阔区域数据库中的数据拟合线性混合模型，从而估计固定效应、随机效应和共变量的参数。在一个或多个指定空间位置中的每一个，在一个或多个指定期间中的每一个，使用参数估计结果，预测农作物品种的平均性能。“预测平均性能”仅仅表明预测一个品种的一个农艺学特征的一个数值以及业内周知的该词组的正常用途。

在其它方面，估计这些参数的过程可能包括约束最大似然方法。数据库可能也包括共变量数据。共变量数据可能包括一个反应变量。预测平均性能的过程可能包括通用协克里金方法。共变量数据可能包括固定效应。预测平均性能的过程可能包括通用克里金方法。预测平均性能的过程可能包括预测该农作物品种与另一个农作物品种之间的平均性能差异。本方法可能也包括在参数估计之前，使用一种倍率或Cook距离的方法排除数据。本方法可能也包括计算与预测平均性能相关联的一个标准差。本方法可能也包括形成预测平均性能的一个输出。

另一方面，本发明是一种在广阔区域中评价农作物品种性能的方法，它使用的线性混合模型包含地质统计部分，还包括固定效应、随机效应和共变量的参数。构建一个广阔区域的数据库，其中包括一个或多个农作物品种的试验点空间坐标、试验点所处的地理区域以及一个或多个农作物品种的性能特征值。根据广阔区域数据库中的数据拟合线性混合模型，从而估计固定效应、随机效应和共变量的参数。在一个或多个指定地理区域中的每一个，在一个或多个指定期间中的每一个，使用参数估计结果，预测农作物品种的平均性能。

另一方面，本发明是一种杂交品种培育方法，培育一个杂交品种，获得该杂交品种和对照杂交品种的性能数据。采用广义最小平方法使一个三次多项式曲面拟合到每个杂交品种的性能数据上，使用球面变量图模拟残差方差。新的杂交品种与对照杂交品种的性能进行对比。

另一方面，本发明是一种系统，包括一台计算机和一种程序。该程序在该计算机上执行并包括以下程序代码：采用广义最小平方法使一个三次多项式曲面拟合到每个杂交品种的性能数据上；使用球面变量图模拟残差方差；以及新的杂交品种与对照杂交品种的性能进行对比。

本公开文件的优点应当理解为，预测、区块预测、估计和点估计能够以任何数目的各种排列进行组合，以获得有价值的性能评价。所有这些组合都不脱离本发明的范围。

附图简要说明

本发明采用举例说明，但不受附图的限制。本专利的文件包含至少一张彩色图。专利和商标事务所在收到请求和必需的费用之后，将会提供本专利的拷贝，包括彩色附图。

图1是一个等值线图，显示1996年两个商业化杂交品种——DK512和DK527——产量差异的预测曲面。图中显示了三条产量差异等值线(-5、0和+5蒲式耳/英亩)和四条显著级别等值线(产量差异±5％和±20％)，也显示了两个杂交品种的试验点以及市场区的边界和名称。

图2是一个等值区域图，显示1996年在若干地理区域中的每一个区域内，两个商业化杂交品种——DK512和DK527——产量差异的区块预测。等值区域图上的图例既表明了产量差异的级别(-3、0和+3蒲式耳/英亩)，又表明了差异是否达到5％的显著级别。

图3是一个等值区域图，显示1996年在若干地理区域中的每一个区域内，两个商业化杂交品种——DK512和DK527——产量差异的简单均值。等值区域图上的图例既表明了产量差异的级别(-3、0和+3蒲式耳/英亩)，又表明了差异是否达到5％的显著级别。

图4、图5、图6表示的信息分别与图1、图2、图3相同，只是时间为1997年。

图7、图8、图9表示的信息分别与图1、图2、图3相同，只是时间为1998年。

图10是一个等值线图，显示1996年、1997年和1998年两个商业化杂交品种——DK512和DK527——产量差异的跨年度估计曲面。图中显示了三条产量差异等值线(-5、0和+5蒲式耳/英亩)和四条显著级别等值线(产量差异±5％和±20％)，也显示了两个杂交品种的试验点以及市场区的边界和名称。

图11是一个等值区域图，显示1996年、1997年和1998年在若干地理区域中的每一个区域内，两个商业化杂交品种——DK512和DK527——产量差异的跨年度区块估计。等值区域图上的图例既表明了产量差异的级别(-3、0和+3蒲式耳/英亩)，又表明了差异是否达到5％的显著级别。

图12是一个等值区域图，显示1996年、1997年和1998年在若干地理区域中的每一个区域内，两个商业化杂交品种——DK512和DK527——产量差异的跨年度简单均值。等值区域图上的图例既表明了产量差异的级别(-3、0和+3蒲式耳/英亩)，又表明了差异是否达到5％的显著级别。

示范实施例的说明

现在将要详细介绍在此公开之发明的优选实施例，它的实例在附图中展示。

本发明涉及以上在背景技术中论述的传统统计分析方法中的所有五个限制。第一，它并不要求位置匹配的数据，而是使用全部试验点的数据，只要该点至少具有试验中的品种之一。第二，为了进行地区性能对比，不限制它只使用本地区内的数据，而是也使用周围地区的试验点上的数据。第三，它易于使用共变量。第四，它使用的模型加入年度作为一个模型因素，因此可以容纳多年度的数据。第五，它并不假设不同试验点的观测结果具有独立性，而是利用不同试验点之间的空间相关性以提供改善的统计结论。

期望品种性能数据具有空间相关性是相当自然的，因为品种具有对环境条件的敏感反应。然而，为了对品种性能数据尤其是产量数据中空间相关的程度作出定量评价，对来自Monsanto(DEKALB)玉米研究机构的Zea mays产量数据进行了一项研究。分析了三年中(1994-1996)包括许多Zea mays品种的历史产量数据。分析结果表明，大约80％品种的数据具有很强的空间自相关性，其证明是拟合的协方差模型中很大的正范围(Bhatti，Mulla，Kooehler and Gurmani，1991)。这一结果受到一个1994年和1995年Zea mays产量数据混合模型分析的支持，那是一项似然比试验，分别对每个品种试验空间自相关的程度。平均说来，75％的品种在其数据中表现出显著的自相关性。

本发明使用一种统计模型，它既能进行空间估计，又能进行空间预测，这使它在品种性能评价中具有超过传统方法的某些明显的优点。空间估计能够评价品种的长期平均性能，也能够评价不同品种的长期平均性能差异。这种信息的价值在于确定一个指定的品种相对于其它竞争品种是否具有某种长期的性能优势。可以在特定的位置确定，也可以在更大的区域中确定。这一点对于评价该品种在市场上的商业价值是必需的，因为度量一个品种的价值是要跨越它在市场上的整个寿命。空间预测能够在指定的时间点(年度)和/或指定的空间点(地理位置)或聚集在更大区域的空间点上，评价品种的性能和性能差异。在确定一个品种跨越空间和/或时间情况下性能的一致性时，这种特定空间和/或特定时间的性能评价是有用的。对于一个商业化品种来说，性能的一致性——在作物培育术语中称为品种稳定性——是一个极为必要的属性。采用品种性能的混合模型分析，可以解决这种稳定性问题(Littell等人，1996，250页)。

本发明在对比品种性能的方法中引入了地质统计模型。地质统计模型使用从一种地貌的不同位置获取的样本点，并由这些样本点产生(插值出)一个连续曲面。这些样本点为某种现象的测量结果，比如空气中的污染物质含量或地质岩芯的矿物级别。据信地质统计学的起源有数个来源。使用协方差的插值——也称为最优线性无偏预测(BLUP)——来自Wold(1938)、Kolmogorov(1941)、Wiener(1949)、Gandin(1959)和Goldberger(1962)。使用变量图的插值归功于Gandin(1959，1963)和Matheron(1962，1969)。Cressie(1990)对地质统计学的起源给出了更多细节。地质统计学最初应用于矿业领域，随后应用的领域诸如气象学和环境科学。近年来它已经应用于农业有关的领域，比如土壤性质和农作物产量的模拟。

混合效应模型：

“混合效应模型”形成了理解本发明的背景，所以应当进行更为详细的讨论，以向读者提供与在此讨论的实施例有关的背景信息。考虑以下普通模型：

y＝Xθ+Dr+ε，                [1]

式中y为反应向量，θ为未知固定效应参数的向量，r为未知随机效应参数的向量，X和D为已知结构的矩阵，ε为未知随机误差分量的向量。因为这种模型具有固定和随机两种效应参数，所以称为混合效应模型，或者更简单地称为“混合模型”。假设r和ε为正态分布，其参数为

$E\left(\begin{array}{c}r\\ \mathrm{\ϵ}\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),---\left[2\right]$

以及

$\mathrm{var}\left(\begin{array}{c}r\\ \mathrm{\ϵ}\end{array}\right)=\left(\begin{array}{cc}G& 0\\ 0& \mathrm{\Σ}\end{array}\right).---\left[3\right]$

在以下段落中，将要说明地质统计模型与这种混合模型的关系。

地质统计学中的混合模型：

地质统计学假设的一个空间线性模型由下式给定

Z(s)＝μ(s)+ε(s)，              [4]

式中向量s包含空间坐标(如经纬度)，表示随机变量Z(s)的空间位置。[4]式由两部分组成：一个确定部分μ(s)和一个随机部分ε(s)。对于两个不同的点s和u，随机变量ε(s)和ε(u)是空间相关的。因为这种相关性存在于同一随机过程ε(·)内的各项之间，在术语中称为“自相关”。为了在统计模拟中进行复制，假设这种空间自相关仅仅依赖于分离s和u的距离和方向，而与它们的确切位置无关。这种性质在术语中称为“平稳性”。

在地质统计学中最普通的模型是[4]式的一种特殊情况，确定曲面μ(s)不随位置s而变化。这种模型由下式给定

Z(s)＝μ+ε(s).                [5]

[4]式和[5]式的模型都可以用于预测，术语“预测”对应于随机量的推断(Cressie，1993，15页)。这里，推断一词是指随机量的预测以及与该预测相关联的不确定性。在[4]式的模型中，确定曲面μ(s)随位置变化，对于某个位置s₀，Z(s₀)的预测在术语中称为“通用克里金”。在[5]式的模型中，模型的确定部分是一个常数，不随位置变化，Z(s₀)的预测在术语中称为“普通克里金”。

[4]式和[5]式的模型也可以用于估计，术语“估计”是指固定或确定量的推断(Cressie，1993，13页)。这里，推断是指固定量的估计以及与该估计相关联的不确定性。注意，在[4]式的模型中估计指的是确定μ(s)，而在[5]式的模型中估计指的是确定μ。

以上关于预测和估计的介绍涉及在若干点上的推断；也就是在特定点上某些数量的预测或估计。预测和估计也可以应用于空间聚集或称为区域的级别；例如在整个州。在这种情况下，在某个区域A中，我们或者预测随机变量Z(s)的平均值，或者估计确定趋势曲面μ(s)的平均值。这两种运算在术语中分别称为“区块预测”Z(A)≡∫_AZ(s)ds/|A|和“区块估计”μ(A)≡∫_Aμ(s)ds/|A|，式中|A|≡∫_A1ds为A的面积。

[4]式的模型可以写为一个普通线性模型，

z＝Xθ+ε                       [6]

表示一个品种所有随机变量的联合期望。这里，随机误差的方差-协方差矩阵给定为var(ε)＝∑。注意，[6]式的模型可视为[1]式模型的一种特殊情况。它可以推广为容纳多变量的情况(例如两个品种)，写为

$\left(\begin{array}{c}{z}_1\\ z_2\end{array}\right)=\left(\begin{array}{cc}{x}_1& 0_1\\ 0_2& x_2\end{array}\right)\left(\begin{array}{c}{\mathrm{\θ}}_1\\ {\mathrm{\θ}}_2\end{array}\right)+\left(\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\end{array}\right),---\left[7\right]$

式中

$\mathrm{var}\left(\mathrm{\ϵ}\right)=\mathrm{var}\left(\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\end{array}\right)\≡\mathrm{\Σ}=\left(\begin{array}{cc}{\mathrm{\Σ}}_{11}& {\mathrm{\Σ}}_{12}\\ {\mathrm{\Σ}}_{12}& {\mathrm{\Σ}}_{22}\end{array}\right),---\left[8\right]$

它仍然是[1]式模型的一种特殊情况。例如，若是第i个品种有n_i个观测结果，i＝1，2，那么[7]式中的向量z_i和ε_i为n_i维的向量，i＝1，2。矩阵0_i，i＝1，2是元素为零的矩阵。0_i和X_i的行数是n_i，i＝1，2，而这两个矩阵的列数和向量θ_i，i＝1，2，的长度取决于模型中使用的参数数目。

在[7]式的模型中，对第i个品种预测Z_i(s₀)，i＝1，2时，这个过程称为简单情况下的协克里金

$\mathrm{X\θ}=\left(\begin{array}{cc}{1}_{n_1}& 0_{n_1}\\ 0_{n_2}& 1_{n_2}\end{array}\right)\left(\begin{array}{c}{\mathrm{\μ}}_1\\ {\mathrm{\μ}}_2\end{array}\right),---\left[9\right]$

式中1_ni和0_ni分别为n_i个1和0的向量，i＝1，2。对于比[7]式的模型更一般的情况，这个过程称为通用协克里金。当联合预测多个品种时——由向量z(s₀)表示，这个过程称为多变量空间预测(Ver Hoef和Cressie，1993)。

协方差的地质统计模型

到目前为止，我们还没有指定随机误差协方差矩阵var(ε)＝∑的结构。确实，使本方法成为“地质统计学”的，是使用空间坐标指定矩阵∑的结构。令一个位置的空间坐标包含在2×1的向量s中。令第i个品种(如第一个品种)的第j个观测结果位于s＝(x_s，y_s)，并令第k个品种(如第二个品种)的第m个观测结果位于u＝(x_u，y_u)。那么第i个品种的第j个观测结果与第k个品种的第m个观测结果之间的协方差可以由一个地质统计协方差模型来模拟。这种协方差模型的一个实例是球面模型(Cressie，1993，61页)，由下式给出

$\mathrm{cov}(Z_{\mathrm{ij}},Z_{\mathrm{km}})=\left\{\begin{array}{ccc}{\mathrm{\φ}}_{\mathrm{ik}}\[1-(\frac{3\left|\right|h\left|\right|}{{2\mathrm{\ρ}}_{\mathrm{ik}}}-\frac{{\left|\right|h\left|\right|}^3}{{2\mathrm{\ρ}}_{\mathrm{ik}}^3})]& \mathrm{for}& \left|\right|h\left|\right|\≤{\mathrm{\ρ}}_{\mathrm{ik}},\\ 0& \mathrm{for}& \left|\right|h\left|\right|>{\mathrm{\ρ}}_{\mathrm{ik}},\end{array}---\left[10\right]\right.$

式中h＝u-s，‖h‖是向量h的长度。本质上‖h‖是从s到u的欧几里得距离，等于 $\sqrt{{(x_s-x_u)}^2+{(y_s-y_u)}^2}.$ 在地质统计学中，参数_ik称为门槛，而ρ_ik称为射程。如果i＝k，[10]式就定义了自相关，否则它就定义了互相关。从物理上说，门槛参数表示方差，而射程表示一个距离，超过后相关性下降至零。

必须注意确保方差-协方差矩阵var(ε)＝∑是正定的。在单一品种的情况下，使用[10]式保证了∑是正定的。例如，注意到对于自相关，每个i的φ_ii和ρ_ii都大于零，足以保证正定性。然而对于互相关，要保证正定性，对φ_ii、φ_kk和φ_ik就有所限制(Goovaerts，1977，108页)。传统方法中使用了一类称为共区域化法的模型，以保证矩阵∑是正定的。Ver Hoef和Barry(1998)已经给出了模拟互相关的其它方式。注意到本实施例中∑₁₂＝∑₂₁＝0，当使用[10]式时，∑₁₁和∑₂₂中每个i的φ_ii和ρ_ii都大于零，显而易见[8]式中的方差-协方差矩阵是正定的。

现在我们要介绍一类共区域化法的模型，作为本发明和本实例的细节。对于一个指定的整数m，m＝1，2，……，M，让我们开始为第i个变量构建一个过程如下，i＝1，……，K，

$Y_{\mathrm{im}}\left(s\right)=\sqrt{1-r_{\mathrm{im}}^2}{W}_{\mathrm{im}}\left(s\right)+r_{\mathrm{im}}{W}_{0m}\left(s\right),---\left[11\right]$

式中-1≤r_im≤1，对于m＝1，2，……，M，W_im是另一个空间过程，对于所有i≠k(i，k＝0，1，……，K)，具有以下协方差性质：

cov[W_im(s)，W_km(s+h)]＝0，

以及对于i＝0，1，……，K，

f_m(h；γ_m)≡cov[W_im(s)，W_im(s+h)]，

在这里，函数f_m(h；γ_m)是某种无尺度(f_m(0；γ_m)＝1)的模型，作为协方差模型是有效的，也就是所有可能的协方差矩阵都是正定的。例如，按照[10]式定义、但是没有尺度参数φ_ik的球面模型，作为协方差模型是有效的，如下式给出

式中γ_m＝ρ_m。对于f_m(h；γ_m)，另一个显而易见的例子是无尺度的“熔核效应”，对于所有h，都有f_m(h；γ_m)＝1。熔核效应通常是由于误差和/或可变性的采样尺度小于观测数据集中最小的试验点距离造成的。大多数地质统计模型对于某种其它的模型都包括一种熔核效应，比如以上[12]式中的球面模型，当h＝0时，这种效应造成原点处的不连续。

使用[11]式和f_m(h；γ_m)的定义，我们有

$\mathrm{cov}[Y_{\mathrm{im}}\left(s\right),Y_{\mathrm{km}}\left(s\right)]=\left\{\begin{array}{cc}{f}_m(h;{\mathrm{\γ}}_m)& \mathrm{ifi}=k,\\ r_{\mathrm{im}}{r}_{\mathrm{km}}{f}_m(h;{\mathrm{\γ}}_m)& \mathrm{ifi}\≠k.\end{array},---\left[13\right]\right.$

下面，对于i＝1，2，……，K，我们通过定义以下求和，构建共区域化法的线性模型，

${\mathrm{\ϵ}}_i\left(s\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\φ}}_{\mathrm{im}}{Y}_{\mathrm{im}}\left(s\right),---\left[14\right]$

注意，[14]式的构建提供了一个有效的空间随机过程。换句话说，对于i＝1，2，……，K；m＝1，2，……，M，只要分量过程Y_im(s)存在，对于i＝1，2，……，K，过程ε_i(s)也就存在。还要注意，对于m＝1，2，……，M，过程Y_im(s)可能具有不同的基本协方差模型(例如球面的、指数的等等)。

[14]式的构建为构建ε_i(s)提供了以下协方差模型(共区域化法的模型)：

$C_{\mathrm{ik}}(h;\mathrm{\γ})\≡\mathrm{cov}[{\mathrm{\ϵ}}_i\left(s\right),{\mathrm{\ϵ}}_k\left(s\right)]=\left\{\begin{array}{cc}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\φ}}_{\mathrm{im}}^2{f}_m(h;{\mathrm{\γ}}_m)& \mathrm{ifi}=k,\\ \underset{m=1}{\overset{M}{\mathrm{\Σ}}}{r}_{\mathrm{im}}{r}_{\mathrm{im}}{\mathrm{\φ}}_{\mathrm{im}}{\mathrm{\φ}}_{\mathrm{km}}{f}_m(h;{\mathrm{\γ}}_m)& \mathrm{ifi}\≠k.\end{array}\right.---\left[15\right]$

注意，当K＝2时，参数r_1m和r_2m并不能各自确定，但是它们的积是可以确定的，我们简单记为r_m≡r_1m r_2m。现在只剩下估计以上[15]式协方差模型中的参数，以及固定效应θ。我们在下一节讨论这些问题。

约束最大似然：

对于[1]式中论述的混合模型，[6]式中的地质统计模型是一种特殊情况，必须使用数据来估计模型的参数。一般说来，参数可以分为两部分：[1]式中给出的固定效应θ，以及[8]式给出的∑中包含的协方差参数。

让我们把∑对参数向量的依赖记为∑()。例如，可能是[12]式中球面协方差模型的参数，或者[15]式中共区域化法模型中的参数。一个普通的估计方法是使用约束(也称为残差)最大似然(REML)方法估计θ和。

我们现在简要介绍REML估计方法，它是由Patterson和Thompson(1971，1974)研究的。假设数据服从一个多变量正态分布N(Xθ，∑())，传统的最大似然估计方法能够用于估计θ和，只要在θ和的参数空间中，使双倍负对数似然，

L(θ，)＝Nlog(2π)+log|∑()|+(y-Xθ)′∑()^-1(y-Xθ)，[16]

同时对于θ和达到最小。在的条件下，θ的一个解由下式给出

它是一个众所周知的、θ的广义最小平方估计量。把[17]式代入[16]式，协方差参数的似然函数变为

众所周知，使[18]式中的L()最小化而得出的协方差参数的最大似然估计是有偏的(Mardia和Marshall，1984)。REML估计是最大似然估计的一种改进，它校正了的最大似然估计结果的偏置。对于最小化

就获得了REML估计量，式中d为X的秩。[19]式需要一个数值解法。关于空间协方差结构的更多REML细节可以参见Cressie(1993，91页)。[19]式最小化的计算细节可以参见Harville(1977)，Zimmerman(1989)和Wolfinger爹人(1994)。虽然研究REML时假设数据服从一个多变量正态分布，Heyde(1994)和Cressie和Lahiri(1996)说明了[19]式也是一个估计方程(Godambe，1991)；所以即使数据并不服从一个多变量正态分布，的REML估计也是无偏的。

参考了以上讨论，现在就可以详细讨论本发明的实施例了。

模型：

对于单一品种或者品种对比，为了在集合的点级别或某种区块级别进行估计和预测，研究了一种统计模型。在品种培育的过程中，品种的对比是最重要的目的，几乎不可能发生实验者一次只考虑一个品种或者拟合单品种模型的情况。注意，如果对杂交品种各自的数据拟合单品种模型，就不可能进行品种对比的统计试验。相反，使用一个双品种或多品种的模型是最普通的情况，而且在品种对比的同时，也能够得出对单一品种的推断(估计结果和预测结果)。所以，详细考虑这种一般模型。此处考虑的模型是一种统计模型，既包含固定效应，又包含随机效应，称为混合效应模型。该模型为

y＝Xθ+Dr+z+ε.                     [20]

[20]式模型中的向量y包含了所关注的主要变量(如粮食产量)，可能还有其它伴随变量。作为这些伴随变量的一个例子，考虑生长季节中的平均降水量。提供这些降水量数据的气象台，处于测量所关注变量的试验点所在大区域。向y引入这些其它变量的观测结果的理由，是使用这些变量与预测和估计中所关注的变量之间的协方差。在预测的情况下，读者可以参考前面有关协克里金的讨论。

在一个实施例中，y可能包含恰好两个品种中一个变量的反应，后面将要如此假设。注意，这个实施例仍然表示向量y的一般形式，因为它能够使用两个品种的该变量反应之间的协方差。在这个实施例中，向量y写为y＝(y₁′，y₂′)′，式中 $y_i^{\′}=(Y_{i1},Y_{i2},...,Y_{{\mathrm{in}}_i}),i=\mathrm{1,2},$ 包括第一个品种的n₁个观测结果和第二个品种的n₂个观测结果。换句话说，一共有来自两个品种的N≡n_l+n₂个观测结果。

[20]式模型中的矩阵X是固定效应的设计矩阵，

$X=\left(\begin{array}{cc}{X}_1& 0\\ 0& X_2\end{array}\right),---\left[21\right]$

式中

$X_i=\left(\begin{array}{cc}{v}_{i1}^{\′}& x_{i1}^{\′}\\ v_{i2}^{\′}& x_{i2}^{\′}\\ \·& \·\\ \·& \·\\ \·& \·\\ v_{m_i}^{\′}& x_{{\mathrm{in}}_i}^{\′}\end{array}\right),i=\mathrm{1,2},---\left[22\right]$

i＝1，2。在[22]式中，列向量V_ij，i＝1，2，j＝1，2，……，n_i包含用于模拟趋势曲面的固定效应。这些向量包括第j个观测结果的空间坐标函数，j＝1，2，……，n_i。例如，第j个观测结果的x坐标可能是是纬度，而y坐标可能是该观测结果的经度。也可以使用其它坐标系统，比如通用横墨卡托(UTM)。如果x坐标和y坐标记为x_j和y_j，i＝1，2，那么对于第i个品种的第j个观测结果，V_ij的实施例包含一个三次多项式内的项，

$v_{\mathrm{if}}={(1,x_j,y_j,x_j^2,y_j^2,x_j,y_j,x_j^3,y_j^3,x_j{y}_j^2,x_j^2{y}_j)}^{\′}.---\left[23\right]$

在[22]式中，对于第i个品种的第j个观测结果，列向量x_ij，i＝1，2，j＝1，2，……，n_i包含其它固定的共变量效应。这些共变量可能包括某些因素，比如有没有灌溉、年度的影响等等。所以固定效应参数向量θ可以写为

$\mathrm{\θ}=\left(\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\β}}_1\\ {\mathrm{\α}}_2\\ {\mathrm{\β}}_2\end{array}\right),---\left[24\right]$

式中子向量α_i是一个a_i×1的固定效应向量，用于多项式v_ij′a_i

形成的、反映第i个品种的空间趋势曲面，β_i是一个b_i×1的固定效应向量，反映第i个品种的共变量，i＝1，2。令θ的维数为d×1，也就是d＝a₁+b₁+a₂+b₂。在[23]式给出的实施例中，对于i＝1，2，α_i是一个10×1的向量。

在[20]式的模型中，r、z和ε是随机效应的向量。向量r包含位置的随机效应。如果有n_c个位置两个品种都有观测结果，那么就只有N-n_c个不同的位置，r就变成(N-n_c)×1的向量。假设E(r)＝0以及

$\mathrm{var}\left(r\right)={\mathrm{\σ}}_0^2{I}_r,---\left[25\right]$

式中I_r为(N-n_c)×(N-n_c)的单位矩阵。在[20]式的模型中，D为一个N×(N-n_c)的关联矩阵，反映位置的随机效应。关联矩阵D的行中所有元素都置为0，每行只有一个1表明一个特定位置与该观测结果的关联(例如，参见Searle et al.，1992，138-139页)。此外，一列中出现两个1表明在一个位置两个品种都有。

在[20]式的模型中，z表示一个空间随机过程。向量z可能包含两个品种的空间相关随机变量；z＝(z₁′，z₂′)′式中 $Z_i^{\′}=(Z_{i1},Z_{i2},...,Z_{{\mathrm{in}}_i}).$ ，i＝1，2。假设E(z)＝0以及

$\mathrm{cov}\left(z\right)=\left(\begin{array}{cc}{\mathrm{\Σ}}_{11}& {\mathrm{\Σ}}_{12}\\ {\mathrm{\Σ}}_{21}& {\mathrm{\Σ}}_{22}\end{array}\right),---\left[26\right]$

具有某种空间自相关和互相关的变量图/协方差模式。令一个位置的空间坐标包含在2×1的向量s中。令第i个品种的第j个观测结果位于s＝(x_s，y_s)，并令第k个品种的第m个观测结果位于u＝(x_u，y_u)。那么第i个品种的第j个观测结果与第k个品种的第m个观测结果之间协方差的一个实施例是一个球面协方差模型(Cressie，1993，61页)，由下式给出

$\mathrm{cov}(Z_{\mathrm{ij}},Z_{\mathrm{km}})=\left\{\begin{array}{ccc}{\mathrm{\φ}}_{\mathrm{ik}}[1-(\frac{3\left|\right|h\left|\right|}{{2\mathrm{\ρ}}_{\mathrm{ik}}}-\frac{{\left|\right|h\left|\right|}^3}{{2\mathrm{\ρ}}_{\mathrm{ik}}^3})]& \mathrm{for}& \left|\right|h\left|\right|\≤{\mathrm{\ρ}}_{\mathrm{ik}},\\ 0& \mathrm{for}& \left|\right|h\left|\right|>{\mathrm{\ρ}}_{\mathrm{ik}},\end{array}---\left[27\right]\right.$

式中h＝u-s，‖h‖是向量h的长度。本质上‖h‖是从s到u的欧几里得距离， $\sqrt{{(x_s-x_u)}^2+{(y_s-y_u)}^2.}$ 在地质统计学中，参数φ_ik称为门槛，而ρ_ik称为射程。如果i＝k，[10]式就定义了自相关，否则它就定义了互相关。

必须注意确保方差-协方差矩阵[26]是正定的。传统方法中使用了一类称为共区域化法的模型，以保证这个矩阵是正定的，但是这些方法相当受限。注意到对于自相关，每个i的φ_ii和ρ_ii都大于零，足以保证正定性。然而对于互相关，要保证正定性，对φ_ii、φ_kk和φ_ik就有所限制(Goovaerts，1977，108页)。Ver Hoef and Barry(1998)已经给出了模拟互相关的其它方式。注意到一个实施例中∑₁₂＝∑₂₁＝0，当使用[10]式时，∑₁₁和∑₂₂中每个i的φ_ii和ρ_ii都大于零，显而易见[26]式中的方差-协方差矩阵是正定的。

在[20]式的模型中，ε是一个包含独立随机误差的向量。向量ε可能包含两个品种的随机变量；ε＝(ε₁′，ε₂′)′式中 ${\mathrm{\ϵ}}_i^{\′}=[{\mathrm{\ϵ}}_{i1},]{\mathrm{\ϵ}}_{i2},...,{\mathrm{\ϵ}}_{{\mathrm{in}}_i]}$ ，i＝1，2。假设E(ε)＝0以及

$\mathrm{cov}\left(\mathrm{\ϵ}\right)\left(\begin{array}{cc}{\mathrm{\σ}}_1^2{I}_1& 0\\ 0& {\mathrm{\σ}}_2^2{I}_2\end{array}\right),---\left[28\right]$

式中I_i为n_i×n_i的单位矩阵，i＝1，2。还假设r、z和ε相互独立。

根据[20]、[25]、[26]和[28]式，以及r、z和ε相互独立的事实，可知反应向量y的方差-协方差矩阵由下式给出

$\mathrm{\Σ}={\mathrm{\σ}}_0^2{\mathrm{DD}}^{\′}+\left(\begin{array}{cc}{\mathrm{\Σ}}_{11}& {\mathrm{\Σ}}_{12}\\ {\mathrm{\Σ}}_{21}& {\mathrm{\Σ}}_{22}\end{array}\right)+\left(\begin{array}{cc}{\mathrm{\σ}}_1^2{I}_1& 0\\ 0& {\mathrm{\σ}}_2^2{I}_2\end{array}\right).---\left[29\right]$

注意，矩阵DD’仅包含0和1，1在对角线上。在DD’中偏离对角线的位置出现1表明在某些位置两个品种一起进行了试验。

干扰诊断：

在拟合[20]式的模型以便得出推断之前，很重要的是定位有干扰的数据点以及评价它们对模型的冲击。在进行模型拟合之前，最好排除这些数据点。注意，对于数据的空间分布进行一次目测，即可排除一个点或一簇点，因为在空间上脱离其它观测结果。在我们的过程中，除了对数据的目测，我们还使用干扰的以下两种目标测度，以甄别异常或有干扰的观测结果。第一种干扰测度称为倍率。对于混合效应模型最一般的情况，Chrestensen et al.(1992)定义了所谓的“帽子”矩阵，

H＝X(X′∑^-1X)^-1X′                   [30]

作为倍率的一种测度。如果h_jj＞2t/N，那么第j个数据点就被视为一个高倍率点，式中h_jj为H的第j个对角线元素，t为H的迹，N为数据点的总数目。在一个实施例中，当∑是一个单位矩阵时，帽子矩阵退化为以下形式(Hoaglin and Welsh，1978)：

H＝X(X′X)^-1X′.                            [31]

在这种情况下，如果h_jj＞2d/N，那么第j个数据点就被视为一个高倍率点(Montgomery and Peck，1992，182页)，正如前面的定义，式中d为θ中参数的数目。在我们的过程中使用的第二种干扰测度为对于第j个观测结果定义的Cook距离D_j(Cook，1977)。对于混合效应模型最一般的情况，这个距离由下式给出，

$D_j=\frac{{({\widehat{\mathrm{\θ}}}_{\left(j\right)}-\widehat{\mathrm{\θ}})}^{\′}{X}^{\′}{\mathrm{\Σ}}^{-1}X({\widehat{\mathrm{\θ}}}_{\left(j\right)}-\widehat{\mathrm{\θ}})}{d},---\left[32\right]$

式中 为排除第j个观测结果后θ的一个估计结果，∑为估计的方差-协方差矩阵。D_i＞1的观测结果通常视为有干扰的(Montgomery and Peck，1992，182页)。在一个实施例中，当∑是一个单位矩阵时，为了甄别有干扰的观测结果，采用以下表达式计算Cook距离：

$D_j=\frac{{({\widehat{\mathrm{\θ}}}_{\left(j\right)}-\widehat{\mathrm{\θ}})}^{\′}{X}^{\′}X({\widehat{\mathrm{\θ}}}_{\left(j\right)}-\widehat{\mathrm{\θ}})}{d{\widehat{\mathrm{\σ}}}^2},---\left[33\right]$

如前面一样，D_j＞1的观测结果视为有干扰的。在我们的过程中为了进行干扰诊断，我们使用[31]式和[33]式，以及一个没有任何共变量效应的模型。

为了筛选出有干扰的观测结果，既使用H又使用D_j的理由如下。H的计算完全根据数据点的空间位置(而不是数值)。所以H对空间的异常敏感，可以用于筛选空间异常。相反，Cook距离D_j则是部分根据空间位置，部分根据观测值。所以，D_j可以甄别的数据点具有异常的高低值但是不必是空间异常。

混合模型参数的估计：

介绍了[20]式给出的数据模型之后，下一步应当涉及参数估计的问题。一般说来，这些参数可以划分为两个范畴：[24]式给出的固定效应θ和[29]式给出的∑中包含的协方差参数。让我们把∑对于参数

的依赖记为∑()。在一个优选实施例中，约束最大似然(REML)方法用于估计参数θ和。采用REML估计出∑之后，利用广义最小平方估计θ，由Searle(1971，87页)给出如下

$\widehat{\mathrm{\θ}}={\left(X^{\′}{\mathrm{\Σ}}^{-1}X\right)}^{-1}{X}^{\′}{\mathrm{\Σ}}^{-1}y,---\left[34\right]$

正如前面的定义，式中y为全部选定数据点的向量。

参数估计结果用于两个不同的过程：估计和预测。在两个过程中，所关注的都是或者在单一的位置，或者是一个地理区域的平均。前者称为点估计(预测)，后者称为区块估计(预测)。注意，点估计结果也用于构建估计趋势曲面，也就是对于覆盖着数据的网格，在每一个网格点上进行点估计。同样，对于覆盖着数据空间范围的网格，在所有网格点上进行点预测，就可以构建预测曲面。

在以下两个小节中，我们将要考虑通过拟合[20]式的模型以进行估计和预测。注意，在各个区域进行区块估计或区块预测之前，使用所有选定的数据点(原始数据集减去有干扰的观测结果)拟合[20]式的模型只进行了一次。由这个模型估计出的参数然后用于各个点和区域均值的估计和预测。

区块和点估计：

假定目标是估计固定空间效应趋势曲面的均值。对于一个位置s＝(x，y)’和第i个品种，i＝1，2，让我们定义一个空间坐标的向量值函数

$w_i\left(s\right)={(w_{i1}\left(s\right),w_{i2}\left(s\right),......w_{i{a}_i}\left(s\right))}^{\′}.---\left[35\right]$

注意，在[23]式给出的实施例中，我们有a₁＝a₂＝10和w_i(s)＝(1，x，y，x²，y²，xy，x³，y³，xy²，x²y)′。所以，在某个位置s的趋势曲面值由w_i′(s)α_i给出。不过，往往受到关注的是某个地理区域A的趋势曲面平均值。如上所述，估计A上的这个平均值称为区块估计，对于i＝1，2，区块估计结果由下式给出

${\mathrm{\μ}}_i\left(A\right)=x_{i0}^{\′}{\mathrm{\β}}_i+\frac{1}{\left|A\right|}{\∫}_A{w}_i^{\′}\left(s\right){\mathrm{\α}}_i\mathrm{ds},---\left[36\right]$

式中|A|为区域A的面积，对于i＝1，2，向量x₁₀包含一组对应于区域A的共变量数值。假设共变量x_i0与趋势曲面没有相互作用；指定x_i0的数值仅仅在相应的区域内使趋势曲面升降而不影响其形状。

在[20]式中模型的一个实施例中，因为某些理由，比如没有共变量数据或缺少一个适当的共变量，[22]式中的共变量效应x_ij可以完全从模型中消除。在这种情况下，对应于共变量的x_i0′B_i，项就会从[36]式中消失。注意，从[36]式的积分中提出x_i0′β_i项有助于区分共变量与趋势曲面分量。这些共变量可能是空间的，也可能是非空间的，可能是观测的，也可能是用户指定的。它们的数值可能是连续的，也可能是离散的。这些共变量类型的某些实例是有规律的。作为离散非空间共变量的一个实例，考虑年度效应，它不是空间可积的。当一个共变量是空间的，并且数值是连续的，向量x_i0中的相应分量就表示一个积分平均结果，与[36]式中的积分相似。对位于s的第i个品种，一个空间共变量的数值记为g_i(s)，并令此共变量对应于[22]式中向量x_ij的第k个分量。那么x_i0的第k个分量会是 $\frac{1}{\left|A\right|}{\∫}_A{g}_i\left(u\right)\mathrm{du},$ ，式中g_i(u)假设是黎曼可积的。作为空间共变量观测值的一个实例，考虑土系。在所有点上此共变量都必须已知。幸运的是所有主要的农业区域都有土系地图，并对任何区域(例如县)都按多种土系进行了分类。这意味着在所有可能的空间位置，都能获得作为一个共变量的土系数据。作为用户指定的空间共变量的一个实例，考虑灌溉，它可以作为二进制值函数(1≡灌溉的，2≡非灌溉的)对待。例如，用户可能关注在100％灌溉情况下的一个区域中，对两个品种之间的差异进行估计(预测)。如果作为一个共变量，灌溉对应于向量x_ij的第l个分量，那么这就意味着x_i0的第l个分量指定的数值为1表示灌溉条件。

在一个实施例中有恰好两个品种的一个反应变量，对于区域A中两个品种的区块平均，让我们将其差异定义为μ₀(A)。那么μ₀(A)由下式给出，

${\mathrm{\μ}}_0\left(A\right)\≡{\mathrm{\μ}}_1\left(A\right)-{\mathrm{\μ}}_2\left(A\right)=x_{10}^{\′}{\mathrm{\β}}_1-x_{20}^{\′}{\mathrm{\β}}_2+\frac{1}{\left|A\right|}{\∫}_A[w_1^{\′}\left(s\right){\mathrm{\α}}_1-w_2^{\′}\left(s\right){\mathrm{\α}}_2]\mathrm{ds}.---\left[37\right]$

现在，让我们定义一个集合S_A＝{p：sp∈A}，其中的点{s_p}或者是位于覆盖在A上的规则网格上，或者是在A之内随机选取的，S_A中有n_A个元素。让我们再定义

${\mathrm{\λ}}_1\left(s_p\right)\≡\left(\begin{array}{c}{w}_1\left(s_p\right)\\ x_{10}\\ 0\\ 0\end{array}\right),{\mathrm{\λ}}_2\left(s_p\right)\≡\left(\begin{array}{c}0\\ 0\\ w_2\left(s_p\right)\\ x_{20}\end{array}\right),\mathrm{and}{\mathrm{\λ}}_0\left(s_p\right)\≡\left(\begin{array}{c}{w}_1\left(s_p\right)\\ x_{10}\\ -w_2\left(s_p\right)\\ -x_{20}\end{array}\right).---\left[38\right]$

那么，近似[36]式和[37]式中的积分，并使用[24]式中θ的定义，我们有

${\mathrm{\μ}}_i\left(A\right)\≈x_{10}^{\′}{\mathrm{\β}}_i+\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{w}_i^{\′}\left(s_p\right){\mathrm{\α}}_1=\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{\left[{\mathrm{\λ}}_i\left(s_p\right)\right]}^{\′}\mathrm{\θ},---\left[39\right]$

i＝1，2，以及

${\mathrm{\μ}}_0\left(A\right)\≈x_{10}^{\′}{\mathrm{\β}}_1-x_{20}^{\′}{\mathrm{\β}}_2+\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}[w_1^{\′}\left(s_p\right){\mathrm{\α}}_1-w_2^{\′}\left(s_p\right){\mathrm{\α}}_2]=\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{\left[{\mathrm{\λ}}_0\left(s_p\right)\right]}^{\′}\mathrm{\θ}.---\left[40\right]$

在[39]式和[40]式中代入[34]式中θ的广义最小平方估计量，我们就得到了μ_i(A)的估计量

${\widehat{\mathrm{\μ}}}_1\left(A\right)\≈\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{\left[{\mathrm{\λ}}_1\left(s_p\right)\right]}^{\′}{\left(X^{\′}{\mathrm{\Σ}}^{-1}X\right)}^{-1}{X}^{\′}{\mathrm{\Σ}}^{-1}y,---\left[41\right]$

i＝0，1，2。

的方差为

$\mathrm{var}\left({\widehat{\mathrm{\μ}}}_1\left(A\right)\right)\≈\frac{1}{n_A^2}\underset{{p\∈S}_A}{\mathrm{\Σ}}\underset{{q\∈S}_A}{\mathrm{\Σ}}{\left[{\mathrm{\λ}}_i\left(s_p\right)\right]}^{\′}{\left(X^{\′}{\mathrm{\Σ}}^{-1}X\right)}^{-1}\left[{\mathrm{\λ}}_i\left(s_q\right)\right],---\left[42\right]$

其标准差为 $\mathrm{se}\left({\widehat{\mathrm{\μ}}}_i\left(A\right)\right)=\sqrt{\mathrm{var}\left[{\widehat{\mathrm{\μ}}}_i\left(A\right)\right]}$ i＝0，1，2。

现在让我们考虑点估计的情况。目的是估计以下函数，它是[36]式的空间情况，其中区域A现在退化为一个位置(点)。在s₀的点估计为

μ_i(s₀)＝x_i0′β_i+w_i′(s₀)α_i.          [43]

正如区块估计时的情况，共变量x_i0的数值仅仅通过使点估计结果或增或减，在[43]式中引入了一个附加的效应。然而，由于这种附加效应的幅度和方向在不同的点可能有变化，模型中有了这些共变量之后，采用点估计结果构建的整个趋势曲面的形状可能发生实质性的变化。

作为[41]式中给出的区块估计量的一种特殊情况，μ_i(s₀)的估计量可以由下式获得

${\widehat{\mathrm{\μ}}}_i\left(s_0\right)={\left[{\mathrm{\λ}}_i\left(s_0\right)\right]}^{\′}{\left(X^{\′}{\mathrm{\Σ}}^{-1}X\right)}^{-1}{X}^{\′}{\mathrm{\Σ}}^{-1}y,---\left[44\right]$

i＝0，1，2。类似于[42]式中的区块估计结果，可以获得 标的估计方差，由下式给出，

$\mathrm{var}\left[{\widehat{\mathrm{\μ}}}_i\left(s_0\right)\right]={\left[{\mathrm{\λ}}_i\left(s_0\right)\right]}^{\′}{\left(X^{\′}{\mathrm{\Σ}}^{-1}X\right)}^{-1}\left[{\mathrm{\λ}}_i\left(s_0\right)\right].---\left[45\right]$

点估计结果的标准差为 $\mathrm{se}\left[{\widehat{\mathrm{\μ}}}_i\left(s_0\right)\right]=\sqrt{\mathrm{var}\left[{\widehat{\mathrm{\μ}}}_i\left(s_0\right)\right]}$ i＝0，1，2。

通用区块和点协克里金：

让我们首先考虑区块预测，它也称为区块协克里金，正如在本发明的背景技术中的讨论。这里的目标是在某个指定的地理区域，估计反应曲面的平均值，它是以下三个分量的和：固定的空间效应的趋势曲面、固定的协方差效应和自相关的随机曲面。定义随机函数Z_i(s)，其中集合{Z_i(s)；s∈A}；i＝1，2形成了一个随机曲面，可以进行Cressie(1993，107页)定义的积分。令Z₀(s)≡Z₁(s)-Z₂(s)为一个随机函数，其中{Z₀(s)；s∈A}形成了一个随机曲面，如果Z₁(s)和Z₂(s)可积，它就可积。现在，对于指定的、x_i0，i＝1，2给出的共变量数值集合，让我们定义函数

$S_i\left(A\right)\≡x_{i0}^{\′}{\mathrm{\β}}_i+\frac{1}{\left|A\right|}{\∫}_A[w_i^{\′}\left(s\right){\mathrm{\α}}_i+Z_i\left(s\right)]\mathrm{ds},---\left[46\right]$

i＝1，2，为第i个品种在区域A上平均实现值。正如区块估计时的情况(参见[36]式后的讨论)，[46]式中的x_i0′β_i项可以容纳所有类型的共变量，并表示对应于共变量x_i0的一个平均的结果。

这里的目的是根据数据的某种线性组合，比如说ω_i′y预测区块均值S₁(A)、S₂(A)和S₀(A)≡S₁(A)-S₂(A)。向量ω_i，i＝0，1，2的估计是由一种称为最优线性无偏预测的方法完成的，预测量是数据的一个线性函数，是无偏的，并且在所有线性无偏预测量中具有最小的均方预测误差。

让我们首先假设[20]式给出的线性模型。为了预测区块均值S_i(A)，i＝0，1，2，通过最小化

E[S_i(A)-ω_i′y]²，                 [47]

可以获得最优线性无偏预测量(BLUP)，此时需要无偏条件

E[S_i(A)]＝E(ω_i′y)，              [48]

i＝0，1，2。

现在，使用[20]式的模型和[46]式，i＝1，2时，[48]式可以写为

$\frac{1}{\left|A\right|}{\∫}_A{w}_i^{\′}\left(s\right){\mathrm{\α}}_i\mathrm{ds}+x_{0i}^{\′}{\mathrm{\β}}_i={\mathrm{\ω}}_i^{\′}\mathrm{X\θ},---\left[49\right]$

而i＝0时，[48]式可以写为

$\frac{1}{\left|A\right|}{\∫}_A{w}_1^{\′}\left(s\right){\mathrm{\α}}_1\mathrm{ds}+x_{01}^{\′}{\mathrm{\β}}_1-\frac{1}{\left|A\right|}{\∫}_A{w}_2^{\′}\left(s\right){\mathrm{\α}}_2\mathrm{ds}-x_{02}^{\′}{\mathrm{\β}}_2={\mathrm{\ω}}_0^{\′}\mathrm{X\θ}.---\left[50\right]$

一般说来，使用[35]式中定义的向量值函数w_i(s)，i＝0，1，2，[49]式和[50]式可以结合成单一的方程

u_i′θ＝ω_i′Xθ，                         [51]

i＝0，1，2，式中

$u_1\≡\left(\begin{array}{c}\frac{1}{\left|A\right|}{\∫}_A{w}_{11}\left(s\right)\mathrm{ds}\\ .\\ .\\ .\\ \frac{1}{\left|A\right|}{\∫}_A{w}_{{1a}_1}\left(s\right)\mathrm{ds}\\ x_{01}\\ 0\\ .\\ .\\ .\\ 0\\ 0\end{array}\right),u_2\≡\left(\begin{array}{c}0\\ .\\ .\\ .\\ 0\\ 0\\ \frac{1}{\left|A\right|}{\∫}_A{w}_{21}\left(s\right)\mathrm{ds}\\ .\\ .\\ .\\ \frac{1}{\left|A\right|}{\∫}_A{w}_{{2\mathrm{a\α}}_2}\left(s\right)\mathrm{ds}\\ x_{02}\end{array}\right),\mathrm{and}{u}_0\≡\left(\begin{array}{c}\frac{1}{\left|A\right|}{\∫}_A{w}_{11}\left(s\right)\mathrm{ds}\\ .\\ .\\ .\\ \frac{1}{\left|A\right|}{\∫}_A{w}_{{1a}_1}\left(s\right)\mathrm{ds}\\ x_{01}\\ -\frac{1}{\left|A\right|}{\∫}_A{w}_{21}\left(s\right)\mathrm{ds}\\ .\\ .\\ .\\ -\frac{1}{\left|A\right|}{\∫}_A{w}_{{2a}_2}\left(s\right)\mathrm{ds}\\ -x_{02}\end{array}\right).---\left[52\right]$

为了满足无偏的条件，对于每一个可能的θ，[51]式都必须成立，因此有

u_i′＝ω_i′X，                             [53]

i＝0，1，2。

[52]式中的向量可以近似为

$u_1\≈\left(\begin{array}{c}\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{w}_{11}\left(s_p\right)\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{w}_{a11}\left(s_p\right)\\ x_{01}\\ 0\\ .\\ .\\ .\\ 0\\ 0\end{array}\right),u_2\≈\begin{array}{c}\left(\begin{array}{c}0\\ .\\ .\\ .\\ 0\\ 0\\ \frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{w}_{21}\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{w}_{2a_2}\left(s_p\right)\\ x_{02}\end{array}\left(s_p\right)\right)\end{array},\mathrm{and}{u}_0\≈\left(\begin{array}{c}\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{w}_{11}\left(s_p\right)\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{w}_{{1\mathrm{a\α}}_1}\left(s_p\right)\\ x_{01}\\ -\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{w}_{21}\left(s_p\right)\\ .\\ .\\ .\\ -\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{w}_{{2\mathrm{a\α}}_2}\left(s_p\right)\\ -x_{02}\end{array}\right),---\left[54\right]$

式中S_A已经在前面定义过。

现在，使用[20]式、[46]式和[52]式，均方预测误差[47]式可以写为

$E{[S_i\left(A\right)-{\mathrm{\ω}}_i^{\′}y]}^2=E{[u_i^{\′}\mathrm{\θ}+\frac{1}{\left|A\right|}{\∫}_A{Z}_i\left(s\right)\mathrm{ds}-{\mathrm{\ω}}_i^{\′}(\mathrm{X\θ}+\mathrm{Dr}+z+\mathrm{\ϵ})]}^2---\left[55\right]$

i＝0，1，2。

让我们定义N向量 $\mathrm{\δ}=({\mathrm{\δ}}_{11},{\mathrm{\δ}}_{12},......,{\mathrm{\δ}}_{1n_1},{\mathrm{\δ}}_{21},{\mathrm{\δ}}_{22},.....{\mathrm{\δ}}_{2n_2})$ δ＝Dr+z+ε，并令 $Z_i\left(A\right)\≡\frac{1}{\left|A\right|}{\∫}_A{Z}_i\left(s\right)\mathrm{ds},$ i＝0，1，2。然后，使用[51]式，[55]式可以写为

E[S_i(A)-ω_i′y]²＝E[Z_i(A)-ω_i′δ]²，         [56]

i＝0，1，2。

求取[56]式中的数学期望需要求取cov[Z_i(A)，δ]，它可以近似为

$\left(\begin{array}{c}\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}\mathrm{cov}[Z_i\left(s_p\right),{\mathrm{\δ}}_{11}]\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}\mathrm{cov}[Z_i\left(s_p\right),{\mathrm{\δ}}_{{1n}_1}]\\ \frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}\mathrm{cov}[Z_i\left(s_p\right),{\mathrm{\δ}}_{21}]\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}\mathrm{cov}[Z_i\left(s_p\right),{\mathrm{\δ}}_{{2n}_2}]\end{array}\right),---\left[57\right]$

i＝0，1，2。

现在，由于cov(z，ε)＝0和cov(z，r)＝0，[57]式中的向量可以写为

$c_1\≡\mathrm{cov}[Z_1\left(A\right),\mathrm{\δ}]\≈\left(\begin{array}{c}\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}\mathrm{cov}[Z_i\left(s_p\right),Z_1\left(s_1\right)]\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}\mathrm{cov}[Z_i\left(s_p\right),Z_1\left(s_{n_1}\right)]\\ \frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}\mathrm{cov}[Z_i\left(s_p\right),Z_2\left(s_1\right)]\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}\mathrm{cov}[Z_i\left(s_p\right),Z_2\left(s_{n_2}\right)]\end{array}\right),---\left[58\right]$

i＝0，1，2

求取[56]式中的数学期望还需要求取C_i(A，A)≡var[Z_i(A)]，它可以近似为

$C_i(A,A)\≈\frac{1}{n_A^2}\underset{{p\∈S}_A}{\mathrm{\Σ}}\underset{{q\∈S}_A}{\mathrm{\Σ}}\mathrm{cov}[Z_i\left(s_p\right),Z_i\left(s_q\right)],---\left[59\right]$

i＝0，1，2。现在，由于 $Z_i\left(A\right)\≡\frac{1}{\left|A\right|}{\∫}_A{Z}_i\left(s\right)\mathrm{ds},$ 它可以近似为 $\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{Z}_i\left(s_p\right),$ 我们注意到

$E\left[Z_i^2\left(A\right)\right]=E{\left[\frac{1}{n_A}\underset{{p\∈S}_A}{\mathrm{\Σ}}{Z}_i\left(s_p\right)\right]}^2$

$=\frac{1}{n_A^2}\underset{{p\∈S}_A}{\mathrm{\Σ}}\underset{{q\∈S}_A}{\mathrm{\Σ}}E\left[Z_i\left(s_p\right)Z\left(s_q\right)\right]$

$\frac{1}{n_A^2}\underset{{p\∈S}_A}{\mathrm{\Σ}}\underset{{q\∈S}_A}{\mathrm{\Σ}}\mathrm{cov}[Z_i\left(s_p\right),Z_i\left(s_q\right)]$

$=C_i(A,A)$

因为对于任何p∈S_A以及i＝0，1，2，E[Z_i(s_p)]＝0。所以，[56]式中的均方预测误差可以写为

$E{[Z_i\left(A\right)-{\mathrm{\ω}}_i^{\′}\mathrm{\δ}]}^2=E[Z_i^2\left(A\right)-2Z_i\left(A\right){\mathrm{\ω}}_i^{\′}\mathrm{\δ}+{\mathrm{\ω}}_i^{\′}\mathrm{\δ}{\mathrm{\δ}}^{\′}{\mathrm{\ω}}_i]$

$=E\left[Z_i^2\left(A\right)\right]-2{\mathrm{\ω}}_i^{\′}E\left[Z_i\left(A\right)\mathrm{\δ}\right]+{\mathrm{\ω}}_i^{\′}E\left(\mathrm{\δ}{\mathrm{\δ}}^{\′}\right){\mathrm{\ω}}_i$

$=C_i(A,A)-2{\mathrm{\ω}}_i^{\′}{c}_i+{\mathrm{\ω}}_i^{\′}\mathrm{\Σ}{\mathrm{\ω}}_i,$

[60]

i＝0，1，2。现在，为了满足[53]式中给出的无偏条件，增加一个拉格朗日乘子m_i，[60]式可以写为

L(ω_i，m_i)≡C_i(A，A)-2ω_i′c_i+ω_i′∑ω_i-2(u_i′-ω_i′X)m_i，   [61]

i＝0，1，2。[61]式对ω_i和m_i求偏导，有

$\frac{\∂L({\mathrm{\ω}}_i,m_i)}{\∂{\mathrm{\ω}}_i^{\′}}=-2c_i+2\mathrm{\Σ}{\mathrm{\ω}}_i+2X{m}_i,---\left[62\right]$

以及

$\frac{\∂L({\mathrm{\ω}}_i,m_i)}{\∂m_i^{\′}}=-2(u_i-X^{\′}{\mathrm{\ω}}_i),---\left[63\right]$

i＝0，1，2。

[62]式和[63]式中的这两个偏导数置为0，然后联立求解，使[60]式中的均方预测误差在[53]式中给出的无偏条件下最小化，并产生区域A中品种性能(或性能差异)的最优线性无偏预测量(BLUP)如下。区块预测(协克里金)方程为

$\left(\begin{array}{cc}\mathrm{\Σ}& X\\ X^{\′}& 0\end{array}\right)\left(\begin{array}{c}{\mathrm{\ω}}_i\\ m_i\end{array}\right)=\left(\begin{array}{c}{c}_i\\ u_i\end{array}\right),---\left[64\right]$

i＝0，1，2。对方程[64]求解ω_i由下式给出，

    ∑ω_i+Xm_i＝c_i

    ω_i＝∑^-1(c_i-Xm_i)，                          [65]

i＝0，1，2。所以，区域A中品种性能(或性能差异)的预测量由下式给出，

${\hat{S}}_i\left(A\right)={\mathrm{\ω}}_i^{\′}y,---\left[66\right]$

i＝0，1，2。把从方程[65]解出的ω_i代入[60]式的表达式中，就提供了均方预测误差或称为预测方差的表达式如下，

$\mathrm{var}\left[{\hat{S}}_i\left(A\right)\right]=C_i(A,A)-2{\mathrm{\ω}}_i^{\′}{c}_i+{\mathrm{\ω}}_i^{\′}\mathrm{\Σ}\left[{\mathrm{\Σ}}^{-1}(c_i-X_{m_i})\right]$

$=C_i(A,A)-2{\mathrm{\ω}}_i^{\′}{c}_i+{\mathrm{\ω}}_i^{\′}{c}_i-{\mathrm{\ω}}_i^{\′}X{m}_i$

$=C_i(A,A)-{\mathrm{\ω}}_i^{\′}{c}_i-{\mathrm{\ω}}_i^{\′}X{m}_i$

$=C_i(A,A)-{\mathrm{\ω}}_i^{\′}(c_i+X{m}_i),$

[67]

i＝0，1，2。把从方程[65]解出的ω_i代入从矩阵方程[64]式下半部获得的方程X’ω_i＝u_i中，就提供了拉格朗日乘子m_i的表达式如下，

X′ω_i＝u_i

X′∑^-1(c_i-Xm_i)＝u_i

X′^-1c_i-X′∑^-1Xm_i＝u_i                [68]

m_i＝(X′∑^-1X)^-1(X′∑^-1c_i-u_i)，

i＝0，1，2。

我们现在考虑点预测或称为点协克里金的情况。注意，方程[46]的一种特殊情况是区域A是单一的点，比如说s₀，也就是，S_i(s₀)≡x_i0′β_i+w_i′(S₀)α_i+Z_i(s₀)，i＝1，2。现在的目标是根据数据的某种线性组合，比如说ω_i*′y，在点s₀进行点预测，也就是，预测S₁(s₀)、S₂(s₀)和S₀(s₀)≡S₁(s₀)-S₂(s₀)。换句话说，点预测量由下式给出，

${\hat{S}}_i\left(s_0\right)={{\mathrm{\ω}}_i^{*}}^{\′}y,---\left[69\right]$

i＝0，1，2，式中协克里金权向量，ω_i ^*需要确定。与方程[65]类似，权向量ω_i ^*由下式给出，

${\mathrm{\ω}}_i^{*}={\mathrm{\Σ}}^{-1}(c_i^{*}-X{m}_i^{*}),---\left[70\right]$

式中向量c_i ^*类似于为了[58]式中的区块预测而定义的向量c_i，拉格朗日乘子向量m_i ^*类似于为了[68]式中的区块预测而定义的向量m_i。所以向量m_i ^*由下式给出，

$c_i^{*}\≡\mathrm{cov}[Z_i\left(s_0\right),\mathrm{\δ}]\≈\left(\begin{array}{c}\mathrm{cov}[Z_i\left(s_0\right),Z_1\left(s_1\right)]\\ .\\ .\\ .\\ \mathrm{cov}[Z_i\left(s_0\right),Z_1\left(s_{n_1}\right)]\\ \mathrm{cov}[Z_i\left(s_0\right),Z_2\left(s_1\right)]\\ .\\ .\\ .\\ \mathrm{cov}[Z_i\left(s_0\right),Z_2\left(s_{n_2}\right)]\end{array}\right),---\left[71\right]$

以及

$m_i^{*}={\left(X^{\′}{\mathrm{\Σ}}^{-1}X\right)}^{-1}(X^{\′}{\mathrm{\Σ}}^{-1}{c}_i^{*}-u_i^{*}),---\left[72\right]$

i＝0，1，2。在方程[72]中，向量u_i ^*，i＝0，1，2类似于为了区块预测的情况而在方程[52]中给出的向量u_i，i＝0，1，2，现在定义为

$u_1^{*}\≡\left(\begin{array}{c}{w}_{11}\left(s\right)\\ .\\ .\\ .\\ w_{{1\mathrm{a\α}}_1}\left(s\right)\\ x_{01}\\ 0\\ .\\ .\\ .\\ 0\\ 0\end{array}\right),u_2^{*}\≡\left(\begin{array}{c}0\\ .\\ .\\ .\\ 0\\ 0\\ w_{21}\left(s\right)\\ .\\ .\\ .\\ w_{{2\mathrm{a\α}}_2}\left(s\right)\\ x_{02}\end{array}\right),\mathrm{and}{u}_0^{*}\≡\begin{array}{c}\left(\begin{array}{c}{w}_{11}\left(s\right)\\ .\\ .\\ .\\ w_{{1\mathrm{\α}}_1}\left(s\right)\\ x_{01}\\ -w_{21}\left(s\right)\\ .\\ .\\ .\\ -w_{{2a}_2}\left(s\right)\\ -x_{02}\end{array}\right).---\left[73\right]\end{array}$

最后，正如在方程[67]中一样，获得点预测量

的预测方差为

$\mathrm{var}\left[{\hat{S}}_i\left(s_0\right)\right]=C_i(s_0,s_0)-{{\mathrm{\ω}}_i^{*}}^{\′}(c_i^{*}+X{m}_i^{*}).---\left[74\right]$

以下包括的实例是为了说明本发明的特定实施例。本领域的技术人员应当理解，在以下的实例中公开的技术是发明人发现的，在本发明的实践中表现良好，因此可视为构成了其实践的特定模式。不过，利用本公开文件的利益，本领域的技术人员将会理解，在这些公开的特定实施例中可以进行许多改变，依然可以获得相象或相似的结果，而不脱离本发明的实质和范围。阅读以下实例时应当参考图1-图12，在以上的附图简要说明中已有小结。

实例一

两个Zea mays品种在广阔区域试验中的产量对比

我们首先说明本发明如何用于两个Zea mays品种在广阔区域的产量对比。这两个品种称为DK512和DK527，使用从1996年到1998年三年的广阔区域试验数据进行对比。

本过程的第一步包括使用帽子矩阵[31]式和Cook距离[33]式整理数据。对位置也进行了目测，以排除任何异常的空间簇。整理之后，对于DK512和DK527这两个杂交品种，跨年度的样本数目分别为n₁＝445，n₂＝389。两个杂交品种的数据合并之后，样本数目为N＝834。按年度的样本数目由下表给出。

表1

年度	杂交品种	样本数目
			1996	DK512	141
		113
			1997	DK512	150
		155
			1998	DK512	154
		121

本过程的下一步包括拟合[20]式中给出的混合模型。对于两个杂交品种，[20]式中的固定效应Xθ都是由[23]式中给出，但是没有截距项1。因此，[24]式中向量θ的β₁和β₂分量，每个的维数都是9×1。在向量θ中，α₁和α₂分量是3×1的子向量，表示1996年到1998年三年中每一年的截距参数。所以，[22]式中的X_ii＝1，2为n_i’12，[21]式中的X为N’24矩阵。对于[26]式中的方差协方差矩阵，我们使用[15]式中给出的共区域化法模型，而且M＝2，其中f₁(h；γ_m)和f₂(h；γ_m)都是[12]式中给出的球面模型。[28]式中给出了对每个变量带有一种单独的熔核效应的模型。使用的共区域化法模型带有熔核效应，以及如方程[25]之前的段落所介绍的一种位置匹配的随机效应，所有数据的方差-协方差模型由[29]式给出。

下一步是估计固定效应参数和协方差参数。使[19]式中给出的REML方程最小化，可以估计方差一协方差参数。估计出的参数由下表给出。

表2

参数	估计结果	参考方程
			σ1 2	59.6	[28]，[29]
σ2 2	57.6	[28]，[29]
			φ11 2	144.6	[15]，[29]
φ21 2	108.3	[15]，[29]
			r1	0.941	[15]，[29]
ρ1	0.050	[12]，[15]，[29]
			φ12 2	0.307	[15]，[29]
φ22 2	36.1	[15]，[29]
			r2	0.999	[15]，[29]
ρ2	0.012	[12]，[15]，[29]
			σ0 2	367.6	[25]，[29]

获得了参数估计结果之后，使用该模型进行两个杂交品种之间产量差异的空间预测和空间估计。对于单独的地理位置(点预测和估计)和地理区域(区块预测和估计)都要进行产量差异的空间预测和空间估计。除了这些空间计算之外，还使用传统的简单均值法估计地理区域的平均产量差异。(在这个实例中选用的地理区域表示推销商业化品种所使用的区域。)

回顾一下，估计的目标是确定两个品种产量差异的长期平均。预测的目标是在特定的年度和/或在特定的位置确定两个品种的平均产量差异。因为这个理由，空间估计是以跨年度的方式进行的，而空间预测是以按年度的方式进行的。为了提供与传统均值方法对比的基础，对按年度和跨年度的分析，都进行了简单均值分析。(跨年度简单均值方法是基于所有三年的汇总数据。)

使用表2中的参数估计结果和[38]式中的λ₀(s_p)，以及x₁₀＝x₀＝(1/3，1/3，1/3)￠，从[41]式获得市场区的跨年度平均产量差异的区块估计结果。(简单地设置x₁₀＝x₀＝(1/3，1/3，1/3)￠，在估计中为所有三个年度提供了相同的加权系数。)使用[41]式获得区块估计结果时，为每个市场区{s_p}∈A产生一个致密、规则的点网格，其中A是市场区。另外，从[41]式中估计出的数值，以其标准差——[42]式的平方根——划分时，假设具有标准正态分布；称这为z统计。如果z统计小于-1.96或者大于1.96，对于两个品种之间无差异的虚假设，就假定其显著程度在α＝0.05的级别。

使用表2中的参数估计结果和[54]式中的u₀，从[66]式获得多个市场区按年度平均产量差异的区块预测。在使用[54]式时，对于1996、1997和1998年，x₁₀和x₂₀都分别设置为(1，0，0)￠、(0，1，0)￠和(0，0，1)￠。对于每个市场区，在使用致密、规则的点网格方面，以及为了确定统计显著程度而计算和使用z统计方面，区块预测过程都类似于区块估计过程。在区块预测中计算z统计时，从[66]式中预测出的数值，以其标准差——[67]式的平方根——划分时，假设具有标准正态分布。

下一步我们将考虑两个品种之间产量差异的点估计和点预测。点估计和点预测都是在25×25的矩形网格上进行的，该网格覆盖着采集每个杂交品种数据的区域。在25×25网格上的点估计结果和点预测结果共同形成了产量差异的趋势曲面，它可以用于研究杂交品种试验的地理区域中产量差异的变化。

在以上讨论的25×25网格上，跨年度的点估计使用[44]式计算产量差异μ₀(s₀)。为了确定统计显著程度，已经发现，拒绝一个真实的、两个品种之间无差异的虚假设的概率为

$2\{1-\mathrm{\Φ}\left(\right|{\widehat{\mathrm{\μ}}}_0\left(s_0\right)/\sqrt{\mathrm{var}\left[{\widehat{\mathrm{\μ}}}_0\left(s_0\right)\right]}\left|\right)\},$

式中

由[45]式给出，Φ(·)为一个标准正态随机变量的累积分布函数。

在以上讨论的25×25网格上，按年度的点预测使用[69]式计算产量差异S₀(s₀)。为了确定统计显著程度，已经发现，拒绝一个真实的、两个品种之间无差异的虚假设的概率为

$2\{1-\mathrm{\Φ}\left(\right|{\hat{S}}_0\left(s_0\right)/\sqrt{\mathrm{var}\left[{\hat{S}}_0\left(s_0\right)\right]}\left|\right)\},$

式中

由[74]式给出，如同前面的情况，Φ(·)为一个标准正态随机变量的累积分布函数。

以上分析的全部结果以图形方式展示在图1至图12中。这些图分为两种类型，一种是等值线图，另一种是等值区域图。等值区域图是一种二维彩色地图，表示一个变量数值在一个区域的变化。

在等值线图中，三个级别(-5、0和+5蒲式耳/英亩)的产量差异和四个级别(产量差异±5％和±20％)的统计显著程度都用等值线表示。在等值线图中还显示了两个杂交品种的试验点和市场区的边界和名称。

在等值线图中，实等值线表示产量差异，虚等值线表示显著程度级别。如同任何等值线图，各条线的关系表示同一变量(如产量差异)的不同级别，用于确定曲面的大体形态。例如，对于产量差异，一个跨越红色、蓝色和绿色实线的横切面(以此顺序)，就是在产量差异曲面上向上移动，分别是从-5至0至+5蒲式耳/英亩。

对于产量差异和显著程度级别而言，等值线之间的关系用于在产量差异和统计显著程度的多种级别，确定产量差异发生显著变化的地理区域。例如，一条绿虚线围绕的区域落在一条绿实线围绕的区域中，后者又落在一条蓝实线围绕的区域中，表明该区域具有大于+5蒲式耳/英亩的、统计意义上显著的产量差异(在显著程度级别0.05)。

另一种类型的图形，也就是等值区域图，对于若干地理区域(市场区)中的每一个，显示这两个杂交品种在其中的平均产量差异。等值区域图上的图例既表明产量差异的级别(小于-5、-5和+5之间以及大于+5蒲式耳/英亩)，又表明差异是否达到5％级别的统计显著程度。

对于空间预测方法(仅仅按年度，1996年到1998年)和空间估计方法(仅仅跨年度)产生等值区域图和等值线图。等值线图不适合简单均值方法，所以只产生了等值区域图(既有按年度也有跨年度)。为了保持一致，如果某个市场区没有一个试验点种植了两个杂交品种，该市场区的产量差异就不会展示在等值区域图上(回顾一下，简单均值方法需要位置匹配的数据。)

为了对比传统的简单均值方法和空间方法，简单均值方法的按年度地图可以与按年度的区块预测地图对比，简单均值方法的跨年度地图可以与区块估计方法的跨年度地图对比。

实例二

有效性研究

为了表明本发明在农作物品种性能评价的应用条件下的功效，进行了几项有效性研究。对于本发明的两个组件：估计和预测，进行了分组的有效性研究。回顾以上讨论的实例，这两个组件中的每一个又包括两个部分：(1)点估计(或预测)部分，用于构建带有等值线的估计(或预测)曲面地图；(2)区块估计(或区块预测)部分，在实例中通过映射每个区域的平均值，用于构建市场区性能地图。

我们以对预测的有效性研究开始。首先，我们讨论点预测或点克里金的有效性研究。这项研究的基础是互证有效(Cressie，1993，101页)，使用Monsanto玉米研究系统对于两个DEKALB玉米杂交品种，DK512和DK527，从1996年到1998年的实际数据。对于两个杂交品种估计模型的预测质量，互证有效用作一种诊断检验。在互证有效的过程中，基本的思路是逐个地每次删除一个数据，然后使用其它数据预测已排除的数据。对于估计和预测的地质统计学方法，这是一种常见的评价方式。在进行互证有效时，834个数据中的每一个，都是一次排除一个，然后使用我们的过程预测已排除的数据。注意，这里我们将要预测随机变量Y_i(s₀)；i＝1，2，它是从数据集中排除的观测值之一。我们不想预测光滑值S_i(s₀)；i＝1，2，它具有与排除的熔核效应相关联的误差项。我们使用了协克里金预测量，使用Cressie(1993，138页)介绍的协方差。

在互证有效的研究中，计算了以下五种统计量。

1.均方预测误差(MSPE)，以量 ${[Y_i\left(s_0\right)-{\hat{Y}}_i\left(s_0\right)]}^2$ 的平均值计算，其中

是已排除数据的预测值，而Y_i(s₀)是观测值。

2.平均估计预测方差(Cressie，1993，140页，3.2.52式)。

3.80％预测区间包含真实值次数的比例。

4.标准化偏置： $[Y_i\left(s_0\right)-{\hat{Y}}_i\left(s_0\right)]/\mathrm{se}\left({\hat{Y}}_i\left(s_0\right)\right),$ ，(Cressie，1993，102页，2.6.15式)。

5.标准化MSPE，以量 ${[Y_i\left(s_0\right)-{\hat{Y}}_i\left(s_0\right)]}^2/\mathrm{var}\left({\hat{Y}}_i\left(s_0\right)\right).$ 的平均值的平方根计算。

在构建预测区间时，预测区间的下限计算为 ${\hat{Y}}_i\left(s_0\right)-1.28*\mathrm{se}\left({\hat{Y}}_i\left(s_0\right)\right),$ ${\hat{Y}}_1\left(s_0\right)+1.28*\mathrm{se}\left({\hat{Y}}_i\left(s_0\right)\right).,$ 上限计算为

这里原文53页19行2是预测标准差，它计算为估计预测方差(Cressie，1993，140页，3.2.52式)的平方根。

跨越所有三个年度归纳了以上的统计量。我们现在按照支持本发明的需要和实际服从的规律，展示互证有效研究的结果。

a.MSPE(统计量1)和平均估计预测方差(统计量2)应当相互接近。这两个统计量的互证有效研究结果分别为269.1和297.9。

b.预测区间(统计量3)应当接近名义(80％)覆盖。研究中的平均覆盖为82.7％。

c.标准化偏置(统计量4)应当接近0。研究中的标准化偏置计算结果为-0.021。

d.标准化MSPE(统计量5)应当接近1。研究中的标准化MSPE为1.049。

我们可以看到，互证有效的结果非常接近理论上需要的值。换句话说，该有效性结果支持以下说法：作为点预测工具，在[20]式中提出的模型符合我们的期望。

对于本发明中区块克里金组件的有效性研究，也是基于玉米杂交品种产量的实际数据分析。这种有效性研究的基础是把这个组件与传统的方法对比。回顾一下，在传统的方法中，来自一个区域的所有位置匹配的、并排种植的品种差异数据的简单均值，用于预测该区域的均值。我们将称这种方法为均值方法。有效性研究包括的步骤如下：

1.选择在国内有充足重叠数据的一对玉米杂交品种。

2.对于这对杂交品种，搜集其全部位置匹配的、并排种植的差异数据。这将被称为“原始种群数据”。

3.在原始种群数据的空间范围之内构建一个矩形。我们将称这个矩形为“预测矩形”，我们将在其中构建平均预测值。

4.使用预测矩形内部的所有原始种群数据的均值表示该矩形的“真实均值”。

5.随机选择原始种群数据的50％。这表示可用于分析的“样本数据”。

6.使用落在预测矩形内部的样本数据的均值表示用均值方法预测的、该矩形的均值。

7.使用(预测矩形内外的)全部样本数据，应用本发明的通用区块克里金组件预测该预测矩形的均值。这表示区块预测的数值。

8.在步骤4的真实均值与步骤6和步骤7的两个预测均值之中的每一个之间，计算平方差异。

9.对于几对杂交品种重复以上所有步骤。

10.对于所有的杂交品种对，使用步骤8的平方差异，对均值方法和区块克里金方法计算均方预测误差(MSPE)。

在这项有效性研究中，使用了Monsanto(DEKALB)研究系统的12对玉米杂交品种的数据。研究结果表明，区块克里金方法的MSPE比均值方法的MSPE低36％。注意，MSPE是预测误差的一种直接度量，所以，较低的MSPE对应于较高的预测精确度。

下一步，我们讨论对本发明的估计组件进行的有效性研究。这项分析是基于一项广泛的仿真研究，它以一种独特的方式，既使用真实数据，又使用仿真数据，正如我们下面的介绍。首先回顾一下，本发明的估计组件是基于广义最小平方(GLS)过程，它是对更传统的普通最小平方(OLS)过程的一种改进。与GLS过程不同，OLS忽略数据中的空间自相关，仅仅模拟大尺度的趋势。这项有效性研究的计划，是要对比三种截然不同的方法：(a)GLS方法、(b)OLS方法和(c)传统的均值方法，在品种对比的情况下，后者仅仅包括有关的位置匹配的、并排种植的两个品种的差异均值的计算。这项有效性研究的主要目的之一，是要了解先从均值方法到OLS方法，再从OLS方法到GLS方法，是否有增加的利益。

下面提供了这项有效性研究过程各步骤的介绍。

1.(从1994年、1995年或1996年中)选择一年，再选择相对成熟程度(RM)接近的一对玉米杂交品种。不考虑位置匹配，搜集这两个杂交品种的全部研究试验数据。

2.整理数据，排除空间异常值，并用一个矩形框围绕整理后的数据。我们将称之为数据矩形。这将表示可用于分析的样本数据。

3.在数据矩形中构建另一个内部矩形，在东西方向和南北方向都覆盖1/3到2/3的数据范围。我们将命名这个矩形为估计矩形。这就是我们将要进行估计的矩形。

4.使用数据矩形内的全部观测结果，对每个杂交品种分别拟合一个三次多项式趋势曲面。拟合这个模型是通过GLS完成的，先拟合一个球面协方差模型(参见背景技术一节)到每个杂交品种的残差协方差，然后使用两个品种共同的估计熔核效应。注意，在这两个杂交品种之间，并未假设相互的协方差。

5.在估计矩形中对两个品种及其差异(也就是两个拟合趋势曲面的差异)，计算趋势曲面均值的估计结果。我们将称它们为原始估计结果。

6.对于欧几里得距离0.01和0.05(分别对应于大约40和200英里)，计算每个杂交品种的自相关和两个杂交品种之间的互相关。

7.使用对应趋势曲面的估计参数和协方差参数以及步骤4的熔核效应，在数据矩形的0.01网格上，对每个杂交品种模拟数据。

8.在x方向和y方向都采用随机的振幅和频率产生两种正弦波。对于第一种正弦波，随机产生之频率的均值和标准差分别为100和5，而对于第二种正弦波它们分别是17和3。同样，对于第一种正弦波，随机产生之振幅的均值和标准差分别为4和2，而对于第二种正弦波它们分别是3和1。(注意，在模拟的数据加入这些正弦波，意味着引入某种不易被三次模型拟合的系统性变化。现在的数据结果等于三次曲面加正弦波加自相关后的随机误差。)

9.在估计矩形中，使用(通过步骤4中的GLS拟合出的)三次曲面参数和(步骤8中的)正弦波参数，计算两个杂交品种均值的“真实”值及其差异。这些真实的均值将用于在以下步骤中计算估计的偏置和置信区间的覆盖范围。

10.独立地对每个杂交品种随机地去除40％的数据点，在网格点上产生两个杂交品种的人工不匹配。

11.在估计矩形内部，仅在两种模拟值都有的网格点上，使用模拟值计算杂交品种的平均值及其差异。这将称为均值估计结果。

12.对于估计矩形内外的全部样本值(也就是使用整个数据矩形)，通过拟合三次多项式(两个杂交品种为不同的多项式)，以及假设两个杂交品种的残差相互独立但是残差的方差不同，使用OLS估计两个杂交品种(在估计矩形上)的均值及其差异。这就提供了OLS估计结果。

13.使用估计矩形内外的全部样本值(也就是使用整个数据矩形)，通过拟合三次多项式(两个杂交品种的不同)和一个带有熔核效应和残差互相关参数的球面协方差结构(两个杂交品种的不同)，使用GLS估计两个杂交品种(在估计矩形上)的均值及其差异。这就提供了GLS估计结果。

14.对于使用(i)均值估计结果、(ii)OLS估计结果和(iii)GLS估计结果的真实均值及其相应的标准差构建80％置信区间。标注置信区间是否覆盖步骤9中获得的真实均值。对以上三个估计结果中的每一个，还要计算对于真实均值的偏置。

15.重复步骤7到步骤14一百次。这是模拟的数目。

16.使用所有100次模拟的估计结果，对于(a)均值的均值估计结果、(b)均值的OLS估计结果和(c)均值的GLS估计结果计算标准差。通过采用步骤14中获得的偏置(包括所有100次模拟)的平均值，还要计算三个估计结果中每一个的平均偏置。

17.对于均值、OLS和GLS方法中的每一种，计算(在100次模拟中)置信区间覆盖真实均值次数的百分比。

18.对于覆盖宽广的RM范围的27对杂交品种和对于1994年、1995年和1996年三年的所有数据，重复以上的所有步骤。

注：如上所述，注意到这种有效性研究以一个独特的方式结合了实际数据和模拟结果，首先使用每对玉米杂交品种的实际数据获得趋势和相关参数，然后使用这些参数产生模拟研究所用的数据。注意，在这项有效性研究中不能仅仅使用实际数据的理由是，这个过程需要“真实”均值来计算覆盖范围的百分比和偏置，如步骤14中的介绍。

我们现在归纳这项有效性研究的结果。

1.以标准差度量提高了精度，从均值方法到GLS方法为37％。这种提高在从均值到OLS和从OLS再到GLS之间等分。

2.GLS方法的80％置信区间覆盖真实均值次数的百分比(77％)非常接近名义值，而对于均值和OLS方法则比80％低很多(均值为46％，OLS方法为56％)。

3.对于均值和OLS方法，这项研究揭示了真实标准差的总体欠估计，这意味着对无差异的虚假设有太多的假拒绝，也就是品种之间确实没有差异时，结论却是品种差异存在。

结论：

我们在这里介绍过的所有这些有效性研究，共同表明了本发明在农作物品种性能的空间估计和空间预测方面的功效。

参考文献

这里专门引用下列参考文献和说明书中提及的所有参考文献作为参考。

Besag，J.，and Kempton，R.(1986)，Statistical analysis of field experiments usingneighboring plots，Biometrics，42，231-251.

Bhatti，A.U.，Mulla，D.J.，and Frazier，B.E.(1991)，″Estimation of soil propertiesand wheat yields on complex eroded hills using geostatistics and thematicmapper images，″Remote Sens.Environ.，37，181-191.

Bhatti，A.U.，Mulla，D.J.，Koehler，F.E.，and Gurmani，A.H.(1991)，″Identifying andremoving spatial correlation from yield experiments，″Soil Sci.Soc.Am.J.，55，1523-1528.

Bradley，J.P.，Knittle，K.H.，and Troyer，A.F.(1988)，″Statistical methods in seedcom product selection，″J.Prod.Agric.，1，34-38.

Brownie，C.，Bowman，D.T.，and Burton，J.W.(1993)，″Estimating spatial variationin analysis of data from yield trials：a comparison of methods，″Agron.J.，85，1244-1253.

Cook，R.D.(1977)，“Detection of influential observation in linear regression，”Technometrics，19，15-18.

Cressie，N.(1990)，″The origins of kriging，″Mathematical Geology，22，239-252.

Cressie，N.(1993)，Statistics for spatial data，Revised Edition，John Wiley andSons，New York，USA.

Cressie，N.，and Lahiri，S.N.(1996)，“Asymptotics for REML estimation of spatialcovariance parameters，”Journal of Statistical Planning and Inference，50，327-341.

Christensen，R.，Pearson，L.M.，and Johnson，W.(1992)，“Case-deletion diagnosticsfor mixed models，”Technometrics，34，38-45.

Cullis，B.，Gogel，B.，Verbyla，A.，and Thompson，R.(1998)，“Spatial analysis ofmulti-environment early generation variety trials，”Biometrics，54，1-18.

Gandin，L.S.(1959)，“The problem of optimal interpolation，”Trudy GGO 99，67-75.

Gandin，L.S.(1963)，Objective Analysis of Meteorological Fields，Gidrometeorologicheskoe Izdatel′stvo(GIMIZ)，Leningrad(translated by IsraelProgram for Scientific Translations，Jerusalem，1965).

Godambe，V.P.(editor)(1991)，Estimatingfunctions，Clarendon Press，Oxford.

Goldberger，A.S.(1962)，“Best linear unbiased prediction in the generalized linearregression model，”Journal of the American Statistical Association，57，369-375.

Goovaerts，P.(1997)，Geostatistics for Natural Resources Evaluation，OxfordUniversity Press，New York.

Gotway，C.A.，and Hartford，A.H.(1996)，“Geostatistical methods for incorporatingauxiliary information in the prediction of spatial variables，”Journal ofAgrricultural Biological and Environmental Statistics，1，1，17-39.

Gotway，C.A.，and Stroup，W.W.(1997)，“A generalized linear model approach tospatial data anlysis and prediction，”Journal of Agricultural Biological andEnvironmental Statistics，2，2，157-178.

Harville，D.A.(1977)，“Maximum likelihood approaches to variance componentestimation and to related problems”Journal of the American StattisticalAssociation，72，320-340.

Heyde，C.C.(1994)，“A quasi-likelihood approach to the REML estimatingequations，”Statistics and Probability Letters，21，381-384.

Hoaglin，D.C.，and Welsch，R.E.(1978)，“The hat matrix in regression andANOVA，”American Statistician，32，17-22.

Kolmogorov，A.N.(1941)，“Interpolation and extrapolation of stationary randomsequences，”Isvestiia Akademii Nauk SSSR，Seriia Matematicheskiia，5，3-14.(Translation，1926，Memo RM-3090-PR，Rand Corp.，Santa Monica，CA).

Littell，R.C.，Milliken，G.A.，Stroup，W.W.，and Wolfinger，R.D.(1996)，SASSystem for Linear Models，SAS Institute Inc.，Cary，North Carolina，USA.

Mardia，K.V.，and R.J.Marshall.(1984)，“Maximum likelihood estimation of modelsfor residual covariance in spatial regression，”Biometrika，71，135-146.

Matheron，G.(1962)，Traitéde Géostatistique Appliquée，Tome I.Mémoires duBureau deRecherches Géologlogiques et Minières，14，Editions Technip，Paris.

Matheron，G.(1969)，“Le Krigeage Universel，”Cahiers du Centre de MorphologieMathematique，1，Fontainebleau，France.

Montgomery，D.C.，and Peck，E.A.(1992)，Introduction to Linear RegressionAnalysis Second Edition，Wiley，New York.

Ovalles，F.A.，and Collins，M.E.(1988).“Evaluation of soil variability in northwestFlorida using geostatistics，”Soil Sci.Soc.Am.J.，52，1702-1708.

Patterson，H.D.，and R.Thompson.(1971)，“Recovery of interblock informationwhen block sizes are unequal，”Biometrika，58，545-554.

Patterson，H.D.，and R.Thompson.(1974)，“Maximum likelihood estimation ofcomponents of variance，”Pages 197-207in Proceedings of the 8thInternational Biometric Conference，Biometric Society，Washington，DC，USA.

Searle，S.R.(1971)，Linear Models，Wiley，New York.

Searle，S.R.，Casella，G.，and McCulloch，C.E.(1992)，Variance Components，Wiley，New York.

Sprague，G.F.，and Eberhart，S.A.(1976)，“Corn breeding，”Corn and CornImprovement，J.A.Dudley and G.F.Sprague(eds)，Iowa State Univ.Press.

Stroup，W.W.，Baenziger，P.S.，and Mulitze，D.K.，(1994)“Removing spatialvariation from wheat yield trials：a comparison of methods，”Crop Science，86，62-66.

Ver Hoef，J.M.，and Barry，R.D.(1998)，“Modeling crossvariograms for cokrigingand multivariable spatial prediction，”Journal of Statistical Planning andInference，69，275-294.

Ver Hoef，J.M.，and Cressie，N.(1993)，“Multivariable spatial prediction，”Mathematical Geology，25，219-240.

Weiner，N.(1949)，Extrapolation，Interpolation，and Smoothing of Stationary TimeSeries，MIT Press，Cambridge，MA.

Wold，H.(1938)，A Study in the Analysis of Stationary Time Series，Almqvist andWiksells，Uppsala.

Wolfinger，R.D.，Tobias，R.D.，and Sall，J.(1994)，“Computing Gaussian likelihoodsand their derivatives for general linear mixed models，”SL4M Journal onScientifc Computing，15，1294-1310.

Wu，T.，Mather，D.E.，and Dutilleul，P.(1998)，“Application of geostatistical andneighbor analyses to data from plant breeding trials，”Crop Science，38，1545-1553.

Yost，R.S.，Uehara，G.，and Fox，R.L.(1982a)，“Geostatistical analysis of soilchemical properties of large land areas.I.Semi-variograms，”Soil Sci.Soc.Am.J.，46，1028-1032.

Yost，R.S.，Uehara，G.，and Fox，R.L.(1982b)，“Geostatistical analysis of soilchemical properties of large land areas.II.Kriging，”SoilSci.Soc.Am.J.，46，1033-1037.

Zimmerman，D.L.(1989)，“Computationally efficient restricted maximum likelihoodestimation of generalized covariance functions，”Mathematical Geology，21，655-672.

Zimmerman，D.L.，and Harville，D.A.(1991)，“A random field approach to theanalysis of field-plot experiments and other spatial experiments，”Biometrics，47，1，223-239.

Claims (87)

1.使用线性混合模型评价两个或多个农作物品种广阔区域性能的一种方法，该模型内含地质统计部分并包括固定效应、随机效应和协方差的参数，本方法包括：

构建一个广阔区域数据库，包括两个或多个农作物品种的试验点空间坐标，两个或多个农作物品种在其各自试验点的性能特征值，以及用于评价两个或多个农作物品种的性能的其它地理位置的空间坐标；

为广阔区域数据库中的数据构建一个线性混合模型，包括通过固定效应趋势表面对两个或多个农作物品种中的每一个建立大规模性能模型，建立两个或多个农作物品种中的每一个的随机误差分量之间的空间自协方差的模型，以及在多变量核心引力场量子化模型中的两个或多个农物作品种之间建立随机误差的空间交叉协方差模型；

利用广阔区域数据库中的数据拟合线性混合模型，估计固定效应、随机效应和协方差的参数；以及

使用参数估计结果，对于一个或多个地理位置中的任意一个，估计所述两个或多个农作物品种中任意一个的长期期望性能。

2.根据权利要求1的方法，其特征在于，数据库进一步包括共变量数据。

3.根据权利要求1的方法，其特征在于，估计长期期望性能包括估计两个或多个农作物品种之间的长期期望性能差异。

4.根据权利要求1的方法，其特征在于，估计协方差参数包括一种约束最大似然方法。

5.根据权利要求1的方法，其特征在于，估计固定效应参数包括一种广义最小平方方法。

6.根据权利要求1的方法，其特征在于，数据库进一步包括用于评价两个或多个农作物品种的性能的地理区域，每个地理区域都包含多个地理位置。

7.根据权利要求6的方法，进一步包括：

使用参数估计结果，对于一个或多个地理区域中的任意一个，估计两个或多个农作物品种中的任意一个的长期期望性能。

8.根据权利要求1的方法，进一步包括形成长期期望性能的一种输出。

9.根据权利要求8的方法，其特征在于，该输出包括文本输出。

10.根据权利要求8的方法，其特征在于，该输出包括图形输出。

11.根据权利要求10的方法，其特征在于，该输出包括一个等值线图，表示长期期望性能的连续曲面。

12.使用线性混合模型评价两个或多个农作物品种广阔区域性能的一种方法，该模型内含地质统计部分并包括固定效应、随机效应和协方差的参数，本方法包括：

构建一个广阔区域数据库，包括两个或多个农作物品种的试验点空间坐标、用于评价两个或多个农作物品种的性能的地理区域，每个地理区域都包含多个地理位置，和两个或多个农作物品种在其相应试验点的性能特征值；

为广阔区域数据库中的数据构建一个线性混合模型，包括通过固定效应趋势表面对两个或多个农作物品种中的每一个建立大规模性能模型，建立两个或多个农作物品种中的每一个的随机误差分量之间的空间自协方差的模型，以及在多变量核心引力场量子化模型中的两个或多个农物作品种之间建立随机误差的空间交叉协方差模型；

利用广阔区域数据库中的数据拟合线性混合模型，估计固定效应、随机效应和协方差的参数；以及

使用参数估计结果，对于一个或多个地理区域中的任意一个，估计所述两个或多个农作物品种的长期期望性能。

13.根据权利要求12的方法，其特征在于，数据库进一步包括共变量数据。

14.根据权利要求12的方法，其特征在于，估计长期期望性能包括估计两个或多个农作物品种之间的长期期望性能差异。

15.根据权利要求12的方法，其特征在于，估计协方差参数包括一种约束最大似然方法。

16.根据权利要求12的方法，其特征在于，估计固定效应参数包括一种广义最小平方方法。

17.根据权利要求12的方法，进一步包括形成长期期望性能的一种输出。

18.根据权利要求17的方法，其特征在于，该输出包括文本输出。

19.根据权利要求17的方法，其特征在于，该输出包括图形输出。

20.根据权利要求19的方法，其特征在于，该输出包括一个等值区域图，对于一个或多个指定地理区域中的任意一个，表示长期期望性能。

21.使用线性混合模型评价两个或多个农作物品种广阔区域性能的一种方法，该模型内含地质统计组件并包括固定效应、随机效应和协方差的参数，本方法包括：

构建一个广阔区域数据库，包括两个或多个农作物品种的试验点空间坐标，两个或多个农作物品种在其各自试验点的来自一个或多个时间期间的性能特征值，以及用于评价两个或多个农作物品种的性能的其它地理位置的空间坐标；

为广阔区域数据库中的数据构建一个线性混合模型，包括通过固定效应趋势表面对两个或多个农作物品种中的每一个建立大规模性能模型，建立两个或多个农作物品种中的每一个的随机误差分量之间的空间自协方差的模型，以及在多变量核心引力场量子化模型中的两个或多个农物作品种之间建立随机误差的空间交叉协方差模型；

利用广阔区域数据库中的数据拟合线性混合模型，估计固定效应、随机效应和协方差的参数；以及

使用参数估计结果，对于一个或多个空间位置中的任意一个，以及对于一个或多个时间期间的任意一个，预测两个或多个农作物品种的性能。

22.根据权利要求21的方法，其特征在于，估计协方差参数包括一种约束最大似然方法。

23.根据权利要求21的方法，其特征在于，数据库进一步包括共变量数据。

24.根据权利要求23的方法，其特征在于，共变量数据包括一个或多个固定效应。

25.根据权利要求21的方法，其特征在于，预测性能包括预测两个或多个农作物品种之间的性能差异以及相关的偏差。

26.根据权利要求21的方法，其特征在于，数据库进一步包括用于评价两个或多个农作物品种的性能的地理区域，每个地理区域都包含多个地理位置。

27.根据权利要求26的方法，进一步包括：

使用参数估计结果，对于一个或多个指定地理区域中的任意一个，以及对于一个或多个指定时间期间的任意一个，预测两个或多个农作物品种中的任意一个的性能。

28.根据权利要求21的方法，进一步包括形成预测性能的一种输出。

29.根据权利要求28的方法，其特征在于，该输出包括文本输出。

30.根据权利要求28的方法，其特征在于，该输出包括图形输出。

31.根据权利要求30的方法，其特征在于，该输出包括一个等值线图，表示预测性能的连续曲面。

32.使用线性混合模型评价两个或多个农作物品种广阔区域性能的一种方法，该模型内含地质统计部分并包括固定效应、随机效应和协方差的参数，本方法包括：

构建一个广阔区域数据库，包括两个或多个农作物品种的试验点空间坐标、用于评价两个或多个农作物品种中的性能的地理区域，每个地理区域都包含多个地理位置，和两个或多个农作物品种在其相应试验点的来自一个或多个时间期间的性能特征值；

为广阔区域数据库中的数据构建一个线性混合模型，包括通过固定效应趋势表面对两个或多个农作物品种中的每一个建立大规模性能模型，建立两个或多个农作物品种中的每一个的随机误差分量之间的空间自协方差的模型，以及在多变量核心引力场量子化模型中的两个或多个农物作品种之间建立随机误差的空间交叉协方差模型；

利用广阔区域数据库中的数据拟合线性混合模型，估计固定效应、随机效应和协方差的参数；以及

使用参数估计结果，对于一个或多个地理区域中的任意一个，以及对于一个或多个时间期间的任意一个，预测两个或多个农作物品种的性能。

33.根据权利要求32的方法，其特征在于，估计协方差参数包括一种约束最大似然方法。

34.根据权利要求32的方法，其特征在于，数据库进一步包括共变量数据。

35.根据权利要求34的方法，其特征在于，共变量数据包括一个或多个固定效应。

36.根据权利要求32的方法，其特征在于，预测性能包括预测该两个或多个农作物品种之间的性能差异以及相关的偏差。

37.根据权利要求32的方法，进一步包括形成预测性能的一种输出。

38.根据权利要求37的方法，其特征在于，该输出包括文本输出。

39.根据权利要求37的方法，其特征在于，该输出包括图形输出。

40.根据权利要求39的方法，其特征在于，该输出包括一个等值区域图，对于一个或多个指定地理区域中的任意一个，表示预测性能。

41.使用线性混合模型评价两个或多个农作物品种广阔区域性能的一种方法，该模型内含地质统计部分并包括固定效应、随机效应和协方差的参数，本方法包括：

构建一个广阔区域数据库，包括两个或多个农作物品种的试验点空间坐标、用于评价两个或多个农作物品种的性能的地理区域，每个地理区域包含多个地理位置，用于评价两个或多个农物作品种的性能的其它地理位置的空间坐标和两个或多个农作物品种在其相应试验点的性能特征值；

为广阔区域数据库中的数据构建一个线性混合模型，包括通过固定效应趋势表面对两个或多个农作物品种中的每一个建立大规模性能模型，建立两个或多个农作物品种中的每一个的随机误差分量之间的空间自协方差的模型，以及在多变量核心引力场量子化模型中的两个或多个农物作品种之间建立随机误差的空间交叉协方差模型；

利用广阔区域数据库中的数据拟合线性混合模型，估计固定效应、随机效应和协方差的参数；

使用参数估计结果，对于一个或多个地理位置中的任意一个，估计两个或多个农作物品种中的任意一个的长期期望性能；

使用参数估计结果，对于一个或多个地理区域中的任意一个，估计两个或多个农作物品种的长期期望性能；

使用参数估计结果，对于一个或多个地理位置中的任意一个，以及对于一个或多个时间期间的任意一个，预测两个或多个农作物品种的性能；以及

使用参数估计结果，对于一个或多个地理区域中的任意一个，以及对于一个或多个时间期间的任意一个，预测两个或多个农作物品种的性能。

42.根据权利要求41的方法，其特征在于，预测性能包括预测两个或多个农作物品种之间的性能差异以及相关的偏差。

43.根据权利要求41的方法，其特征在于，估计固定效应的参数包括一种广义最小平方方法。

44.根据权利要求41的方法，其特征在于，预测性能包括通用协克里金或通用区块协克里金方法。

45.根据权利要求41的方法，其特征在于，估计协方差参数包括一种约束最大似然方法。

46.根据权利要求41的方法，其特征在于，估计长期期望性能包括估计两个或多个农作物品种之间的长期期望性能差异。

47.根据权利要求41的方法，进一步包括形成估计长期期望性能或预测平均性能的一种输出。

48.根据权利要求47的方法，其特征在于，该输出包括文本输出。

49.根据权利要求47的方法，其特征在于，该输出包括图形输出。

50.根据权利要求1的方法，其特征在于，构建线性混合模型进一步包括建立试验点的模型作为随机效应。

51.根据权利要求2的方法，其特征在于共变量数据包括一个或多个固定效应。

52.根据权利要求1的方法，其特征在于估计长期期望性能包括估计两个或多个农作物品种中任意一个的固定效应趋势表面的值及其相关的偏差。

53.根据权利要求3的方法，其特征在于估计两个或多个农作物品种之间的长期期望性能差异包括估计两个或多个农作物品种的固定效应趋势表面中的差异以及所述差异的相关的偏差。

54.根据权利要求53的方法，进一步包括形成长期期望性能差异的输出。

55.根据权利要求54的方法，其特征在于输出包括文本输出。

56.根据权利要求54的方法，其特征在于输出包括图形输出。

57.根据权利要求56的方法，其特征在于输出包括一个等值线图，表示长期期望性能差异及其统计显著性。

58.根据权利要求1的方法，其特征在于共变量数据包括年的固定效应，并进一步包括：使用参数估计结果，对于一个或多个地理位置中的任意一个以及对于一个或多个时间期间中的任意一个，预测两个或多个农作物品种中任意一个的性能。

59.根据权利要求12的方法，其特征在于，构建线性混合模型进一步包括建立试验点的模型作为随机效应。

60.根据权利要求13的方法，其特征在于共变量数据包括一个或多个固定效应。

61.根据权利要求12的方法，其特征在于估计长期期望性能包括估计跨过一个或多个地理区域中的每个地理区域中的所有地理位置的固定效应趋势表面的平均值，以及两个或多个农作物品种中任意一个的所述平均值的偏差。

62.根据权利要求14的方法，其特征在于估计两个或多个农作物品种之间的长期期望性能差异包括估计两个或多个农作物品种的固定效应趋势表面中的差异以及所述差异的相关的偏差。

63.根据权利要求62的方法，进一步包括形成长期期望性能差异的输出。

64.根据权利要求63的方法，其特征在于输出包括文本输出。

65.根据权利要求63的方法，其特征在于输出包括图形输出。

66.根据权利要求65的方法，其特征在于输出包括一个地区分布图，表示长期期望性能差异及其统计显著性。

67.根据权利要求12的方法，其特征在于共变量数据包括年的固定效应，并进一步包括：使用参数估计结果，对于一个或多个地理位置中的任意一个以及对于一个或多个时间期间中的任意一个，预测两个或多个农作物品种中任意一个的性能。

68.根据权利要求21的方法，其特征在于，构建线性混合模型进一步包括建立试验点的模型作为随机效应。

69.根据权利要求21的方法，其特征在于估计固定效应参数包括广义最小平方方法。

70.根据权利要求21的方法，其特征在于预测性能包括应用于两个或多个农作物品种中的任意一个的通用克里金方法。

71.根据权利要求70的方法，其特征在于所述通用克里金方法包括：使用多变量核心引力场量子化模型，利用每个品种内的自相关以及各品种之间的互相关建立两个或多个农作物品种中任意一个的固定效应趋势表面和相应的随机表面的模型，并且估计预测偏差。

72.根据权利要求25的方法，其特征在于预测性能差异包括应用于两个或多个农作物品种的差异的通用克里金方法。

73.根据权利要求72的方法，进一步包括形成预测性能差异的输出。

74.根据权利要求73的方法，其特征在于输出包括文本输出。

75.根据权利要求73的方法，其特征在于输出包括图形输出。

76.根据权利要求75的方法，其特征在于输出包括一个等值线图，表示预测的性能差异的连续表面及其统计显著性。

77.根据权利要求32的方法，其特征在于，构建线性混合模型进一步包括建立试验点的模型作为随机效应。

78.根据权利要求32的方法，其特征在于估计固定效应参数包括广义最小平方方法。

79.根据权利要求32的方法，其特征在于预测性能包括应用于两个或多个农作物品种中的任意一个的通用克里金方法。

80.根据权利要求79的方法，其特征在于所述通用克里金方法包括：使用多变量核心引力场量子化模型，利用每个品种内的自相关以及各品种之间的互相关建立两个或多个农作物品种中任意一个的固定效应趋势表面和相应的随机表面的模型，并且估计预测偏差。

81.根据权利要求36的方法，其特征在于预测性能差异包括应用于两个或多个农作物品种的差异的通用克里金方法。

82.根据权利要求81的方法，进一步包括形成预测性能差异的输出。

83.根据权利要求82的方法，其特征在于输出包括文本输出。

84.根据权利要求82的方法，其特征在于输出包括图形输出。

85.根据权利要求56的方法，其特征在于输出包括一个地区分布图，表示预测性能差异及其统计显著性。

86.根据权利要求41的方法，其特征在于，构建线性混合模型进一步包括建立试验点的模型作为随机效应。

87.根据权利要求41的方法，其特征在于估计长期期望性能包括估计固定效应趋势表面的值及其相关的偏差。

Images