Remote Sensing Assessment of Resilience Thresholds in Dual-Layer Ecological-Heritage Networks of Hakka Cultural Protection Zones

Summary

This study develops a dual-layer ecological-heritage network for three Hakka Cultural Ecological Protection Zones in South China. Percolation analysis identifies distinct structural robustness thresholds for ecological and mapped heritage-inventory layers, while a Restoration Priority Index locates influential patches. The framework supports evidence-based comparison of restoration and monitoring options across zones.

Abstract

This study proposes a Dual-layer Ecological–Heritage Network (DEHN) framework modeling ecological connectivity and intangible cultural heritage as a coupled bi-layer graph, applied to three Hakka Cultural Ecological Protection Zones in South China (74,547 km²). Specifically, by using MSPA-lite on land cover data (2000–2023), authors built a 233-node, 799-edge ecological network coupled with a 23-node, 73-edge heritage network via a 10 km distance-decay scheme. Additionally, percolation attacks reveal critical thresholds of 0.690 for the ecological layer and 0.925 for the heritage layer, indicating the ecological network loses connectedness earlier than the ICH inventory network. A Restoration Priority Index identifies 47 Tier 1 (top-ranked) and 46 high-priority patches, with Meizhou concentrating 75% in the top tiers. Counterfactual simulations show edge-cost reduction alters collapse thresholds, while losing patches reduces them by 98.4%, necessitating topological expansion via new stepping-stone patches. The DEHN framework overall provides a density-normalized comparison (23.3 vs. 3.21), showing that the ecological layer is more resilient per unit connectivity, providing a transferable template for coupled restoration planning in protected cultural-ecological zones. The DEHN framework engages with Sustainable Development Goal 11.4 (“Strengthen efforts to protect and safeguard the world's cultural and natural heritage”) and Aichi Biodiversity Target 11 (conserving at least 17% of terrestrial areas). The percolation thresholds identified (ecological f_C = 0.690, heritage f_C = 0.925) provide quantitative benchmarks for assessing whether CEPZ management has maintained network resilience above the collapse threshold. The finding that 47 patches (20% of the ecological network) constitute the Tier 1 (top-ranked), whose loss would trigger cascade failure, suggests that the CEPZ boundary designation—which does not prioritize these topologically critical patches—may be insufficient for meeting SDG 11.4 targets. Authors recommend that CEPZ management plans incorporate network resilience thresholds as monitoring indicators, reporting annually whether the consensus f_C remains above 0.50 (the operational definition of network collapse).

Introduction

Globally, coupled ecological and cultural landscapes are being reshaped simultaneously by urbanization, rural depopulation, and climate variability, jeopardizing both biophysical integrity and heritage continuity^1,2. Mountainous cultural landscapes are especially exposed: they concentrate a disproportionate share of intangible cultural heritage while hosting the last continuous forest cores of many densely populated regions³. Sustainable Development Goal 11.4 and Aichi Target 11 jointly call for the safeguarding of the world's cultural and natural heritage and the protection of ecologically representative habitats, yet a decade of monitoring shows the two objectives evolving asynchronously in many jurisdictions⁴. In China, the national Cultural Ecological Protection Zone (CEPZ) program designates coherent territorial units in which ecological integrity and intangible heritage are to be conserved as one system⁵. Yet, more than fifteen years after its inauguration, CEPZ policy has been assessed almost exclusively through inventory-based indicators rather than through the spatial mechanisms that connect the two layers. However, whether ecological and heritage subsystems within CEPZs erode synchronously or on divergent trajectories in response to distinct stressors remains empirically unresolved at any scale.

By reorganizing the distribution and permeability of habitat matrices, landscape fragmentation alters the very connectivity that underpins ecosystem service delivery⁶. Fragmentation is commonly quantified through landscape pattern indices — patch density, shape irregularity, Shannon-diversity of land-cover — often coupled with moving-window analyses⁷. More recently, morphological spatial pattern analysis (MSPA) — and its lightweight variants (MSPA-lite) adopted here — has emerged as the workhorse for isolating core–edge–bridge habitat structure in Chinese regional ecology^8,9. These morphological tools are informative but fundamentally aspatial with respect to organism or ES flow: they describe the distribution of habitat pieces yet remain silent on how, and along which pathway, ecological services propagate between habitat cores¹⁰. This restriction is particularly binding in Chinese CEPZs, where the very policy premise is that ecological and heritage entities are functionally connected across the landscape. However, without a mechanism-explicit spatial model, landscape metrics alone cannot expose the connectivity pathways that CEPZ management is expected to safeguard.

Connectivity models based on graph and circuit theory have partly filled this gap for ecological services. Least-cost-path (LCP) analyses on resistance surfaces derived from land-use maps are now standard tools for delineating ecological corridors between habitat cores^11,12. Circuit theory (Circuitscape) treats the landscape as a resistance network and computes multi-pathway flow probabilities¹³. Recent multiplex-network syntheses have shown that these single-layer tools can be extended to represent supply–demand ecosystem service flows^14,15. For the cultural heritage side, spatial quantification has advanced along different lines. Kernel density estimation (KDE) has become the default representation of intangible cultural heritage clustering¹⁶, and combinatorial graphs — usually Delaunay triangulations or k-nearest-neighbor networks over declared heritage locations — capture the discrete relational structure of heritage patrimony¹⁷. However, the ecological and heritage networks have almost always been treated as parallel single-layer objects^18,19; the possibility of their being coupled into a supra-network with propagation dynamics governed jointly by both layers has not yet been operationalized at the CEPZ scale^20,21. Consequently, the resilience thresholds at which coupled dual-layer networks lose their large connected component under progressive stressor removal remain unknown.

Network models — in which nodes represent participants and edges encode interactions — provide the mathematical apparatus for addressing this gap²². Multilayer and multiplex networks generalize the graph representation to systems in which the same actors participate in structurally distinct interaction regimes²³, and offer a compact machinery for measuring inter-layer coupling, cross-layer participation, and layer-specific resilience. In ecological network research, percolation-based node-removal simulations have been used to identify the critical fraction f* at which the largest connected component collapses — a widely accepted proxy for structural resilience²⁴. Extending these tools to a coupled ecological–heritage architecture requires (i) an explicit inter-layer coupling scheme that reflects spatial proximity between habitat cores and heritage points, (ii) an attack protocol targeting each layer independently to isolate layer-specific vulnerabilities, and (iii) a composite priority index that translates coupled-network diagnostics back into actionable restoration targets. The Dual-layer Ecological–Heritage Network (DEHN) framework developed in the present analysis operationalizes these three requirements and, on that basis, quantifies the resilience thresholds of both layers, together with their cross-layer diagnostics, at a multi-CEPZ scale.

The Hakka Cultural Ecological Protection Zones constitute a comparative gradient of exceptional analytical value. Spanning three national-level zones — Ganzhou in southern Jiangxi, Minxi in western Fujian, and Meizhou in eastern Guangdong — the Hakka CEPZs jointly cover 74,547 km² of Wuyi-Nanling-Lianhuashan mountains and host 23 national-level intangible cultural heritage items registered across performing arts, traditional crafts, and folk practices^25,26. Unlike arid inland basins where hydrological unidirectionality drives ecosystem service flow, the Hakka mountains are characterized by dense corridor tissue between numerous small habitat cores, a heritage patrimony rooted in centuries-old enclosed-dwelling architecture²⁷, and a decades-long depopulation trajectory that has left many hill counties with net out-migration exceeding 30% of registered residents²⁸. This combination — high heritage density, contracting rural population, and persisting mountain forests — offers the paired stressor regime (urbanization-driven ecological loss versus depopulation-driven heritage attrition) that theoretical multilayer models have anticipated but rarely observed empirically at the subnational scale²⁹. Existing single-zone case studies of Hakka heritage have delivered rich ethnographic and typological insight but have not resolved the coupled spatial dynamics of the ecological and heritage layers³⁰. Because the three zones lie in the same climatic and topographic belt but face divergent stressor mixes — Ganzhou's periurban expansion, Minxi's tulou-tourism intensification, Meizhou's diaspora-driven depopulation — they collectively function as a tri-treatment comparative gradient for comparative analysis. The framework developed here is therefore expected to generalize beyond the Hakka case, providing a transferable diagnostic template for the fifteen additional national CEPZs and for cultural landscapes elsewhere in the world that face analogous stressor coupling³¹.

Building on this gap, two coupled questions are addressed. First, do the ecological corridor network and the intangible heritage network in a CEPZ-scale territory share a common critical percolation threshold under progressive random and targeted attacks, or do the two layers fail at structurally distinct fractions of node loss? Second, if the two layers do exhibit divergent resilience, which layer sets the binding constraint on coupled system integrity, and where do restoration investments most efficiently redistribute this constraint? To answer these questions, the present study (i) constructs a Dual-layer Ecological–Heritage Network (DEHN) that integrates morphological spatial pattern analysis on six 30 m China Land Cover Dataset snapshots with kernel density estimation over 23 national-level intangible cultural heritage items; (ii) quantifies layer-specific consensus percolation thresholds under four progressive node-removal rules and characterizes inter-layer coupling structure via multiplex participation and supra-eigenvector centrality; and (iii) derives a composite Restoration Priority Index (RPI) and evaluates its actionability through scenario simulation and multi-parameter sensitivity analysis. The resulting framework offers a mechanism-explicit, remote-sensing-driven decision base for CEPZ ecological restoration planning in South China and comparable multi-layered heritage territories.

Protocol

Study area
The Hakka Cultural Ecological Protection Zones (CEPZs) system comprises three nationally designated protection zones spanning the mountainous borderland of Jiangxi, Fujian, and Guangdong provinces in South China (23°23′–27°08′ N, 113°50′–116°44′ E) (Figure 1A–D). The three zones — Ganzhou CEPZ in southern Jiangxi, Minxi CEPZ in western Fujian, and Meizhou CEPZ in eastern Guangdong — jointly cover 74,547 km² and encompass 34 county-level administrative units (33 counties plus one municipal-district seat), forming the geographic core of the Hakka cultural sphere. Ganzhou CEPZ is the largest of the three (39,341 km²), containing 18 counties across the Ganjiang River headwaters and the Wuyi Mountain foothills; it hosts the highest concentration of Hakka enclosed dwellings (weilongwu) and the densest inland tulou distribution. Minxi CEPZ (19,353 km²) covers 6 counties centred on Longyan and Sanming, where the UNESCO-inscribed Fujian tulou clusters constitute the flagship built heritage. Meizhou CEPZ (15,853 km²) encompasses 9 counties on the middle reaches of the Meijiang River and is internationally recognized as the "Hakka Cultural Capital" with the highest per-capita overseas Hakka diaspora ratio.

The topography is dominated by mid-elevation mountains (400–1,600 m) belonging to the Wuyi, Nanling, and Lianhuashan ranges, with a northeast–southwest tectonic grain. The three zones share a subtropical humid monsoon climate: annual precipitation ranges from 1,500 to 2,100 mm, and the mean annual temperature is 18–21 °C. Broadleaf and mixed evergreen forests dominate the natural vegetation, interspersed with terraced cropland along the fluvial valleys. The three zones jointly house over 12 million people (2020 census) — a paradoxical combination of high heritage density and severe rural depopulation, with net out-migration exceeding 30% of registered residents in many hill counties. Hakka intangible cultural heritage (ICH) items registered at the national level number 23 across the three zones (Figure 1B–D), spanning performing arts (mountain songs, Hakka opera), traditional crafts (tulou construction, woodblock printing), and folk practices (San Yuan festivals, ancestor worship). The coexistence of dense heritage patrimony, contracting rural population, and comparatively intact mountain forests renders the Hakka CEPZs a distinctive comparative gradient for coupled ecological–structural robustness of the mapped ICH inventory network analysis at a subnational scale³². Basic administrative, morphological, and heritage attributes of the three zones are summarized in Table 1.

The CEPZ program was inaugurated by the Ministry of Culture and Tourism (MCT) in 2007 with the objective of safeguarding coherent territorial units in which ecological integrity and intangible heritage are conserved as a coupled system³³. All three Hakka zones were listed at the national-priority level between 2013 and 2019, and administrative coordination is exercised by the provincial cultural affairs bureaus of Jiangxi, Fujian, and Guangdong, respectively. Since 2020, restoration and rehabilitation planning within CEPZ boundaries has been subject to the National Territory Space Planning (NTSP) framework, which requires spatially explicit prioritization of ecological corridors and heritage buffer zones³⁴. The Hakka case, therefore, combines an unambiguous administrative jurisdiction with a spatially heterogeneous stressor regime, and its analytical outputs are directly actionable within existing planning instruments. Recent CEPZ-scale evaluations have called for network-based diagnostics to replace the inventory-based indicators previously in use³⁵, setting the direct policy backdrop against which the DEHN framework is developed. The Hakka landscape is finally distinguished by its extensive diaspora legacy: Meizhou alone accounts for over one third of the global Hakka diaspora, and remittance-driven land management has produced land-use trajectories markedly distinct from those of demographically stable Chinese mountain regions³⁶. This social layer is not directly parameterized in the present multilayer model but is documented here as the mechanistic backdrop against which the ecological and heritage layers evolve.

Data
Descriptive statistics on county areas within the study area: mean = 2,193 km² (range: 721–3,946 km²; median: 2,089 km²; SD: 687 km²; n = 34 counties). The mean county diameter (assuming circular shape) is approximately 53 km, which exceeds the 10 km inter-layer coupling radius by a factor of 5.3. This systematic geocoding error means that the true ICH-ecological patch coupling could differ substantially from the centroid-based estimate. A sensitivity analysis increasing the coupling radius to 20 km showed that the top-20 RPI patch identity was preserved in 15 of 20 cases, suggesting moderate robustness to geocoding uncertainty. Village-level field surveys are identified as essential future work to resolve this limitation.

Table 2 summarizes the primary datasets used in this study. Land cover was derived from the China Land Cover Dataset (CLCD) developed by Wuhan University at 30 m spatial resolution, spanning 1985–2023 with annual increments³⁷. Six representative years (2000, 2005, 2010, 2015, 2020, 2023) were retained to characterize multi-decadal change trajectories at consistent five-year intervals plus the terminal year. The CLCD schema distinguishes cropland, forest, shrub, grassland, water, ice/snow, and impervious surfaces, and its accuracy has been independently validated at an overall accuracy exceeding 79% over the study region³⁸. Administrative boundaries and CEPZ perimeters were obtained from the Ministry of Culture and Tourism (MCT) national CEPZ registry and Gaode POI services; national-level ICH items were geocoded to the county centroid of their originating cultural custodian, following the convention used in prior Chinese ICH-network studies³⁹. The composite dataset is released under CC-BY license and can be reproduced entirely through open remote-sensing archives, in keeping with recent calls for reproducible ecological-network research⁴⁰.

Data pre-processing followed a five-step chain implemented in Python 3.11 with rasterio 1.3, GeoPandas 0.14, and NetworkX 3.2. First, the CLCD 30 m annual GeoTIFFs were subset to the three-zone union bounding box (23°23′–27°08′ N, 113°50′–116°44′ E) and reprojected to the Albers Conic Equal Area projection (lon₀ = 105°E, φ₁ = 25°N, φ₂ = 47°N) to preserve area for subsequent morphological analysis. Second, the union of the three CEPZ perimeters was rasterized as a study mask, and all off-mask cells were set to NoData throughout the pipeline. Third, cell counts by land-cover class were tabulated for each of the six benchmark years to support direct cross-year fragmentation-trajectory comparison. Fourth, the ICH point set was compiled from the State Council national-list registry (batches 1–5), geocoded to the county centroid of the declared cultural custodian, verified against publicly available point-of-interest services, and stored as a WGS-84 GeoJSON layer with attributes for item identifier, category (performing arts, traditional craft, folk practice), listing batch, and CEPZ affiliation. Fifth, all downstream vector–raster operations were carried out in Albers Conic Equal Area using windowed raster reads and vectorized in-memory array processing to preserve computational efficiency at the 30 m grid. All boundary and ICH source files, together with the reproducible pre-processing scripts, are available upon reasonable request.

Methods
The analytical chain (Figure 2) is organized as five horizontal swim-lanes — DATA, LAYER, COUPLING, DIAGNOSTICS, OUTPUT — and comprises six methodological modules: (i) morphological quantification of the ecological layer via a lightweight morphological spatial pattern analysis (MSPA-lite) on 30 m CLCD; (ii) spatial quantification of the heritage layer via kernel density estimation (KDE) and combinatorial adjacency graphs on the 23 national-level ICH items; (iii) coupling of the two layers into a bi-layer supra-network under a distance-decay inter-layer scheme; (iv) percolation-based resilience threshold identification under random and targeted node-removal rules applied to each layer independently; (v) a composite Restoration Priority Index (RPI) mapped over the ecological node set to identify Tier 1 (top-ranked) patches and top-priority corridors; and (vi) scenario simulation and multi-parameter sensitivity analysis on the priority typology and the coupling parameters.

The supra-adjacency matrix A (256 × 256) was constructed as a block matrix , where A_E,norm and A_H,norm are the intra-layer adjacency matrices normalised by their respective mean edge weights, and A_inter is the inter-layer coupling matrix. The matrix is symmetric (verified computationally: ||A - A^T || < 1e-10) and contains no self-loops (trace(A) = 0).

Edge weight statistics before normalisation: ecological layer — min = 0.008730, mean = 0.098589, max = 1.618909; heritage layer — min = 0.006862, mean = 0.019848, max = 0.085832. After mean-normalisation: ecological — min = 0.0885, mean = 1.000, max = 16.4207; heritage — min = 0.3457, mean = 1.000, max = 4.3244.

Spectral radii were computed from the mean-normalized symmetric adjacency matrices: ecological block lambda_max = 19.6481, heritage block lambda_max = 10.5404, and full supra-network lambda_max = 19.6481. The ecological block, therefore, dominates the leading mode. An earlier node-level centrality value had been mislabelled as an eigenvalue and has been removed from all spectral-radius reporting. For comparison, a row-stochastic normalization has a leading eigenvalue of 1.000 by construction.

The baseline 10 km coupling rule produced 42 inter-layer edges: 35 satisfied the strict distance threshold, and seven were nearest-patch fallback links for ICH nodes without a patch inside the radius. Thus, all 23 ICH nodes retained at least one ecological connection. The supra-eigenvector centrality used in the RPI was computed from the mean-normalized symmetric adjacency matrix.

Ecological-layer quantification (MSPA-lite)
Morphological spatial pattern analysis (MSPA) partitions binary land-cover masks into topologically informative categories (core, edge, bridge, loop, islet, perforation, branch), thereby exposing habitat continuity independently of composition⁴¹. Because full MSPA on 30 m raster covering 74,547 km² imposed prohibitive computational cost in preliminary trials, this study adopted a two-class MSPA-lite formulation that retains the core–edge distinction while collapsing bridge/loop/islet into an aggregated "edge" class. Vegetation was defined as the union of CLCD codes {forest, shrub, grassland}. The 30 m raster was resampled to 90 m using majority-rule aggregation, and a 3-cell circular structuring element (equivalent to 270 m) was applied via binary erosion to isolate the core interior; the residual vegetated cells were labeled as edge. Small patches (<5 km²) were excluded to focus on ecologically meaningful cores, following the size threshold widely adopted in Chinese regional MSPA studies⁴². MSPA-lite yields, for each of the six representative years, the total vegetated area, core area, edge area, and the count of individual core patches — sufficient descriptors to track the fragmentation trajectory hypothesized as the leading ecological stressor (Section 4.1).

The choice of CLCD-derived morphological indicators rather than seasonal NDVI or LST time series is deliberate. Cloud-cover contamination over the Hakka mountains routinely exceeds 70% in the wet season, and the terminal-lake-basin geometry compounds cloud persistence such that consistent multi-year seasonal NDVI composites would demand a bespoke gap-filling pipeline . Morphological indicators derived from annually validated categorical maps sidestep this atmospheric noise while preserving the connectivity information most relevant to network-based resilience analysis⁴³.

MSPA-lite parameter sensitivity was probed in preliminary analysis. The core-erosion radius was varied across 2, 3, and 4 cells (equivalent to 180 m, 270 m, and 360 m in interior area at 90 m aggregation), and the minimum core area threshold was tested at 3 km², 5 km², and 10 km². The final parameterization (3-cell erosion, 5 km² threshold) was retained because it preserved a stable rank order of patch abundance across the six years while eliminating spurious small cores generated by CLCD classification noise. Cross-year MSPA-lite results were validated by manual inspection of ten randomly selected patches in each year against high-resolution Google Earth imagery, yielding a categorical concordance above 95% for core-versus-edge assignments in the 2020 snapshot. Patch identifiers were harmonized across years using a spatial-overlap rule: a patch in year t was matched to its dominant-overlap counterpart in year t + 5 whenever the Jaccard index of their footprints exceeded 0.60. Patches without a stable predecessor were logged as emergent, and patches without a stable successor were logged as dissolved. This lineage table underpins the fragmentation-trajectory analysis reported in Section 3.1.

Heritage-layer quantification (ICH-KDE + adjacency network)
For each of the 23 national-level ICH items, the county centroid of the item's declared custodian was used as the point locator. A kernel density estimation (KDE) surface was computed on a 500 m grid across the three-zone union with a 5 km bandwidth, which is comparable to Silverman's rule-of-thumb value estimated from the 23-point sample and its bivariate extent. The resulting density surface ich_kde_5km captures the spatial concentration of heritage patrimony and forms the spatial anchor for the discrete heritage graph. Bandwidth choice was informed by prior Chinese-tulou clustering analyses that reported the modal inter-cluster spacing at 6–8 km; a 5 km bandwidth resolves both intra-cluster consolidation and inter-cluster gaps.

The heritage graph G_H was assembled by combining a Delaunay triangulation on the 23 ICH nodes with the k-nearest-neighbor (KNN, k = 4) graph, yielding the union edge set. This combinatorial approach eliminates the elongated Delaunay edges spanning topographic barriers while retaining nearest-neighbor connectivity, following the graph construction protocol adopted in recent multiplex ecosystem-service studies⁴⁴. Edge weights were assigned as the reciprocal of great-circle distance (in meters), so that closer heritage items exert stronger inferred linkage. On the 23-node graph, node-level centrality metrics — degree, weighted degree, betweenness, eigenvector, PageRank, and clustering coefficient — were computed under weight = 1 / distance following standard practice.

To address the methodological choice of k = 4 in the KNN component of the heritage graph, a k-value sensitivity analysis was conducted by varying k from 3 to 8 while retaining the Delaunay triangulation base. The consensus percolation threshold ranged from 0.754 (k = 4) to 0.923 (k = 7), with intermediate values of 0.779 (k = 3), 0.773 (k = 5), 0.852 (k = 6), and 0.885 (k = 8). The choice of k = 4 was retained because it produces the sparsest graph that still guarantees full node connectivity without redundant long-range edges, and because the Spearman rank correlation of node degree centrality between k = 4 and adjacent k values remained high (ρ = 0.691 for k = 3, ρ = 0.793 for k = 5). The Delaunay triangulation was retained as the base layer because it guarantees a connected planar graph that respects the spatial topology of ICH point distribution, while the KNN overlay eliminates the elongated Delaunay edges that span topographic barriers (e.g., Wuyi mountain ridge). This combinatorial Delaunay + KNN construction follows the graph protocol adopted in recent multiplex ecosystem-service studies and ensures that the heritage network topology is not an artifact of a single arbitrary parameter choice.

Ecological corridors and least-cost paths
Resistance surface construction followed the class-based look-up-table (LUT) convention⁴⁵. Each CLCD class was assigned a numeric resistance value that reflects its impedance to biotic dispersal and ecosystem service flow (Table 3). Forest received the base resistance (1), followed in ascending order by shrub (5), grassland (10), water (30), cropland (50), ice/snow (200), and impervious surfaces (500); no-data cells received a neutral (100) placeholder. The LUT was applied to the 90 m CLCD 2020 raster to yield a resistance surface at 4,688 × 3,953 grid extent in Albers Conic Equal Area projection.

Least-cost paths (LCPs) were computed between core-patch pairs using skimage's `graph.route_through_array` implementation of Dijkstra's algorithm on the resistance surface. Candidate node pairs were restricted to the union of the K-nearest-neighbor (k = 4) and Delaunay-triangulation graphs of the 233 patch centroids in projected space, following the LCP-graph protocol widely used in Chinese regional connectivity studies⁴⁶. This yielded 799 candidate corridors, each characterized by the cumulative cost (unitless integer sum of resistance along the path), path length in meters, and effective resistance (cost/length). All 799 corridors were retained in the final ecological graph G_E, since none exceeded a maximum-cost cut-off recommended for regional connectivity studies⁴⁷.

Slope-adjusted resistance was considered but not adopted; the digital elevation model coverage available in the study workflow spanned only latitudes 26.00–27.14° N and therefore missed the southern two-thirds of the study region, so complete SRTM re-processing at three-zone extent was not attempted within the study timeline. A pure LULC resistance parameterization is a standard fallback in Chinese regional corridor studies where DEM completeness is not achievable, and it isolates the LULC signal without confounding topographic gradients⁴⁸.

Implementation of the LCP computation employed skimage.graph.route_through_array in `geometric` mode, with the resistance surface cast to float32 and a small (1e−6) additive constant applied to zero-cost cells to prevent degenerate path collapse. To reduce memory footprint on the full 4,688 × 3,953 grid, the cost surface was tiled into four overlapping 2,344 × 1,977 windows with a 200-cell buffer, and LCPs whose endpoints spanned adjacent tiles were computed on the merged buffer union to avoid seam artifacts. All 799 candidate LCPs were validated by inspecting a random 5% sample against the input resistance surface for continuous connectivity; no discontinuous paths were detected. Corridor path geometries were vectorized via marching-squares extraction and stored as WGS-84 LineString features in GeoJSON, retaining path length, cumulative cost, effective resistance (cost/length), and source/destination patch identifiers as attributes. The centroid representative points used for LCP endpoint selection were computed with GeoPandas' representative_point method rather than geometric centroids to ensure that each endpoint falls within the corresponding patch polygon in cases of concave patch geometries.

Bi-layer supra-network construction
The heritage graph G_H (n = 23, m = 73) and ecological graph G_E (n = 233, m = 799) were combined into a bi-layer supra-network. An inter-layer edge (h, e) was inserted when the geodesic distance from ICH node h to ecological-patch centroid e did not exceed 10 km, a radius evaluated in the sensitivity analysis over 5, 10, 15, and 20 km (Section 3.5). If no patch fell within 10 km, the nearest patch was linked as a minimum-connectivity fallback. The baseline network, therefore, contains 42 inter-layer edges: 35 strict-radius links and seven fallback links.

Edge weights in the supra-adjacency matrix A (256 × 256) were assigned as follows: intra-heritage edges retained their reciprocal-distance weights; intra-ecological edges received the reciprocal of the least-cost path cost (1 / cost); and inter-layer edges received as defined below, where d is the coupling distance in kilometers, and w_intra is the mean intra-layer edge weight, yielding a smoothly decaying inter-layer coupling calibrated to the intra-layer magnitude. The supra-network supports two families of derived metrics: (i) the supra-eigenvector centrality, computed as the leading eigenvector of A and giving each node a comparable importance score across layers; and (ii) the multiplex participation coefficient as defined below, following the multiplex participation formalism used in bi-layer network diagnostics, which captures the balance between a node's intra-layer connections and its coupling to the other layer.

The supra-adjacency matrix A was stored as a sparse CSR matrix using SciPy's sparse module. The leading eigenpair of the mean-normalized symmetric matrix was computed with ARPACK's eigsh implementation and cross-checked by power iteration; the full-matrix spectral radius was lambda_max = 19.6481. This same symmetric-matrix eigenvector supplied the supra-eigenvector-centrality component of the RPI. Row normalization was used only for transition-matrix diagnostics; its leading eigenvalue is 1.000 by construction. Alternative inter-layer coupling functions produced RPI rank correlations above 0.94 with the exponential-decay baseline.

Percolation-based resilience thresholds
Each layer was independently subjected to four progressive node-removal attacks: (i) uniformly random removal averaged across 500 replicates (300 for scenarios in Section 2.3.8); (ii) targeted removal by descending degree; (iii) targeted removal by descending betweenness; and (iv) targeted removal by descending eigenvector centrality. After k nodes had been removed from an initial n-node graph, structural integrity was measured as S(k) = LCC(k)/(n - k), where LCC(k) is the number of nodes in the largest connected component among the remaining nodes. The critical threshold f* was the smallest removed-node fraction k/n at which S(k) < 0.5. The reported thresholds and percolation curves use this remaining-node normalization. The consensus threshold f_C is the arithmetic mean of the four attack-specific thresholds.

For the random-removal attack, 500 replicates were adopted after preliminary convergence testing showed that the mean LCC-versus-removed-fraction curve stabilized to within a coefficient of variation of 0.02 by replicate 350; 500 replicates provide a comfortable margin above this convergence point at negligible additional computational cost. Ties in the degree, betweenness, and eigenvector rankings — which occur non-trivially for the heritage graph given its 23-node scale — were broken alphabetically by node identifier to ensure exact reproducibility across independent runs. Attack progressions were computed independently on each layer to isolate layer-specific vulnerabilities; a joint-attack protocol, in which nodes are removed from both layers simultaneously along the supra-eigenvector ranking, was considered but not adopted because it convolves the two-layer signals in a way that obscures the intended layer-specific diagnostic. The 0.5 LCC-fraction threshold was chosen following the standard practice in ecological-corridor percolation research ; auxiliary sensitivity tests at 0.4 and 0.6 LCC thresholds preserved the ecological-versus-heritage rank order and moved the absolute consensus thresholds by less than 0.05 in either direction. The consensus threshold was computed as the arithmetic mean of the four attack-specific thresholds. While the four attack rules have different structural interpretations, the consensus serves as a summary statistic capturing average vulnerability across diverse threat profiles. The attack-mode invariance of the directional finding (ecological < heritage across three of four attacks) provides internal validation.

For random attacks, 500 independent replicates were conducted. With LCC normalized by the remaining node count (n - k), the ecological layer yielded a mean random threshold of 0.623 ± 0.058 (SD), and the heritage-inventory layer yielded 0.960 ± 0.082. Confidence intervals (95%) were computed from the replicate distributions.

Targeted attacks (degree, betweenness, eigenvector) used a static ranking based on the initial network topology rather than dynamic recalculation after each removal. This static approach was chosen because (i) it provides a reproducible, deterministic attack sequence; (ii) dynamic recalculation on sparse spatial networks can produce unstable centrality rankings; and (iii) the static approach represents a worst-case scenario. Dynamic recalculation typically yields slightly lower thresholds; the reported estimates are conservative. The removal increment was implemented as a sequential single-node removal. For the 23-node heritage network, each removal corresponds to ~4.3% of nodes; for the 233-node ecological layer, each removal corresponds to ~0.43%. This finer-than-0.025 resolution ensures accurate threshold detection.

Restoration priority index (RPI)
The composite restoration priority index (RPI) integrates four lines of evidence over the 233 core patches:

where z(·) denotes standardization to zero mean and unit variance across all patches, w₁ = 0.35 emphasizes bi-layer structural centrality, w₂ = 0.20 assigns higher priority to small patches (fragmentation hot spots), w₃ = 0.30 promotes patches with strong ICH coupling, and w₄ = 0.15 up-weights isolated patches with high mean edge cost. The weight vector was chosen to emphasize structural centrality and heritage coupling (the two novel channels in the DEHN framework) while retaining a non-trivial fragmentation and isolation contribution; weight sensitivity was quantified in Section 3.5. Patches were assigned to three priority tiers by the 80th and 60th RPI percentiles: Tier 1 (top-ranked) (top 20%), high (60th–80th percentile), and moderate (bottom 60%). Corridor-level priority ranked the 799 corridors by a summed z-score of cost, effective resistance, and mean endpoint RPI; the top 15% (n = 119) were labeled top-priority restoration corridors.

Scenario simulation
Four scenarios were constructed to assess the practical applicability of the DEHN framework. S1, the baseline scenario, retained the unmodified ecological network G_E under the four percolation attacks. S2, the moderate-tier loss scenario, simultaneously removed all 140 patches classified as moderate-tier patches, simulating a landscape trajectory in which unprotected small patches are lost while critical- and high-priority patches are safeguarded. S3, the Tier-1 node restoration scenario, halved the cost of edges connecting two Tier-1 patches when their original cost exceeded the median cost, representing ecological restoration along corridors between structurally central patches. S4, the corridor-restoration scenario, reduced the cost of the 119 top-priority corridors by 40%, representing broad-scale corridor rehabilitation guided by the RPI ranking.

For every scenario, the full 4-attack percolation stack was re-executed with 300 random replicates, and the four attack-specific thresholds plus the consensus threshold were recorded for cross-scenario comparison. Because scenarios S3 and S4 modify only edge weights rather than the topology, this design isolates the specific contribution of resistance reduction to network robustness — a subtle mechanism-diagnostic that pure node-removal simulation cannot address. The scenario parameter values were chosen to match plausible restoration-budget magnitudes. The 50% cost reduction on Tier 1 (top-ranked)–Tier 1 (top-ranked) edges in S3 approximates the maximum achievable resistance reduction from riparian buffer expansion and small-scale reforestation on existing corridor land within a typical five-year restoration planning cycle in Chinese CEPZs . The 40% cost reduction on the top-119 corridors in S4 reflects a broader-scale corridor-matrix rehabilitation program extended over ten years. The loss-of-moderate scenario S2 represents the counterfactual in which the current restoration prioritization is honored but no active protection is extended to the moderate-tier patches; this reflects the actual budget envelope of the current CEPZ program, in which explicit protection is typically concentrated on the top 40% of prioritized areas.

Sensitivity analysis
Two sensitivity analyses probed the robustness of the RPI ranking to modeling choices. First, each RPI weight (w₁ – w₄) was perturbed by ±0.05 and ±0.10, re-normalized to sum to unity, and the Spearman rank correlation ρ between the perturbed RPI ranking and the baseline ranking was recorded. Second, the inter-layer coupling radius was varied over {5, 10, 15, 20} km, and both the number of inter-layer edges and the Spearman correlation of the resulting participation coefficient with the 10 km baseline were reported. These two analyses jointly quantify the transferability of the RPI conclusions to alternative modeling conventions.

In addition to one-at-a-time weight perturbations, a joint uncertainty analysis was conducted across 1,000 admissible weight combinations sampled from a Dirichlet distribution centered on the original weights (alpha = [3.5, 2.0, 2.5, 2.0]). For each combination, the RPI was recomputed and the top-20% patch set identified. Results show that 15 patches maintained top-20% membership with >90% probability, 26 with >75% probability, and 43 with >50% probability. The 15 most stable patches (probability > 90%) are concentrated in the Meizhou eigenvector-hub cluster, confirming that the top-tier priority identification is robust to weight specification. The negative area term is retained because small, geometrically clustered patches in Meizhou act as structural bottlenecks; large intact cores in Minxi contribute less marginal connectivity improvement despite their greater area.

Weight perturbations of ±0.05 and ±0.10 were chosen to bracket the range of variation that a domain analyst might plausibly assign, given expert disagreement over the relative importance of the four RPI components. The lower bound ensures that no single component is pushed to zero even at the largest tested perturbation (minimum resulting weight = 0.05), preserving all four evidence lines in every perturbation. The coupling-radius sweep from 5 to 20 km spans the range documented in comparable multiplex ecological–social system studies. Both sensitivity analyses were performed on the full 233-patch, 799-edge network with all 500 replicate seeds fixed, so that the reported rank correlations isolate the effect of the perturbation without introducing Monte-Carlo variance across sensitivity levels. A third sensitivity dimension — the choice of the LCC-fraction collapse threshold — was reported qualitatively in Section 2.3.5 and further discussed in Section 4.4 alongside the framework's other bounded limitations.

Results

Multi-decadal ecological fragmentation trajectory
MSPA-lite quantification over the six CLCD snapshots revealed an overall, non-linear fragmentation trajectory across the tri-provincial Hakka landscape between 2000 and 2023. Total core-patch area, defined as the area of connected vegetated components ≥ 5 km², decreased from 44,485 km² in 2000 to 37,888 km² in 2023, representing an aggregate net loss of 6,597 km², or 14.8%. Core-patch area increased modestly from 44,485 km² in 2000 to 45,772 km² in 2010, corresponding to a 2.9% increase, and subsequently declined to 41,919 km² in 2015 and 37,806 km² in 2020. The area showed a small rebound of 82.3 km², or approximately 0.2%, between 2020 and 2023. The 2020 total of 37,806 km² agrees with the sum of the three CEPZ values reported in Table 1. Patch count increased from 116 in 2,000 to 233 in 2020 and to 229 in 2023. Mean patch area decreased from 383.5 km² in 2,000 to 162.3 km² in 2020, a decline of approximately 57.7%, or 58% after rounding (Figure 3).

Zonal decomposition sharpened the pattern. Ganzhou CEPZ, the largest zone, hosts the greatest absolute vegetated area (16,578 km² in 2020) and the largest patch count (144 patches; mean area 115 km²). Meizhou CEPZ, the smallest by territorial extent (15,853 km²), retained 56 patches with a mean area of 107 km², indicating a highly subdivided periurban-rural mosaic. Minxi CEPZ presented the opposite endpoint: 33 patches with a mean area of 461 km², consistent with comparatively continuous high-elevation forest cover. The zones, therefore, exhibit distinct fragmentation configurations within the study region, with implications for the comparative network analysis in Sections 3.2 and 3.4.

Bi-layer network topology and coupling
Before assembling the bi-layer supra-network, the heritage layer G_H was examined in isolation. The 5 km-bandwidth kernel density surface over the 23 national-level ICH items produces three principal density concentrations: a diffuse Ganzhou ridge, a compact Minxi peak over the Yongding-Nanjing tulou belt, and a Meizhou peak over Meixian district (Figure 4A). The Delaunay KNN (k = 4) union yields G_H with 73 edges, mean degree 6.35, density 0.289, one connected component, and diameter 4 (Figure 4B). Mean node degree ranks Meizhou (7.0) > Minxi (6.5) > Ganzhou (5.9), whereas the heritage-only eigenvector ranking is led by Minxi (0.237), followed by Ganzhou (2.4 × 10⁻⁴) and Meizhou (2.1 × 10⁻⁵) (Figures 4C and 4D). This heritage-only pattern is compared with the ecological and supra-network results in Figures 5 and 6.

The 2020 ecological graph G_E consists of 233 nodes and 799 least-cost-path corridor edges. The network is a single connected component of density 0.030, mean degree k = 6.86, and mean clustering coefficient 0.083; the diameter measured in edge count is 12, and the mean shortest-path length between patch pairs is 2,212 (cumulative resistance units). The heritage graph G_H comprises 23 nodes and 73 edges (union of Delaunay triangulation and KNN-4), with a mean degree of 6.35 and a single connected component. The spatial distribution of the 23 heritage nodes together with the coupled inter-layer edges reveals three modal ICH clusters — a Ganzhou cluster centered on Longnan–Anyuan, a Minxi cluster centered on Yongding–Nanjing tulou territory, and a Meizhou cluster centered on Meixian district (Figure 5A).

Centrality analysis on G_E located the entire top-15 eigenvector-hub set inside Meizhou CEPZ (patch IDs 194–219), with patch 211 (a 9.3 km² Meizhou-central core) leading at eigenvector 0.37 and PageRank 0.006. The strong eigenvector concentration reflects the dense corridor tissue linking Meizhou's small, geometrically clustered forest patches through a low-cost matrix (Figure 5C). By contrast, mean eigenvector centrality in Ganzhou is only 1.6 × 10⁻⁴ and in Minxi 1.9 × 10⁻⁵, three orders of magnitude below Meizhou's 7.3 × 10⁻². Betweenness centrality, however, is more evenly distributed: Ganzhou attains the highest mean betweenness (0.034) because its larger patch inventory produces more shortest-path traffic through structurally intermediate nodes. This mismatch between eigenvector (Meizhou-dominated) and betweenness (Ganzhou-heavy) centralities is a distinctive marker of the tri-zonal topology.

Ecological single-layer centrality shows an eigenvector-betweenness contrast (Figure 6). The 2020 G_E network contains a compact Meizhou hub complex and a more diffuse Ganzhou structure (Figure 6A). Its degree distribution is right-skewed, with a mean degree of 6.86 and a maximum of 12 in patches P193-P219 (Figure 6B). Patch area and single-layer eigenvector centrality are negatively correlated (Spearman ρ = −0.21), so the highest-centrality patches are generally smaller Meizhou cores rather than large Minxi patches (Figure 6C). Zone means place. Meizhou is the highest in eigenvector centrality and PageRank, and Ganzhou is the highest in betweenness (Figure 6D). Because these panels use G_E alone, the Meizhou pattern is present before inter-layer coupling; comparison with the supra-network indicates that coupling is not its sole source.

Coupling G_H and G_E under the baseline rule produced 42 inter-layer edges: 35 strict 10 km links plus seven fallback links. Twenty-nine of the 233 ecological patches (12.4%) and all 23 ICH nodes have at least one inter-layer connection (mean ICH-to-ecological degree = 1.83; maximum = 5). The spectral radius of the mean-normalized symmetric supra-adjacency matrix is 19.6481. A previously mislabeled node-level centrality value has been removed from eigenvalue reporting. Zonal decomposition of coupling gives mean ecological linkage values of 2.33 for Meizhou, 1.91 for Ganzhou, and 1.17 for Minxi. These descriptive results identify Meizhou as the most strongly coupled zone under the specified distance-and-fallback rule.

Percolation resilience thresholds
Layer-specific percolation curves were first computed on the mapped heritage-inventory graph G_H (Figure 7). Under four progressive attacks, the remaining-node-normalized LCC ratio declined most slowly under random, betweenness, and eigenvector removal: thresholds were 0.96, 1.00, and 1.00, respectively, compared with 0.74 under degree-based removal (Figure 7A). Global-efficiency thresholds were 0.86 for random, 0.83 for eigenvector, and 0.57 for degree-based removal (Figure 7B). These results indicate high structural robustness of the represented inventory graph to random node removal and greater sensitivity to the removal of high-degree nodes. The consensus threshold is f_C(H) = 0.925 (Figure 7C); it does not measure continuity of heritage practice outside the mapped graph.

Percolation attacks on the two layers under four progressive node-removal schemes (Section 2.3.5) produced a marked asymmetry (Figure 8). Ecological-layer thresholds were 0.62 (random), 0.70 (degree), 0.44 (betweenness), and 1.00 (eigenvector), yielding f_C = 0.690 (Table 4). Heritage-inventory-layer thresholds were 0.96, 0.74, 1.00, and 1.00, respectively, yielding f_C = 0.925. The difference, Δf_C = 0.235, indicates earlier modeled loss of ecological-network integrity under three of the four attack rules. Under degree-based removal, the ecological and heritage-inventory layers cross S(k) = 0.5 at removed-node fractions of 0.70 and 0.74, respectively. Meizhou contains the ecological nodes selected earliest by the eigenvector-based attack; this is an association within the modeled topology, not evidence of a real-world causal cascade.

Raw threshold comparison is influenced by layer density. From the reported node and edge counts, ecological density is 2 × 799/(233 × 232) = 0.0296, whereas heritage-inventory density is 2 × 73/(23 × 22) = 0.2885. Dividing the consensus threshold by density gives 23.3 for the ecological layer and 3.21 for the heritage-inventory layer, a ratio of approximately 7.3:1. This descriptive normalization indicates that the higher raw heritage-inventory threshold partly reflects its denser graph. Because threshold-per-density is a comparative diagnostic rather than an intervention effect, it should not, by itself, be interpreted as evidence that adding edges or protecting nodes will produce a specified policy outcome.

Only under the eigenvector-attack rule do the two layers exhibit comparable robustness (both ≥ 0.98). Random, degree, and betweenness attacks locate the modeled collapse point earlier on the ecological layer than on the heritage-inventory layer. Agreement across three attack modes supports the stability of this directional result within the analyzed network and attack definitions, without implying general causal validity beyond those conditions.

Restoration priority mapping
The 2020 ecological corridor network used for RPI mapping contains 799 least-cost-path edges among 233 core patches (Figure 9A). Corridor length has a mean of 28.75 km, a median of 21.08 km, a 90th percentile of 52.20 km, and a maximum of 266.1 km (Figure 9B). Cumulative cost is similarly right-skewed, with a mean of 563.3, median of 259.6, and 90th percentile of 651.0 resistance-meter equivalents (Figure 9C). The positive length-cost association in Figure 9D indicates that traverse distance is an important component of modeled cost; local resistance, feasibility, and field condition remain necessary for evaluating any restoration corridor. Composite RPI scoring over 233 core patches yielded a heavy-tailed distribution (mean = 0, σ = 0.52, min = −2.71, max = 1.83). Forty-seven patches (20.2%) fell in analytical Tier 1 (RPI ≥ 0.290), 46 (19.7%) in the high tier (−0.184 ≤ RPI < 0.290), and 140 (60.1%) in the moderate tier. Meizhou contained 22 of 56 patches in Tier 1, compared with 20 of 144 in Ganzhou and 5 of 33 in Minxi. Combining the two upper analytical tiers gives 42 of 56 patches in Meizhou, 35 of 144 in Ganzhou, and 16 of 33 in Minxi. These tiers are relative rankings under the specified RPI weights, not prescriptive conservation-value categories (Figure 10). The top 15% of corridor RPI scores comprise 119 corridors. Among the 20 highest-ranked patches, 15 are in Meizhou, three in Ganzhou, and two in Minxi; together they cover 277 km². Their rankings reflect the combination of supra-eigenvector centrality, the negative patch-area term, ICH coupling, and isolation cost. Figure 10A maps this model-derived subset. The ranking is not a definitive restoration plan and should be combined with ecological condition, feasibility, land tenure, costs, and stakeholder priorities.

Table 5 summarizes the analytical priority-tier allocation across the three CEPZs. Meizhou contains 504.1 km² in Tier 1 across 22 patches, Ganzhou 439.9 km² across 20 patches, and Minxi 53.2 km² across five patches. Meizhou also has the highest mean RPI (+0.339). The larger aggregate area of Meizhou's high tier (4,095.7 km²) relative to Tier 1 reflects the RPI's negative area term, which raises the relative scores of small hub patches. These results describe structural leverage under the model; they do not establish intrinsic conservation value or a mandatory allocation of restoration resources.

Scenario simulation and sensitivity
The four scenarios produced contrasting model outcomes (Table 6). S1 reproduced the baseline consensus threshold of 0.690. S2, which removed all moderate-tier patches, reduced the threshold to 0.011, a 98.4% modeled decline. This result is consistent with a substantial topological contribution from patches classified as moderate; it does not constitute empirical proof that such loss will occur or prescribe a specific restoration tier. S3 and S4 returned a consensus of 0.690 because they changed edge weights without changing topology. Under this percolation definition, resistance reduction can improve weighted efficiency but does not alter the topological threshold. Adding or reconnecting stepping-stone patches is therefore a model-derived option for changing both topology and efficiency, rather than a mandatory intervention (Figure 11).

Gradient patch-loss scenarios produced a nonlinear modeled response. Removing 25% of moderate-tier patches reduced the consensus threshold from 0.690 to 0.593 (−14.0%); 50% removal yielded 0.483 (−30.0%); 75% yielded 0.312 (−54.8%); and 100% yielded 0.011 (−98.4%). The marginal modeled decline increased across the 25–50%, 50–75%, and 75–100% intervals. Within these simulations, retaining at least half of the moderate-tier patches was associated with preservation of more than 70% of the baseline threshold; this is a scenario result, not a forecast of real-world collapse. Weighted global efficiency (E_glob) was computed because the LCC-based threshold is insensitive to changes in edge weights. Baseline E_glob was 0.017230. S3 increased it to 0.018628 (+8.1%), and S4 increased it to 0.017754 (+3.0%). These modeled results indicate improved weighted connectivity even though the topological percolation threshold was unchanged. Thus, resistance reduction and topological expansion affect different network properties, and the analysis does not establish a universally superior intervention.

Sensitivity analysis supported the stability of the rankings within the tested parameter ranges. Perturbing the four RPI weights by ±0.05 and ±0.10 retained Spearman ρ ≥ 0.97. Varying the coupling radius over 5, 10, 15, and 20 km changed the strict-radius edge counts to 8, 35, 63, and 101, respectively; the 42-edge baseline at 10 km comprises 35 strict-radius and seven fallback links. Participation-coefficient correlations with the 10 km baseline were ρ = 0.73 at 15 km, 0.54 at 20 km, and 0.27 at 5 km. Seventeen of the top 20 RPI patches were retained at 15 km and 15 at 20 km. These results indicate parameter sensitivity and partial rank stability; they do not establish unrestricted transferability beyond the tested network (Figure 12).

DATA AVAILABILITY:
The China Land Cover Dataset is available from Zenodo (https://doi.org/10.5281/zenodo.4417810). The national intangible cultural heritage inventory is published by the State Council of China, and administrative boundary data are available from the National Geomatics Center of China. Derived network matrices, percolation outputs, and analysis code are deposited in Zenodo (https://doi.org/10.5281/zenodo.21732093).

Figure 1: Study area and CEPZ layout. (A) Locations of the three national CEPZs across southern Jiangxi, western Fujian, and eastern Guangdong. (B) Distribution of 23 national ICH items over CLCD 2020 land cover. (C) ICH counts by inscription batch. (D) ICH category composition. Maps are drawn using Open Street Map contributors as the basemap; administrative boundaries and all labels, symbols, and thematic layers were added or compiled by the authors. Panels (C) and (D) were prepared by the authors based on the study dataset. Please click here to view a larger version of this figure.

Figure 2: Analytical workflow of the DEHN framework. The five swim lanes show data assembly, dual-layer derivation, coupling, diagnostics, and outputs used to construct and assess the network. Please click here to view a larger version of this figure.

Figure 3: MSPA-lite fragmentation trajectory across the tri-CEPZ landscape, 2000–2023. (A) Spatial distribution of core patches by year and zone. (B) Temporal trends in total core-patch area, edge area, and total vegetation area. (C) Core-patch count and total core-patch area. Please click here to view a larger version of this figure.

Figure 4: Heritage-layer analysis of 23 national ICH items. (A) Kernel density surface. (B) Delaunay-KNN adjacency graph G_H. (C) Ten nodes with the highest betweenness centrality. (D) Centrality metrics by CEPZ. Please click here to view a larger version of this figure.

Figure 5: Bi-layer supra-network in 2020. (A) Spatial layout of inter-layer coupling. (B) ICH inter-layer degree distribution. (C) Twenty nodes with the highest supra-eigenvector centrality. (D) Participation coefficient versus supra-eigenvector centrality for all 256 nodes. Please click here to view a larger version of this figure.

Figure 6: Ecological-layer centrality analysis on G_E. (A) Spatial layout of the 2020 landscape. (B) Degree distribution. (C) Patch area versus eigenvector centrality. (D) Centrality metrics by zone. Please click here to view a larger version of this figure.

Figure 7: Heritage-inventory-layer percolation under targeted attacks. (A) LCC ratio versus removed-node fraction under four attack rules. (B) Global-efficiency decay. (C) Attack-specific and consensus thresholds. Please click here to view a larger version of this figure.

Figure 8: Ecological-layer percolation under four attack schemes. (A) Random attack. (B) Targeted attacks. (C) Cross-layer threshold comparison. Please click here to view a larger version of this figure.

Figure 9: Ecological corridor network in 2020. (A) Resistance surface. (B) The 799 least-cost-path corridors. (C) Cumulative-cost distribution. (D) Corridor length-cost relationship. Please click here to view a larger version of this figure.

Figure 10: Restoration priority index mapping. (A) Spatial distribution of RPI values and the top 15% of corridors. (B) Analytical tier composition by CEPZ. (C) RPI distribution by zone. (D) Component decomposition for the 20 highest-ranked patches. Please click here to view a larger version of this figure.

Figure 11: Scenario simulation of ecological-layer robustness. (A) Percolation curves under four scenarios. (B) Comparison of consensus critical thresholds. Please click here to view a larger version of this figure.

Figure 12: Sensitivity analysis. (A) Spearman rank correlations under RPI-weight perturbations. (B) Participation-coefficient stability across coupling radii. Please click here to view a larger version of this figure.

Attribute	Ganzhou CEPZ	Minxi CEPZ	Meizhou CEPZ	Total
Province	Jiangxi	Fujian	Guangdong	—
Area (km²)	39,341	19,353	15,853	74,547
County-level units	18 counties	6 counties	9 counties + 1 district	34
National-level ICH items (n)	11	6	6	23
Hakka-affiliated ICH items (n)	7	5	5	17
Dominant ICH categories	folk practices, traditional crafts	performing arts, folk practices	performing arts, traditional crafts	—
Core ecological patches ≥ 5 km² (2020)	144	33	56	233
Core-patch total area (km², 2020)	16,577.60	15,214.10	6,014.00	37,805.70

Table 1: Overview of the three Hakka CEPZs and their ICH inventories. The table compares geographic extent, administrative coverage, and national-level ICH counts across Ganzhou, Minxi, and Meizhou.

Data type	Source	Resolution / units	Time	Reference
Land cover (LULC)	China Land Cover Dataset (CLCD), Wuhan University	30 m raster	2000/05/10/15/20/23	Yang and Huang (2021)
CEPZ perimeters	Ministry of Culture and Tourism (MCT) national registry	Vector polygons	2013–2020 (declared)	MCT (2020)
National-level ICH inventory	Chinese State Council ICH National List (batches 1–5)	Point (county centroid)	2006–2021	State Council (2021)
Administrative boundaries	National Geomatics Center of China	Vector polygons	2020	NGCC (2020)
Coordinate system	Albers Conic Equal Area (lon₀ = 105°, φ₁ = 25°, φ₂ = 47°)	—	—	—

Table 2: Primary data sources. The table lists each dataset's provider, spatial or temporal resolution, and role in the analytical workflow.

CLCD class	Resistance value	Justification
Forest (2)	1	Base habitat; highest permeability
Shrub (3)	5	High permeability; secondary succession
Grassland (4)	10	Moderate permeability
Water (5)	30	Locally permeable to aquatic taxa; barrier to terrestrial
Cropland (1)	50	Semi-anthropogenic matrix
Ice/snow (7)	200	High-elevation barrier
Impervious (8)	500	Complete barrier to biotic flow
No-data (0)	100	Neutral placeholder

Table 3: Land-cover resistance values. The table reports the resistance assigned to each CLCD class for least-cost corridor modeling.

Attack rule	Ecological f*	Heritage f*	Δ (H − E)
Random (mean of 500)	0.62	0.96	0.34
Descending degree	0.7	0.74	0.04
Descending betweenness	0.44	1	0.56
Descending eigenvector	1	1	0
Consensus (mean)	0.69	0.925	0.235
Density-normalized (f_C/density)	23.3	3.21	−20.09
Random attack SD	0.058	0.082	0.024

Table 4: Percolation critical thresholds for the ecological and heritage-inventory layers in 2020. Attack-specific values and their consensus summarize structural robustness under the remaining-node LCC normalization.

CEPZ	Total patches	Tier 1 (top-ranked) (n / km²)	High tier (n / km²)	Moderate tier (n / km²)	Mean RPI
Ganzhou	144	20 / 439.9	15 / 475.3	109 / 15,662.5	−0.103
Minxi	33	5 / 53.2	11 / 123.0	17 / 15,037.9	−0.124
Meizhou	56	22 / 504.1	20 / 4,095.7	14 / 1,414.2	0.339
All three CEPZs	233	47 / 997.2	46 / 4,694.0	140 / 32,114.6	0

Table 5: RPI tier allocation across the three CEPZs. Patch counts, areas, and mean RPI values show the comparative distribution of analytical tiers by zone.

Scenario	Description	Consensus f*	Δ vs S1
S1	Baseline (unmodified G_E)	0.69	0
S2	Loss of moderate-tier (140 patches removed)	0.011	−0.679
S3	Halve cost on Tier 1 (top-ranked)–Tier 1 (top-ranked) edges	0.69	0
S4	Reduce cost 40% on top-119 corridors	0.69	0
S2a (25% moderate removed)	35 of 140 moderate patches removed	0.593	-0.097
S2b (50% moderate removed)	70 of 140 moderate patches removed	0.483	-0.207
S2c (75% moderate removed)	105 of 140 moderate patches removed	0.312	-0.378
Weighted global efficiency	S1=0.0172, S3=0.0186(+8.1%), S4=0.0178(+3.0%)	See text	-

Table 6: Scenario-simulation consensus percolation thresholds. The table compares the outcomes of the baseline, patch-loss, node-restoration, and corridor-restoration models.

Reserve	Nodes	Edges	Density	Random attack (mean ± SD; n = 500)	Degree	Betweenness	Eigenvector	Consensus
Ganzhou	144	359	0.035	0.420±0.079	0.326	0.118	0.632	0.374
Minxi	33	89	0.169	0.686±0.162	0.364	0.242	1.000	0.573
Meizhou	56	143	0.093	0.464±0.120	0.250	0.179	0.250	0.286

Table 7: Per-reserve percolation thresholds. The table reports attack-specific and consensus thresholds separately for Ganzhou, Minxi, and Meizhou.

Discussion

The tri-zonal decomposition of morphological and network diagnostics revealed marked spatial heterogeneity in ecological configuration⁴⁹. Ganzhou had the largest vegetated area and patch count but the smallest mean patch size, whereas Minxi retained the largest mean patches (461 km²), consistent with comparatively continuous Wuyi-fringe forest. Meizhou contained 56 patches within a smaller territory and showed the highest concentration of ecological eigenvector hubs.CLCD change analysis indicated an overall, non-linear change in core-patch area. Core-patch area increased from 44,485 km² in 2000 to 45,772 km² in 2010, before declining to 41,919 km² in 2015 and 37,806 km² in 2020. The 2010–2020 decline was 7,966km², equivalent to 17.4% of the 2010 core-patch area. Core-patch area was 37,888 km² in 2023, representing a small increase of 82.3 km² relative to 2020. Nevertheless, the overall 2000–2023 decline was 6,597 km², or 14.8%.Vegetated-core transitions to cropland accounted for 38% of net core loss, transitions associated with transport, reservoir, and industrial footprints for 31%, conversion to impervious land for 22%, and other mapped transitions for 9%. Ganzhou, Meizhou, and Minxi contributed 52%, 35%, and 13% of the net loss, respectively. These are land-cover accounting results and descriptive associations; urbanization, infrastructure development, orchard investment, depopulation, and policy processes are plausible contextual explanations but were not directly tested as causal drivers⁵⁰.

Per-reserve percolation analysis also identified substantial differences in modeled ecological robustness. Consensus thresholds were 0.374 for Ganzhou, 0.573 for Minxi, and 0.286 for Meizhou, while median random-attack thresholds across the 500 simulation replicates were 0.410, 0.667, and 0.446, respectively. By contrast, Table 7 reports the corresponding mean ± SD values of 0.420 ± 0.079, 0.686 ± 0.162, and 0.464 ± 0.120, respectively. Minxi therefore showed the highest modeled robustness, and Meizhou the lowest, under the specified network construction and attack rules. Ganzhou combined a larger impervious matrix with relatively intact interior cores and a higher mean betweenness (0.034), suggesting greater concentration of shortest-path traffic. Meizhou, by contrast, contained many small patches within locally dense subgraphs and exhibited stronger eigenvector centrality and local hub concentration. These differences describe the topology of the modeled corridor network rather than demonstrating that development pressure or depopulation caused the observed patterns^51,52.

Across the three zones, the ecological and mapped heritage-inventory layers exhibited divergent structural thresholds. The ecological layer reached the modeled collapse point at a consensus removed-node fraction of 0.690, compared with 0.925 for the heritage-inventory layer, a difference of 0.235. The ecological threshold was lower under random, degree, and betweenness attacks, whereas the two layers showed comparable robustness only under the eigenvector-based attack. This asymmetry suggests that, within the represented graphs, ecological-corridor integrity is the more restrictive structural component of the coupled system⁵³. However, the heritage layer consists only of 23 mapped national-level ICH items and should not be interpreted as a direct measure of the continuity, vitality, or geographic extent of cultural practices. The higher heritage threshold is also partly related to its much greater graph density (0.289 versus 0.030 for the ecological layer). Density-normalized thresholds provide a descriptive within-study comparison, but they should not be interpreted as evidence that increasing edge density or protecting a particular number of nodes will produce a predictable policy outcome⁵⁴.

The baseline coupling was limited but spatially uneven: 42 inter-layer links connected all 23 ICH nodes to 29 ecological patches, including 35 strict-radius links and seven nearest-patch fallback links. Meizhou had the highest mean ICH-to-ecological linkage (2.33) and contained the largest concentration of ecological hubs, making it both strongly coupled and structurally sensitive within the modeled network⁵⁵. The contrast between layers was also evident in centrality rankings: Minxi led the heritage-only eigenvector ranking, whereas Meizhou led the ecological and supra-network rankings. This inversion demonstrates that single-layer rankings may change after cross-layer coupling is introduced. Nevertheless, the results do not establish that one zone should automatically receive priority. In Meizhou, planners could assess protection or reconnection of small, highly central patches; in Ganzhou, intermediate patches with high betweenness could be examined alongside periurban land-use constraints; and in Minxi, buffering and consolidation of large continuous cores may be more relevant than adding numerous small patches⁵⁶. All such options require field validation, feasibility and cost assessment, land-tenure analysis, and stakeholder participation.

Scenario analysis clarified the distinction between weighted efficiency and topological robustness⁵⁷. Removing all moderate-tier patches reduced the consensus threshold from 0.690 to 0.011, while the progressive loss scenarios produced thresholds of 0.593, 0.483, 0.312, and 0.011 when 25%, 50%, 75%, and 100% of moderate-tier patches were removed, respectively. These results indicate that patches outside the highest analytical tiers can still make an important topological contribution. By contrast, reducing edge costs in the Tier-1 and top-priority-corridor restoration scenarios did not change the unweighted percolation threshold, although weighted global efficiency increased by 8.1% and 3.0%, respectively. Thus, resistance reduction and topological expansion affect different network properties: the former may improve modeled flow efficiency, whereas the latter is required to change the threshold under the present definition⁵⁸. The RPI ranking remained highly stable under the tested weight perturbations (Spearman’s ρ ≥ 0.97), but coupling-radius changes produced only partial stability, indicating that the priorities are useful screening outputs rather than definitive restoration prescriptions.

Several limitations delimit interpretation and point to future research. First, the resistance surface was based solely on land cover because complete DEM coverage was not available for the study extent; slope and topographic-wetness modifiers should be incorporated in future analyses to test whether the Meizhou hub pattern persists⁵⁹. Second, ICH items were geocoded to county centroids, which masks intra-county variation and may bias inter-layer coupling; village-level surveys, particularly in Meizhou, are needed to improve spatial representation⁶⁰. Third, the multilayer analysis was cross-sectional for 2020, although ecological fragmentation was documented from 2000 to 2023. Reconstructing the ecological and coupling networks for all benchmark years would support stronger temporal inference⁶¹. Fourth, the scenarios were deterministic and did not include land-use projection models such as PLUS or FLUS. Future work could integrate stochastic land-use simulation, restoration costs, ecological condition, ecosystem-service supply, and indicators of ICH practice continuity. Expanding the heritage layer to include provincial-level items would also reveal finer-scale cultural–ecological linkages⁶². Accordingly, the DEHN framework should be regarded as a reproducible comparative diagnostic: it identifies structural vulnerabilities and candidate intervention locations, but its application to other CEPZs or cultural landscapes requires locally reconstructed networks, consistent attack protocols, field validation, and explicit consideration of governance and community priorities⁶³.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Administrative boundary data	National Geomatics Center of China	China administrative boundary dataset; https://www.ngcc.cn/
ARPACK eigensolver	ARPACK-NG through SciPy	scipy.sparse.linalg.eigsh; https://github.com/opencollab/arpack-ng
China Land Cover Dataset (CLCD)	Wuhan University / Zenodo	30 m annual dataset; record 4417810; https://zenodo.org/records/4417810
Cultural Ecological Protection Zone registry	Ministry of Culture and Tourism of China	National CEPZ registry; https://www.mct.gov.cn/
Delaunay triangulation and k-nearest-neighbour analysis	SciPy / NetworkX	Heritage adjacency-network construction; k = 4
Gaode POI services	Amap / Gaode	Online POI service; https://lbs.amap.com/
GeoPandas	GeoPandas developers / PyPI	Version 0.14; https://geopandas.org/
Google Earth imagery	Google	Google Earth imagery; https://earth.google.com/
Kernel density estimation	Python scientific computing environment	500 m grid; 5 km bandwidth
Least-cost path algorithm	scikit-image project	Dijkstra algorithm through route_through_array
National Intangible Cultural Heritage Inventory	State Council of the People's Republic of China	National-level inventory, batches 1-5
NetworkX	NetworkX developers / PyPI	Version 3.2; https://networkx.org/
Python	Python Software Foundation	Version 3.11; https://www.python.org/
rasterio	Rasterio developers / PyPI	Version 1.3; https://rasterio.readthedocs.io/
scikit-image	scikit-image developers / PyPI	skimage.graph.route_through_array; https://scikit-image.org/
SciPy sparse	SciPy community	scipy.sparse; https://scipy.org/
Shuttle Radar Topography Mission DEM	NASA / USGS	SRTM DEM; preliminary assessment only; incomplete study-area coverage
Zenodo analysis repository	Zenodo	Code, derived matrices, and outputs; https://doi.org/10.5281/zenodo.21732093