The impact of urban expansion in Beijing and Metropolitan Area urban heat Island from 1999 to 2019

Land use cover change (LUCC) has an important impact on the surface-atmosphere system. Urban expansion is one of the most prominent hotspots in China’s current research on LUCC temporal and spatial pattern changes. And the heat island effect caused by urban expansion is also a typical representative of regional climate changes caused by land-use changes (Lu et al., 2020). The current global urbanisation process is becoming more and more intense. While the process of urbanisation brings economic and social benefits to mankind, it also produces a series of ecological and environmental problems. The urban heat island effect is one of the most important problems. With the rapid development of China’s economy, Beijing, the capital of politics, culture, and economy, has seen significant urban expansion and steady population growth, so the thermal environment is bound to be affected (Cui et al., 2016). This paper is a quantitative study of the extent of urban expansion in Beijing in the past two decades and its impact on the city’s heat island effect.

The overall methodology of this paper is following a description of basic facts, the revelation of phenomena and laws, basic cause analysis, mechanism analysis and model simulation. This paper utilises Bigmaper Beijing Sixth Ring Road remote sensing image data as the basic data, with ArcGIS to retrieve the surface temperature of Beijing. The heat island effect intensity index is also calculated. Then it can be used to quantitatively analyse the temporal and spatial distribution characteristics and development trend of the urban heat island effect in Beijing, combining with some factors, including land surface temperature (LST), normalised difference vegetation index (NDVI) and normalised difference building index (NDBI), to explore the impact of surface vegetation distribution and building density on the urban heat island effect. The results show that from 1999 to 2014, the area of urban construction land in Beijing continued to grow, and the growth rate was relatively fast from 2000 to 2003. The rate of urban expansion slowed down significantly from 2004 to 2012. Based on the core functional area, we can say the city gradually expands to the surroundings, showing that the external development is centred on the main urban area. The temperature of urban land gradually started to dominate the temperature level of the whole region, and the difference between the temperature of urban land and the average temperature of the region has decreased significantly. The urban heat island effect of Beijing has gradually increased from 2014 to 2019, and the spread of the heat island is mainly the extension of the urban new development zone and the transportation network. The urban surface temperature has an obvious negative correlation with the normalised vegetation index and an obvious positive correlation with the normalised building index. Further, the correlation between the surface temperature and the normalised vegetation index is particularly prominent. It is recommended that in the process of urban development and construction, increasing urban greening investment, rationally increasing urban vegetation coverage and moderately reducing urban building density can alleviate the urban heat island effect.

In this study, the model that analyse the impact of urban expansion on the urban environment can be determined using literature. From a quantitative point of view, the urban expansion process and basic characteristics of climate change in Beijing for the past 20 years can be analysed through the establishment of multi-period LUCC and meteorological data. Then, the remote sensing data is used to analyse the spatial distribution pattern of urban heat islands and the underlying surface, by mainly focusing on the correlation analysis of LST, NDVI and NDBI. In this paper, we have proposed to focus on analysing the relationship between urban building density distribution and the urban heat island effect. This paper takes Beijing as the research object. Based on high-resolution remote sensing images and through artificial visual interpretation, combined with the surface temperature data obtained from remote sensing inversion, we correlate the internal relationship between the urban building density distribution and the urban heat island effect. This paper discusses the impact of urban building density distribution on the urban heat island effect. It can provide a scientific basis for the planning and layout of urban construction and urban thermal environmental management.

Encroachment of urban lands causes constrained evaporation which arises from decreased thermal inertia and the vegetation index, and thus ultimately leading to reduced heat loss by latent heat flux; it is the flux of heat from the surface to the atmosphere that is associated with evaporation of water and subsequent condensation of water vapour in the troposphere (Goswami & Singh, 2008). UHI is further aggravated by reduction in speed of the environmental winds, congested transportation networks, increased energy demands, higher anthropogenic heat release and other sky-view factors (Chang Cao et al., 2016). Since urban development, occasioned by increased migration into the metropolitan areas, is relatively difficult to control, the only good urban form is the only way which can help sustain the metropolitan area’s environmental condition.

Beijing is one of the leading growth drivers in the country getting some of the country’s huge investments. Besides the manufacturing, processing and manufacturing industries, Beijing’s growth also draws from the real estate segment that forms an integral economic zone, whose growth started from the mid-1990s. The city forms part of the coastal economic zone even though it is geographically located on the North China Plain (Chang Cao et al., 2016). The Ring and Radial Highway System stands as the city’s ideal transportation mechanism model. The system has received a great structural boost since its creation in the 1950s through new legislations like the 1982 and 1993 comprehensive plans for Beijing’s spatial structure. The intercity highway connecting the city’s 14 satellite towns is known as the 6th ring, while the 4th ring road acts as the edge of the city centre and 5th ring road as the linkage between 10 other scattered districts within the metropolitan area and the lower the number of rings, the closer it is to the centre, which also represents major buildings in these areas. Following the adoption of economic reform policies, China’s land use/cover has undergone radical changes that also accelerated economic growth.

As reported by Qiao et al. (2014), from 1989 to 2010, Beijing’s urban land grew by 775.82 km² at a rate of 184.31%. The mounting intensity and increasing UHI index grew at a high pace from 1989 to 2000 but at a declining pace from 2000 to 2010. The study also reports that ‘the LSTs and urban heat island ratio index (URIs) of urban land in Beijing increased between 1989 and 2010, heat island areas expanded, and the UHI resulting from urban expansion increased’ (Qiao et al., 2014).

China’s land legalities and policies have changed rapidly consistent with the changes in the social and economic conditions. As discussed by Liu, Dunford, Song and Chen (2016), most of these policies are ‘reflected in Hukou reform, urban and rural rights, the Annual Land Use Quota system, public assets, and the local government finance vehicles.’

