一种从被动微波数据估算瞬时地表发射率的方法【专利摘要】本发明公开了一种从被动微波数据估算瞬时地表发射率的方法，该方法所使用的模型依赖于被动微波的通道亮温数据和大气总可降水量数据，具体包括以下四个步骤：估算模型的推演、方程系数的率定、数据采集、被动微波地表发射率估算。本发明由辐射传输方程推导而来，具有坚实的物理意义和较好的普适性，可以应用于不同的地面情况和不同的季节；其输入数据不多且容易获得，而且不受云的限制，可在包括有云情况下的几乎所有天气情况下获得，保证了模型的全天候应用和瞬时地表发射率的获得，填补了传统辐射传输方法只能生产月时间尺度地表发射率的空白。【专利说明】一种从被动微波数据估算瞬时地表发射率的方法
技术领域：
[0001]本发明涉及一种估算瞬时地表发射率的方法，尤其涉及一种从被动微波数据估算瞬时地表发射率的方法。【
背景技术：
】[0002]微波地表发射率是表征地表辐射能力的重要参数，是反演地表温度、土壤含水量、降水量、云中液态水量等一系列重要参数的关键因子。然而，由于微波地表发射率受多种因素影响，如土壤含水量、地表粗糙度、土壤类型、观测频率、入射角等，而这些参数本身随时间和空间的变化而变化，因此地表发射率也具有时间和空间上的多变性，对其的准确测量和估算一直是一个难题。遥感估算法被认为是获得全球大尺度时间连续性微波发射率估算的唯一方法，受到了越来越多的重视。[0003]目前，比较常见的微波地表发射率遥感估算法主要有统计回归法、神经网络法和辐射传输法三类。统计回归法是将卫星的通道亮温看做地表发射率的直接体现，依据其二者的强相关性，利用统计的方法来建立卫星的通道亮温与对应地表发射率之间的经验方程，从而依据亮温数据推算出地表发射率;该方法的输入参数只有亮温值，很难表达复杂的非线性问题，其结果难以准确区分陆表和大气参数的变化对地表发射率的影响;且回归方程本身不具备明确的物理意义，只是一种统计经验的表述，其精度在不同地区、不同季节都可能会不同；针对不同情况，模型的形式(线性、非线性）、选择的通道(频率)等都存在诸多变数，这对于模型的广泛和业务化应用造成了困难，严重影响了模型的适应性。[0004]神经网络法是模拟大量相互联系的神经元来逼近任意复杂的非线性关系，能够通过大量的训练数据发现地表发射率与其影响因素之间的某种联系，不需要求解复杂的数学方程，从数学层面弱化了非线性求解问题。但是神经网络方法解决问题的过程是不可见的，也无法给出详细的数学表达式和计算过程，其本质更接近"黑箱"，因此该方法本身不具备物理意义；同时，神经网络估算法首先需要大量的训练数据才能更好的逼近最优解，对训练数据的较高要求增加了研究区筛选的难度;此外，发射率估算的精度也会与训练数据的质量直接相关，这也限制了神经网络估算法的推广使用。[0005]辐射传输法的核心思想是基于辐射传输模型，模拟从地表发射辐射能量开始，到卫星接收辐射能量为止的整个过程中，辐射能量与地表、大气、卫星之间的相互作用，用物理的方法对整个过程加以描述，并用数学模型表达出来;这种估算法具有坚实的物理基础，解决问题的过程明确，相比于统计回归法和神经网络法，辐射传输估算法具有更广阔的使用空间。但是传统的辐射传输方法存在两个主要缺点：1)依赖大量的辅助数据作为输入，导致了该方法的使用很大程度上受制于辅助数据的能否获得和数据质量;2)现有的发射率估算模型在有云情况下无法使用，使得同一个地区不能保证每天都有发射率数据，因此现有的发射率产品多是月尺度产品；月尺度发射率产品会平均掉发射率每日的变化情况，对于很多要求瞬时发射率的应用，比如日尺度降水量的反演、天气预测预报等，都会失去意义。【
发明内容】[0006]为了解决上述技术所存在的不足之处，本发明提供了一种从被动微波数据估算瞬时地表发射率的方法。[0007]为了解决以上技术问题，本发明采用的技术方案是：一种从被动微波数据估算瞬时地表发射率的方法，包括以下四个步骤：[0008]I、估算模型的推演：一般星载被动微波传感器所接收到的辐射信息包含三个部分:大气向上辐射，用TTa表示；由地面向上发射并被大气衰减的辐射，用T?e?Ts表示;大气向下与宇宙背景辐射被地面反射进入大气，又被大气衰减，并最终被微波传感器接收的辐射，用（1-小(汍+1^)表示，这三个过程组成了完整的辐射传输方程，如公式(1)所示：[0009]TB=TTa+T?e?Ts+a-ehUk+Tsky)(1)[0010]其中，e为某一通道的地表发射率，无量纲;Tb为通道亮温，单位为Kj是大气透过率，无量纲;TTa是大气上行亮温、Tk是大气下行亮温，单位均为K;Ts为地表温度，单位为K;Tsky是宇宙背景辐射温度，单位为K，可视为定值2.73K;[0011]由于TTa和相差很小，因此假设二者相等，将公式(1)进一步表示为公式(2):[0012]TB=Ta+x?e?Ts+(l-e)x(Ta+2.73)(2)[0013]由于不同的通道具有独立的发射特性以及大气衰减特性，因此1和1以及e对于不同通道而言都是不同的，因此通道亮温Tb对不同频率而言也是不同的；假设有两个不同的频率，将它们的福射传输方程表不为：[0014]TB,i=Ta,i+Ti?￡i?Ts+(l-ei)Ti(Ta,i+2.73)(3)[0015]TBj=Ta>j+Tj?ej?Ts+(l-ej)Tj(Taj+2.73)(4)[0016]其中，i和j分别代表两个不同的频率；根据前人的经验，可以假设18.7GHz和23.8GHz两个频率的发射率相等，由此将公式(3)和公式(4)分别进一步表示为：[0017]TB,18.7=Ta,18.7+Tl8.7*￡*Ts+(1~￡)^18.7(Ta,18.7+2.73)(5)[0018]TB,23.8=Ta,23.8+T23.8?e?Ts+(l-e)T23.8(Ta,23.8+2?73)(6)[0019]将公式(5)和公是(6)右边的第一项和第三项均移到左边，它们的左右两边再同除以各自频率的t得到公式(8)和公式(9):[0022]将公式(7)和公式(8)联立，得到公式(9):[0024]由此，地表温度这个变量被约去;进一步对公式(9)变形，将发射率e单独放到公式左侧，其他项放到公式右侧，得到公式(10):[0026]先前学者已经发现了t和Ta与大气总可降水量(precipitablewatervapor，PWV)呈非常好的线性关系，即：[0027]x=a?PWV+b(11)[0028]Ta=c?PWV+d(12)[0029]其中，a，b，c，d均是线性方程的系数；由此可知，公式（10)中的1^,23.8*丁18.7-Ta,18.7?T23.8和T18.7?123.8(!3,18.7-1'￡1,23.8)都能表不成只有?1￥-个未知数的方程;根据热动力初估检索数据集提供的全球大气温湿度廓线数据，利用大气辐射传输模型计算了18.76泡和23.86泡的了18.7、123.8、1^18.7、1^23.8以及对应的13評，用线性方程描述1^23.8.Tl8.7-Ta，18.7?丁23.8、丁18.7?T23.8(Ta，18.7-Ta，23.8)^PIPWV|^**，_P:[0030]Ta,23.8*Tl8.7-Ta,18.7*X23.8=e*PWV+f(13)[0031]118.7*T23.8(Ta,18.7-Ta>23.8)=g*PWV+h(14)[0032]其中，e，f，g，h均是线性方程的系数；[0033]根据以上分析，将公式（10)进一步变形为一个只包含通道亮温和PWV的方程，即只包含被动微波成像辐射仪的18.7Hz和23.8GHz的H极化通道亮温和大气总可降水量的被动微波地表发射率估算模型，如公式(15)所示：[0035]其中，e为地表发射率，无量纲；Tb,18.7和Tb,23.8分别为被动微波成像辐射仪的18.7GHz和23.8GHz两个通道的H极化亮温，单位均为K;PWV为大气总可降水量，单位为g/0112;1]11、1]12、1]13、1114、1]15、1]16、1]17、1]18均为方程系数；[0036]11、方程系数的率定:针对研究区的实际情况，对模型中的八个系数1111、1112、1113、1114、m5、m6、m7、m8进行率定，确保上述估算模型在特定研究区拥有较高的精度;方程系数率定的具体方法如下：[0037]a、采集研究区的以下三类历史数据:被动微波成像辐射仪的18.7GHz和23.8GHz的H极化通道亮温数据，该类数据的获得途径包括带有被动微波成像辐射仪的卫星接收站点、数据共享网站;与亮温数据时间一致的大气总可降水量的地面测量数据，该类数据的获得途径包括地面观测网点的测量数据，测量方法包括地面气象站点释放的无线电探空气球、地基GPS、地基微波辐射计、太阳光度计；与亮温数据时间一致的被动微波地表发射率测量数据，该类数据从地面测量获得；[0038]b、同一时间的18.7GHz和23.8GHz的H极化亮温数据、大气总可降水量数据、被动微波地表发射率数据，这三类数据均齐全算作一组，采集的历史数据应保证不少于八组；[0039]c、采用最小二乘法率定系数:将采集的历史数据代入公式（15)，建立不少于八个方程的方程组，采用最小二乘法计算出模型的八个系数值。[0040]m、数据采集:采集研究区内待研究时间段的被动微波成像辐射仪的18.7GHz和23.8GHz的H极化通道亮温数据，该类数据的获得途径包括带有被动微波成像辐射仪的卫星的接收站点、数据共享网站;采集研究区内待研究时间段的与亮温数据时间一致的大气总可降水量的测量数据，该类数据的获得途径包括地面观测网点的测量数据，测量方法包括地面气象站点释放的无线电探空气球、地基GPS、地基微波福射计、太阳光度计；同一时间内的卫星亮温和大气总可降水量这两类数据齐全算作一组；[0041]IV、被动微波地表发射率估算:将步骤n中采集的数据代入经过步骤I率定好参数的估算模型中，直接计算出与输入数据的采集时间同时的瞬时被动微波地表发射率。[0042]本发明建立了可以在任何天气情况下估算被动微波瞬时地表发射率的方法，具体表现为：[0043](1)本发明只需要借助被动微波的通道亮温数据和大气总可降水量数据，即可估算瞬时地表发射率，不需要地表温度作为输入，避免了传统方法在有云情况下无法获得地表温度而造成模型无法使用的情况；[0044](2)本发明能够获得的发射率是瞬时的，而不是当前常见的月尺度的，对许多应用和生产具有更加重要的意义；[0045](3)本发明可以应用于不同的地区和不同的季节。【附图说明】[0046]图1为本发明的整体流程图。【具体实施方式】[0047]下面结合附图和【具体实施方式】对本发明作进一步详细的说明。[0048]如图1所示，本发明包括以下四个步骤：[0049]I、估算模型的推演：一般星载被动微波传感器所接收到的辐射信息包含三个部分:大气向上辐射，用TTa表示；由地面向上发射并被大气衰减的辐射，用t?e?Ts表示;大气向下与宇宙背景辐射被地面反射进入大气，又被大气衰减，并最终被微波传感器接收的辐射，用（1-小(汍+1^)表示，这三个过程组成了完整的辐射传输方程，如公式(1)所示：[0050]TB=TTa+T?e?Ts+a-ehUk+Tsky)(1)[0051]其中，e为某一通道的地表发射率，无量纲;Tb为通道亮温，单位为K厂是大气透过率，无量纲;TTa是大气上行亮温、Tk是大气下行亮温，单位均为K;Ts为地表温度，单位为K;Tsky是宇宙背景辐射温度，单位为K，可视为定值2.73K;[0052]由于TTa和T、相差很小，可以假设二者相等，因此公式（1)可以进一步表示为公式(2)：[0053]TB=Ta+x?e?Ts+(l-e)x(Ta+2.73)(2)[0054]虽然地表温度作为土壤的固有属性是唯一的，但是不同的通道具有独立的发射特性以及大气衰减特性，因此1和1以及e对于不同通道而言都是不同的，因此通道亮温Tb对不同频率而言也是不同的;假设有两个不同的频率，它们的辐射传输方程可表示为：[0055]TB,i=Ta,i+Ti?￡i?Ts+(l-ei)Ti(Ta,i+2.73)(3)[0056]TBj=Ta>j+Tj?ej?Ts+(l-ej)Tj(Taj+2.73)(4)[0057]其中，i和j分别代表两个不同的频率；[0058]Matzler等（1994)通过观测数据发现22GHz附近频率的发射率基本相等，Weng和Grody(1998)也曾假设19GHz和22GHz两个频率的发射率相等。本发明根据前人的经验，假设18.7GHz和23.8GHz两个频率的发射率相等，因此公式(3)和(4)可分别进一步表示为：[0059]TB，18.7=Ta，18.7+Tl8.7?e?Ts+(l-e)Tl8.7(Ta，18.7+2.73)(5)[0060]TB，23.8=Ta，23.8+T23.8?e?Ts+(l-e)T23.8(Ta，23.8+2?73)(6)[0061]将公式(5)和(6)右边的第一项和第三项均移到左边，它们的左右两边再同除以各自频率的T可得：[0064]将公式(7)和(8)联立，得到公式(9):[0066]由此，地表温度这个变量被约去，因此本发明不需要地表温度产品作为模型的输入参数，可以在减少输入数据的同时有效避免传统方法在有云情况下不能获得地表温度进而不能获得地表发射率的缺点;对公式(9)进行进一步变形，将发射率e单独放到公式左侧，其他项放到公式右侧，得到公式(10):[0068]先前学者（Liu等，2013;Westwater等，1990)已经发现了t和Ta与大气总可降水量(precipitablewatervapor,PWV)呈非常好的线性关系，即：[0069]x=a?PWV+b(11)[0070]Ta=c?PWV+d(12)[0071]其中，a，b，c，d均是线性方程的系数；由此可知，公式（10)中的1^23.8*丁18.7-Ta,18.7?T23.8和T18.7?123.8(1'￡1,18.7-1'￡1,23.8)都可以表不成只有?1￥-个未知数的方程;根据热动力初估检索数据集(ThermodynamicInitialGuessRetrieval,TIGR)提供的全球大气温湿度廓线数据，使用大气辐射传输模型计算了18.7GHdP23.8GH4h18.7、T23.8、Ta,18.7、Ta,23.8以及对应的PWV，发现可以用线性方程来描述Ta,23.8?Ti8.7-Ta,18.7?T23.8、Tl8.7?T23.8(Ta,18.7-Ia,23.8)和PWV的关系，即：[0072]Ta,23.8?Tl8.7-Ta,18.7?T23.8=e?PWV+f(13)[0073]118.7?T23.8(Ta,18.7-Ta>23.8)=g*PWV+h(14)[0074]其中，e，f，g，h均是线性方程的系数；[0075]根据以上分析，公式（10)可进一步变形为一个只包含通道亮温和PWV的方程，即只包含被动微波成像辐射仪的18.7Hz和23.8GHz的H极化通道亮温和大气总可降水量的被动微波地表发射率估算模型，如公式(15)所示：[0077]其中，e为地表发射率，无量纲；Tb,18.7和Tb,23.8分别为被动微波成像辐射仪的18.7GHz和23.8GHz两个通道的H极化亮温，单位均为K;PWV为大气总可降水量，单位为g/0112;1]11、1]12、1]13、1114、1]15、1]16、1]17、1]18均为方程系数。[0078]II、方程系数率定:首先要针对研究区的实际情况，对模型中的八个系数nu、m2、m3、m4、m5、m6、m7、m8进行率定，以确保上述估算模型能够拥有较高的精度;方程系数率定的具体操作方法包括以下三个步骤：[0079]a、采集研究区的以下三类历史数据:被动微波成像辐射仪的18.7和23.8GHz的H极化通道亮温数据，该类数据可以从带有被动微波成像辐射仪的卫星接收站点和数据共享网站直接获得;与亮温数据时间一致的大气总可降水量的地面测量数据，该类数据可以从地面观测网点的测量数据直接获得，测量方法有地面气象站点释放的无线电探空气球、地基GPS、地基微波辐射计、太阳光度计;与亮温数据时间一致的被动微波地表发射率测量数据，这类数据可以从地面测量获得；[0080]b、同一时间的18.7GHz和23.8GHz的H极化亮温数据、大气总可降水量数据、被动微波地表发射率数据，三类数据均齐全算作一组，采集的历史数据应保证不少于八组（因为有八个系数要率定，至少要建立八个方程才能求解，这里采集的历史数据越多则模型精度越尚）；[0081]c、采用最小二乘法率定系数:将采集的历史数据代入公式（15)，建立不少于八个方程的方程组，采用最小二乘法计算出模型的八个系数值。由于最小二乘法有很多种，本发明采用麦夸特法(Levenberg_Marquardt，LM)作为最优选的方法，因为该方法的计算时间短、精度高。[0082]m、数据采集:采集研究区内待研究时间段的被动微波成像辐射仪的18.7GHz和23.8GHz的H极化通道亮温数据，该类数据可以从带有被动微波成像辐射仪的卫星的接收站点和数据共享网站直接获得;采集研究区内待研究时间段的与亮温数据时间一致的大气总可降水量的测量数据，该类数据可以从地面观测网点的测量数据直接获得，测量方法有地面气象站点释放的无线电探空气球、地基GPS、地基微波福射计、太阳光度计；同一时间内的卫星亮温和大气总可降水量这两类数据齐全算作一组。[0083]IV、被动微波地表发射率估算:将步骤n中采集的数据代入经过步骤I中已经率定好参数的估算模型中，直接计算出与输入数据的采集时间同时的瞬时被动微波地表发射率。[0084]传统的基于辐射传输方程的被动微波地表发射率估算方法，需要大量的辅助数据，其中部分辅助数据并不容易获得，使得传统方法的应用受到限制，也进一步造成了传统方法只有月尺度产品而无瞬时产品的现状。本发明跟现有的发射率估算方法相比，其优点主要为：1)输入参数少，只需要输入通道亮温和大气总可降水量，而且这两个参数容易获得，即使是在有云的天气情况下也能获得，因此本发明可以全天候使用，可以保证获得瞬时地表发射率;2)本发明基于辐射传输模型，具有坚实的物理基础，普适性更强，可以应用于不同的地面情况和不同的季节。[0085]上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本
技术领域：
的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或替换，也均属于本发明的保护范围。【主权项】1.一种从被动微波数据估算瞬时地表发射率的方法，其特征在于:所述方法包括以下四个步骤：I、估算模型的推演：一般星载被动微波传感器所接收到的辐射信息包含三个部分:大气向上辐射，用TTa表示；由地面向上发射并被大气衰减的辐射，用τ·ε·Ts表示;大气向下与宇宙背景辐射被地面反射进入大气，又被大气衰减，并最终被微波传感器接收的辐射，用(1-ε)τ(?ν?^)表示，这三个过程组成了完整的辐射传输方程，如公式(1)所示：Τβ=ΤΤ3+τ·ε.Ts+(l-e)T(Tia+Tsky)(1)其中，ε为某一通道的地表发射率，无量纲;TB为该通道亮温，单位为K;τ是大气透过率，无量纲;TTa是大气上行亮温、1^是大气下行亮温，单位均为Κ;Ts为地表温度，单位为Κ;Tsky是宇宙背景辐射温度，单位为K，可视为定值2.73K;由于TTa和1^相差很小，因此假设二者相等，将公式（1)进一步表示为公式(2):Τβ=Τ3+τ·ε.Ts+(l-e)T(Ta+2.73)(2)由于不同的通道具有独立的发射特性以及大气衰减特性，1和1以及ε对于不同通道而言都是不同的，因此通道亮温Τβ对不同频率而言也是不同的；假设有两个不同的频率，将它们的福射传输方程表不为：TB，i=Ta,i+Ti·ε?·Ts+(l-￡i)Ti(Ta,i+2.73)(3)Tb,j=Ta,j+Tj·εj·Ts+(1_εj)Tj(Ta,j+2·73)(4)其中，i和j分别代表两个不同的频率;根据前人的经验，假设18.7GHz和23.8GHz两个频率的发射率相等，由此将公式(3)和公式(4)分别进一步表示为：TB，18.7=Ta，18.7+Tl8.7·ε·Ts+(l-e)Ti8.7(Ta，18.7+2.73)(5)TB，23.8=Ta,23.8+T23.8·ε·Ts+(l_￡)T23.8(Ta,23.8+2.73)(6)将公式(5)和公式(6)右边的第一项和第三项均移到左边，它们的左右两边再同除以各自频率的τ得到公式(8)和公式(9):将公式(7)和公式(8)联立，得到公式(9):由此，地表温度这个变量被约去；进一步对公式（9)变形，将发射率ε单独放到公式左侦U，其他项放到公式右侧，得到公式(10):先前学者已经发现了τ和Ta与大气总可降水量(precipitablewatervapor，PWV)呈非常好的线性关系，即：i=a·PWV+b(11)Ta=c·PWV+d(12)其中，a，b，c，d均是线性方程的系数；由此可知，公式（10)中的Ta,23.8·T18.7-Ta,18.7·丁23.8和118.7*123.8(1'￡1,18.7-1'￡1,23.8)都能表不成只有?1￥-个未知数的方程;根据热动力初估检索数据集提供的全球大气温湿度廓线数据，利用大气辐射传输模型计算了18.7GHz和23.86他的118.7、123.8、!^18.7、!^23.8以及对应的13評，用线性方程描述1^23.8.118.7-!^18.7·Τ23.8、τ?8.7·Τ23.8(Τ3，18.7-Τ3,23.8)^ΡΙΡ￥ν|^]**，·Ρ:Ta,23.8*Ti8.7~Ta,18.7*^23.8=θ·PWV+f(13)^18.7*T23.8(Ta,18.7-Ta,23.8)=g*PWV+h(14)其中，e，f，g，h均是线性方程的系数；根据以上分析，将公式（10)进一步变形为一个只包含通道亮温和PWV的方程，即只包含被动微波成像辐射仪的18.7Hz和23.8GHz的Η极化通道亮温和大气总可降水量的被动微波地表发射率估算模型，如公式(15)所示：其中，ε为地表发射率，无量纲;ΤΒ,18.7和ΤΒ,23.8分别为被动微波成像辐射仪的18.7GHz和23.86抱两个通道的!1极化亮温，单位均为1(;?￥￥为大气总可降水量，单位为8/〇]12;1]11、1]12、1113、1114、1115、1116、1117、1118均为方程系数；Π、方程系数的率定:针对研究区的实际情况，对模型中的八个系数mi、m2、m3、m4、m5、m6、m7、ni8进行率定，确保上述估算模型在特定研究区拥有较高的精度；ΙΠ、数据采集：采集研究区内待研究时间段的被动微波成像辐射仪的18.7GΗζ和23.8GHz的Η极化通道亮温数据，该类数据的获得途径包括带有被动微波成像辐射仪的卫星的接收站点、数据共享网站;采集研究区内待研究时间段的与亮温数据时间一致的大气总可降水量的测量数据，该类数据的获得途径包括地面观测网点的测量数据，测量方法包括地面气象站点释放的无线电探空气球、地基GPS、地基微波福射计、太阳光度计；同一时间内的卫星通道亮温和大气总可降水量这两类数据齐全算作一组；iv、被动微波地表发射率估算:将步骤m中采集的数据代入经过步骤π率定好参数的估算模型中，直接计算出与输入数据的采集时间同时的瞬时被动微波地表发射率。2.根据权利要求1所述的从被动微波数据估算瞬时地表发射率的方法，其特征在于:所述方程系数率定的具体方法如下：a、采集研究区的以下三类历史数据:被动微波成像辐射仪的18.7GHz和23.8GHz的Η极化通道亮温数据，该类数据的获得途径包括带有被动微波成像辐射仪的卫星接收站点、数据共享网站;与亮温数据时间一致的大气总可降水量的地面测量数据，该类数据的获得途径包括地面观测网点的测量数据，测量方法包括地面气象站点释放的无线电探空气球、地基GPS、地基微波辐射计、太阳光度计；与亮温数据时间一致的被动微波地表发射率测量数据，该类数据从地面测量获得；b、同一时间的18.7GHz和23.8GHz的Η极化亮温数据、大气总可降水量数据、被动微波地表发射率数据，这三类数据均齐全算作一组，采集的历史数据应保证不少于八组；C、采用最小二乘法率定系数:将采集的历史数据代入公式（15)，建立不少于八个方程的方程组，采用最小二乘法计算出模型的八个系数值。【文档编号】G06F17/12GK105930664SQ201610266664【公开日】2016年9月7日【申请日】2016年4月26日【发明人】周芳成,吴骅,李召良,段四波,唐伯惠,唐荣林【申请人】中国科学院地理科学与资源研究所,中国农业科学院农业资源与农业区划研究所