Integrating Vegetation Indices Models and Phenological Classification with Composite SAR and Optical Data for Cereal Yield Estimation in Finland (Part I)

Abstract

During 1996–2006 the Ministry of Agriculture and Forestry in Finland, MTT Agrifood Research Finland and the Finnish Geodetic Institute carried out a joint remote sensing satellite research project. It evaluated the applicability of composite multispectral SAR and optical satellite data for cereal yield estimations in the annual crop inventory program. Three Vegetation Indices models (VGI, Infrared polynomial, NDVI and Composite multispectral SAR and NDVI) were validated to estimate cereal yield levels using solely optical and SAR satellite data (Composite Minimum Dataset). The average R² for cereal yield (y_b) was 0.627. The averaged composite SAR modeled grain yield level was 3,750 kg/ha (RMSE = 10.3%, 387 kg/ha) for high latitude spring cereals (4,018 kg/ha for spring wheat, 4,037 kg/ha for barley and 3,151 kg/ha for oats).

1. Introduction

Remote sensing has been extensively applied in world crop production estimations by the European Union, the United States of America and the Food and Agriculture Organization (FAO) of the United Nations. The advantage of applying satellite based remote sensing data is global coverage and equally calibrated data, which enables temporal and spatial comparison over the years in monitoring areas.

New precision farming techniques in agriculture can be used more accurately to assess crop growth and soil conditions during growing seasons [1,2,3,4,5]. New high accuracy GPS satellites used for position control, e.g., the new EU-funded GALILEO system [6], and new hyperspectral spectrometers operating over a wide range of wavelengths (λ = 400–2,400 nm) provide new possibilities and tools for cereal cultivation [1,7,8,9]. Major climatic events, such as global drought periods during the forthcoming climate change can be monitored with new generation satellites of the enhanced spectral resolution [10,11] . The Joint Research Centre of EU has applied JRC-FAPAR index (Fraction of Absorbed Photosynthetically Active Radiation) for monitoring large area vegetation stress states and biomass productivity during Pan-European drought periods [12,13,14,15,16,17]. In addition, JRC-FAPAR model with Radiative Transfer models (RT) and using SeaWiFS data from the OrbView-2 (AKA SeaStar) satellite has bee applied to identify key vegetation phenology events (start, end, length) during growing season.

Recently Harrison et al. [18] used the AFRCWHEAT2 crop model for scaling-up wheat phenological development in Europe. Moulin et al. [19] applied SPOT satellite data to validate the AFRCWHEAT2 and SAIL reflectance models to estimate wheat yield and biomass NPP (Net Primary Production).

In northern latitudes the cloudiness during the growing season seriously reduces the applicability of optical spectrum satellite data . However, new generation microwave based satellite systems (e.g., COSMO-SkyMed/Italy, ENVISAT/ESA/EU, ALOS/PALSAR, Radarsat-2 /Canada, TerraSAR-X/ Germany can map images through cloud coverage and during the night-time. The Japanese ADEOS II (Miradori II) and GOSAT/IBUKI (Greenhouse gases Observing SATellite) satellites can also measure atmospheric aerosol columns, vertical greenhouse gas profiles (CO₂, CH_4, NO₂, N₂O_, O_3, CFC-11, CFC-12) and water vapour . Global greenhouse gas estimates can be used in large area climate change studies to assess changes in biomass NPP production and cereal yield production [19, 20]. Recently multispectral remote sensing studies have focused on integrating both optical and microwave SAR data [1,21]. Recently McNairn et al. [1] integrated multispectral classification techniques with SAR and optical data in Canada. Cereal, grass and oil crop fields in four prairie testing sites were classified with Neural Network and Decision trees methodologies.

In Finland Kuittinen [22], Kuittinen et al. [23] and Karvonen et al. [24] have reviewed the yield estimation of spring wheat and barley (Hordeum vulgare L.) combined with Landsat and SPOT data for Nordic high latitude agricultural regions. Yara Co. (formerly Kemira GrowHow Finland) [25] has established in Finland a commercial Kemira Loris™ (LOcal Resource Information System) integrated expert system for farmers. Loris utilizes GPS based precision farming techniques with optimum fertilization and cultivation practises combined with infrared aerial photographs.

During 1996–2006 the Ministry of Agriculture and Forestry in Finland (MAFF), MTT Agrifood Research Finland (MTT), the Finnish Geodetic Institute (FGI) and the Technical Research Centre of Finland (VTT, Dept. of Space Research) performed a joint remote sensing satellite research project combined with ground truth measurements to evaluate the applicability of multispectral SAR and optical satellite data for cereal yield estimations in the annual crop inventory program. Both optical and composite SAR microwave Vegetation Indices models were validated to estimate spring cereal yield levels using solely optical and SAR satellite data (Composite Minimum Dataset).

In summary, new composite multispectral optical and SAR classification and modeling techniques provide new integrated tools for cereal yield estimations in national inventory programs [1,17,26,27,28]. In Finland the Ministry of Agriculture and Forestry performs the annual ground based crop inventory sampling to estimate cereal yield production and cultivation areas on national level. New classification and modeling techniques integrated with digitized Finnish Land Parcel Identification System (FLPIS) and Integrated Administration and Control System (IACS) can be used to support Finnish national crop inventory programs [29,30]. More accurate cereal inventory estimates will also support Finnish national CAP-policy (Common Agricultural Policy) goals and objectives on EU scale. These aspects are reviewed in Part II of this publication.

The aim of the present study was to estimate actual non-potential grain yield levels for high latitude spring cereals (spring wheat, barley and oats, Avena Sativa L.) in large area field conditions in southern Finland. The cereal theoretical maximum yielding capacity is limited by environmental and vegetation stresses (e.g., drought periods, nutrient deficiencies, pathogen epidemics) during growing season in actual field growing conditions. These stress factors result to reduced non-potential baseline yield levels (y_b, kg/ha) on field parcel level.

The objectives of the present study were: (i) to construct a dynamic SatPhenClass phenological classification model, which classifies both optical and SAR satellite data based on cereal actual phenological development in both vegetative and generative phases (ii) to calibrate and validate multispectral Composite Vegetation Indices (VGI) models, which integrate both phenologically pre-classified optical (Models I–II) and microwave SAR data (Composite SAR and NDVI Model III), and finally (iii) VGI models were used to estimate cereal non-potential baseline yield (y_b) levels in growing zones (I–IV) in southern Finland during 1996–2006.

2. Methodology and Study Area

2.1. System Analysis and Strategy

The research applied in this study focused on utilizing multispectral optical and microwave SAR data extracted from experimental sites in southern Finland. A dynamic SatPhenClass phenological algorithm was developed to classify measured satellite data based on phenological development of different spring cereals. The pre-classified SAR and optical satellite data were used as input for VGI models. The system analysis [31,32] of this study is depicted in Figure 1.

Figure 1. System Analysis and design diagram by satellite systems and VGI models [31,38].

Figure 1. System Analysis and design diagram by satellite systems and VGI models [31,38].

The analysis strategy applied in this study is depicted in phases I–XI. Phases I–II (Section 2.2) illustrate satellite and ground truth data pool consisting of solely SAR and optical satellite data (Composite SAR and Optical Minimum Datasets 1996–2006). Correspondingly phase III (Section 2.4) explains geometric and radiometric calibration of SAR and optical data. Phases IV–V (Section 2.3.1, Section 3.2.1) illustrate SatPhenClass phenological classification model and LAI and ETS ground truth sampling methodology. Phases VI–VII (Section 2.3.2, Section 3.2.2) explain calibration of VGI models using Minimum Datasets without extensive agrometeorological ground truth data required by conventional dynamic crop models [33,34,35,36]. Phases VIII–IX (Section 2.5, Section 3.3) explain the validation of VGI models. Phases X–XI(a,b) (Section 3.2, Section 3.3) explain the estimation of baseline grain yield levels (y_b) for different spring cereals resulting outputs for SAR & NDVI models (Phase XI(a), Part I) and for detailed optical VGI models (Phase XI(b), Part II).

In addition, Phase XI(c) scheme presents an alternative cereal specific LAI-bridge coupling system [23] reviewed in Part II. LAI-bridge coupling system utilizes optical satellite data as input data for dynamic crop models. Hodges [33] and Bowen [36] have reviewed the principles of dynamic crop models explaining phenological development and crop physiological processes (root and soil systems, photosynthesis, respiration and translocation of assimilates to heads). Dynamic crop models require extensive agrometeorological and phenological ground truth data for daily integration of biomass and yield accumulation [34,35].

In LAI-bridge coupling system, the GEMI-index (Global Environment Monitoring Index, [37]) is applied to estimate cereal LAI development (LAI_satellite) from optical reflectance data during growing season. A similar model, VGI GEMI (model IV) is used in this study in Part II. The LAI-bridge coupling method applies iterative Kalman filter (KF) algorithm [38] commonly applied in System Control theory, Markov chains and Bayesian estimations. The GEMI based LAI estimates (LAI_Satellite) and LAI values estimated by crop models (LAI_CropModel) are corrected in Kalman filter iterations with LAI_Satellite and LAI_CropModel variances to obtain optimized LAI (_Final) values for non-potential baseline yield (y_b) calculations. LAI-bridge coupling method was originally developed by Karvonen et al. [24] using SPOT and Landsat data as input for PotCropF [24] dynamic crop model used in Finnish cereal yield estimations. Later on, Kuittinen et al. [23] applied LAI-bridge coupling system with GEMI index by using SPOT and NOAA optical data as input for CropWatN [39] and WOFOST [40] crop models.

In this study (Part I) three Vegetation Indices (VGI) models were calibrated and validated with optical (NDVI) and microwave SAR satellite and ground truth (1996–2006) data. The satellite data were measured in five different experimental locations in southern Finland (Table 1, Figure 2). The detailed analysis strategy (Figure 1, Figure 3a,b) consisted of following specific procedures:

(a) Pre-classification of SAR and optical data with the SatPhenClass phenological classification algorithm for spring cereals (Figure 3a, b). The SatPhenClass algorithm exploits BBCH [41,42] and Zadoks [43] growth scales (Table 2) with four phenological classes for spring cereals (a_p: BBCH 0–12, b_p: BBCH 12–50, c_p: BBCH 50–90, d_p: BBCH > 90). Class a_p corresponds to development period between sowing and two leaf stage with double ridge formation. Class b_p corresponds to period between two leaf stage and ear emergence with maximum Leaf Area Index (LAI_max) exposure with fully closed canopy. Class c_p corresponds to period between ear emergence and anthesis with grain filling until full maturity. Finally, class d_p corresponds to senescence and post-harvest phases.

(b) Calibration of Composite Vegetation Indices (VGI) models (models I–III in Part I, models IV–V in Part II) with pre-classified C-band SAR and optical data. Both linear and non-linear polynomial response functions in VGI models were applied for cereal baseline yield (y_b) estimations.

(c) Validation of each calibrated VGI model using two independent data sources separately: (i) the MTT Agrifood Research Finland official variety trial data (1996–2006) and (ii) the MAFF (Ministry of Agriculture and Forestry in Finland) annual crop yield inventory data (1996–2006). MTT and MAFF validation datasets take into account temporal and spatial variation between years and locations.

(d) Estimating non-potential baseline grain yield levels (y_b) for different spring cereals in actual field growing conditions in Southern Finland with calibrated VGI models and optical and SAR data.

Table 1. Sensors, datasets, growing zones, locations and soil classification [44].

Sensor Type (1)	Publ. part	Data set	Growing Zone (2)	MAFF Agriculture Advisory Centre (Figure 2 and Figure 4)	Location (Figure 2)	Coarse soils (%)					Clay soils (%)		Organic & mould (%)
						Gravel with coarse sand	Fine sand	Coarse sand	Silt	Sandy clay	Silt clay	Gyttja clay 3)
SAR & Optical	I/II	1.1	III (IV)	Etelä-Pohjanmaa	Lapua 23° 10' E, 62° 50' N		11.2	19.8		23.5
	I/II	1.2	III (IV)	Etelä- Pohjanmaa	Seinäjoki 23° 10' E, 62° 50' N		15.5	10.1		20.5			25.3
	I/II	1.3	III (IV)	Etelä- Pohjanmaa	Ilmajoki 23° 10' E, 62° 50' N		13.2	12.0	3.9	21.4
	I/II	2.1	I	Nylands Svenska	Porvoo 25° 50' E, 60° 50' N	6.2	17.8	8.5				38.0	10.9
Optical	II	3.1	I	Nylands Svenska	Kirkkonummi 24° 30'E, 60° 10’N		9.4			27.5	8.5	41.4	5.4
	II	3.2	II	Häme	Jokioinen 23° 50' E, 60° 50' N (Kuuma Exp. Area)					56.0	7.0	15.1	7.7 (70–80)
	II	3.3	II	Häme	Mellilä 22° 20'E, 60° 50' N	7.4	8.3	13.2		28.6		36.0

(e) ⁽¹⁾ For abbreviations refer to

Table 7. Abbreviations used.