Committed to Habitat II in 1996 principle, the government or china implemented radical measures associated urban and rural housing. The goal of the policy was to ensure all households have self-contained rooms and to achieve a 9-metre per-capita construction space (Assembly, 2015). The policies have led to an increase in the overall level of urban housing in China. The yearly and persistent large-scale housing construction has reached an unparalleled level in China such that houses constructed increased from 1.22 billion in 1996 to 1.93 billion in 2013. During the same time housing space grew from 395 million to 1.07 billion m² in urban centres (Assembly, 2015).

Internationally recognised l-use regulations such as ‘minimum lot sizes, minimum parking requirements, maximum floor-to-area ratios or floor space index (FSI)’ have had a significant impact on UHI in Chinese cities such as Beijing (Menon et al., 2019). For instance, while FSI limits reduce the ‘density of buildings”, they cannot limit the ‘density of people’ since most individuals opt for smaller and especially informal spaces. Prohibitions on lower height structures inflated prices causing urban sprawl such as the case of Beijing, which caused a 12% expansion of its city boundaries (Deng & Huang, 2004).

China’s land use/cover changes have led to the subdivision of surface structures into three layers based on the specific area’s land use and climate (Siddique et al., 2020). Vegetation cover often acts as a sieve that purifies the air from dust particles and heating thereby improving evapotranspiration which in turn promotes precipitation. In Beijing, after the air passes through the green belt, it is intercepted by the city’s built-up areas and experience solar radiation and absorption, and as a result greenhouse gas emission occurs at the cost of anthropogenic disturbances which makes it warmer, hence urban heat island (Yuan, 2020).

With the intensification of urbanisation, the building area continues to expand, and human activity has increased significantly. These factors may lead to an increase in the intensity of the heat island.

Jing (2019) studied the current situation of the heat island effect in the main urban area of Beijing’s Fifth Ring Road in the summer of 2017, and selected 65 urban green spaces of different sizes and types within Beijing’s Fifth Ring Road through a large amount of data linear fitting and controlled variable method on the influencing factors of the internal temperature of urban green space and the quantitative and influencing factors of the green space cooling range. The study suggests that there is a relatively obvious heat island effect in the main urban area of Beijing. The proportion of high-temperature areas from the Second Ring to the Fifth Ring gradually decreases, and the proportion of low-temperature areas gradually increases. The southwest and southeast quadrants of the study area have the highest proportion of high temperature, and the northeast and the northwest quadrant have the highest proportion of low-temperature areas. Large areas of high-temperature heat island areas are mostly old bungalow residential areas, factories, wholesale markets, train stations, etc. The analysis of the factors affecting the cooling range of green space concludes that the smaller the green area, the smaller the cooling range. When the green area is between 0 ha and 7 ha or the perimeter of the green area is between 0 km and 1.5 km, the cooling range of the green area will follow the increase of green area increases. The percentage of impervious surface in the surrounding land of urban green space also has a significant impact on the cooling range of urban green space. In the case of the same green area, whenever the percentage of the impervious surface area surrounding the green area increases by 1%, the green area cooling range will be decreased by 4.01 m.

Li et al. (2008) used the temperature data from the meteorological observatory in Beijing from 1990 to October 2004 to analyse the characteristics of the urban heat island in autumn in Beijing for the past 15 years. The results showed that the urban heat island in the autumn night in Beijing is stronger than that in the daytime. Besides, a comparative analysis of the characteristics of a strong heat island, a weak heat island and their meteorological influence factors show that there is a strong heat island under certain conditions at night in Beijing in autumn. The formation and maintenance of a strong heat island is the result of the combined effects of multiple factors. During the day and sunny nights, the surface wind field in the suburbs of Beijing is very weak and the atmospheric wind field <47 m in the vertical direction of the urban area continues to be very strong. The rate and amplitude of the atmospheric temperature drop in the suburbs after sunset are much greater than that in the urban area, which promotes the formation and maintenance of a strong heat island at night. The weakening of the urban atmospheric stability and the disappearance of the urban atmospheric inversion are the main reasons for the weakening of strong heat island and its final disappearance at night.

Based on discussing the concept and characteristics of urbanisation, Zhao et al. (2014) analysed the interaction between urbanisation and climate change. Urbanisation has caused significant effects on different climatic elements, but related research is mainly concentrated in terms of local microclimate changes. Climate change is a very complex process.

Assuming the experimental data (x_i,y_i),(i = 1,2,...,n), find the function f (x,ŷ), and by guaranteeing the position of the function at the point XI, where I = 1, 2... N, the sum of squares of the deviation between the value of the function and the value of the observation in this region should be minimised, then the function conforming to the following conditions should be calculated: ∑i=1n(f(Xi,y^)−yi)2 \sum\limits_{i = 1}^n {\left({f\left({{X_i},\hat y} \right) - {y_i}} \right)^2}

Therefore, the nonlinear fitting method should be used to solve the impact of the urbanisation process in Beijing on the urban warming. The following steps should be followed: First, the scatter diagram should be defined and the function category should be analysed; Second, the initial value of the undetermined parameters is obtained and the best parameters are calculated by using MATLAB software. Third, the determination coefficient is used to compare the effect. The formula of determination coefficient is as follows: R2=1−∑i=1n(yi−y^i)2∑i=1n(yi−y¯)2 {R^2} = 1 - {{\sum\limits_{i = 1}^n {{\left({{y_i} - {{\hat y}_i}} \right)}^2}} \over {\sum\limits_{i = 1}^n {{\left({{y_i} - \bar y} \right)}^2}}} where y¯=1n∑i=1nyi \bar y = {1 \over n}\sum\limits_{i = 1}^n {y_i} , and the closer R² is to 1, the higher the actual fitting effect is. In MATLAB software to achieve the command of the coefficient of determination: R2=1−sum((y−y1).Λ2)/sum((y−mean(y)).Λ2) R2 = 1 - sum(\left({y - y1} \right){.^\Lambda}2)/sum((y - mean(y)){.^\Lambda}2)