Definition, abbreviation	Unit, (range)	Description
X		Mean of sample
Sd		Standard deviation of sample (n)
Cv		Coefficient of variation (%) = Sd/X
Sx		Standard error of mean = $\frac{\mathrm{Sd}}{\sqrt{\text{n}}}$
MSE		Mean squared error $\mathrm{MSE}\left(\overline{X}\right) = \text{E}\left({(\overline{X}-\mu )}^2\right)={\left(\frac{\sigma }{\sqrt{n}}\right)}^2$
RMSE		Root Mean Square Error , square root of MSE
R2		Coefficient of determination
LSE		Least-Square Estimation-algorithm
VGI with submodels (I-III)		Vegetation Indices submodels (I-III: I—Infrared polynomial, II—NDVI, III—composite NDVI and backscattering model (Table 6)
ap, bp, cp, dp		DVS (Phenological Development stage) four classification values in SatphenlClass algorithm (Table 3 a,b), used in VGI models: ap—vegetation stage class from emergence until 2 leaf and double ridge stages, bp—generative stage class until heading, cp—grain filling stage in generative phase between anthesis and full maturity, dp—senescence phase (Used only in microwave polynomial model III)
Rf (Ch,month)		Landsat or SPOT calibrated reflectance values with index denoting channel and month during growing season (a–May, b–June, c–July); used in VGI Infrared model (I) as independent variables for crop*cultivar*soil covariance interaction and yield estimations (Table 3, Models. 1.1–7.3, Table 12 App. C), * - general notation for covariance effects
DYMTT		Yield Difference (kg/ha) modelled (VGI)—observed MTT average (Table 11,14)
DRTMTT		Yield Difference Ratio, Modelled (VGI)/MTT Observed; over/underestimation (%) from the reference (100%) (Table 11, 14)
DYMAFF		Yield Difference (kg/ha) modelled (VGI)—obs. MAFF average (Table 10, 13).
DRTMAFF		Yield Difference Ratio, modelled (VGI)/MAFF Observed; over/underestimation (%) from the reference (100%) (Table 10,13)
Oe / ue		Overestimated / Underestimated by the corresponding VGI model vs. observed (Table 12, Table 13 and Table 14).
spc., Cv.		species, cultivar
spc.		Spring sown cereal species (spring wheat, barley, oats)
Swh		Spring wheat (Triticum aestivum L.) including cv. Heta, Kadett, Manu, Reno, Ruso, Satu, Tjalve
Brl		Barley (Hordeum vulgare L.) including cv. Arra, Arttu, Artturi, Arve, Eero, Ida, Inari, Kustaa, Kymppi, Loviisa, Mette, Pohto, Pokko
Oats		Oats (Avena Sativa L.) including cv. Aarre, Salo
loc.		Loc—Location: J-Jokioinen, K—Kirkkonummi, M-Mellilä, P—Porvoo, L—Lapua
Tb	degree (°C)	Threshold temperature
Dd	degree days
ETS(Tb)	dd—degree days	Cumulative temperature sum over threshold temperature (Tb = 5°)
PAR	MJ/d/m2 (10–20)	Photosynthetically Active Radiation (λ = 400–700 nm)
IR NIR—Near IR Mid IR Thermal IR	MJ/d/m2	infrared radiation (IR), λ = 630–690 nm near infra, λ = 760–900 nm mid infra, λ = 1.55–1.75 μm thermal IR, λ = 10.4–12.5 μm
Rf	(0.0– 1.0) Optical (λ = 400–700 nm) and infrared sensors (λ = 630–12.5 μm).	Reflectance; reflected radiation from soil and vegetation canopies and measured by optical satellites [22,23,57,58,68,69]
σ0 (sigma zero)	(–20–10 dB). Calibrated SAR (Synthetic Aperture radar) backscattering signal with microwave 5.4 GHz (C-band, λ = 5.7 cm) and 9.8 GHz (X-band) and polarization levels (HH, VV, VH, HV).	Backscatter coefficient (sigma zero) for microwave backscattering signal, which is a combined signal reflected from soil and vegetation canopies [60,61,62,63,26,27,28].
Potential, non-limited yield, yield potential	kg/ha	Modelled maximum yield capacity (kg/ha) for a specific cultivar without limiting environmental stress factors during growing season (vegetation water stress, nutrient deficiencies, pathogen epidemics etc.)
Non-potential, limited yield	kg/ha	VGI modeled yield level (kg/ha) for a specific cultivar with limiting environmental stress factors during growing season reducing maximum yield capacity, see potential yield.
yb(spc, cv, soil type.) spc=swh,brl,oats	kg/ha, 15% moisture content	Baseline yield for spring cereal species (swh, brl, oats).VGI (I-III) modeled cereal yield level (kg/ha) using time series for a specific cultivar under field conditions, see non-potential yield.Index denotes crop, cultivar and soil type; used as dependent variable in VGI models (Models. 1.1–7.3, Table 12, App. C).
NDVI	%	Normalized Difference Vegetation Index
SatphenlClass	BBCH and Zadok’s scaling	Satellite data classification algorithm based on cereal phenology
Minimum dataset		Experimental dataset without ground truth or meteorological data, containing only optical or microwave satellite data
VGI with submodels (I-II)		Vegetation Indices submodels (I-IIII: I—Infrared polynomial, II—NDVI, III—composite NDVI and backscattering model.
MAFF		Ministry of Agriculture and Forestry in Finland
IIASA		The International Institute for Applied Systems Analysis [32]

, Appendix B. ⁽²⁾ Growing zones (I–IV) are depicted in Figure 4a. ⁽³⁾ Gyttja clay contains peat and mud fractions.

Table 1. Sensors, datasets, growing zones, locations and soil classification [44].

Sensor Type (1)	Publ. part	Data set	Growing Zone (2)	MAFF Agriculture Advisory Centre (Figure 2 and Figure 4)	Location (Figure 2)	Coarse soils (%)					Clay soils (%)		Organic & mould (%)
						Gravel with coarse sand	Fine sand	Coarse sand	Silt	Sandy clay	Silt clay	Gyttja clay 3)
SAR & Optical	I/II	1.1	III (IV)	Etelä-Pohjanmaa	Lapua 23° 10' E, 62° 50' N		11.2	19.8		23.5
	I/II	1.2	III (IV)	Etelä- Pohjanmaa	Seinäjoki 23° 10' E, 62° 50' N		15.5	10.1		20.5			25.3
	I/II	1.3	III (IV)	Etelä- Pohjanmaa	Ilmajoki 23° 10' E, 62° 50' N		13.2	12.0	3.9	21.4
	I/II	2.1	I	Nylands Svenska	Porvoo 25° 50' E, 60° 50' N	6.2	17.8	8.5				38.0	10.9
Optical	II	3.1	I	Nylands Svenska	Kirkkonummi 24° 30'E, 60° 10’N		9.4			27.5	8.5	41.4	5.4
	II	3.2	II	Häme	Jokioinen 23° 50' E, 60° 50' N (Kuuma Exp. Area)					56.0	7.0	15.1	7.7 (70–80)
	II	3.3	II	Häme	Mellilä 22° 20'E, 60° 50' N	7.4	8.3	13.2		28.6		36.0

(e) ⁽¹⁾ For abbreviations refer to Table 7, Appendix B. ⁽²⁾ Growing zones (I–IV) are depicted in Figure 4a. ⁽³⁾ Gyttja clay contains peat and mud fractions.

Figure 2. Satellite and ground truth measurement locations in Finland (Original Data (C) NASA, visibleearth.nasa.gov/).

Figure 2. Satellite and ground truth measurement locations in Finland (Original Data (C) NASA, visibleearth.nasa.gov/).

2.2. Overview of Satellite and Ground Truth Sites

A general overview of satellite and ground truth experimental sites is displayed in Figure 2. A detailed description of the experimental sites is given by Kuittinen et al. [23]. Table 1 depicts the microwave SAR and optical data sources with publication division (I–II), growing zones and Agriculture Advisory Centres in Finland. A general scheme of the measuring satellite systems and calibration parameters used in this study is displayed in

Table 12. SAR and optical VGI (I-IV) linear (REG/Stepwise) and non-linear (RSREG) yield models for spring cereals^(1),(2).

Model category	Factor : Crop /Species 1)		Model (I-IV)	R2	RMSE kg/ha	Model equation
I Optical	Swh		1.1 (I)	0.764	282.3	yb(swh) = 4941.9 – 5455.9 × rf3ap – 1351.4 × rf4ap + 957.1 × rf3bp + 656.2 × rf4bp + 4742.1 × rf3cp – 4983.5 × rf4cp
	Swh		1.2 (I)	0.794	42.46	yb(swh) = 985.93 + 13337.46 × rf3ap + 8355.88 × rf4ap – 387.19 × rf4bp + 255.58 × rf3cp
	Swh		1.3 (II, NDVI)	0.737	297.6	yb(swh,NDVI) = 4659.2 + 175.4 × nap – 25.6 × nbp – 3215.5 × ncp + 93.7 × nap2 – 19.1 × nbp × nap – 178.5 × nbp2 – 560.4 × ncp × nap + 3250.4 × ncp × nbp – 864.7 × ncp2
	Swh		1.4 (II, NDVI)	0.732	300.1	yb(swh,NDVI) = 4692.90 + 109.84 × nap – 210.38 × nbp – 969.22 × ncp
	Swh		1.5(5) (III,GEMI)	0.704	316.1	Y(swh,GEMI) = 5074.1 – 766.2 × gap – 143.1 × gbp – 3378.2 × gcp + 4254.1 × gap2 + 1079.8 × gbp × gap – 1326.2 × gbp2 – 7264.8 × gcp × gap + 7717.2 × gcp × gbp – 2622.1 × gcp2
	Swh		1.6(5) (III,GEMI)	0.570	536.8	yb(swh,GEMI) = 5506.31 + 7780.47 × gap – 12478.3 × gbp + 761.96 × gcp
	Swh		1.7(5) (IV,PARND)	0.712	311.6	yb (swh, PARND)=4397.9+2736.9*pap-379.4*pbp-493.8*pcp-4235.1*pap2+ 2701.1*pbp*pap+605.8*pbp2 -2161.6*pcp*pap+883.1*pcp*pbp+493.9pcp2
	Swh		1.8(5) (IV,PARND)	0.509	406.3	Y(swh, PARND)=4502.31+529.36*pap+2377.05*pbp-779.31*pcp
	Brl		2.1 (I)	0.615	449.3.	Y(brl) = 5348.8 – 600.8 × rf3ap – 184.5 × rf4ap – 8562.4 × rf3bp – 2105.6 × rf4bp – 2766.9 × rf3cp – 1556.2 × rf4cp
	Brl		2.2 (II, NDVI)	0.611	449.6	Y(brl,NDVI) = 3481.5 –2 05.5 × nap–312.1 × nbp + 396.5 × ncp – 60.5 × nap2 – 584.5 × nbp × nap + 864.6 × nbp2 + 1081.4 × ncp × nap – 710.5 × ncp × nbp – 232.4 × ncp2
	Brl		2.3(5) (III, GEMI)	0.614	448.6	Y(brl,GEMI)=5184.0 – 2802.9 × gap – 187.4 × gbp – 3457.9 × gcp – 374.8 × gap2 + 91.5 × gbp × gap + 1161.4 × gbp2 + 4817.1 × gcp × gap – 1932.1 × gcp × gbp + 2255.6 × gcp2
	Brl		2.4(5) (IV,PARND)	0.587	463.7	Y(brl,PARND) = 4621.7 + 1852.8 × pap-6418.5 × pbp – 5850.4 × pcp + 48.9 × pap2 – 5952.3 × pbp × pap + 10853.0 × pbp2 –2040.4 × pcp × pap + 12879.0 × pcp × pbp + 6002.9 × pcp2
	Oats		2.5 (I)	0.760	55.0	yb(oats) = 3457.4 – 3762.5 × rf3ap– 2135.8 × rf4ap+6643.1 × rf3bp+ 1566.5 × rf4bp+ 571.4 × rf3cp– 179.1 × rf4cp
II Optical Species*soil Covariance	Crop /Species 1)		Model	R2	RMSE	Model equation
	Swh * sandy clay		2.6 (I)	0.764	282.3	yb(swheat,clay) = 6765.5 – 38407.0 × rf3ap – 14979.0 × rf4ap– 23698.0 × rf3bp + 9552.7 × rf4bp + 7261.4 × rf3cp – 22022.0 × rf4cp
	Brl * sandy clay		2.7 (I)	0.166	1382	yb(barley,clay) = 7141.4-3842.1 × rf3ap-2612.4 × rf4ap – 3032.1 × rf3bp-11708.0 × rf4bp – 28637.0 × rf3cp – 1636.2 × rf4cp
III Optical Species* Cultivar Cov.	Swh* cv. Manu		3.1 (I)	0.089	1292	yb (cv. Manu*clay) = 5696.5 + 2293.8 × rf3bp + 5387.9 × rf3cp – 736396.0 × (rf3bp)2
	Swh * cv. Satu		3.2 (I)	0.046	1031	yb (cv. Satu*clay) = –9798.2 + 736440.0 × rf2bp – 9925014.0 × (rf2bp)2
	Brl * cv. Inari		3.3 (I)	0.144	1220	yb (cv. Inari*clay) = 8336.8 – 71506.0 × rf3bp – 45627.0 × (rf3bp)2
IV Microwave SAR	Sensor	Cereal specie	Model4)	R2	RMSE	Model equation
	ERS SAR(4)	Brl	4.1(III)	0.448	482.7	yb (brl, ERS2) = 4345.7 + 109.4 × NDVIap – 211.6 × NDVIbp – 983.3 × NDVIcp–0.57 × VV(5GHz,cp) + 5.61 × VV(5GHz,dp)
		Oats	4.2(III)	0.417	584.2	yb (oats, ERS2) = 3739.5 + 108.9 × NDVIap – 212.1 × NDVIbp – 938.8 × NDVIcp – 0.47 × VV(5GHz,cp) + 4.28*VV(5GHz,dp)
	Radarsat SAR4)	Swh	5.1(III)	0.731	300.8	yb (swh, Radarsat) = 4690.7 + 111.8 × NDVIap – 213.4 × NDVIbp – 982.6 × NDVIcp – 2.69 × HH(5GHz,cp) + 3.9 × HH(5GHz,dp)
		Brl	5.2(III)	0.702	322.8	yb(brl,Radarsat) = 4430.1 + 109.4 × NDVIap – 211.6 × NDVIbp – 983.3 × NDVIcp – 0.52 × HH(5GHz,cp) + 5.07 × HH(5GHz,dp)
		Oats	5.3(III)	0.624	483.6	yb (oats, Radarsat) = 3843.3 + 108.9 × NDVIap – 212.1 × NDVIbp – 983.8 × NDVIcp – 0.47 × HH(5GHz,cp) + 4.03 × HH(5GHz,dp)
	Envisat ASAR(4)	Swh	6.1(III)	0.723	302.1	yb(swh,Envisat) = 4701.1 + 108.2 × NDVIap – 208.8 × NDVIbp – 983.1 × NDVIcp – 3.9 × VH(5GHz,cp) + 17.4 × VV(5GHz,cp) – 3.1 × VH(5GHz,dp)+ 5.2 × VV(5GHz,dp)
		Brl	6.2(III)	0.694	349.8	yb (brl, Envisat) = 4261.4 + 109.4 × NDVIap – 211.6 × NDVIbp – 983.3 × NDVIcp – 4.59 × VH(5GHz,cp) + 18.24 × VV(5GHz,cp) – 4.04 × VH(5GHz,dp) + 6.15 × VV(5GHz,dp)
		Oats	6.3(III)	0.617	389.7	yb (oats, Envisat) = 3635.7 + 108.9 × NDVIap – 212.8 × NDVIbp– 983.8 × NDVIcp – 2.59 × VH(5GHz,cp) + 16.46 × VV(5GHz,cp) – 2.03 × VH(5GHz,dp) + 4.05 × VV(5GHz,dp)
	HUTSCAT Scattero meter(3),(4)	Swh	7.1(III)	0.582	416.8	yb(swh,HUTSCAT) = 4258.4 + 109.4 × NDVIap – 198.2 × NDVIbp – 937.4 × NDVIcp + 5.2 × VV(5GHz,cp) + 18.4 × HH(5GHz,cp) – 2.9 × VH(5GHz,cp) – 16.4 × HV(5GHz,cp) + 4.4 × VV(5GHz,dp) + 12.4 × HH(5GHz,dp) – 2.3 × VH(5GHz,dp) – 14.4 × HV(5GHz,dp)
		Brl	7.2(III)	0.518	490.1	yb (brl, HUTSCAT) = 4294.2 + 107.2 × NDVIap– 209.2 × NDVIbp – 928.2 × NDVIcp + 3.2 × VV(5GHz,cp) + 17.4 × HH(5GHz,cp) – 3.9 × VH(5GHz,cp) – 15.4 × HV(5GHz,cp) + 5.4 × VV(5GHz,dp) + 11.4 × HH(5GHz,dp) – 4.3 × VH(5GHz,dp) – 15.4 × HV(5GHz,dp)
		Oats	7.3(III)	0.424	544.2	yb(oats,HUTSCAT) = 3782.5 + 106.8 × NDVIap – 207.2 × NDVIbp – 942.5 × NDVIcp + 4.2 × VV(5GHz,cp) + 15.2 × HH(5GHz,cp) – 5.9 × VH(5GHz,cp) – 17.1 × HV(5GHz,cp) + 4.8 × VV(5GHz,dp) + 11.8 × HH(5GHz,dp) – 4.3 × VH(5GHz,dp) – 12.2 × HV(5GHz,dp)

⁽¹⁾ For abbreviations refer to Table 9 ⁽²⁾ Independent variables classified with SatPhenClass-algorithm (Figure 3a,b). ⁽³⁾ HUTSCATused for calibration verification purposes only (helicopter mounted) ⁽⁴⁾ Only in Part I (SAR+Optical models) ⁽⁵⁾ Only in Part II (Optical models)

and 13 (Appendix C). The abbreviations applied in this study are given in Table 7 (Appendix B).

The satellite measurement and ground truth experiment sites (Lapua, Kirkkonummi , Jokioinen , Mellilä and Porvoo) with corresponding soil classifications [44] are depicted in Table 1 and

Table 9. Satellite measurements in Finnish experimental locations [23,26,28].

Location	Date (Satellite measurement)	Satellite type, sensor and Image no.	Incidence angle (φ) –1)	Soil type4)
Lapua, Seinäjoki 23° 10 ' E, 62° 50' N	6.7.1994, 24.6.1995, 1996: 21.5., 3.7.	Optical: SPOT XS 65 220	–21.3, –4.3, 10, 0.2	Clay (Sandy clay 56%), coarse (Coarse sand 33.5%)
	1997: 17.5.,7.6.,1.7	Optical: SPOT XS 65 220–221	9.3,5.9,–23.8
	1996: 17.6, 30.7, 25.8, 25.9	Microwave SAR: ERS1 SAR5), f = 5.3 GHz, C band, (λ = 5.7 cm), VV polarization
	2001: 6.5,16.5, 30.5, 6.6, 16.6, 23.6.,10.7, 17.7, 24.7, 3.8, 17.8, 27.8, 10.9, 17.9, 20.9, 27.9, 28.10.	Microwave SAR: Radarsat1 SAR, f = 5.3 GHz, C band, HH polarization
	2003–2004: 15.6,18.6,21.6,28.6,4.7,7.7,14.7,23.7,02.8,11.8,24.8,15.9 2005–2006: Only ground truth observations 5)	Microwave ASAR: ENVISAT ASAR , f = 5.3 GHz, C-band), VV, HH, VV/HH, HV/HH, or VH/VV polarizations
Kirkkonummi 24° 30' E, 60° 10' N	1994: 6.6., 7.7	Optical: SPOT XS 73 227	24.5, –21.6	Clay (2) (gyttja clay 41.4%)
	25.7.1994, 3.6.1995	Optical: Optical: LANDSAT TM 189 18	23
	15.6.1995, 13.7.1994	Optical: SPOT XS 73 227, 69 225	–2.2
Jokioinen 23° 50' E, 60° 50' N	25.7.1994	Optical: LANDSAT TM 189 18	–23.4	Clay (Sandy clay 56%), Kuuma exp. area (70–80% organic top layer, Table 3)
Mellilä 22° 20' E, 60°, 50' N	1989:24.5., 25.6., 27.7.	Optical: LANDSAT 5/TM	48.28, 50.39, 46.05	Clay (Sandy clay 29%, gyttja clay 36%)(2)
Porvoo 25° 50 ' E, 60° 50' N	1990: 13.5., 21.6 29.7.	Optical: SPOT/HRV2/XS	50.01, 48.1	Clay (2) (gyttja clay 38%)
	1990: 21.6	Optical: LANDSAT 5/TM	-
	1990: 25.7, 24.8	Microwave SAR: HUTSCAT Scatterometer (f = 5.4 GHz, C-band, 9.8 GHz, X-band), VV, HH, VH, HV polarizations(3)
	1995: 25.6., 15.7., 1996:25.7.	Optical: SPOT XS 69 225	–1.6
	14.6.1994	Optical: SPOT XS 65 220	5.2

⁽¹⁾ Used in Equation 2 (Appendix C) ⁽²⁾ Gyttja clay contains peat and mud fractions) ⁽³⁾ Helicopter mounted [61,62,63] ⁽⁴⁾ Main soil type classification with soil sub fractions ⁽⁵⁾ No Envisat data obtained in 2005–2006 because of re-programming of ASAR sensors and technical problems in ESA. Only ground truth phenological and LAI measurements measured in experimental areas.

(Appendix B). The Lapua experimental site, consisting of Lapua, Ilmajoki and Seinäjoki experimental sites (Figure 2), was located near the Gulf of Bothnia on sandy clay type soils. Respectively Porvoo and Kirkkonummi experimental sites were located close to the Baltic Sea with humid marine climatic conditions. Jokioinen and Mellilä sites were located on inland plateaus with mainly clay type soils.

2.3. Calibration

2.3.1. Phenological classification algorithm (SatPhenClass) for satellite data