If the coefficients are polynomial, it is called polynomial regression and the parameters are polynomial coefficients. But if it is of the following forms namely exponential function, logarithmic function, power function, trigonometric function, etc., then it is called nonlinear fitting. Equations with an S-shaped curve include the Rogersti model: y=a1+βe−γx y = {a \over {1 + \beta {e^{- \gamma x}}}} , Gompertz model: y = aexp(−β e^−kx), Richards model: y = a/[1 + exp(β − γx)]^1/δ, and Weibull model: y = a − β exp(−γt^δ).

In order to better realise nonlinear fitting, the function should be defined first: the definition function of inline can be used for numerical calculation of curve fitting. On the one hand, the M file should be built first, and on the other hand, f un = inline(′ f (x)′,′ DependingOnAParameter′,′ x′) should be obtained. In practical problems, in many cases, the scatter diagram is not necessarily the graph of the polynomial, but also may be other curves, which requires the combination of [beta,r,J] = nl inf it(x,y, f un,beta0) command to make guess the type. In the above equation, x and y represent the original data, fun represents the defined function in the M file, beta0 refers to the initial value of the parameter in the function, beta represents the optimal value of the parameter, r refers to the fitting residuals of each point and J represents the value of the Jacobian matrix.

After combining the factors that influence the urbanisation process of Beijing on the urban temperature increase, the scatter plot is constructed by combining the above nonlinear fitting method, and an appropriate nonlinear regression equation can be selected, as shown below: y=μ(x¯,a,b…) y = \mu (\bar x,a,b \ldots) where a and b represent unknown parameters, and to calculate the estimated values of both, on the one hand, a linear model can be used for analysis, such as y = b₀ + b₁t + b₂t² and x₁ = t,x₂ = t², to obtain y = b₀ + b₁x₁ + b₂x₂. In this case, x₁ and x₂ can be regarded as independent variables, and y can be regarded as a linear function of the two variables. On the other hand, it cannot be directly converted into a linear model to calculate parameters, such that Y = B₁ ln(1 + B₂V), and B1 and B2 are undetermined parameters, V represents urban warming and Y represents the urbanisation process and so the formula can be obtained as follows: {b=∑i=1N(xi−x¯)(yi−y¯)∑i=1N(xi−x¯)2=∑i=1Nxiyi−1N∑i=1Nxi∑i=1Nyi∑i=1Nxi2−1N(∑i=1Nxi)2a=y¯−bx¯ \left\{{\matrix{{b = {{\sum\limits_{i = 1}^N \left({{x_i} - \bar x} \right)\left({{y_i} - \bar y} \right)} \over {\sum\limits_{i = 1}^N {{\left({{x_i} - \bar x} \right)}^2}}} = {{\sum\limits_{i = 1}^N {x_i}{y_i} - {1 \over N}\sum\limits_{i = 1}^N {x_i}\sum\limits_{i = 1}^N {y_i}} \over {\sum\limits_{i = 1}^N x_i^2 - {1 \over N}{{(\sum\limits_{i = 1}^N {x_i})}^2}}}} \hfill \cr {a = \bar y - b\bar x} \hfill \cr}} \right.

In this case, x¯=1N∑i=1Nxi \bar x = {1 \over N}\sum\limits_{i = 1}^N {x_i} , y¯=1N∑i=1Nyi \bar y = {1 \over N}\sum\limits_{i = 1}^N {y_i} . After the values of A and B are determined, the regression equation y = a + bx can be obtained, and B refers to the regression coefficient.

When calculating the least square method of data fitting, the concept of this theory should be clarified first. Assume that is a set of data given (x_i,y_i)(i = 0,1,...,m), and ω_i (i = 0,1,...,m) represents the weight coefficient of all points, which in general needs to be >0, then in the function class: φ=span{ϕ0(x),ϕ1(x),...,ϕn(x)} \varphi = span\left\{{{\phi _0}\left(x \right),{\phi _1}\left(x \right),...,{\phi _n}\left(x \right)} \right\}

And the calculation function can be obtained: S*(x)=∑j=0naj*ϕj(x),(n≤m) {S^*}\left(x \right) = \sum\limits_{j = 0}^n a_j^*{\phi _j}\left(x \right),\left({n \le m} \right)

More in line with the following requirements: ‖δ‖22=∑i=0mωi(S*(xi)−yi)2=minS∈φ∑i=0m(S(xi)−yi)2 \left\| \delta \right\|_2^2 = \sum\limits_{i = 0}^m {\omega _i}{\left({{S^*}\left({{x_i}} \right) - {y_i}} \right)^2} = \mathop {\min}\limits_{S \in \varphi} \sum\limits_{i = 0}^m {\left({S\left({{x_i}} \right) - {y_i}} \right)^2}

The Φ in the S(x)=∑j=0najϕj(x) S\left(x \right) = \sum\limits_{j = 0}^n {a_j}{\phi _j}\left(x \right) formula represents all functions. Based on the above analysis, represents the least square solution of S* (x), and S (x) refers to the fitting function.

Assuming Ψ(a0,a1,...,an)=∑i=0mωi(S(xi)−yi)2=∑i=0mωi(∑j=0najϕj(xi)−yi)2 \Psi \left({{a_0},{a_1},...,{a_n}} \right) = \sum\limits_{i = 0}^m {\omega _i}{\left({S\left({{x_i}} \right) - {y_i}} \right)^2} = \sum\limits_{i = 0}^m {\omega _i}{\left({\sum\limits_{j = 0}^n {a_j}{\phi _j}\left({{x_i}} \right) - {y_i}} \right)^2} , then the point of calculating the least squares solution is to specify the undetermined parameter aj*(j=0,1,...,n) a_j^*\left( {j = 0,1,...,n} \right) . The important conditions to be observed by the extreme value of the function are: ∂Ψ∂ak=0,(k=0,1,…,n) {{\partial \Psi} \over {\partial {a_k}}} = 0,(k = 0,1, \ldots,n)

Make ϕr=(ϕr(x0),ϕr(x1),...,ϕr(xm)),r=0,1,...,nf=(y0,y1,...,ym)∑j=0n(ϕj,ϕk)aj=(f,ϕk),(k=0,1,...,n) \matrix{{{\phi _r} = \left({{\phi _r}\left({{x_0}} \right),{\phi _r}\left({{x_1}} \right),...,{\phi _r}\left({{x_m}} \right)} \right),r = 0,1,...,n} \hfill \cr {f = \left({{y_0},{y_1},...,{y_m}} \right)} \hfill \cr {\sum\limits_{j = 0}^n \left({{\phi _j},{\phi _k}} \right){a_j} = \left({f,{\phi _k}} \right),\left({k = 0,1,...,n} \right)} \hfill \cr}

And define the inner product: (ϕj,ϕk)=∑i=0mωiωj(xi)ϕk(xi)(f,ϕk)=∑i=0mωiyiϕk(xi) \matrix{{\left({{\phi _j},{\phi _k}} \right) = \sum\limits_{i = 0}^m {\omega _i}{\omega _j}\left({{x_i}} \right){\phi _k}\left({{x_i}} \right)} \hfill \cr {\left({f,{\phi _k}} \right) = \sum\limits_{i = 0}^m {\omega _i}{y_i}{\phi _k}\left({{x_i}} \right)} \hfill \cr}

The system can be called: ∑j=0n(ϕj,ϕk)aj=(f,ϕk),(k=0,1,...,n) \sum\limits_{j = 0}^n \left({{\phi _j},{\phi _k}} \right){a_j} = \left({f,{\phi _k}} \right),\left({k = 0,1,...,n} \right)

Then the normal equations of the function system at discrete points are as follows: [(ϕ0,ϕ0)(ϕ0,ϕ1)…(ϕ0,ϕn)(ϕ1,ϕ0)(ϕ1,ϕ1)…(ϕ1,ϕn)••••••••••••(ϕn,ϕ0)(ϕn,ϕ1)…(ϕn,ϕn)][a0a1•••an]=[(f,ϕ0)(f,ϕ1)•••(f,ϕn)] \left[ {\matrix{{({\phi _0},{\phi _0})} & {({\phi _0},{\phi _1})} & \ldots & {({\phi _0},{\phi _n})} \cr {({\phi _1},{\phi _0})} & {({\phi _1},{\phi _1})} & \ldots & {({\phi _1},{\phi _n})} \cr \matrix{\bullet \hfill \cr \bullet \hfill \cr \bullet \hfill \cr} & \matrix{\bullet \hfill \cr \bullet \hfill \cr \bullet \hfill \cr} & \matrix{\bullet \hfill \cr \,\,\,\, \bullet \hfill \cr \,\,\,\,\,\,\,\, \bullet \hfill \cr} & \matrix{\bullet \hfill \cr \bullet \hfill \cr \bullet \hfill \cr} \cr {({\phi _n},{\phi _0})} & {({\phi _n},{\phi _1})} & \ldots & {({\phi _n},{\phi _n})} \cr}} \right]\left[ {\matrix{{{a_0}} \cr {{a_1}} \cr \matrix{\bullet \hfill \cr \bullet \hfill \cr \bullet \hfill \cr} \cr {{a_n}} \cr}} \right] = \left[ {\matrix{{(f,{\phi _0})} \cr {(f,{\phi _1})} \cr \matrix{\bullet \hfill \cr \bullet \hfill \cr \bullet \hfill \cr} \cr {(f,{\phi _n})} \cr}} \right]

The determinant formed according to the basis function ϕ_r = (ϕ_r (x₀),ϕ_r (x₁),...,ϕ_r (x_m)),r = 0,1,...,n= is: ∑j=0n(ϕj,ϕk)aj=(f,ϕk),(k=0,1,...,n) \sum\limits_{j = 0}^n \left({{\phi _j},{\phi _k}} \right){a_j} = \left({f,{\phi _k}} \right),\left({k = 0,1,...,n} \right)

Therefore, the solution of the system is aj=aj*,(j=0,1,...,n) {a_j} = a_j^*,\left({j = 0,1,...,n} \right)

And S*(x)=∑j=0naj*ϕj(x) {S^*}\left(x \right) = \sum\limits_{j = 0}^n a_j^*{\phi _j}\left(x \right) refers to the least square solution, ‖δ‖22=∑i=0mωi(S*(xi)−yi)2 \left\| \delta \right\|_2^2 = \sum\limits_{i = 0}^m {\omega _i}{\left({{S^*}\left({{x_i}} \right) - {y_i}} \right)^2} refers to the square error of the multiplication solution and ‖δ‖2=∑i=0mωi(S*(xi)−yi)2 {\left\| \delta \right\|_2} = \sqrt {\sum\limits_{i = 0}^m {\omega _i}{{\left({{S^*}\left({{x_i}} \right) - {y_i}} \right)}^2}} represents the mean square error. The specific formula is as follows: ‖δ‖22=|(f,f)−(S*,f)|=|∑i=0mωiyi2−∑k=0nak*(ϕk,f)| \left\| \delta \right\|_2^2 = \left| {\left({f,f} \right) - \left({{S^*},f} \right)} \right| = \left| {\sum\limits_{i = 0}^m {\omega _i}y_i^2 - \sum\limits_{k = 0}^n a_k^*\left({{\phi _k},f} \right)} \right|

In the case that the basis is ϕj (x) = x ^j (j = 0,1,...,n), the relevant normal equations are as follows: [∑i=0mωi∑i=0mωixi…∑i=0mωixin∑i=0mωixi∑i=0mωixi2…∑i=0mωixin+1••••••••••••∑i=0mωixin∑i=0mωixin+1…∑i=0mωixi2n][a0a1an•••]=[∑i=0mωiyi∑i=0mωixiyi•••∑i=0mωixinyi] \left[ {\matrix{{\sum\limits_{i = 0}^m {\omega _i}} & {\sum\limits_{i = 0}^m {\omega _i}{x_i}} & \ldots & {\sum\limits_{i = 0}^m {\omega _i}x_i^n} \cr {\sum\limits_{i = 0}^m {\omega _i}{x_i}} & {\sum\limits_{i = 0}^m {\omega _i}x_i^2} & {\ldots} & {\sum\limits_{i = 0}^m {\omega _i}x_i^{n + 1}} \cr \bullet & \bullet & \bullet & \bullet \cr \bullet & \bullet & \bullet & \bullet \cr \bullet & \bullet & \bullet & \bullet \cr {\sum\limits_{i = 0}^m {\omega _i}x_i^n} & {\sum\limits_{i = 0}^m {\omega _i}x_i^{n + 1}} & \ldots & {\sum\limits_{i = 0}^m {\omega _i}x_i^{2n}} \cr}} \right]\left[ {\matrix{{{a_0}} \cr {{a_1}{a_n}} \cr \bullet \cr \bullet \cr \bullet \cr}} \right] = \left[ {\matrix{{\sum\limits_{i = 0}^m {\omega _i}{y_i}} \cr {\sum\limits_{i = 0}^m {\omega _i}{x_i}{y_i}} \cr \matrix{\bullet \hfill \cr \bullet \hfill \cr \bullet \hfill \cr} \cr {\sum\limits_{i = 0}^m {\omega _i}x_i^n{y_i}} \cr}} \right]

When calculating the data fitting method of least square integration, it is necessary to make clear its concrete concept. First, it is clear that n groups of observed data are (x_i,y_i)i = 1,2,...,n, and they are all obtained by the same model: y(x)=g(x)+ε(x),x1=a<x<b=xn y\left(x \right) = g\left(x \right) + \varepsilon \left(x \right),{x_1} = a < x < b = {x_n}

Wherein, G (x) refers to the error function of function ɛ (X), Y (x) refers to the function obtained by interference of G (x) when function ɛ (X) and y(x_i) = y_i,ɛ_i = ɛ (X_i),X_i ∈ [a,b](i = 1,2,..,n). Here, yi represents the observed value, and ɛ_j represents the random error of observation.

Assuming that ĝ(X) is a data fitting function with p parameters, it can be calculated as: e(X)=y(X)−g^(X)ek=yk−g^(Xk) \matrix{{e\left(X \right) = y\left(X \right) - \hat g\left(X \right)} \hfill \cr {{e_k} = {y_k} - \hat g\left({{X_k}} \right)} \hfill \cr}

Then in the fitting problem of the traditional least square method, the following formula can reach the minimum value if all the parameters contained therein are calculated according to the above formula. This method is called the least square sum method. S=∑k=1nek2=∑k=1n(yk−g^(Xk))2 S = \sum\limits_{k = 1}^n e_k^2 = \sum\limits_{k = 1}^n {\left({{y_k} - \hat g\left({{X_k}} \right)} \right)^2}

By taking the observed data (x_i,y_i)i = 1,2,...,n as a point on the y (x) function, and assuming that ɛ (X), y (x) and g (x) are all continuous functions in combination with model analysis, the fitting function can be transformed into: Q=∫aek2(x)dx=∫a(y(x)−g^(x))2dx Q = \int_a e_k^2\left(x \right)dx = \int_a {\left({y\left(x \right) - \hat g\left(x \right)} \right)^2}dx

In general, the fitting function G (x) represents: g(x)=F(c,d,x)=c1φ1(d,x)+c2φ2(d,x)+...+ciφi(d,x) g\left(x \right) = F\left({c,d,x} \right) = {c_1}{\varphi _1}\left({d,x} \right) + {c_2}{\varphi _2}\left({d,x} \right) +... + {c_i}{\varphi _i}\left({d,x} \right)

Where φ_i (d,x),i = 1,2,...,l, belongs to a group of linearly independent functions, also known as the basis function, which refers to the nonlinear function of X, while c = (c₁,c₂,...,c_i) represents the linear fitting parameter and d = (d₁,d₂,...,d_q) refers to the nonlinear fitting parameter. In this case, if you satisfy the condition φi(d,x) = e^d_ix,i = 1,2,...,l, then you get g(x)=∑i=11cie−dix g\left(x \right) = \sum\limits_{i = 1}^1 {c_i}{e^{- {d_i}x}} .

Here, we first analyse the specific situation when g (x) is used as a linear fitting function, and then we can get: g(x)=c1φ1(x)+c2φ2(x)+...+ciφi(x) g\left(x \right) = {c_1}{\varphi _1}\left(x \right) + {c_2}{\varphi _2}\left(x \right) +... + {c_i}{\varphi _i}\left(x \right)

Where, c_i,i = 1,2,...,l refers to the regression coefficient, also known as the linear regression model. According to the following formula: Q=∫abe2(x)dx=∫ab(y(x)−g(x))2dx=∫ab(y(x)−c1φ1(x)−c2φ2(x)−...−ciφi(x))2dx Q = \int_a^b {e^2}\left(x \right)dx = \int_a^b {\left({y\left(x \right) - g\left(x \right)} \right)^2}dx = \int_a^b {\left({y\left(x \right) - {c_1}{\varphi _1}\left(x \right) - {c_2}{\varphi _2}\left(x \right) -... - {c_i}{\varphi _i}\left(x \right)} \right)^2}dx

Calculate the minimum value, which is ∂Q∂ci=0 {{\partial Q} \over {\partial {c_i}}} = 0 , i = 1, 2, … l, and you can get: ∑i=11[∫abφi(x)φk(x)dx]ci=∫aby(x)φk(x)dx,k=1,2,...l \sum\limits_{i = 1}^1 \left[ {\int_a^b {\varphi _i}\left(x \right){\varphi _k}\left(x \right)dx} \right]{c_i} = \int_a^b y\left(x \right){\varphi _k}\left(x \right)dx,k = 1,2,...l

Assuming that (f,g)=∫abf(x)g(x)dx \left({f,g} \right) = \int_a^b f\left(x \right)g\left(x \right)dx == is regarded as the inner product of an integral, then the above formula can be solved: [c1c2•••cl]=(X′X)−1X′Y \left[ {\matrix{{{c_1}} \hfill \cr {{c_2}} \hfill \cr \bullet \hfill \cr \bullet \hfill \cr \bullet \hfill \cr {{c_l}} \hfill \cr}} \right] = {\left({X'X} \right)^{- 1}}X'Y

The X=((φ1,•••φ1)..•••.(φ1,•••φl)(φl,φ1)...(φl,φl)) X = \left({\matrix{{({\varphi _1}\mathop,\limits_{\bullet \matrix{{} \cr {\bullet \matrix{{} \cr \bullet \cr}} \cr}} {\varphi _1}).\mathop.\limits_{\bullet \matrix{{} \cr {\bullet \matrix{{} \cr \bullet \cr}} \cr}}.({\varphi _1}\mathop,\limits_{\bullet \matrix{{} \cr {\bullet \matrix{{} \cr \bullet \cr}} \cr}} {\varphi _l})} \hfill \cr {\left({{\varphi _l},{\varphi _1}} \right)...\left({{\varphi _l},{\varphi _l}} \right)} \hfill \cr}} \right) , Y=((y,φ1)(y,•••φ2)(y,φl)) Y = \left({\matrix{{\left({y,{\varphi _1}} \right)} \hfill \cr {(y\mathop,\limits_{\bullet \matrix{{} \cr {\bullet \matrix{{} \cr \bullet \cr}} \cr}} {\varphi _2})} \hfill \cr {\left({y,{\varphi _l}} \right)} \hfill \cr}} \right)

In data fitting analysis, two of the most common fitting function is: φ_i (d,x) = h_i (x)(i = 1,2,...l) on the one hand, the basis function is orthogonal polynomial fitting function, according to the basis function analysis, the h_i (x) refers to the weight function of 1 integral inner product of orthogonal polynomials, and as a result of the (x) function hi conform to the requirements of the ∫abhi(x)hj(x)dx(i≠j) \int_a^b {h_i}\left(x \right){h_j}\left(x \right)dx\left({i \ne j} \right) , so can be combined with least square integral: [∫abh12(x)dx∫abh22(x)dx...∫abhl2(x)dx][c1c2•••cl]=[∫abh1(x)y(x)dx∫abh2(x)y(x)dx•••∫abhl(x)y(x)dx] \left[ {\matrix{{\int_a^b h_1^2\left(x \right)dx} \hfill \cr {\int_a^b h_2^2\left(x \right)d{x_{...}}} \hfill \cr {\int_a^b h_l^2\left(x \right)dx} \hfill \cr}} \right]\left[ {\matrix{{{c_1}} \hfill \cr {{c_2}} \hfill \cr \bullet \hfill \cr \bullet \hfill \cr \bullet \hfill \cr {{c_l}} \hfill \cr}} \right] = \left[ {\matrix{{\int_a^b {h_1}\left(x \right)y\left(x \right)dx} \hfill \cr {\int_a^b {h_2}\left(x \right)y\left(x \right)dx} \hfill \cr \bullet \hfill \cr \bullet \hfill \cr \bullet \hfill \cr {\int_a^b {h_l}\left(x \right)y\left(x \right)dx} \hfill \cr}} \right]

The corresponding regression coefficient is: ci=∫abhi(x)y(x)dx∫abhi2dx,i=1,2,...l {c_i} = {{\int_a^b {h_i}\left(x \right)y\left(x \right)dx} \over {\int_a^b h_i^2dx}},i = 1,2,...l

On the other hand, the basis function is an exponential fitting function. First observe the data (t,y_i)(t = 1,2,...,n), and then confirm the m-th algebraic equation with the coefficient of the first term 1 as follows: P(λ)=λm+a1λm−1+...+am−1λ+am=0 P\left(\lambda \right) = {\lambda ^m} + {a_1}{\lambda ^{m - 1}} +... + {a_{m - 1}}\lambda + {a_m} = 0

In addition, the above coefficient Ai is defined. At the same time, the root λ₁,λ₂,...,λ_m of the above formula is regarded as the basis function of X return trip. Then, the fitting function can be obtained as: g(x)=cλ1x+cλ2x+...+cλmx g\left(x \right) = c\lambda _1^x + c\lambda _2^x +... + c\lambda _m^x

According to different characteristic root analysis, clear λix \lambda _i^x . In the above formula, it is clear that x=r+k, and both r and k are positive integers, and multiply both sides by am-r, and then sum r from 1 to m. In this process, g(r + k) = y_r+k,a₁ = 1 so we can get: ∑r=0myr+kam−k=0,k=0,1,...,m \sum\limits_{r = 0}^m {y_{r + k}}{a_{m - k}} = 0,k = 0,1,...,m

And a₁,a₂,...,a_m is calculated according to the following linear equation: {ym−1a1+ym−2a2+...+y0am=−ymyma1+ym−1a2+...+y1am=−ym+1...................................................yn−1a1+yn−2a2+...+yn−mam=−yn \left\{{\matrix{{{y_{m - 1}}{a_1} + {y_{m - 2}}{a_2} +... + {y_0}{a_m} = - {y_m}} \hfill \cr {{y_m}{a_1} + {y_{m - 1}}{a_2} +... + {y_1}{a_m} = - {y_{m + 1}}} \hfill \cr {...................................} \hfill \cr {{y_{n - 1}}{a_1} + {y_{n - 2}}{a_2} +... + {y_{n - m}}{a_m} = - {y_n}} \hfill \cr}} \right.

Assuming n = 2m − 1, and det A ≠ 0, the Gauss elimination method can be used to calculate the solution of x. But if n > 2m − 1, then you need to use the least square method to calculate the solution. Since λ_i can be a negative number, after specifying λ₁,λ₂,...,λ_m, we change the basis function to the real function. In other words, if this value is negative, consider λ_i and λ_n+1 as the complex conjugate roots of the equation, and we can get: λi=eai(cosβi+isinβi),λi+1=eai(cosβi−isinβi) {\lambda _i} = {e^{ai}}(\cos {\beta _i} + i\sin {\beta _i}),{\lambda _{i + 1}} = {e^{ai}}(\cos {\beta _i} - i\sin {\beta _i})

The corresponding basis function is: φi(x)=eaixcosβix,φi+1(x)=eaixsinβix {\varphi _i}(x) = {e^{{a_i}x}}\cos {\beta _i}x,{\varphi _{i + 1}}(x) = {e^{{a_i}x}}\sin {\beta _i}x

Assuming that there are k real roots and t complex conjugate roots in the characteristic equation of the fitting curve, then the new fitting function will be: λi=eai(cosβi±isinβi)(i=1,2,...,t) {\lambda _i} = {e^{{a_i}}}\left({\cos {\beta _i} \pm i\sin {\beta _i}} \right)\left({i = 1,2,...,t} \right) g(x)=c1φ1(x)+c2φ2(x)+...+cmφm(x)=c1λ1x+c2λ2x+...+ckλkx+ck+1ea1xcosβ1x+ck+2ea1xsinβ1x+...+cm−1ea1xcosβtx+cmeatxsinβtx \matrix{{g\left(x \right) = {c_1}{\varphi _1}\left(x \right) + {c_2}{\varphi _2}\left(x \right) +... + {c_m}{\varphi _m}\left(x \right)} \hfill \cr {= {c_1}\lambda _1^x + {c_2}\lambda _2^x +... + {c_k}\lambda _k^x + {c_{k + 1}}{e^{{a_1}x}}\cos {\beta _1}x + {c_{k + 2}}{e^{{a_1}x}}\sin {\beta _1}x} \hfill \cr {+... + {c_{m - 1}}{e^{{a_1}x}}\cos {\beta _t}x + {c_m}{e^{{a_t}x}}\sin {\beta _t}x} \hfill \cr}

The in the formula needs to be calculated with the least square integration method proposed above.

According to the above steps, we can not only master more influencing factors, but also put forward better solutions based on the numerical changes obtained from calculation and analysis.

The idea of this paper follows the path of problem discovery, subject research and practical reflection, in which both qualitative and quantitative analyses are applied. Figure 1 shows the workflow of this main methodology. One of the main concerns of qualitative researchers is their emphasis on the description of the urbanisation process, and the research analysis model. In this study, literature research can give the key influencing factors of the urban warming effect. Quantitative research is guided by a positivist methodology to explore the different indexes of urbanisation and the heat island effect. At the same time, the model that analyses the impact of urban expansion on the urban environment can be determined from the literature. From a quantitative point of view, the urban expansion process and basic characteristics of climate change in Beijing can be analysed for the past 20 years through the establishment of multi-period LUCC and meteorological data. Then the remote sensing data is used to analyse the spatial distribution pattern of urban heat islands and the underlying surface, mainly focusing on the correlation analysis of LST, NDVI and NDBI. Through this paper, we are hopeful to focus on analysing the relationship between urban building density distribution and the urban heat island effect. This paper takes Beijing as the research object; using high-resolution remote sensing images and through artificial visual interpretation in combination with the surface temperature data obtained from remote sensing inversion, correlation of the internal relationship between the urban building density distribution and the urban heat island effect is obtained. This paper discusses the impact of urban building density distribution on the urban heat island effect. It can provide a scientific basis for the planning and layout of urban construction and urban thermal environmental management.

Fig. 1

The measurement system in this paper uses WGS 1984 and the projection uses equal area projection. The landsat5 and landsat8 data are established through the geospatial data cloud platform, including five dates and specific satellites. The cloud coverage is <10%.

Through ENVI, LST and LUCC are established by data processing, sample establishment, classifier selection and post-processing.

Urban heat island remains one of the biggest attention-drawing environmental problems in China (Zheng & Li, 2016). The data of Beijing districts and counties can be downloaded from the National Basic Information Centre. The datasets used for the study of urban heat island, as an impact of urban expansion, can be obtained through periodic data collection at three different times to get a comparative basis for the changes that occur over a given period. Landsat TM images are then chosen and corrected radiantly to retrieve the land cover types followed by visual interpretation to differentiate false-colour composites of TM/ETM.

Beijing city covers approximately 16410 km², longitudinally and latitudinally lying at about 2”05′ and 1”37′, respectively (Zhi et al., 2014). The city is graced by a sub-humid warm temperate continental monsoon climate that gives rise to four distinct seasons, humid, hot summer, cold and windy winter. Beijing has continued to experience rapid urbanisation, which has consequently led to increased environmental problems like UHI, pollution haze and dust and sandstorms in the urban areas (UA) (Siddiue et al., 2020). The study area is mapped into specific distinctive portions as UA and non-urban areas (NUA). The calculation of urban area land can be performed using various ways depending on the precision and accuracy of the method, using satellite images and change detection. In this study special nets with specific dimensional measurements that match the Moderate-resolution Imaging Spectroradiometer are created for the calculation of the land area percentage.

Fig. 1

NDVI is a common parameter in exploring the influencing factors of ground temperature. It can be used to detect vegetation coverage, growth status and seasonal changes. The formula is as follows, NDVI=(RNIR−R)/(NIR+RNIR) NDVI = (RNIR - R)/(NIR + RNIR)

In the formula, NIR is the near-infrared band reflectivity, R is the red-light band reflectivity, which is represented by the 5th and 4th bands respectively in Landsat 8 TIRS. The value of NDVI ranges from −1 to 1. NDVI>0 means all vegetation. The higher the NDVI, the higher the vegetation coverage, which means good growth and the rich types of vegetation in the area.

LST is a good indicator of the energy balance of the earth’s surface and the greenhouse effect. It is a key factor in the physical process of ground strikes on a regional and global scale. It can reflect the energy flow and material exchange of the soil-vegetation-atmosphere system. It is very necessary for many fields such as climate, hydrology, ecology and biogeochemistry. The ground temperature retrieved from satellite data can reflect in more detail the average temperature of the underlying surface of each pixel and can reflect the spatial distribution characteristics of the underlying surface temperature field. Although the error of using remote sensing data to retrieve the surface temperature is inevitable, it is still the most effective and easiest way to obtain the surface temperature of a large area at present. This is also an advantage of LANDSAT images to reflect the difference of the surface temperature distribution in a large area.

Landsat bands can be used in land detection with the wavelength of 433 nm to 2300 nm and the spatial resolution is 30 m. LST product images, with particular desired properties and resolution levels, are used for analytical purposes. For instance, the Landsat product, which is capable of averaging the LST after every 16 days can be used to obtain the daily LST produced by the split-window algorithm, which is the method for LST retrieval from satellite data. On the other hand, the Landsat Vegetation Index product can be used to obtain data for the EVI data by determining the vegetation index for a given period at a specific spatial resolution. The mean EVI, which is commonly used to determine the vegetation abundance, is calculated for the whole period of the study i.e. 1999 to 2019. The NDBI data is generated using a specific Landsat Surface Reflectance product at a particular spatial resolution. This determination provides the estimated surface spectral reflectance of the study bands, which is then corrected for atmospheric conditions. The NDBI is defined by the equation: NDBI=(RSWIR−RNIR)/(RSWIR+RNIR), NDBI = (RSWIR - RNIR)/(RSWIR + RNIR),

NDBI can be used to characterise the density of urban land?where RSWIR is the surface spectral reflectance of the second last band and RNIR is the surface spectral reflectance of the second band used in the study process.

Using the difference in averaged LST value between the urban and rural areas of study, followed by a determination of the urban area’s pixel levels as the difference in the averaged nonurban LST and each pixel’s LST in the urban area, the SUHII levels of the study areas can be determined by calculation. Investigation of the data obtained from temporal trends of variables like SUHII, EVI, LST and NDBI from 1999 to 2019 is performed at every pixel by the use of linear regression models and ordinary least squares (Minghong et al., 2011). Using the LST as the dependent variable and the study year as the independent variable, multiple ordinary least squares models are run at each pixel to examine the regression coefficient that is used in measuring the annual change in the LSTs magnitude and trend. The significance is then tested at a specific level e.g. at 0.05. Thus, any positive trend with p <0.05 indicates a significantly increased LST at that specific pixel from 1999 to 2019.

Spatiotemporal variation of UHI using LST data using remote sensing images formed at different periods of the study is the major focus of this method. Quantification of the UHI is then performed by introducing URI based on bright temperature normalisation (Siddiue et al., 2020). However, the surrounding conditions of the rural areas affect the magnitude of a UHI and the fact that UHI focuses on the LST’s relative spatial intensity. Normalisation of LSTs is done to allow comparison of their spatial distribution to provide a spatiotemporal pattern variation of the HUI because the remote sensing images only change LST values at different study periods as opposed to changing the spatial LST distribution. The normalised LSTs are then classified into five different thermodynamic levels by density segmentation method to characterise their distribution levels and the area of each level calculated. URI is then introduced for the quantification of the rate of distribution of the urban land to UHI using the following equation: Ts=K1/ln(K2/B(Ts)+1) Ts = K1/\ln (K2/B(Ts) + 1)

Where: T_S is the surface temperature; B (T_S) is the blackbody at the same temperature. K1 and K2 can be obtained from the image header file.

Based on Landsat TM data The URI often reflects the degree of UHI development in built-up land; thus, the more extensive the URI, the more severe the UHI effects and vice versa (Jia Wang et al., 2019). Computational analysis of the UHI area to urban land ratio while considering each temperature level’s weighted values provides the URI value.