The calibrated optical reflectance and microwave SAR backscattering data (

Table 11. The optical calibration parameters of Landsat and SPOT—satellites [57].

Satellite/ sensor (1),(2)	Sensor channel/λ (wave length)/type	Date	Sun elevation Angle (deg.)	S(6)	α(7)	β(8)
TM1(1)	450–520 nm/PAR(3)	24.5.1989	48.28	620	0.602	–1.5
TM2(1)	520–600 nm/PAR(3)	25.6.1989	50.39	577	1.17	–2.8
TM3(1)	630-690nm/infrared(3)	27.7.1989	46.05	493	0.806	–1.2
TM4(1)	760–900 nm/near infra	21.6.1990	50.01	332	0.815	–1.5
TM5(1)	1.55–1.75 μm /mid infra			67.1	0.108	–0.37
TM6(1)	10.4–12.5 μm/thermal IR(4)			n.a.(3)	-	-
TM7(1)	2.08–2.35 μm /mid infra			24.5	0.057	–0.15
HRV2/S1(2)	500–590 nm / PAR(3)	13.5.1990	n.a.(5)	587.0	1.22181	0
HRV2/S2(2)	610–680 nm / infrared(3)	29.7.1990	48.1	502.0	1.22545	0
HRV2/S3(2)	790–890 nm/near infrared			331.0	1.29753	0

⁽¹⁾ Landsat5 / TM: Landsat (1987). Sensor: Thematic Mapper (λ = 450 nm–2.35 μm, resolution 30 × 30 m^{2 (2)} SPOT / HRV2/XS: SPOT (1986). High Resolution Visible 2 (λ=500 nm–890 nm, resolution 20 × 20 m²). ⁽³⁾ PAR – Photosynthetically Active Radiation (λ = 400–700 nm) ⁽⁴⁾ Thermal infrared channel with ground resolution of 120 m. ⁽⁵⁾ n.a.—parameter not available ⁽⁶⁾ S—Equivalent solar radiance ⁽⁷⁾ α—Calibration gain coefficient ⁽⁸⁾ β—Calibration offset coefficient.

, Table 12 and Table 13, Appendix B.) measured from different satellite systems was classified with a phenological classification algorithm (SatPhenClass) for spring cereals. The detailed SatPhenClass classification model and VGI models applied in Parts I and II (detailed optical) are depicted in Figure 3a,b. The SatPhenClass classified satellite data were used in calibration of the VGI yield models (Figure 1). Figure 3a depicts the simplified ’Black Box’ diagram and Figure 3b depicts corresponding components on system analysis level. The pseudocode for the SatPhenClass-algorithm is presented in Appendix D.

Figure 3. (a) Simplified Black Box diagram for SatPhenClass satellite data classification and VGI models. (b) Satellite data classification diagram for spring cereals, ETS with modifications after Kontturi and Mukula & Rantanen [45,46].

Figure 3. (a) Simplified Black Box diagram for SatPhenClass satellite data classification and VGI models. (b) Satellite data classification diagram for spring cereals, ETS with modifications after Kontturi and Mukula & Rantanen [45,46].

The radiometric and geometric corrected surface reflectance (rf) and backscattering (σ ⁰) data were used as input data for the dynamic SatPhenClass classification model using Julian date (DOY) as a driving variable. Figure 3a,b depict SatPhenClass-model diagram for spring cereal (wheat, barley and oats) classification using optical and microwave satellite data.

2.3.1.1. SatPhenClass Classification Categories

SatPhenClass algorithm yields four classification classes (a_p, b_p c_p d_p, Table 2) based on spring cereal DVS (Development stage) and corresponding satellite measurement Julian date [39,45]. The SatPhenClass classification model developed in this project utilizes BBCH numerical growth scale [41,42]. In addition, Zadoks, Feekes and Haun growth scales were used as scaling references [43,47]. The SatPhenClass model contains four BBCH threshold values for phenological classification: 0, 12, 50, and 90 (Figure 3b, phases VII–X, [43]).

Table 2. SatPhenClass Effective classifiers: Effective temperature sum (ETS) and LAI_max values for high latitude cereals on growing zones I–IV [23,45].

Development period [35,50,51]	BBCH DVS Class (Range) (1)	Estimated mean Julian DOY (JDay) (2)	ETS(Tb) minimum requirements (dd) in cultivation zones I–IV (3)	Observed LAImax (Li-Cor 2000) X ± sd
Sowing-two-leaf & double ridge stages	ap (0–12)	165, 175	Swh Tb 5°: 130 (e)–140 (l) Brl Tb 5°: 130 (e)–140 (l) Oats Tb 5°: 130 (e)–140(l)	Swh 2.12 ± 0.46 Brl 1.85 ± 0.89 Oats 1.54 ±0.65
Two leaf-ear emergence, LAImax exposure with fully closed canopy structures	bp (12–50)	225	480 (Tb 4°) Swh Tb 5°: 450 (e)–460 (l) Brl Tb 5°: 800 (e)–950 (l) Oats Tb 5°: 370 (e)–400(l)	Swh 4.27 ± 0.84 Brl 4.05 ± 0.58 Oats 3.44 ± 0.78
Ear emergence, anthesis-maturity, grain filling	cp (50–90)	255	399 (Tb 8°)	Swh 3.87 ± 0.24 Brl 3.24 ± 0.87 Oats 2.14 ± 0.25
Senescence, post-harvest phase	dp (>90)	> 255
Sowing-maturity	ap, bp, cp	165–255	1050 ± 30, Tb 5° Swh Tb 5°: 970 (e)–1040 (l) Brl Tb 5°: 800 (e)–950 (l) Oats Tb 5°: 900 (e)–990(l)

⁽¹⁾ BBCH phenological development scale (Figure 3a, b) [41,42,43] ⁽²⁾ Mean values for JDay_ETS and JDay_LAI (Figure 3b, sections IV–V) ⁽³⁾ T_b [°C] threshold baseline temperature ^(e)—Finnish early cultivars, ^(l)—Finnish late cultivars.

Table 2. SatPhenClass Effective classifiers: Effective temperature sum (ETS) and LAI_max values for high latitude cereals on growing zones I–IV [23,45].

Development period [35,50,51]	BBCH DVS Class (Range) (1)	Estimated mean Julian DOY (JDay) (2)	ETS(Tb) minimum requirements (dd) in cultivation zones I–IV (3)	Observed LAImax (Li-Cor 2000) X ± sd
Sowing-two-leaf & double ridge stages	ap (0–12)	165, 175	Swh Tb 5°: 130 (e)–140 (l) Brl Tb 5°: 130 (e)–140 (l) Oats Tb 5°: 130 (e)–140(l)	Swh 2.12 ± 0.46 Brl 1.85 ± 0.89 Oats 1.54 ±0.65
Two leaf-ear emergence, LAImax exposure with fully closed canopy structures	bp (12–50)	225	480 (Tb 4°) Swh Tb 5°: 450 (e)–460 (l) Brl Tb 5°: 800 (e)–950 (l) Oats Tb 5°: 370 (e)–400(l)	Swh 4.27 ± 0.84 Brl 4.05 ± 0.58 Oats 3.44 ± 0.78
Ear emergence, anthesis-maturity, grain filling	cp (50–90)	255	399 (Tb 8°)	Swh 3.87 ± 0.24 Brl 3.24 ± 0.87 Oats 2.14 ± 0.25
Senescence, post-harvest phase	dp (>90)	> 255
Sowing-maturity	ap, bp, cp	165–255	1050 ± 30, Tb 5° Swh Tb 5°: 970 (e)–1040 (l) Brl Tb 5°: 800 (e)–950 (l) Oats Tb 5°: 900 (e)–990(l)

⁽¹⁾ BBCH phenological development scale (Figure 3a, b) [41,42,43] ⁽²⁾ Mean values for JDay_ETS and JDay_LAI (Figure 3b, sections IV–V) ⁽³⁾ T_b [°C] threshold baseline temperature ^(e)—Finnish early cultivars, ^(l)—Finnish late cultivars.

In class a_p (Figure 3b, phase VII), when the BBCH value is below 12, the cereal plant is in vegetative phase between germination, two-leaf and double ridge stages.

In class b_p (Figure 3b, phase VIII), the BBCH value is between 12 and 50. When BBCH value is between 12 and 40, the cereal leaf development and stem elongation occur. The booting of the uppermost flag leaf and the initialization of heading occur, when BBCH value is between 40 and 50. In addition, the Leaf Area Index maximum is usually observed in class b_p.

In class c_p (Figure 3b, phase IX), the BBCH value is between 50 and 90. The cereal inflorescence emergence occurs, when BBCH values are between 50 and 60. When the BBCH value is above 60, the cereal transition to generative grain filling phase begins with the initialization of anthesis and flowering. When the BBCH value is between 70 and 80, the development of grains in the head occurs with early milk, hard dough and fully ripened phases. Finally in class d_p (Figure 3b, phase X), when the BBCH value is above 90, the cereal is in senescence phase. On the average, in Finnish long day growing conditions, the phenological class a_p corresponds for May, class b_p for June, class c_p for July and d_p for August and September.

The SatPhenClass phenological classification model enabled the temporal synchronization of different optical and microwave satellite data along with spring cereal phenological development. The classified optical and microwave SAR data can in turn be used as input for cereal VGI (Vegetation Indices) or dynamic crops models, like CERES-Wheat model [33,34,35,36,48]. In CERES-Wheat phenological submodel, the genotype x cultivar variation is taken into account by using genetic coefficients controlling cultivar photoperiodism, vernalization and phyllochron development. With SatPhenClass submodel, the cultivar and genotype level effects for phenological development were assumed to be constant. The genetic coefficients for high latitude Nordic spring cereals have been previously calibrated and assessed by Laurila [35] and later reviewed by Saarikko [49]. In Sweden Åfors et al. [50] have reviewed the use phenological growth scales for cereals grown in high latitude Nordic growing conditions. Recently Peltonen-Sainio et al. [51] revised the Feekes and Zadoks cereal growth scaling [43,47] for Finnish long day growing conditions.

2.3.1.2 LAI and ETS Ground Truth Sampling Methodology for the SatPhenClass Algorithm

During growing season crop phenological observations (Julian Day of Year, DOY), LAI measurements (Leaf Area Index), biomass, growing density and other ground truth parameters were recorded simultaneously with satellite measurement dates (JDay_SAT, Table 9, Appendix B) in the experimental sites (Table 1). Kuittinen et al. [23] presented a detailed description for the ground truth sampling methodology. Portable Li-Cor 2000 Plant Canopy Analyzer [52] was used to measure spring cereal LAI values with corresponding Julian dates (JDay_LAI) during the growing period.

The sampling methodology on experimental sites was applied after Kuittinen et al. [23] and is displayed in Figure 6 (Appendix A). Sampling plots were randomly selected from the field parcels in experimental areas (Table 1). Sampling plots were adjusted to Landsat pixel size (30 × 30 m). In each plot inside selected parcels, the LAI and biomass were measured using 10 points. The distance between the points varied between 20 and 30 m depending on the size of the parcel. The LAI value of each point was the mean of 9 LAI measurements around the original point integrated automatically by the Li-Cor 2000 Plant Canopy Analyzer. The average plant biomass on the field parcel was estimated by selecting randomly 20 plant samples. Late on, the dry-weights of the dried plant samples were measured. The growing density was determined by calculating the number of plants from the randomly selected 1 meter long sowing row with two replicates. After the harvest in the autumn, the final cereal grain yield levels (kg/ha corrected to 15% moisture content) were estimated by the farmers in the experimental areas. Yield samples were measured from the granary silos.

The Finnish Meteorological Institute provided Cumulative Effective Temperature Sum (ETS, T_b = 5 °C) measurements with corresponding Julian dates (JDay_ETS, Figure 4b). Kriging interpolation method was applied to estimate ETS values in experimental sites from the nearest weather stations (Figure 4c).

Figure 4. (a) Growing zones (I–V) and MTT Experimental Station locations in Finland, (b) ETS (T_b 5 °C) cumulative isolines and (c) Meteorological Weather Stations in Finland (Original Data (C) MTT Agrifood Research Finland, Finnish Meteorological Institute and Finfood Finland).

Figure 4. (a) Growing zones (I–V) and MTT Experimental Station locations in Finland, (b) ETS (T_b 5 °C) cumulative isolines and (c) Meteorological Weather Stations in Finland (Original Data (C) MTT Agrifood Research Finland, Finnish Meteorological Institute and Finfood Finland).

The cumulative ETS(T_b =5 °C) and measured LAI values were used as categorical classifiers in the SatPhenClass algorithm [23,24]. The Julian dates of the measured satellite data (JDay_Sat) were compared with the estimated Julian dates for spring cereals using the SatPhenClass algorithm (Figure 3b, sections II–VI: LAI development and ETS accumulation). The LAI_max-values for high latitude spring cereals (Table 2) with corresponding Julian dates (JDay_LAI, Figure 3b, section V) were used to estimate the transition date from generative (class b_p) to grain filling phase (class c_p). The LAI_max for spring cereals varied between 3.44 and 4.27 in b_p phase and between 2.14 and 3.87 in c_p phase. Table 2 depicts the cereal ETS temperature sum requirements in Finnish growing conditions for different phenological classes (a_p–d_p) with corresponding threshold temperatures (T_b) and typical Julian dates (JDay_ETS, Figure 3b, section IV). In growing zones I–IV, the ETS requirements for a spring cereals [53,54] vary over years and locations (Figure 4a). In Nordic high latitude long day growing conditions the growing season is generally defined as the period when mean air temperature exceeds +5 °C (T_b 5 °C, Table 2 [45]). The average thermal time requirement of 1050 ± 30° degree-days (dd, T_b 5 °C) from sowing to yellow ripening stage is considered adequate for spring cereal ideotypes grown in growing zones I–IV [23,55,56]. The corresponding averaged Finnish ETS isolines (T_b 5 °C) are presented in Figure 4b. According to Hakala [53], Hakala et al. [54] typical Julian dates for spring cereals in southern Finland are: Sowing on average on 14–16th of May (DOY 164–166), emergence 25–30th of May (DOY 175–180), anthesis on 5–15th of July (DOY 220–225), full maturity on 5–15th of August (250–255).

In the final stage the SatPhenClass algorithm (Figure 3b, sections IV–V: ETS and LAI Classifiers) compares the ETS and LAI classification results. The classification result is accepted only, if both submodels reach concord (Figure 3b, section VI: ETS and LAI Evaluation). The measured optical or microwave satellite data with corresponding Julian date (JDay_Sat) is finally set to corresponding cereal phenological class a_p-d_p shown in Figure 3b (graphical plot in left lower corner). The phenologically classified optical and SAR satellite data can be used as input for optical VGI infrared and NDVI models (Table 3, Figure 3b, XI: Output VGI_I–II) and for Composite NDVI and SAR model (Figure 3b, XII: Output VGI_III).

Table 3. Composite Vegetation Indices Models (I–III) with dependent baseline yield (y_b) and independent variables [23,57,58,59]^(1),(2).

Model (Type)	Model equation, Table 12, Appendix C (1), (2)	Independent variables (2)	Model name, description of derived satellite parameters used in regression equations
I Optical	yb= rf3(ap,bp,cp) + rf4(ap,bp,cp)	rf3(ap,bp,cp) rf4(ap,bp,cp)	Polynomial infrared model Calibrated reflectance (rfchannel) values for infrared (channel = 3) and near infrared (channel=4) during growing season. - ap,bp,cp classes correspond on average to Zadoks crop phenological growth scale with cereals: ap: 0–12, bp: 12–50, cp: 50–90, dp: >90 [43]. (2)
II Optical	yb = NDVI (ap,bp,cp) NDVI = (NIR(ap,bp,cp) RED(ap,bp,cp))/ (NIR(ap,bp,cp) + RED(ap,bp,cp))	NDVI (ap,bp,cp)	Normalized Difference Vegetation Index (NDVI) model [57,58,59] (2). - Ratio between near infrared and infrared channel surface reflectance (rf) values - rf3 infrared (λ = 0.63–0.69 µm) - rf4 near infrared (λ = 0.76–0.90 µm) - Landsat NDVI = (rf4(ap,bp,cp) – rf3(ap,bp,cp))/(rf4(ap,bp,cp) + rf3(ap,bp,cp)) , - SPOT NDVI = (rf3(ap,bp,cp) – rf2(ap,bp,cp))/(rf3(ap,bp,cp) + Rf2(ap,bp,cp))
III SAR + optical	yb = NDVI (ap,bp,cp ) + σ0 HH5GHz,(ap,bp,cp,dp) + σ0 VV5GHz,(ap,bp,cp,dp) + σ0 HV5GHz,(ap,bp,cp,dp) + σ0 VH5GHz,(ap,bp,cp,dp)	NDVI (ap,bp,cp), σ0 5GHz,(ap,bp,cp, dp)	Composite multispectral SAR and NDVI model for spring cereals (swh, oats, barley) using NDVI reflectance and microwave backscattering (σ0, f = 5.3, 5.4 GHz) data with horizontal (HH), vertical (VV), and cross-polarization (HV, VH) levels. Instruments: HUTSCAT Scatterometer, ERS/SAR, Radarsat/SAR, Envisat/ASAR ( Table 10. SAR and optical satellite systems used in remote sensing campaigns. Satellite type Name Sensor Experimental locations & years Reference Microwave ASAR ENVISAT (3) ASAR(6), f = 5.3 GHz, C-band), VV, HH, VV/HH, HV/HH, or VH/VV polarizations Seinäjoki, Lapua (2002–2004, Table 11). http://earth.esa.int, http://envisat.esa.int/object/index.cfm?fobjectid=3772, [27,28] Microwave Scatterometer HUTSCAT(1),(6) Scatterometer (f = 5.4 GHz, C-band, 9.8 GHz, X-band), VV, HH, VH, HV polarizations Porvoo, calibration data (1990, Table 11). www.space.hut.fi/research/equipment/ hutscat.html. [61,62,63] Microwave SAR ERS12) SAR (5), f = 5.3 GHz, C band, (λ = 5.7 cm), VV polarization Seinäjoki, Lapua (1995–1996, Table 11). http://earth.esa.int/ers, earth.esa.int/ers/sar, [26,62] Microwave SAR Radarsat1 SAR(5), f = 5.3 GHz, C band, HH polarization Seinäjoki, Lapua (2001, Table 11). Reference [27], Canadian Space Agency (CSA). ccrs.nrcan.gc.ca/, radar/spaceborne/radarsat1/index_e.php Optical Landsat 5 Thematic Mapper (TM) (λ = 450 nm–2.35 μm) Porvoo, Mellilä, Kirkkonummi, Jokioinen, Lapua (1989–1997, Table 11). www.landsat.org, [76,22,23] Optical SPOT 2 HRV2/XS (λ = 450 nm–890 nm) Porvoo, Mellilä, Kirkkonummi, Jokioinen, Lapua (1989–1997, Table 11). www.spot.com, www.spotimage.fr, [77,22,23] Multi- sensor ADEOS1 Advanced Earth Observing Satellite4) AVNIR, ILAS, RIS, IMG,TOMS: atmospheric greenhouse gas (CO2, O3, CH4) columns (4) (1996-1997, non-operational) Collaboration with Prof. Hiroshi Koizumi, NIAES/ Tsukuba, Japan NASDA/ JAXA, http//home.gna.org/adeos/ http://kuroshio.eorc.jaxa.jp/ADEOS,http://msl.jpl.nasa.gov/QuickLooks/ adeosQL.html (1) HUTSCAT is helicopter mounted, see space.hut.fi/research/equipment/hutscat.html [63](2) European Remote Sensing satellite (ESA, European Space Agency)(3) ENVISAT (Environmental satellite, ESA) AOS (Announce of Opportunity) contract: AOE-488 for ENVISAT(4) ADEOS1 AOS contract: NASDA Contract 1062/Vegetation and Biology, MAFF, 5118/416/94(5) SAR, Synthetic Aperture Radar (6) Advanced Synthetic Aperture Radar(6) HUTSCATused for calibration verification purposes only (helicopter mounted) , [60,61,62,63]).

⁽¹⁾ Model equations presented in Table 12 (Appendix C), for abbreviations refer to Table 7. ⁽²⁾ Independent variables classified with SatPhenClass-algorithm (Figure 3a,b).

Table 3. Composite Vegetation Indices Models (I–III) with dependent baseline yield (y_b) and independent variables [23,57,58,59]^(1),(2).

Model (Type)	Model equation, Table 12, Appendix C (1), (2)	Independent variables (2)	Model name, description of derived satellite parameters used in regression equations
I Optical	yb= rf3(ap,bp,cp) + rf4(ap,bp,cp)	rf3(ap,bp,cp) rf4(ap,bp,cp)	Polynomial infrared model Calibrated reflectance (rfchannel) values for infrared (channel = 3) and near infrared (channel=4) during growing season. - ap,bp,cp classes correspond on average to Zadoks crop phenological growth scale with cereals: ap: 0–12, bp: 12–50, cp: 50–90, dp: >90 [43]. (2)
II Optical	yb = NDVI (ap,bp,cp) NDVI = (NIR(ap,bp,cp) RED(ap,bp,cp))/ (NIR(ap,bp,cp) + RED(ap,bp,cp))	NDVI (ap,bp,cp)	Normalized Difference Vegetation Index (NDVI) model [57,58,59] (2). - Ratio between near infrared and infrared channel surface reflectance (rf) values - rf3 infrared (λ = 0.63–0.69 µm) - rf4 near infrared (λ = 0.76–0.90 µm) - Landsat NDVI = (rf4(ap,bp,cp) – rf3(ap,bp,cp))/(rf4(ap,bp,cp) + rf3(ap,bp,cp)) , - SPOT NDVI = (rf3(ap,bp,cp) – rf2(ap,bp,cp))/(rf3(ap,bp,cp) + Rf2(ap,bp,cp))
III SAR + optical	yb = NDVI (ap,bp,cp ) + σ0 HH5GHz,(ap,bp,cp,dp) + σ0 VV5GHz,(ap,bp,cp,dp) + σ0 HV5GHz,(ap,bp,cp,dp) + σ0 VH5GHz,(ap,bp,cp,dp)	NDVI (ap,bp,cp), σ0 5GHz,(ap,bp,cp, dp)	Composite multispectral SAR and NDVI model for spring cereals (swh, oats, barley) using NDVI reflectance and microwave backscattering (σ0, f = 5.3, 5.4 GHz) data with horizontal (HH), vertical (VV), and cross-polarization (HV, VH) levels. Instruments: HUTSCAT Scatterometer, ERS/SAR, Radarsat/SAR, Envisat/ASAR (Table 10, [60,61,62,63]).

⁽¹⁾ Model equations presented in Table 12 (Appendix C), for abbreviations refer to Table 7. ⁽²⁾ Independent variables classified with SatPhenClass-algorithm (Figure 3a,b).

2.3.2. Composite Vegetation Indices (VGI) models I–III

Three different Composite Vegetation Indices (VGI) models (I–III, Table 3) were calibrated with linear and non-linear polynomial regression models (Equations 1–2, Appendix C). Calibrated optical reflectance data were used for Models I–II (Models I-Infrared, II-NDVI) and composite SAR and optical data for the Model III (Composite NDVI and SAR). The NDVI optical and SAR satellite data were pre-classified with SatPhenClass model into the vegetative pre-anthesis phenological classes (a_p, b_p) and into the generative post-anthesis and senescence classes (c_p, d_p, Figure 3b, phases IV–VI).

The VGI models (I–III) were calibrated with SAS^© linear stepwise REG/Stepwise (STW) regression algorithm (Equation 1, Appendix C, [57,58]) and with non-linear polynomial RSREG regression algorithm (Equation 2, Appendix C) by using SAR and optical satellite data (Table 12, Appendix C). The non-linear RSREG algorithm takes into account (i) the general linear component, (ii) the polynomial quadratic response of non-linear variation (iii) the cross product effect of dependent variables and (iv) the total model component (Table 5). The statistical significance levels are given in

Table 8. Statistical significance levels.

Significance levels (α = reference probability) (1)
N.S.	Statistically non-significant
0	Moderately significant on 10% error level, α(0.10)
*	Significant on 5% error level, α(0.05)
**	Highly significant on 1% error level, α(0.01)
***	Highly significant on 0.1% error level, α (0.001)

⁽¹⁾ Used in testing the p-value (Pr) in F-test

(Appendix B).

The polynomial infrared model (I) is a polynomial linear regression model, which estimates baseline yield production (y_b, kg/ha) with infra red (rf₃) and near infrared (rf₄) channels. Respectively in VGI model II the NDVI Index [57,58,59] and in Model III the composite SAR and optical NDVI components were applied.

2.4. Calibration of Data

2.4.1. Microwave SAR backscattering (σ⁰ ) and optical reflectance data

The primary combined microwave SAR and optical dataset pool (Table 1) consisted of (i) microwave SAR/ASAR (ERS, Radarsat, Envisat) and (ii) optical (Landsat and SPOT) datasets from the Etelä-Pohjanmaa (Lapua, Seinäjoki and Ilmajoki sites), Nylands Svenska (Kirkkonummi) and Häme (Jokioinen) Agricultural Advisory and Rural Development Centres 1996–2006 (Table 1, Figure 2, [27,28]). The satellite measurement dates (JDay_SAT) in ground truth measurement sites (Lapua, Kirkkonummi, Jokioinen, Mellilä and Porvoo) are depicted in Table 9 (Appendix B). The Technical Research Centre of Finland (VTT, Dept. of Space Research) and the Finnish Geodetic Institute digitized and performed the classification of field parcels from the satellite images [23,27,28]. The total cultivation area (Seinäjoki, Ilmajoki and Lapua, Mellilä, Porvoo, Kirkkonummi, Jokioinen) estimated from satellite images was 1,253 ha for all spring cereal crops.

The composite SAR and optical reflectance datasets were used for spring cereal yield modeling under actual non-potential field conditions. Water stress limited significantly cereal yield accumulation on clay type soils in Kuuma and Porvoo experimental areas during growing season (Table 1, Figure 2).

The SAR calibration details for ERS, Radarsat and Envisat data have been previously published by Karjalainen and et al. [27,28] and by Matikainen et al. [26]. The set of Radarsat-1 images consisted of scenes acquired with various beam modes of Radarsat-1 satellite. The Radarsat incidence look angles varied between 39–44° in the experimental sites (Table 1). Envisat ASAR images were acquired in dual-polarization mode, the VV and VH polarization channels were used in calibration of VGI models (Table 11 and Table 12, Appendix B).

In addition, two older auxiliary optical [23] and microwave (HUTSCAT scatterometer, Table 10, Appendix B) [61,63] datasets consisting from Mellilä and Porvoo sites (Nylands Svenska Agric. Advisory Centre, Figure 2) were used in calibration of SatPhenClass phenology and VGI models (I–III, Table 3). The Porvoo experimental site (Kullogård experimental farm, Figure 5) was used for both optical Landsat and SPOT and microwave HUTSCAT measurements for spring cereals.

Figure 5. Landsat/TM satellite image from the Porvoo field experimental area. The image location was used for optical Landsat, SPOT and microwave HUTSCAT measurements. Porvoo River on the right corner of the image.

Figure 5. Landsat/TM satellite image from the Porvoo field experimental area. The image location was used for optical Landsat, SPOT and microwave HUTSCAT measurements. Porvoo River on the right corner of the image.

The helicopter mounted HUTSCAT (Helsinki Univ. of Technology, Lab. Of Space Technology [63]) is a monostatic scatterometer system measuring with C- (f = 5.4 GHz) and X-bands (f = 9.8 GHz) with horizontal (HH), vertical (VV) and cross-polarizations levels (VH, HV). In this study only C-band data from HUTSCAT was used for the calibration of VGI models. Both the transmitter and receiver are mounted in the same flight rack under the helicopter [63]. The technical configuration of HUTSCAT resembles SAR/ASAR spectrometers currently onboard ERS, Radarsat and ENVISAT satellites (Table 10, Appendix B). The technical configuration of the monostatic HUTSCAT scatterometer has been previously published by Hallikainen et al. [61] and Koskinen et al. [62]. The HUTSCAT calibration details have been published Hyyppä et al. [63].

The Finnish Geodetic Institute applied ortho-rectification and radiometric correction for the optical Landsat and SPOT satellite data (Table 1, [23]). The optical calibration procedure for different Vegetation Indices models was originally presented by Price [57] and later revised and utilized for cereals by Flowers et al. [64], Prasad et al. [65], Maas [66,67,68] and Maas and Dunlap [69]. Table 11 (Appendix B) presents calibration parameter values applied in optical calibration equations 3–5 (Appendix C), which convert binary digital count (dc) values from satellite optical sensors into surface reflectance values (rf). The Landsat and SPOT raw data coded as 8 bit binary values (digital counts, dc) were converted to integer values between 0–255 (surface reflectance, r_f ) with higher values indicating higher reflectance radiation to satellite sensors from ground cover and atmosphere. The satellite digital count-value (dc) to spectral radiance (R) transformation was performed after Equation 3. The spectral radiance —values were transformed into surface reflectance (r_f) by means of solar equivalent radiance (S) (Equation 4). In addition, the sun zenith-angle during the observation time was taken into account (Equation 5).

The calibrated SAR backscattering and ortho-rectified and radiometric corrected optical reflectance data were applied in the calibration of VGI models. Total of 18 SAR and cloud free optical satellite images (Landsat5 /TM and SPOT /HRV2) were obtained during the 1996–2006 campaign. Images were temporarily well dispersed along the growing season between May, June, July and August (Table 11 and Table 12, Appendix C). Additional details of the optical calibration procedure and corresponding calibration parameters are presented in Electronic Supplementary Information Appendix (ESI).

2.5. Validation

2.5.1. Validation of Composite Vegetation Indices (VGI) Models

After the calibration procedure for VGI Composite models (I–III)), the validation procedure was applied to test the VGI Composite models for cereal non-potential baseline grain yield (y_b, kg/ha) predictions. Yield differences were calculated as absolute values (kg/ha) and per cent (%) differences vs. observed MAFF annual yield inventory [29,70] and MTT Official Variety Trial grain yield estimates in corresponding sites and years [55,56,71]. MAFF and MTT validation datasets with over 3º latitude variation were used to assess cereal baseline yield variation (y_b) in growing zones I–IV (Figure 4). In addition, the Polynomial Infrared and NDVI models (Models I, II) were analyzed with soil*species and soil*cultivar covariance categories. The MTT Official Variety Trial yield data for one Finnish malting barley cultivar (cv. Inari) and for two Finnish spring wheat cultivars (cv. Manu, Satu) in cultivation zones I–IV were used in analyzing soil*cultivar covariances (Table 1, Figure 4a).

2.5.2. Validation Data

The calibrated VGI yield models (I–III) were validated with two independent cereal yield datasets consisting of (i) the Ministry of Agriculture and Forestry Official Inventory Statistics 1994–2006 data for cereal yields [29,30] and (ii) the MTT Agrifood Research Finland Official Variety trial results 1994–2006 (Figure 4, [56,71]). The MAFF validation dataset consisted yield inventory statistics from Etelä Pohjanmaa (Averaged III–IV Zone) Agricultural Advisory Centre. The MTT Official Variety Trial dataset consisted of yield trials from the Ylistaro (Zone III) and Ruukki (Zone IV) MTT Experimental Stations Averaged to Zone III–IV. Zones III–IV were pooled together because the satellite and ground truth measurement sites in Lapua, Seinäjoki and Ilmajoki were located in the middle of 1,100 ETS_{(Tb = 5 °C)} isoline (Figure 4b).

3. Results

3.1. Cereal Canopy Soil Backscattering Covariance Variation

The SAR (Radarsat, ERS), ASAR (Envisat) and HUTSCAT backscattering variance (σ⁰) between different spring cereals in Lapua, Seinäjoki and Porvoo experimental sites is presented in Table 4. In addition, SAR backscattering horizontal and vertical polarization levels (HH, VV) and the cross-polarization (VH/HV) of the monostatic HUTSCAT scatterometer are presented. The Radarsat horizontal (HH) backscattering (σ⁰) variation between malt (MBrl) and fodder barley (FBrl) is presented in Radarsat sensor column. Table 4 depicts cereal canopy*soil backscattering covariance during generative phenological phases c_p and d_p (anthesis, grain full maturity and canopy senescence, Table 2) grown on clay, coarse and fine sand type soils.

Both the soil type and soil*cereal canopy covariances expressed extensive variation with SAR backscattering signal (σ ⁰) on clay type soils, which were dominant in experimental areas in zones I–IV (Table 1). According to Henderson & Lewis [60] the vertically oriented components (VV) with cereals, especially the stems, interact effectively with vertically polarized signals resulting increased backscattering attenuation. Oat species with more planophile canopy and panicle inflorescence structures have different backscattering polarization properties from corresponding properties of wheat and barley, which have more erectophile head and canopy structures [23,28,72].

On sandy clay soils the HUTSCAT backscattering signal (σ⁰) varied between –6 dB and –29 dB. On fine and coarse sands the Envisat and ERS SAR signal amplitude varied between –17 dB (Envisat in VV, VH dual polarization mode) and –8 dB (ERS) respectively. In the Porvoo experimental area the topsoil (5–10 cm) and the subsoil (<10 cm) consisted of fine and coarse sandy alluvium deposits from the Porvoo river (Figure 5) and gyttja clays containing minor fractions of organic compounds (mould, peat and mud) and silt. In Lapua and Seinäjoki areas the major soil type was sandy clay with minor fractions of organic mould and humus fractions (Table 1).

The SAR backscattering signals indicated extensive topsoil*cereal canopy layer covariance interaction especially after the anthesis (phenological stage c_p) and during grain filling, yellow ripening and canopy senescence (d_p). When analyzing the variation between different SAR sensors the HUTSCAT backscattering signal (f = 5.4 GHz, Table 10, Appendix B) varied with wheat, barley and oat cereals between –29 dB (oats* c_p*VH*coarse sandy clay) and –6 dB (wheat*c_p*VV*coarse and fine sand clay, barley*c_p*VV*fine sand and gyttja clay). The ERS signal (f = 5.3 GHz) varied between –8 dB (barley*d_p*VH*fine and coarse sandy clay) and –10 dB (barley*c_p*VV* sandy clay, oats*c_p*VV*sandy clay). The Radarsat signal (f = 5.3 GHz) varied between –13 dB (malt barley*c_p*HH* fine and coarse sandy clay) and –10 dB (wheat*d_p*HH*sandy clay). The Envisat signal (f = 5.3 GHz) varied between—16 dB (oats*d_p*VH* fine and coarse sandy clay) and –9 dB (barley*d_p*VV* sandy clay).

The helicopter mounted HUTSCAT signal amplitude varied more significantly when compared to other backscattering signals (Envisat, ERS, Radarsat) measured on clay and sandy soils. In addition, the HUTSCAT C-band measurement frequency was slightly higher (f=5.4 GHz) than with Envisat, ERS and Radarsat (f = 5.3 GHz).

Table 4. SAR backscattering variation (σ ⁰, f=5.3–5.4 GHz) on sand and clay type soils with different polarization levels during anthesis (c_p), full grain maturity and canopy senescence (d_p). ^(1),(2)

SARS ensor	Main soil Type	S.wheat [dB] (X±Sd)				Barley 1) [dB] (X±Sd)				Oats [dB] (X ± Sd)
DVS3)		cp (July)		dp (Aug.)		cp (July)		dp (Aug.)		cp (July)		dp (Aug.)
HUT SCAT Scattero Meter	Sandy clay	VV/ HH	VH/HV 3)	VV/ HH	VH/ HV 3)	VV/ HH	VH/ HV 3)	VV/ HH	VH/ HV 3)	VV/ HH	VH/ HV 3)	VV/ HH	VH/ HV 3)
		–6.40 (±0.84)/ –16.21 (±0.32)	–20.26 (±0.35)/ –20.16 (±0.14)	–10.53 (±0.24)/ –10.09 (±0.23)	–20.71 (±0.15)/ –20.05 (±0.13)	–5.78 (±1.79)/ –12.27 (±1.72)	–22.85 (±0.90)/ –22.97 (±0.87)	–11.26 (±1.79)/ –10.21 (±1.72)	–22.54 (±0.67)/ –22.48 (±0.49)	–9.75 (±1.58)/ –18.52 (±1.87)	–21.09 (±0.90)/ –21.31 (±0.73)	–13.39 (±1.47)/ –11.62 (±1.97)	–21.54 (±0.76)/ –21.17 (±0.58)
Envisat ASAR	Fine, coarse sand	VV	VH	VV	VH	VV	VH	VV	VH	VV	VH	VV	VH
		–12.47 (±1.17)	–17.52 (±0.95)	–10.09 (±1.41)	–14.81 (±1.18)	–11.45 (±2.17)	–17.16 (±1.59)	–9.45 (±1.38)	–14.59 (±1.51)	–11.08 (±2.51)	–17.22 (±1.47)	–11.37 (±1.37)	–16.25 (±1.87)
Radarsat SAR (1),(2)	Fine, coarse sand	HH		HH		MBrl HH 2)	FBrl HH 2)	MBrl HH 2)	FBrl HH 2)	HH		HH
		–12.27 (±2.28)		–10.48 (±0.79)		–13.88 (±1.25)	–12.57 (±1.58)	–12.28 (±2.27)	–11.39 (±1.28)	–12.64 (±1.65)		–11.23 (±0.97)
ERS SAR	Fine, coarse sand					VV		VV		VV		VV
						–10.24 (±1.17)		–8.12 (±2.11)		–10.24 (±1.28		–9.28 (±1.02)

⁽¹⁾MBrl- Malt barley, FBrl—feed barley ⁽²⁾ Soil types in Table 1 ⁽³⁾ HUTSCAT σ⁰ VH/HV cross-polarization

Table 4. SAR backscattering variation (σ ⁰, f=5.3–5.4 GHz) on sand and clay type soils with different polarization levels during anthesis (c_p), full grain maturity and canopy senescence (d_p). ^(1),(2)

SARS ensor	Main soil Type	S.wheat [dB] (X±Sd)				Barley 1) [dB] (X±Sd)				Oats [dB] (X ± Sd)
DVS3)		cp (July)		dp (Aug.)		cp (July)		dp (Aug.)		cp (July)		dp (Aug.)
HUT SCAT Scattero Meter	Sandy clay	VV/ HH	VH/HV 3)	VV/ HH	VH/ HV 3)	VV/ HH	VH/ HV 3)	VV/ HH	VH/ HV 3)	VV/ HH	VH/ HV 3)	VV/ HH	VH/ HV 3)
		–6.40 (±0.84)/ –16.21 (±0.32)	–20.26 (±0.35)/ –20.16 (±0.14)	–10.53 (±0.24)/ –10.09 (±0.23)	–20.71 (±0.15)/ –20.05 (±0.13)	–5.78 (±1.79)/ –12.27 (±1.72)	–22.85 (±0.90)/ –22.97 (±0.87)	–11.26 (±1.79)/ –10.21 (±1.72)	–22.54 (±0.67)/ –22.48 (±0.49)	–9.75 (±1.58)/ –18.52 (±1.87)	–21.09 (±0.90)/ –21.31 (±0.73)	–13.39 (±1.47)/ –11.62 (±1.97)	–21.54 (±0.76)/ –21.17 (±0.58)
Envisat ASAR	Fine, coarse sand	VV	VH	VV	VH	VV	VH	VV	VH	VV	VH	VV	VH
		–12.47 (±1.17)	–17.52 (±0.95)	–10.09 (±1.41)	–14.81 (±1.18)	–11.45 (±2.17)	–17.16 (±1.59)	–9.45 (±1.38)	–14.59 (±1.51)	–11.08 (±2.51)	–17.22 (±1.47)	–11.37 (±1.37)	–16.25 (±1.87)
Radarsat SAR (1),(2)	Fine, coarse sand	HH		HH		MBrl HH 2)	FBrl HH 2)	MBrl HH 2)	FBrl HH 2)	HH		HH
		–12.27 (±2.28)		–10.48 (±0.79)		–13.88 (±1.25)	–12.57 (±1.58)	–12.28 (±2.27)	–11.39 (±1.28)	–12.64 (±1.65)		–11.23 (±0.97)
ERS SAR	Fine, coarse sand					VV		VV		VV		VV
						–10.24 (±1.17)		–8.12 (±2.11)		–10.24 (±1.28		–9.28 (±1.02)

⁽¹⁾MBrl- Malt barley, FBrl—feed barley ⁽²⁾ Soil types in Table 1 ⁽³⁾ HUTSCAT σ⁰ VH/HV cross-polarization

When analyzing the cereal variation for the wheat, the SAR signal varied between –25 dB (HUTSCAT *c_p*VH* fine and coarse sandy clay) and –6 dB (HUTSCAT *d_p*VH* sandy clay). The wheat VH cross-polarization amplitude was higher when compared with VV vertical levels both in anthesis (c_p) and full maturity (d_p) on clay and sandy soils. The wheat HH horizontal signal amplitude was higher in stage c_p compared to d_p.

The oat SAR amplitude varied between –29 dB (HUTSCAT *c_p*VH* sandy clay) and –6 dB (HUTSCAT *d_p*VH* sandy clay). The oats VH cross-polarization amplitude was higher when compared with VV vertical levels both in anthesis (c_p) and full maturity (d_p) on clay and sandy soils. Both the oats horizontal and vertical (HH, VV) signal amplitudes were higher in anthesis stage (c_p) compared to maturity stage (d_p).

With barley the SAR amplitude varied between –22 dB (HUTSCAT *c_p*VH* sandy clay) and –5 dB (HUTSCAT *c_p*VV* sandy clay). Especially the barley experimental plots with clay and peat topsoil layers in Porvoo (Table 1, Figure 2) were affected by excessive water from melting snow in the beginning of the growing season (a_p) causing the topsoil to reach the full soil water capacity. Later on during mid-summer, drought periods especially in generative phases (c_p, d_p), caused the soil moisture content to recede close to wilting point. These factors affected the translocation of assimilates to the head with a reduction effect on final grain filling. The barley VH cross-polarization amplitude was higher when compared with VV vertical levels both in anthesis (c_p) and full maturity (d_p) on clay and sandy soils. Both the barley horizontal (HH) and vertical (VV) σ ⁰ signal amplitudes were higher in anthesis stage (c_p) than in full maturity stage (d_p). The malt barley (MBrl) σ ⁰ amplitude with Radarsat varied between –13 dB (malt barley*c_p*HH* fine and coarse sandy clay) and –12 dB (malt barley*d_p*HH*sandy clay). Respectively with fodder barley (FBrl) the Radarsat signal varied between—12 dB (feed barley*c_p*HH* fine and coarse sandy clay) and –11 dB (feed barley*d_p*HH* sandy clay).

3.2. Calibration Results

3.2.1. Testing the SatPhenClass Phenological Model Accuracy

The SatPhenClass classification results for spring cereals are depicted in Table 5. The SatPhenClass classification accuracy for spring cereals varied between 61% and 89% with corresponding classification error ranging between phenological classes (a_p 21–28%, b_p 11–22%, c_p 14–28% and d_p 25–39%) as calculated from the observed mean phenological Julian dates (DOY). The observed phenological Julian dates for spring cereals were averaged from ground truth experimental sites in cultivation zones I–IV (Figure 2 and Figure 4). In addition, the mean observed Julian dates for spring cereals from MTT experimental stations in cultivation zones I–III were used for phenological validation (Zone I-Pernaja and Inkoo, Zone II-Jokioinen, Zone III and IV-Ylistaro and Ruukki).

The early emergence in vegetative phase (a_p, BBCH 0–12) in two leaf stage before double ridge induction and the senescence phase after full maturity and harvest (d_p), BBCH > 90) were difficult to estimate. This increased the classification error in a_p and d_p classes [35]. In addition, the SatPhenClass classification result was accepted only, if both submodels (ETS and LAI comparisons) were in line (Figure 3b).

Table 5. VGI composite linearand non-linear model (I–III) non-potential baseline yield estimates (y_b, kg/ha) for spring cereals [64,65].

Crop

Model

Grow. Zone

Main Soil Type

Dataset, (Model Equation) (1)

R2lin(a)/ non-lin. polyn.(b),(2)

RMSE lin(c)/ non-lin. polyn. (d)(4)

Model mean yield (I–III) kg/ha(2)

Cv (%)

Pr > F linear

Pr > F quadratic

Pr > F cross product

Pr > F total model

Sat Phen Class phonological classification error range (% DOY, Table 2) (3)

Swh

I IR

I–II

Gyttja & Sandy clay

2.1–3.3 (1.1)

0.764(a)

282.3(c)

4219.0

6.66

<.0001 ***

<.0001 ***

<.0001 ***

<.0001 ***

ap 21–29 bp 11–16 cp 14–19 dp 25–32

III–IV

Sandy clay

II (1.2)

0.794(a)

42.46(c)

3768.6

5.89

II NDVI

I–II

Gyttja & Sandy clay

2.1–3.3 (1.3)

0.737(b)

297.6(d)

4219.0

7.02

<.0001 ***

<.0001 ***

<.0001 ***

<.0001 ***

I–II

Gyttja & Sandy clay

2.1–3.3(1.4)

0.732(a)

300.1(c)

3556.7

7.88

III SAR,NDVI(2)(a)

III–IV

Sandy clay

1.1–1.3 (5.1,6.1)

0.723 E 0.731 R

302.1 300.8

4127.0 4213.0

3.78 5.56

<.0001 ***

<.0001 ***

<.0001 ***

<.0001 ***

Brl

I IR

I–II

Gyttja & Sandy clay

2.1–3.3(2.1)

0.615(a)

449.3(c)

4395.0

10.3

<.0001 ***

<.0001 ***

N.S.

<.0001

ap 24–26 bp 15–19 cp 17–21 dp 29–34

IINDVI

I–II

Gyttja & Sandy clay

2.1–3.3(2.2)

0.611(b)

449.6(d)

4298.0

10.3

<.0001 ***

<.0001 ***

0.0014 *

***

IIISAR,NDVI (2)(a)

III–IV

Sandy clay

1.1–1.3 (4.1,5.2,6.2)

0.694 E 0.702 R 0.448 S

349.8 322.8 482.7

3750.0 3909.0 3835.0

3.12 6.92 7.22

<.0001 ***

<.0001 ***

<.0001 ***

<.0001 ***

Oat

I IR

I,II,IV

Sandy clay

1.1,2.1,3.3(2.5)

0.760(a)

55.0. (c)

3740.0

1.58

<.0001 ***

<.0001 ***

<.0001 ***

<.0001

ap 28–34 bp 17–22 cp 19–28 dp 36–39

I,II

Gyttja & Sandy clay

2.1,3.3(e)

0.756

55.1

3488.7

1.58

<.0001 ***

N.S.

N.S.

***

I,II

Gyttja & Sandy clay

2.1,3.3(e)

0.056

994.8

3462.0

28.7

<.0001 ***

<.0001 ***

N.S.

<.0001

III SAR, NDVI(2) (a)

III–IV

Sandy clay

1.1–1.3 (4.2,5.3,6.3) (a)

0.617 E 0.624 R 0.417 S

389.7 483.6 584.2

2826.0 3038.0 2942.0

4.22 7.89 8.47

<.0001 ***

<.0001 ***

<.0001 ***

<.0001 ***

Cereal Mean

0.627

386.9

3752.2

7.47

⁽¹⁾ Table 1, Table 12 (Appendix C) ⁽²⁾ SAR-components: Envisat (E), Radarsat (R), ERS2 (S) SAR backscattering (σ ⁰)⁽³⁾ SatPhenClass phenological classification accuracy error range (% _DOY) ^(a) R² linear ^(b) R² non-linear polynomial equation ^(c) RMSE linear ^(d) RMSE(kg/ha) non-linear polynomial equation ^(e) Equation not shown

Table 5. VGI composite linearand non-linear model (I–III) non-potential baseline yield estimates (y_b, kg/ha) for spring cereals [64,65].

Crop

Model

Grow. Zone

Main Soil Type

Dataset, (Model Equation) (1)

R2lin(a)/ non-lin. polyn.(b),(2)

RMSE lin(c)/ non-lin. polyn. (d)(4)

Model mean yield (I–III) kg/ha(2)

Cv (%)

Pr > F linear

Pr > F quadratic

Pr > F cross product

Pr > F total model

Sat Phen Class phonological classification error range (% DOY, Table 2) (3)

Swh

I IR

I–II

Gyttja & Sandy clay

2.1–3.3 (1.1)

0.764(a)

282.3(c)

4219.0

6.66

<.0001 ***

<.0001 ***

<.0001 ***

<.0001 ***

ap 21–29 bp 11–16 cp 14–19 dp 25–32

III–IV

Sandy clay

II (1.2)

0.794(a)

42.46(c)

3768.6

5.89

II NDVI

I–II

Gyttja & Sandy clay

2.1–3.3 (1.3)

0.737(b)

297.6(d)

4219.0

7.02

<.0001 ***

<.0001 ***

<.0001 ***

<.0001 ***

I–II

Gyttja & Sandy clay

2.1–3.3(1.4)

0.732(a)

300.1(c)

3556.7

7.88

III SAR,NDVI(2)(a)

III–IV

Sandy clay

1.1–1.3 (5.1,6.1)

0.723 E 0.731 R

302.1 300.8

4127.0 4213.0

3.78 5.56

<.0001 ***

<.0001 ***

<.0001 ***

<.0001 ***

Brl

I IR

I–II

Gyttja & Sandy clay

2.1–3.3(2.1)

0.615(a)

449.3(c)

4395.0

10.3

<.0001 ***

<.0001 ***

N.S.

<.0001

ap 24–26 bp 15–19 cp 17–21 dp 29–34

IINDVI

I–II

Gyttja & Sandy clay

2.1–3.3(2.2)

0.611(b)

449.6(d)

4298.0

10.3

<.0001 ***

<.0001 ***

0.0014 *

***

IIISAR,NDVI (2)(a)

III–IV

Sandy clay

1.1–1.3 (4.1,5.2,6.2)

0.694 E 0.702 R 0.448 S

349.8 322.8 482.7

3750.0 3909.0 3835.0

3.12 6.92 7.22

<.0001 ***

<.0001 ***

<.0001 ***

<.0001 ***

Oat

I IR

I,II,IV

Sandy clay

1.1,2.1,3.3(2.5)

0.760(a)

55.0. (c)

3740.0

1.58

<.0001 ***

<.0001 ***

<.0001 ***

<.0001

ap 28–34 bp 17–22 cp 19–28 dp 36–39

I,II

Gyttja & Sandy clay

2.1,3.3(e)

0.756

55.1

3488.7

1.58

<.0001 ***

N.S.

N.S.

***

I,II

Gyttja & Sandy clay

2.1,3.3(e)

0.056

994.8

3462.0

28.7

<.0001 ***

<.0001 ***

N.S.

<.0001

III SAR, NDVI(2) (a)

III–IV

Sandy clay

1.1–1.3 (4.2,5.3,6.3) (a)

0.617 E 0.624 R 0.417 S

389.7 483.6 584.2

2826.0 3038.0 2942.0

4.22 7.89 8.47

<.0001 ***

<.0001 ***

<.0001 ***

<.0001 ***

Cereal Mean

0.627

386.9

3752.2

7.47

⁽¹⁾ Table 1, Table 12 (Appendix C) ⁽²⁾ SAR-components: Envisat (E), Radarsat (R), ERS2 (S) SAR backscattering (σ ⁰)⁽³⁾ SatPhenClass phenological classification accuracy error range (% _DOY) ^(a) R² linear ^(b) R² non-linear polynomial equation ^(c) RMSE linear ^(d) RMSE(kg/ha) non-linear polynomial equation ^(e) Equation not shown

3.2.2. Testing the Composite VGI model (I–III) calibration performance

The spring cereal calibration results derived from phenologically pre-classified (SatPhenClass) SAR and reflectance data, are depicted in Table 5. It contains linear (REG/Stepwise) and non-linear polynomial (RSREG) VGI model (I–III, Table 3) baseline yield estimates (y_b, kg/ha) with corresponding R² and RMSE, C_v (%) and Anova F-test significance level estimates. The calibrated VGI model (I–III) equations are presented in Table 12 (Appendix C).

There was a large variation in VGI Composite model cereal yield responses between years and locations with averaged R² 0.627 and RMSE 387 kg/ha in growing zones I–IV in southern Finland (Table 1). The average modeled yield response varied between 3.5 t/ha and 4.2 t/ha with spring wheat, between 4.1 t/ha and 4.4 t/ha with barley and between 3.4 t/ha and 3.7 t/ha with oats. The coefficient of variation (C_v,%) varied between 1.6% and 28.7% with spring cereals, respectively the VGI model subcomponents (linear, quadratic and cross products) were significant on 0.1% error level.

The calibration results indicate that the R² of the VGI modeled (I–III, Table 3) final grain yield accuracy varied significantly between different spring cereals. The overall R² average for yield (kg/ha) prediction was for spring sown cereals 0.63 (mean response yield 3750 kg/ha, RMSE 387 kg/ha). According to calibration results obtained from the optical VGI models (I–II) and Composite VGI model (III), the R² tends to stabilize on the 0.60–0.70 level. With relatively few observations (n < 500) the R² tends to increase above 0.80, whereas with larger data sets (n > 1,000 observations) the R² varies between 0.50 and 0.70. The overall R² for wheat, barley and oats varied between 0.615 and 0.794 for Infrared model (I), between 0.611 and 0.737 for NDVI (II), and between 0.417 (oats using ERS data) and for 0.730 (wheat using Envisat and Radarsat data) for Composite SAR model (III). The total cultivation area used for the calibration of VGI models was 1253 hectares. The R² varied for wheat ranged between 0.731 (SAR-model (IIII), RMSE 300 kg/ha) and 0.794 (Infrared model (I), RMSE 42 kg/ha), for barley between 0.448 (III) and 0.615 (I) and for oats between 0.417 (III) and 0.760 (I).

The Composite Model III (NDVI * SAR) calibration results for ERS, Radarsat and Envisat are presented in Table 5 (Model III _{SAR, NDVI}). The R² for the Composite model III varied for wheat between 0.723 (Envisat SAR-component, RMSE 302 kg/ha) and 0.731 (Radarsat SAR-component, RMSE 300 kg/ha), for barley between 0.694 (ERS component, RMSE 482 kg/ha) and 0.702 (Radarsat SAR component, RMSE 322 kg/ha) and for oats between 0.417 (ERS SAR component, RMSE 584 kg/ha) and 0.624 (Radarsat, RMSE 483 kg/ha). The corresponding optical calibration results for soil*species and soil _(clay) * cultivar covariances with Infrared Polynomial Model I are reviewed in Part II.

3.3. Validation Results

3.3.1. Soil Species Covariance Validation with Composite Models

The validation results assessing both soil*species (Composite SAR and NDVI Model III) and soil cultivar covariance (Infrared NDVI Models I, II) categories for cereal baseline yield (y_b) are depicted in Table 6 (Model Categories I–II).

In model category I (Table 6), the SAR Composite Model (III) validation results were compared and validated with averaged inventory estimates (kg/ha) from the MAFF Etelä-Pohjanmaa Agricultural Advisory Centre. The SAR and NDVI modeling results depict detailed soil*canopy interactions between cereal canopies and soil top cover (5–10 cm) with fine and coarse sand soil textures. The SAR composite model validation results take into account both the pre-anthesis a_p and b_p phenological phases (Figure 3b, phases I–II) using the NDVI component and post-anthesis and senescence phases (c_p and d_p, phase III) using the SAR component with HH and VV levels.

In model category II (Table 6), the soil cultivar covariance validation results are depicted. Two spring wheat (cv. Manu, Satu) and one malting barley (cv. Inari) cultivars grown on heavy clay type soils and currently cultivated in Finland were analyzed by using optical VGI Infrared and NDVI models [44,56,71]. The response mean yields (kg/ha) were compared with the MTT Official Variety Trial mean yield results correspondingly grown on clay type soils in cultivation zone I (Pernaja and Inkoo Exp. Stations), II (Jokioinen Exp. Station) and III–IV (Ylistaro and Ruukki Exp. Station, Figure 4a).

The Model III soil*species validation results (Model Category I, Table 6) using MAFF validation data, indicated with spring wheat an overestimation (DP_MAFF,%) between +104% (Envisat) and 107% (Radarsat) from the observed MAFF inventory reference (100% [29]). With oat cereals, the underestimation was between 85% (Envisat) and 88% (ERS2). Respectively with barley between 97% (Envisat) and 112% (ERS) and more precisely with malt barley 96% underestimation (Radarsat) vs. fodder barley 102%. Respectively the yield difference calculated as kg/ha (DY_MAFF) ranged between—177 kg/ha (swh, Envisat ASAR) and 484 kg/ha (Oats, Envisat ASAR). The mean composite model baseline yield (y_b, kg/ha) was 4,170 kg/ha for wheat, 3,848 kg/ha for barley (3,909 for malting barley and 3,899 kg/ha for fodder barley) and 3,386 kg/ha for oats.

The Model I and II validation results (Model Category II, Table 6) using MTT validation data, indicated that both observed wheat and barley yield levels from the MTT official variety trials (kg/ha) exceeded corresponding modeled spring wheat and barley yield levels. The soil*species and soil*cultivar covariance validation averaged results with Infrared and NDVI models indicated an underestimated yield difference (DP_MTT, %) by –10% in soil*species category and between –31% and –14% in soil*cultivar category. The modeled results for spring wheat and barley baseline yield levels (y_b) in both categories were lower than the corresponding yield levels in MTT Official Variety trials.

Spring wheat and barley cereals in MTT field trials were grown on more optimal, non-limiting growing conditions. The optical VGI modeling results indicate a general underestimation of modeled vs. MTT observed yield levels. More specifically, the results indicate that yield difference (DP_MTT,%) ranged between 69% [Clay*Swh*(cv. Satu)] and 92% (Barley*clay) from the baseline reference (100%). In both categories (soil*species, soil*cultivar) the VGI models underestimated the baseline yield levels (y_b) when compared with the corresponding observed MTT yield levels. Especially cv. Satu yield level on clay type soils was clearly underestimated (–44.9% DP_MTT). Respectively the yield difference calculated as kg/ha (DY_MTT) varied between –363 kg/ha (barley*clay) and –1,427 kg/ha (Satu*clay).

Table 6. Cereal composite SAR and optical VGI model baseline yield (y_b kg/ha) validation on different soil types vs. MAFF and vs. MTT observed values in the Etelä-Pohjanmaa Agric. Advisory Centre [30,56,71]⁽¹⁾.

Model Category I	Crop	Composite model	Covariance Category	SAR Sensor	Mean yb, kg/ha, X ± Sd	MAFF mean estimate kg/ha (X ± Sd) (2)	DYMAFF Difference kg/ha from MAFF obs.(1) (5)	DPMAFF Difference (%) from MAFF obs. (100 ref.) (4)
I (2) Composite SAR & NDVI (II,III)	Swh	II+III	Fine coarse sandy clay (Table 1)	Envisat ASAR	4127 ± 68	3950 ± 72	–177.0 oe	104.4
		II+III		Radarsat SAR	4213 ± 41	3840 ± 86	–373.0 oe	109.7
		Swh ave.			4170 ± 54	3895 ± 78	–275.0 oe	107.1
	Brl (general)	II+III	Fine coarse sandy clay (Table 1)	Envisat ASAR	3750 ± 91	3880 ± 47	130.0 ue	96.6
		II+III		ERS2 SAR	3835 ± 98	3 420 ± 82	–415.0 oe	112.1
	Brl malt	II+III		Radarsat SAR	3909 ± 24	4050 ± 76	141.0 ue	96.5
	Brl feed & fodder	II+III		Radarsat SAR	3899 ± 32	3820 ± 82	–79.0 oe	102.7
		Brl ave.			3848 ± 74	3792 ± 72	–55.7 oe	101.8
	Oats	II+III	Fine coarse sandy clay (Table 1)	Envisat ASAR	2826 ± 85	3310 ± 54	484.0 ue	85.4
		II+III		ERS2 SAR	2942 ± 49	3 330 ± 54	388.0 ue	88.4
		II+III		Radarsat SAR	3038 ± 23	3520 ± 81	482.0 ue	86.3
		Oats ave.			2935 ± 28	3386 ± 48	451.3 ue	86.7
		Mean tot.			3615 ± 12	3680 ± 49	64.5 ue	98.2
Model Category II	Species	Model	Covariance Category	Soil type/ Cultivar(1)	Averaged model I,II yield yb X ± Sd kg/ha	MTT mean X ± Sd kg/ha (3)	DYMTT Difference kg/ha from MTT obs (7)	DPMTT, Difference (%) from MTT obs (100 ref.) (6)
II (3) Optical Infrared (I), NDVI (II)	Swheat	I+II	Soil* Species	Sandy clay	4240 ± 52	4645 ± 546	–404*ue)	91.3*ue)
	Barley	I+II		Sandy clay	4428 ± 48	4791 ± 29	–363*ue)	92.4*ue)
	Swheat	I+II	Soil*Cultivar	Sand clay*Manu	4015 ± 62	4423 ± 72	–408*ue)	90.8*ue)
	Swheat	I+II		Sandy clay*Satu	3181 ± 31	4608 ± 92	–1427*ue)	69.0*ue)
	Barley	I+II		Sandy clay*Inari	4749 ± 89	5483 ± 44	–733*ue)	86.6*ue)

⁽¹⁾ Yield estimates corrected to 15% moisture content ⁽²⁾ MAFF mean inventory estimate ⁽³⁾ MTT Official Variety Trial averaged results from cultivation zones I–IV (Figure 4a) ⁽⁴⁾ DP_MAFF—difference (%) modeled vs. MAFF inventory estimate ⁽⁵⁾ DY_MAFF - difference modeled vs. MAFF inventory estimate ⁽⁶⁾ DP_MTT —difference (%) modeled vs. MTT estimate ⁽⁷⁾ DY_MTT — difference (kg/ha) modeled vs. MTT estimate ^*(ue) —underestimated by the corresponding model vs. observed ^*(oe) —overestimated by the corresponding model vs. observed.

Table 6. Cereal composite SAR and optical VGI model baseline yield (y_b kg/ha) validation on different soil types vs. MAFF and vs. MTT observed values in the Etelä-Pohjanmaa Agric. Advisory Centre [30,56,71]⁽¹⁾.

Model Category I	Crop	Composite model	Covariance Category	SAR Sensor	Mean yb, kg/ha, X ± Sd	MAFF mean estimate kg/ha (X ± Sd) (2)	DYMAFF Difference kg/ha from MAFF obs.(1) (5)	DPMAFF Difference (%) from MAFF obs. (100 ref.) (4)
I (2) Composite SAR & NDVI (II,III)	Swh	II+III	Fine coarse sandy clay (Table 1)	Envisat ASAR	4127 ± 68	3950 ± 72	–177.0 oe	104.4
		II+III		Radarsat SAR	4213 ± 41	3840 ± 86	–373.0 oe	109.7
		Swh ave.			4170 ± 54	3895 ± 78	–275.0 oe	107.1
	Brl (general)	II+III	Fine coarse sandy clay (Table 1)	Envisat ASAR	3750 ± 91	3880 ± 47	130.0 ue	96.6
		II+III		ERS2 SAR	3835 ± 98	3 420 ± 82	–415.0 oe	112.1
	Brl malt	II+III		Radarsat SAR	3909 ± 24	4050 ± 76	141.0 ue	96.5
	Brl feed & fodder	II+III		Radarsat SAR	3899 ± 32	3820 ± 82	–79.0 oe	102.7
		Brl ave.			3848 ± 74	3792 ± 72	–55.7 oe	101.8
	Oats	II+III	Fine coarse sandy clay (Table 1)	Envisat ASAR	2826 ± 85	3310 ± 54	484.0 ue	85.4
		II+III		ERS2 SAR	2942 ± 49	3 330 ± 54	388.0 ue	88.4
		II+III		Radarsat SAR	3038 ± 23	3520 ± 81	482.0 ue	86.3
		Oats ave.			2935 ± 28	3386 ± 48	451.3 ue	86.7
		Mean tot.			3615 ± 12	3680 ± 49	64.5 ue	98.2
Model Category II	Species	Model	Covariance Category	Soil type/ Cultivar(1)	Averaged model I,II yield yb X ± Sd kg/ha	MTT mean X ± Sd kg/ha (3)	DYMTT Difference kg/ha from MTT obs (7)	DPMTT, Difference (%) from MTT obs (100 ref.) (6)
II (3) Optical Infrared (I), NDVI (II)	Swheat	I+II	Soil* Species	Sandy clay	4240 ± 52	4645 ± 546	–404*ue)	91.3*ue)
	Barley	I+II		Sandy clay	4428 ± 48	4791 ± 29	–363*ue)	92.4*ue)
	Swheat	I+II	Soil*Cultivar	Sand clay*Manu	4015 ± 62	4423 ± 72	–408*ue)	90.8*ue)
	Swheat	I+II		Sandy clay*Satu	3181 ± 31	4608 ± 92	–1427*ue)	69.0*ue)
	Barley	I+II		Sandy clay*Inari	4749 ± 89	5483 ± 44	–733*ue)	86.6*ue)

⁽¹⁾ Yield estimates corrected to 15% moisture content ⁽²⁾ MAFF mean inventory estimate ⁽³⁾ MTT Official Variety Trial averaged results from cultivation zones I–IV (Figure 4a) ⁽⁴⁾ DP_MAFF—difference (%) modeled vs. MAFF inventory estimate ⁽⁵⁾ DY_MAFF - difference modeled vs. MAFF inventory estimate ⁽⁶⁾ DP_MTT —difference (%) modeled vs. MTT estimate ⁽⁷⁾ DY_MTT — difference (kg/ha) modeled vs. MTT estimate ^*(ue) —underestimated by the corresponding model vs. observed ^*(oe) —overestimated by the corresponding model vs. observed.

4. Discussion

4.1. Implications from SAR Soil*Canopy Interactions

Both the SAR backscattering (σ ⁰) signal and SAR and NDVI baseline yield (y_b) levels varied significantly both in cereal soil*species and species*canopy covariance categories. In the species*soil covariance category SAR backscattering signal varied significantly especially on clay type soils with minor fractions of sand, silt and organic mould in the Porvoo, Mellilä, Kirkkonummi, Jokioinen and Lapua experimental areas (Table 1).

With spring wheat the VH cross-polarization amplitude was higher when compared with VV vertical levels both in anthesis (c_p) and full maturity (d_p) on clay and sandy soils. The wheat horizontal signal (HH) amplitude was higher in c_p stage compared to d_p.

The microwave Radarsat VH cross-polarization amplitude in barley plots was higher when compared with VV vertical levels both in anthesis (c_p) and full maturity (d_p) on clay and sandy soils. Both the barley horizontal and vertical (HH, VV) σ ⁰ signal amplitudes were higher in anthesis stage when compared to full maturity stage.

With oat cereals, the VH cross-polarization amplitude was higher when compared with VV vertical levels both in anthesis (c_p) and full maturity (d_p) on clay and sandy soils. Both the horizontal and vertical (HH, VV) signal amplitudes were higher in anthesis stage than in full maturity stage containing grain filling, yellow ripening and canopy senescence stages.

In the species*canopy covariance category level there was a significant variation both in backscattering amplitude and polarization properties between spring wheat, barley and oats. Especially the canopy structure of oat cereals differs morphologically from other spring cereals. Oat cereals with more planophile canopy and panicle inflorescence structures differ in polarization properties from those of wheat and barley with more erectophile head and canopy structures [23,60,72].

4.2. Implications from the SAR Composite Modeling Results for Baseline Yield Levels (Y_b)

The SAR calibration results with Composite NDVI and SAR models indicated that averaged baseline yield (y_b) response varied between 3.5 t/ha and 4.2 t/ha with spring wheat, between 4.1 t/ha and 4.4 t/ha with barley and between 3.4 t/ha and 3.7 t/ha with oats. Correspondingly the SAR validation results compared with MAFF data indicated that with spring wheat the overestimation was between +104% (Envisat) and 107% (Radarsat) from the observed MAFF inventory reference. Respectively with oats the underestimation was between 85% (Envisat) and 88% (ERS2). With barley the underestimation was 97% with Envisat and overestimation 112% with ERS. With malt barley the yield was underestimated (96%) with Radarsat data and overestimated with fodder barley (102%).

The use of Composite Multispectral model (III) combined with the SatPhenClass phenological classification algorithm for spring cereals provides a promising integrating technique for combining both microwave SAR and optical reflectance data. The Composite Multispectral model takes into account both the pre-anthesis phenological phases using the NDVI component and post-anthesis and senescence phases using the SAR component with backscattering polarization levels. In addition, the Composite multispectral model can be used in assessing the soil*canopy covariances between cereal canopies and soil top layers on different soil types.

The Composite model (III) calibration results indicate that the R² tends to stabilize on the 60%–70% level similar to optical VGI models (I–II). The Composite model results indicated that the use of the Envisat ASAR additional cross-polarization component (VH) did not increase the R² level compared to ERS (VV) and Radarsat (HH) with only one polarization level in the 5.3 GHz measurement spectrum. The Radarsat signal with the HH horizontal polarization component yielded highest R² values with wheat, barley and oats. ERS signal with only one VV vertical polarization component produced consistently lower R² values. In Finland Karjalainen et al. [28] reported SAR modeling results with the mean R² of 0.55. The average crop height was used to estimate the amount of biomass using dual-polarization (VV/VH) Envisat SAR data in the Lapua and Seinäjoki experimental sites [28].

Recently McNairn et al. [1] reported the composite VV-VH SAR and optical data to be the most suitable for wheat, maize and soybean classification with over 85% overall accuracy (κ Kappa range 0.47–0.89) in Canadian growing conditions. McNairn et al. [1] applied three primary classification methodologies Neural Networks [73], Gaussian Maximum-Likelihood Classifier [21,74] and Decision trees [75] for crop classification using composite SAR and optical data. Respectively in this study the SatPhenClass classification accuracy for spring cereals varied between 61% and 89% in phenological classes (a_p–d_p). Especially the early vegetative phase before double-ridge induction (a_p) [35,50,51] and post-harvest senescence (d_p) phases were major source for error variation decreasing the overall classification accuracy. Highest classification accuracies (>80%) were obtained during the anthesis (b_p) near the LAI maximum proximity.

The composite model validation results indicated an overestimation with wheat and feed barley cereals. On the contrary, an underestimation was noticed with oats and malting barley cereals by using different SAR/ASAR and NDVI sources. In summary, the results suggest the potential inclusion of horizontal polarization component (HH) in the SAR composite models. The oat canopy and head structures differ from those of wheat and barley potentially changing the backscattering polarization properties [23,28,60].

5. Conclusions

The cereal validation results obtained in this study suggest that spring cereal baseline yield estimates derived from Vegetation Indices Models tend to retain between the observed MAFF yield inventory estimates and MTT Agrifood Official Variety trial data results. The annual MAFF non-potential yield inventory estimates are obtained from suboptimal growing conditions on actual farm level, where as MTT Agrifood Official Variety trial results are obtained from more optimal growing conditions.

Composite multispectral models combined with SatPhenClass phenological classification algorithm provide a new promising integrating technique by combining SAR and optical data with the cereal phenological development both in vegetative and generative phases. This new methodology can be applied in future large area crop monitoring campaigns and also in national crop inventory programs for assessing grain yield production on the Agricultural Advisory Centre level. The use of VGI models combined with satellite based Minimum Datasets would significantly reduce the operational costs and time required for more accurate cereal inventory estimates when compared with the traditional ground based crop inventory systems. The implications of this study for Finnish national inventory system and CAP policy are reviewed in the second part of this publication.